HarpoAnalyseXtalk.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseXtalk.cxx
//
#include "HarpoAnalyseXtalk.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(HarpoAnalyseXtalk);

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

    Double_t qchprev = 0;
    for(Int_t j=0;j<NCHAN;j++){ // Channels
      Double_t qch = 0;
      for(Int_t i=0;i<NADC;i++){ // Time bins
        Double_t q = m->GetData(j,i);
        if(q <= -1000) continue;
        qch += q;
        if(j<=0) continue;
        Double_t qn = m->GetData(j-1,i);
        hQvsQneighbour[ndet]->Fill(q,qn/q);
        if(j<=3) continue;
        Double_t qn3 = m->GetData(j-3,i);
        hQvsQneighbour3[ndet]->Fill(q,qn3/q);
        if(j%2 == 0)
          hQvsQneighbour3e[ndet]->Fill(q,qn3/q);
        else
          hQvsQneighbour3o[ndet]->Fill(q,qn3/q);
      }
      hQchVsQneighbour[ndet]->Fill(qch,qchprev/qch);
      qchprev = qch;
    }
  }
  
  // if(gHDetSet->isExist(PMM2)) { // PMm2 (trigger) data
  //   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
  //   for(imes=0;imes<esize;imes++) {
  //     int ch    = pm->at(imes).getChNum(); // channel
  //     int q     = pm->at(imes).getCharge(); // charge
  //     hPmm2->Fill(ch,q);
  //   } 
  // }

  // //???????????????????????????????
  // if(fEvt->GetRecoEvent()){  // Reconstructed data (clusters, tracks, matched tracks)
  //   HarpoRecoEvent* reco = fEvt->GetRecoEvent();
  //   // Clusters in the event (both X and Y directions)
  //   TClonesArray* clArray = reco->GetClustersArray();
  //   // Tracks in the event (both X and Y directions)
  //   TClonesArray* trArray = reco->GetTracksArray();
  //   // 3D tracks in the event (after X-Y matching)
  //   TClonesArray* tr3dArray = reco->GetReco3DArray();
  //   Int_t NtrX=reco->GetNtracksXEvt();
  //   Int_t NtrY=reco->GetNtracksYEvt();
  //   Int_t nCl = clArray->GetEntries();
  //   Int_t nTr = trArray->GetEntries();
  //   Int_t nTr3d = tr3dArray->GetEntries();
  //   for(Int_t icl = 0; icl<nCl; icl++){
  //     HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
  //     Double_t q = cluster->GetQ();
  //     Int_t plane = cluster->GetPlane();
  //     hQcl[plane]->Fill(q);
  //   }
  //   for(Int_t itr = 0; itr<nTr; itr++){
  //     HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(itr);
  //     Double_t id = track->GetNtrack();
  //     Double_t angle = track->GetAngleTrack();
  //     hAngle->Fill(angle);
  //   }
  //   for(Int_t itr = 0; itr<nTr3d; itr++){
  //     HarpoRecoReco3D* track = (HarpoRecoReco3D*)tr3dArray->At(itr);
  //   }
  // }




  //  if (gHDetSet->isExist(SIMDET)){
  //    HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
  //    Int_t nIonTr = simEvt->GetNtracks();
  //    TH1F* hTmp = new TH1F("hTmp","",512,0,512);
  //    for(Int_t i = 0; i<nIonTr; i++){
  //      TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);
  //      Int_t nPoints = tr->GetNpoints();
  //      for(Int_t j = 0; j<nPoints; j++){
  //     TpcSimIonisationPoint* point = tr->GetPoint(j);
  //     hTmp->Fill(point->GetZ(),point->GetNe());
  //      }
  //    }
  //    hTmp->Delete();
  //  }
}

void   HarpoAnalyseXtalk::Init()
{
  
  // Initialise histograms here
  hQvsQneighbour[0] = new TH2F("hQvsQneighbourX",";index;index;",4096,0,4096,200,0,5);
  hQvsQneighbour[1] = new TH2F("hQvsQneighbourY",";index;index;",4096,0,4096,200,0,5);
  hQvsQneighbour3[0] = new TH2F("hQvsQneighbour3X",";index;index;",4096,0,4096,200,0,5);
  hQvsQneighbour3[1] = new TH2F("hQvsQneighbour3Y",";index;index;",4096,0,4096,200,0,5);
  hQvsQneighbour3e[0] = new TH2F("hQvsQneighbour3eX",";index;index;",4096,0,4096,200,0,5);
  hQvsQneighbour3e[1] = new TH2F("hQvsQneighbour3eY",";index;index;",4096,0,4096,200,0,5);
  hQvsQneighbour3o[0] = new TH2F("hQvsQneighbour3oX",";index;index;",4096,0,4096,200,0,5);
  hQvsQneighbour3o[1] = new TH2F("hQvsQneighbour3oY",";index;index;",4096,0,4096,200,0,5);


   Int_t fNbins = 500;
   Double_t fQmin = 10, fQmax = 1e5;
   
  const Int_t nbinsQcl = fNbins;
  Double_t xminQcl = fQmin;
  Double_t xmaxQcl = fQmax;
  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] = xminQcl + TMath::Exp(logxminQcl+i*binwidthQcl);

  hQchVsQneighbour[0] = new TH2F("hQchVsQneighbourX",";index;index;",nbinsQcl,xbinsQcl,200,0,5);
  hQchVsQneighbour[1] = new TH2F("hQchVsQneighbourY",";index;index;",nbinsQcl,xbinsQcl,200,0,5);
 }

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


   /* TFile *hf = */ OpenHistFile("xtalk");
   // Save histograms here
   hQvsQneighbour[0]->Write();
   hQvsQneighbour[1]->Write();
   hQvsQneighbour3[0]->Write();
   hQvsQneighbour3[1]->Write();
   hQvsQneighbour3e[0]->Write();
   hQvsQneighbour3e[1]->Write();
   hQvsQneighbour3o[0]->Write();
   hQvsQneighbour3o[1]->Write();
   hQchVsQneighbour[0]->Write();
   hQchVsQneighbour[1]->Write();

 }