EMCAL/macros/TestESD.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2
   3 //Root include files
   4 #include <Riostream.h>
   5 #include <TFile.h>
   6 #include <TChain.h>
   7 #include <TParticle.h>
   8 #include <TNtuple.h>
   9 #include <TCanvas.h>
  10 #include <TObjArray.h>
  11 #include <TSystem.h>
  12 #include <TString.h>
  13 #include <TH1F.h>
  14 #include <TVector.h>
  15 #include <TParticle.h>
  16 #include <TRefArray.h>
  17 #include <TArrayS.h>
  18
  19 //AliRoot include files
  20 #include "AliRunLoader.h"
  21 #include "AliStack.h"
  22 #include "AliESDEvent.h"
  23 #include "AliESDVertex.h"
  24 #include "AliESDCaloCluster.h"
  25 #include "AliESDCaloCells.h"
  26 #include "AliPID.h"
  27 #include "AliLog.h"
  28
  29 #endif
  30
  31 void TestESD() {
  32
  33   TFile* f = new TFile("AliESDs.root");
  34   TTree* esdTree = (TTree*)f->Get("esdTree");
  35
  36   AliESDEvent* esd = new AliESDEvent();
  37   esd->ReadFromTree(esdTree);
  38
  39   Int_t nEvt = esdTree->GetEntries();
  40   Float_t pos[3] = {0.,0.,0.};
  41
  42   for(Int_t iev = 0; iev < nEvt; iev++) {
  43     cout << "Event: " << iev+1 << "/" << nEvt << endl;
  44     esdTree->GetEvent(iev);
  45
  46     TRefArray* caloClusters = new TRefArray();
  47     esd->GetEMCALClusters(caloClusters);
  48
  49     //get reconstructed vertex position
  50     Double_t vertex_position[3];
  51     esd->GetVertex()->GetXYZ(vertex_position);
  52
  53     //loop over clusters
  54     Int_t nclus = caloClusters->GetEntries();
  55     for (Int_t icl = 0; icl < nclus; icl++) {
  56
  57       AliESDCaloCluster* clus = (AliESDCaloCluster*)caloClusters->At(icl);
  58       Float_t energy = clus->E();
  59       clus->GetPosition(pos);
  60       TVector3 vpos(pos[0],pos[1],pos[2]);
  61       TLorentzVector p;
  62       clus->GetMomentum(p,vertex_position);
  63       Double_t cphi = vpos.Phi();
  64       Double_t ceta = vpos.Eta();
  65
  66       Int_t nMatched = clus->GetNTracksMatched();
  67       Int_t trackIndex = clus->GetTrackMatched();
  68       Int_t nLabels = clus->GetNLabels();
  69       Int_t labelIndex = clus->GetLabel();
  70
  71       Int_t nCells = clus->GetNCells();
  72
  73       //For later: ADD CHECK THAT CLUSTER IS WITHIN SM FIDUCIAL VOLUME
  74
  75       cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << " Phi: "
  76            << cphi << " Eta: " << ceta << " NCells: " << nCells << " #Matches: " << nMatched
  77            << " Index: " << trackIndex << " #Labels: " << nLabels << " Index: "
  78            << labelIndex << endl;
  79
  80       if(trackIndex >= 0) {
  81         AliESDtrack* track = esd->GetTrack(trackIndex);
  82         Double_t tphi = track->GetOuterParam()->Phi();
  83         Double_t teta = track->GetOuterParam()->Eta();
  84         Double_t tmom = track->GetOuterParam()->P();
  85         cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;
  86
  87         Double_t deta = teta - ceta;
  88         Double_t dphi = tphi - cphi;
  89         if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
  90         if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
  91         Double_t dR = sqrt(dphi*dphi + deta*deta);
  92
  93         Double_t pOverE = tmom/energy;
  94
  95         if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
  96           cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
  97         }
  98
  99       }
 100
 101     }
 102
 103
 104   }
 105
 106
 107
 108 }