e8b3ce964b4a47f7ce47cf1d56781846975ea7ce
[u/mrichter/AliRoot.git] / 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 //Change the bool depending on what information you want to print
32 // when all FALSE, prints minimum cluster information.
33 Bool_t kPrintKine = kTRUE; //Do not use for raw data.
34 Bool_t kPrintCaloCells = kFALSE;
35 Bool_t kPrintTrackMatches = kFALSE;
36 Bool_t kPrintClusterCells = kFALSE;
37
38 void TestESD() {
39
40   TFile* f = new TFile("AliESDs.root");
41   TTree* esdTree = (TTree*)f->Get("esdTree");
42   
43   AliESDEvent* esd = new AliESDEvent();
44   esd->ReadFromTree(esdTree);
45
46   Int_t nEvt = esdTree->GetEntries();
47   Float_t pos[3] = {0.,0.,0.};
48
49   for(Int_t iev = 0; iev < nEvt; iev++) {
50     cout << "<<<< Event: " << iev+1 << "/" << nEvt << " >>>>"<<endl;
51     esdTree->GetEvent(iev);
52           
53         //In case you want to play with MC data
54         AliStack *stack = 0;
55         if(kPrintKine){  
56                 AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),  "read");
57                 rl->LoadKinematics();  
58                 rl->GetEvent(iev);
59                 stack = rl->Stack();
60         }  
61           
62         //get reconstructed vertex position
63         Double_t vertex_position[3];
64         esd->GetVertex()->GetXYZ(vertex_position);
65           
66         //GetCellsClusters Array  
67     AliESDCaloCells &cells= *(esd->GetEMCALCells());
68           // Uncomment to see the full list of digit amplitudes and times.
69         if(kPrintCaloCells){  
70                 Int_t nTotalCells = cells.GetNumberOfCells() ;  
71                 Int_t type        = cells.GetType();
72                 for (Int_t icell=  0; icell <  nTotalCells; icell++) {
73                         cout<<"Cell   : "<<icell<<"/"<<nTotalCells<<" ID: "<<cells.GetCellNumber(icell)<<" Amplitude: "<<cells.GetAmplitude(icell)<<" Time: "<<cells.GetTime(icell)<<endl;        
74                 }// cell loop
75         }
76           
77         //GetCaloClusters Array
78     TRefArray* caloClusters = new TRefArray();
79     esd->GetEMCALClusters(caloClusters);
80
81     //loop over clusters
82     Int_t nclus = caloClusters->GetEntries();
83     for (Int_t icl = 0; icl < nclus; icl++) {
84
85       AliESDCaloCluster* clus = (AliESDCaloCluster*)caloClusters->At(icl);
86       Float_t energy = clus->E();
87       clus->GetPosition(pos);
88       TVector3 vpos(pos[0],pos[1],pos[2]);
89       TLorentzVector p;
90       clus->GetMomentum(p,vertex_position);
91       Double_t cphi = vpos.Phi();
92       Double_t ceta = vpos.Eta();
93
94       Int_t nMatched = clus->GetNTracksMatched();
95       Int_t trackIndex = clus->GetTrackMatched();
96       Int_t nLabels = clus->GetNLabels();
97       Int_t labelIndex = clus->GetLabel();
98
99       Int_t nCells = clus->GetNCells();
100
101       //For later: ADD CHECK THAT CLUSTER IS WITHIN SM FIDUCIAL VOLUME
102
103       cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << " Phi: " 
104            << cphi << " Eta: " << ceta << " NCells: " << nCells << " #Matches: " << nMatched 
105            << " Index: " << trackIndex << " #Labels: " << nLabels << " Index: " 
106            << labelIndex << endl;
107
108         if(kPrintTrackMatches && trackIndex >= 0) {
109                 AliESDtrack* track = esd->GetTrack(trackIndex);
110                 Double_t tphi = track->GetOuterParam()->Phi();
111                 Double_t teta = track->GetOuterParam()->Eta();
112                 Double_t tmom = track->GetOuterParam()->P();
113                 cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;
114
115                 Double_t deta = teta - ceta;
116                 Double_t dphi = tphi - cphi;
117                 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
118                 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
119                 Double_t dR = sqrt(dphi*dphi + deta*deta);
120
121                 Double_t pOverE = tmom/energy;
122
123                 if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
124                         cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
125                 }
126
127         }
128                 
129         //Get CaloCells of cluster and print
130         if(kPrintClusterCells){ 
131                 UShort_t * index = clus->GetCellsAbsId() ;
132                 Double_t * fraction = clus->GetCellsAmplitudeFraction() ;
133                 for(Int_t i = 0; i < nCells ; i++){
134                         Int_t absId =   index[i]; // or clus->GetCellNumber(i) ;
135                         Double_t ampFract =  fraction[i];
136                         Float_t amp       = cells.GetCellAmplitude(absId) ;
137                         Float_t time      = cells.GetCellTime(absId);
138                         cout<<"         Cluster Cell: AbsID : "<< absId << "; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time<<endl;
139                 }
140         }
141                 
142         //Print primary info
143         if(!stack || !kPrintKine) continue;
144         if(labelIndex >= 0 && labelIndex < stack->GetNtrack()){
145                 TParticle * particle = stack->Particle(labelIndex);
146                 //Print primary values
147                 cout<<"         More  contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;   
148                 for(Int_t i = 1; i < nLabels; i++){
149                         particle = stack->Particle(clus->(GetLabels()->At(i)));
150                         cout<<"         Other contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;
151                 }
152         }
153         else if( labelIndex >= stack->GetNtrack()) cout <<"PROBLEM, label is too large : "<<labelIndex<<" >= particles in stack "<< stack->GetNtrack() <<endl;
154         else cout<<"Negative label!!!  : "<<labelIndex<<endl;
155                 
156         }
157
158
159   }
160
161
162
163 }