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