HarpoAnalyseVdrift.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseVdrift.cxx
//
#include "HarpoAnalyseVdrift.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 "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>

ClassImp(HarpoAnalyseVdrift);

void   HarpoAnalyseVdrift::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   HarpoAnalyseVdrift::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  

  if(!gHDetSet->isExist(PMM2)) return;
  
  Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
  //  Int_t triggerType = anaEvt->GetTriggerType(); // Trigger line (with mesh, without, traversing tracks, ...)
  
  int imes;
  int esize = anaEvt->GetHeader()->eventSize;   // number of PMs hit
  Pmm2MesVect * pm = anaEvt->GetMesurements();  // data for each included PM
  Int_t hit[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  for(imes=0;imes<esize;imes++) {
    int ch    = pm->at(imes).getChNum(); // channel
    if(ch<0 || ch>=16) continue;
    hit[ch] = 1;
    // int q    = pm->at(imes).getCharge(); // charge
    // hPmm2->Fill(ch,q);
  } 
  Int_t isTT = 0;
  if((hit[3] || hit[4]) && (hit[5] || hit[6]))
    isTT = 1;
  
  // Process RAW data
  for(Int_t ndet=0;ndet<2;ndet++) { // Detectors 0 and 1 = X and Y planes
    if (!gHDetSet->isExist(ndet)) continue;
    HarpoDccMap* m = fEvt->GetDccMap(ndet);
    if ( m == NULL ) continue;

    Int_t tmin = 512, tmax = 0;
    for(Int_t i=0;i<NADC;i++){ // Time bins
      Double_t qtot = 0;
      for(Int_t j=0;j<NCHAN;j++){ // Channels
        Double_t q = m->GetData(j,i);
        if(q <= -1000) continue;
        qtot += q;
      }
      if(qtot<=0) continue;
      if(isTT)
        hQvT[ndet]->Fill(i,qtot);
      if(qtot<=1000) continue;
      if(i>tmax) tmax = i;
      if(i<tmin) tmin = i;
    }

    hTmin[ndet]->Fill(tmin);
    hTmax[ndet]->Fill(tmax);
    if(isTT)
      hDeltaT[ndet]->Fill(tmax-tmin);

    fEvt->GetRecoEvent()->SetTmin(tmin);
    fEvt->GetRecoEvent()->SetTmax(tmax);
  }
  
}

void   HarpoAnalyseVdrift::Init()
{
  
  // Initialise histograms here
  hTmin[0] = new TH1F("hTminX","",512,0,512);
  hTmin[1] = new TH1F("hTminY","",512,0,512);
  hTmax[0] = new TH1F("hTmaxX","",512,0,512);
  hTmax[1] = new TH1F("hTmaxY","",512,0,512);
  hDeltaT[0] = new TH1F("hDeltaTX","",512,0,512);
  hDeltaT[1] = new TH1F("hDeltaTY","",512,0,512);
  
  hQvT[0]  = new TH1F("hQvTX","",512,0,512);
  hQvT[1]  = new TH1F("hQvTY","",512,0,512);

}

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


   /* TFile *hf = */ OpenHistFile("vdrift");
   // Save histograms here
   hQvT[0]->Write();
   hQvT[1]->Write();
   hTmin[0]->Write();
   hTmin[1]->Write();
   hTmax[0]->Write();
   hTmax[1]->Write();
   hDeltaT[0]->Write();
   hDeltaT[1]->Write();
 
 }