Update testing macro with access to geometry to get cell indices, correct place of...
[u/mrichter/AliRoot.git] / EMCAL / macros / TestESD.C
CommitLineData
cd259404 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
9480984d 31//Change the bool depending on what information you want to print
32// when all FALSE, prints minimum cluster information.
33Bool_t kPrintKine = kTRUE; //Do not use for raw data.
34Bool_t kPrintCaloCells = kFALSE;
35Bool_t kPrintTrackMatches = kFALSE;
36Bool_t kPrintClusterCells = kFALSE;
37
cd259404 38void TestESD() {
39
a6bf15ef 40 TGeoManager::Import("geometry.root");
41 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETE") ;
42
cd259404 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++) {
9480984d 53 cout << "<<<< Event: " << iev+1 << "/" << nEvt << " >>>>"<<endl;
cd259404 54 esdTree->GetEvent(iev);
9480984d 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
cd259404 81 TRefArray* caloClusters = new TRefArray();
82 esd->GetEMCALClusters(caloClusters);
83
cd259404 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: "
a6bf15ef 107 << cphi*TMath::RadToDeg() << " Eta: " << ceta << " NCells: " << nCells << " #Matches: " << nMatched
cd259404 108 << " Index: " << trackIndex << " #Labels: " << nLabels << " Index: "
109 << labelIndex << endl;
110
9480984d 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;
cd259404 117
9480984d 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);
cd259404 123
9480984d 124 Double_t pOverE = tmom/energy;
cd259404 125
9480984d 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 }
cd259404 129
9480984d 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;
a6bf15ef 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 }
9480984d 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
a6bf15ef 166 cout<<" More contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<"; Eta "<<particle->Eta()<<"; Phi "<<particle->Phi()*TMath::RadToDeg()<<endl;
9480984d 167 for(Int_t i = 1; i < nLabels; i++){
a6bf15ef 168 particle = stack->Particle((clus->GetLabels())->At(i));
9480984d 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 }
cd259404 176
177
178 }
179
180
181
182}