HarpoSelectorPileUp.cxx

Go to the documentation of this file.

// 
// File HarpoSelectorPileUp.cxx
//
#include "HarpoSelectorPileUp.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "Pmm2Event.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.h"
#include "HarpoSimEvent.h"
#include "HarpoRecoMonitorGui.h"

#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TF1.h"
#include "TMath.h"
#include "TSystem.h"

#include <cstdlib>
#include <cstring>
#include <cassert>
#include <fstream>
#include <iostream>

using namespace std;


ClassImp(HarpoSelectorPileUp);

void   HarpoSelectorPileUp::print()
 {
   unsigned int nd; // number of detectors
   HarpoEventHeader *hdr = fEvt->GetHeader();

   assert(hdr != NULL);
   hdr->print();
   
   for (nd = 0; nd < gkNDetectors; nd++) {
     //     if (fEvt->isdataExist(nd)) {
       HarpoDccMap *plane = fEvt->GetDccMap(nd);
       if (plane != NULL )plane->print();
   }
 }

void   HarpoSelectorPileUp::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;

  // if(gHDetSet->isExist(PMM2)){
  //   Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
  //   if(anaEvt->GetTriggerType() != 7){ // FIXME
  //     fEvt->GetHeader()->SetEvAnaStatus(-2);
  //     return;
  //   }
  // }

  for(Int_t plane = 0; plane<2; plane++){
    if (!gHDetSet->isExist(plane)) continue;
    // Process RAW data
    HarpoDccMap* m = fEvt->GetDccMap(plane);
    if ( m == NULL ) continue;
    //    Int_t i2 = 0;
    Double_t qTotIn = 0;
    Double_t qTotBefore = 0;
    Double_t qTotAfter = 0;
    for(Int_t j=0;j<NCHAN;j++){ // Channels
      for(Int_t i=0;i<NADC;i++){ //Time bins
        Double_t q = m->GetData(j,i);
        if(q<0) continue;
        //std::cout << i << "-";
        if(i< fTstart            ) qTotBefore += q;
        if(              i >fTend) qTotAfter += q;
        if(i>=fTstart && i<=fTend) qTotIn += q;
      }
    }

    // Info("process","%d: %g %g %g %d %d",plane,qTotBefore,qTotIn,qTotAfter,fTstart,fTend);

    hQbefore[plane]->Fill(qTotBefore);
    hQafter[plane]->Fill(qTotAfter);
    hQin[plane]->Fill(qTotIn);
    // Double_t data[4];
    // data[0] = plane;
    // data[1] = qTotIn;
    // data[2] = qTotBefore;
    // data[3] = qTotAfter;
    HarpoRecoEvent* reco = fEvt->GetRecoEvent();
    Double_t data[7] = {(Double_t) plane,qTotIn,qTotBefore,qTotAfter,
                        (Float_t) reco->GetXstart(),
                        (Float_t) reco->GetYstart(),
                        (Float_t) reco->GetTstart()};
    ntuple->Fill(data);

    if(qTotBefore>fThr || qTotAfter>fThr)
      fEvt->GetHeader()->SetEvAnaStatus(-2);
  }

}

void   HarpoSelectorPileUp::Init()
{
  
  // Initialise histograms here
  fTstart = 90;
  fTend = 430;
   
   Long64_t tstart;
   if( ! gHConfig->Lookup("pileup.tstart",tstart ))
     Info("Constructor","Use default fTstart %d",fTstart);
   else
     fTstart = tstart;

   Long64_t tend;
   if( ! gHConfig->Lookup("pileup.tend",tend ))
     Info("Constructor","Use default fTend %d",fTend);
   else
     fTend = tend;


   fThr = 1e4;
   Double_t thr;
   if( ! gHConfig->Lookup("pileup.thr",thr ))
     Info("Constructor","Use default fThr %g",fThr);
   else
     fThr = thr;

  const Int_t nbinsQcl = 100;
  Double_t xminQcl = 10;
  Double_t xmaxQcl = 1e5;
  Double_t logxminQcl = TMath::Log(xminQcl);
  Double_t logxmaxQcl = TMath::Log(xmaxQcl);
  Double_t binwidthQcl = (logxmaxQcl-logxminQcl)/nbinsQcl;
  Double_t xbinsQcl[nbinsQcl+1];
  xbinsQcl[0] = xminQcl;
  for (Int_t i=1;i<=nbinsQcl;i++)
    xbinsQcl[i] = TMath::Exp(logxminQcl+i*binwidthQcl);
  hQbefore[0] = new TH1F("hQbeforeX",";Q_{before}",nbinsQcl,xbinsQcl);
  hQbefore[1] = new TH1F("hQbeforeY",";Q_{before}",nbinsQcl,xbinsQcl);
  hQafter[0] = new TH1F("hQafterX",";Q_{after}",nbinsQcl,xbinsQcl);
  hQafter[1] = new TH1F("hQafterY",";Q_{after}",nbinsQcl,xbinsQcl);
  hQin[0] = new TH1F("hQinX",";Q_{in}",nbinsQcl,xbinsQcl);
  hQin[1] = new TH1F("hQinY",";Q_{in}",nbinsQcl,xbinsQcl);

  ntuple = new TNtupleD("ntuple","pile up data","plane:qIn:qBefore:qAfter:xS:yS:tS");
  ntuple->SetDirectory(0);

 }

void   HarpoSelectorPileUp::Save(char * /* mode */)
 {


   OpenHistFile("selectorPileup");

   // Save histograms here
   hQbefore[0]->Write();
   hQbefore[1]->Write();
   hQafter[0]->Write();
   hQafter[1]->Write();
   hQin[0]->Write();
   hQin[1]->Write();

   ntuple->Write();

 }