fill the fraction of cell energy in cluster in case of unfolding ON
[u/mrichter/AliRoot.git] / EMCAL / macros / TestESD.C
CommitLineData
c8fe2783 1
cd259404 2#if !defined(__CINT__) || defined(__MAKECINT__)
3
4//Root include files
c8fe2783 5//#include <Riostream.h>
cd259404 6#include <TFile.h>
c8fe2783 7//#include <TSystem.h>
cd259404 8#include <TH1F.h>
cd259404 9#include <TParticle.h>
10#include <TRefArray.h>
cd259404 11
12//AliRoot include files
13#include "AliRunLoader.h"
14#include "AliStack.h"
15#include "AliESDEvent.h"
16#include "AliESDVertex.h"
17#include "AliESDCaloCluster.h"
18#include "AliESDCaloCells.h"
19#include "AliPID.h"
c8fe2783 20#include "AliEMCALGeometry.h"
cd259404 21
22#endif
23
9480984d 24//Change the bool depending on what information you want to print
25// when all FALSE, prints minimum cluster information.
c8fe2783 26Bool_t kPrintKine = kFALSE; //Do not use for raw data.
27Bool_t kPrintCaloCells = kFALSE;
9480984d 28Bool_t kPrintTrackMatches = kFALSE;
29Bool_t kPrintClusterCells = kFALSE;
c8fe2783 30Bool_t kPrintClusterPID = kFALSE;
9480984d 31
cd259404 32void TestESD() {
c8fe2783 33 // Main method to read information stored in AliESDCaloClusters and AliESDCaloCells
34
35 // Init some example histograms
36 // ESD
37 TH1F * hEta = (TH1F*) new TH1F("hEta","reco eta",1000, -0.71,0.71);
38 TH1F * hPhi = (TH1F*) new TH1F("hPhi","reco phi",130, 70,200);
39 TH1F * hE = (TH1F*) new TH1F("hE" ,"reco e",300, 0,30);
40 TH1F * hTime = (TH1F*) new TH1F("hTime" ,"reco time",1000, 0,1000);
41 hEta ->SetXTitle("#eta");
42 hPhi ->SetXTitle("#phi (deg)");
43 hE ->SetXTitle("E (GeV)");
44 hTime->SetXTitle("time (ns)");
45
46 // Monte Carlo
47 TH1F * hMCEta = (TH1F*) new TH1F("hMCEta","MC eta",1000, -0.71,0.71);
48 TH1F * hMCPhi = (TH1F*) new TH1F("hMCPhi","MC phi",130, 70,200);
49 TH1F * hMCE = (TH1F*) new TH1F("hMCE" ,"MC e",300, 0,30);
50 hMCEta->SetXTitle("#eta");
51 hMCPhi->SetXTitle("#phi (deg)");
52 hMCE ->SetXTitle("E (GeV)");
53
54 //ESD - MonteCarlo
55 TH1F * hDEta = (TH1F*) new TH1F("hDEta"," eta cluster - eta MC",500, -0.05,0.05);
56 TH1F * hDPhi = (TH1F*) new TH1F("hDPhi","phi cluster - phi MC",200, -20,20);
57 TH1F * hDE = (TH1F*) new TH1F("hDE" ,"e cluster - e MC",200, -10,10);
58 hDEta->SetXTitle("#eta_{reco}-#eta_{MC}");
59 hDPhi->SetXTitle("#phi_{reco}-#phi_{MC} (deg)");
60 hDE ->SetXTitle("E_{reco}-E_{MC} (GeV)");
61
62 //Open the ESD file, get the tree with events
cd259404 63 TFile* f = new TFile("AliESDs.root");
64 TTree* esdTree = (TTree*)f->Get("esdTree");
cd259404 65 AliESDEvent* esd = new AliESDEvent();
66 esd->ReadFromTree(esdTree);
67
c8fe2783 68 //Init geometry and array that will contain the clusters
69 TRefArray* caloClusters = new TRefArray();
70 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEAR") ;
71 Float_t pos[3];
72
73 //Loop of events
cd259404 74 Int_t nEvt = esdTree->GetEntries();
cd259404 75 for(Int_t iev = 0; iev < nEvt; iev++) {
9480984d 76 cout << "<<<< Event: " << iev+1 << "/" << nEvt << " >>>>"<<endl;
cd259404 77 esdTree->GetEvent(iev);
c8fe2783 78
79 //Pass the geometry transformation matrix from ESDs to geometry
80 for(Int_t mod=0; mod < (geom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
81 //printf("matrix %d\n",mod);
82 if(esd->GetEMCALMatrix(mod)) {
83 //printf("EMCAL: mod %d, matrix %p\n",mod, esd->GetEMCALMatrix(mod));
84 //(esd->GetEMCALMatrix(mod))->Print("");
85 geom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
86 }//matrix
87 }//module
88
89 //In case you want to play with MC data, get stack
90 AliStack *stack = 0;
91 if(kPrintKine){
92 AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(), "read");
93 rl->LoadKinematics();
94 rl->GetEvent(iev);
95 stack = rl->Stack();
96 }
97
98 //get reconstructed vertex position
99 Double_t vertex_position[3];
100 esd->GetVertex()->GetXYZ(vertex_position);
101
102 //------------------------------------------------------
103 //Get Cells Array, all cells in event, print cells info
104 //------------------------------------------------------
105 AliVCaloCells &cells= *(esd->GetEMCALCells());
106 //AliESDCaloCells &cells= *(esd->GetEMCALCells());
cd259404 107
c8fe2783 108 if(kPrintCaloCells){
109 Int_t nTotalCells = cells.GetNumberOfCells() ;
110 //Int_t type = cells.GetType();
111 for (Int_t icell= 0; icell < nTotalCells; icell++) {
112 cout<<"Cell : "<<icell<<"/"<<nTotalCells<<" ID: "<<cells.GetCellNumber(icell)<<" Amplitude: "<<cells.GetAmplitude(icell)<<" Time: "<<cells.GetTime(icell)*1e9<<endl;
113 }// cell loop
114 }
115
116 //------------------------------------------------------
117 // Calo Clusters
118 //------------------------------------------------------
119
120 //Get CaloClusters Array
121 caloClusters->Clear();
122 esd->GetEMCALClusters(caloClusters);
123
cd259404 124 //loop over clusters
125 Int_t nclus = caloClusters->GetEntries();
126 for (Int_t icl = 0; icl < nclus; icl++) {
c8fe2783 127
128 AliVCluster* clus = (AliVCluster*)caloClusters->At(icl);
129 //AliESDCluster* clus = (AliESDCluster*)caloClusters->At(icl);
cd259404 130 Float_t energy = clus->E();
131 clus->GetPosition(pos);
132 TVector3 vpos(pos[0],pos[1],pos[2]);
c8fe2783 133
134 //We can get a momentum TLorentzVector per cluster, corrected by the vertex position
135 //TLorentzVector p;
136 //clus->GetMomentum(p,vertex_position);
137
cd259404 138 Double_t cphi = vpos.Phi();
139 Double_t ceta = vpos.Eta();
c8fe2783 140
141 Int_t nMatched = clus->GetNTracksMatched();
142 Int_t trackIndex = clus->GetTrackMatchedIndex();
143 Int_t nLabels = clus->GetNLabels();
cd259404 144 Int_t labelIndex = clus->GetLabel();
c8fe2783 145 Int_t nCells = clus->GetNCells();
9480984d 146
c8fe2783 147 //Fill some histograms
148 hEta->Fill(ceta);
149 hPhi->Fill(cphi*TMath::RadToDeg());
150 hE ->Fill(energy);
151 hTime->Fill(clus->GetTOF()*1e9);
152
153 //Print basic cluster information
154 cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << "; Phi: "
155 << cphi*TMath::RadToDeg() << "; Eta: " << ceta << "; NCells: " << nCells << " ;#Matches: " << nMatched
156 << "; Index: " << trackIndex << "; #Labels: " << nLabels << " Index: "
157 << labelIndex << "; Time "<<clus->GetTOF()*1e9<<" ns "<<endl;
158
159 //Print primary info
160 if(stack && kPrintKine) {
9480984d 161 if(labelIndex >= 0 && labelIndex < stack->GetNtrack()){
c8fe2783 162 TParticle * particle = stack->Particle(labelIndex);
163 //Fill histograms with primary info
164 hMCEta->Fill(particle->Eta());
165 hMCPhi->Fill(particle->Phi()*TMath::RadToDeg());
166 hMCE ->Fill(particle->Energy());
167 hDEta ->Fill(ceta-particle->Eta());
168 hDPhi ->Fill(cphi*TMath::RadToDeg()-particle->Phi()*TMath::RadToDeg());
169 hDE ->Fill(energy-particle->Energy());
170 //Print primary values
171 cout<<" More contributing primary: "<<particle->GetName()<<"; with kinematics: "<<endl;
172 cout<<" \t Energy: "<<particle->Energy()<<"; Phi: "<<particle->Phi()*TMath::RadToDeg()<<"; Eta: "<<particle->Eta()<<endl;
173 for(Int_t i = 1; i < nLabels; i++){
174 //particle = stack->Particle((((AliESDCaloCluster*)clus)->GetLabelsArray())->At(i));
175 particle = stack->Particle((clus->GetLabels())[i]);
176 //or Int_t *labels = clus->GetLabels();
177 //particle = stack->Particle(labels[i]);
178 cout<<" Other contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;
179 }
9480984d 180 }
181 else if( labelIndex >= stack->GetNtrack()) cout <<"PROBLEM, label is too large : "<<labelIndex<<" >= particles in stack "<< stack->GetNtrack() <<endl;
c8fe2783 182 else cout<<"Negative label!!! : "<<labelIndex<<endl;
183 } // play with stack
184
185 // Matching results
186 if(kPrintTrackMatches && trackIndex >= 0) {
187 AliESDtrack* track = esd->GetTrack(trackIndex);
188 Double_t tphi = track->GetOuterParam()->Phi();
189 Double_t teta = track->GetOuterParam()->Eta();
190 Double_t tmom = track->GetOuterParam()->P();
191 cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl;
192
193 Double_t deta = teta - ceta;
194 Double_t dphi = tphi - cphi;
195 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
196 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
197 Double_t dR = sqrt(dphi*dphi + deta*deta);
198
199 Double_t pOverE = tmom/energy;
200
201 if(dR < 0.02 && pOverE < 1.8 && nCells > 1) {
202 cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl;
9480984d 203 }
c8fe2783 204 }// matching
205
206 //Get PID weights and print them
207 if(kPrintClusterPID){
208 const Double_t *pid = clus->GetPID();
209 printf("PID weights: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
210 pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
211 pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
212 pid[AliVCluster::kPion], pid[AliVCluster::kKaon], pid[AliVCluster::kProton],
213 pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
214 }//PID
215
216 //Get CaloCells of cluster and print their info, position.
217 if(kPrintClusterCells){
218 UShort_t * index = clus->GetCellsAbsId() ;
219 Double_t * fraction = clus->GetCellsAmplitudeFraction() ;
220 Int_t sm = -1;
221 for(Int_t i = 0; i < nCells ; i++){
222 Int_t absId = index[i]; // or clus->GetCellNumber(i) ;
223 Double_t ampFract = fraction[i];
224 Float_t amp = cells.GetCellAmplitude(absId) ;
225 Double_t time = cells.GetCellTime(absId);
226 cout<<"\t Cluster Cell: AbsID : "<< absId << " == "<<clus->GetCellAbsId(i) <<"; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time*1e9<<endl;
227 //Geometry methods
228 Int_t iSupMod = 0 ;
229 Int_t iTower = 0 ;
230 Int_t iIphi = 0 ;
231 Int_t iIeta = 0 ;
232 Int_t iphi = 0 ;
233 Int_t ieta = 0 ;
234 if(geom){
235 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
236 //Gives SuperModule and Tower numbers
237 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
238 iIphi, iIeta,iphi,ieta);
239 //Gives label of cell in eta-phi position per each supermodule
240 Float_t cellPhi = 0;
241 Float_t cellEta = 0;
242 geom->EtaPhiFromIndex(absId,cellEta,cellPhi);
243 cout<< " SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta
244 <<"; Phi "<<iIphi<<"; Index: Cell Eta "<<ieta<<"; Cell Phi "<<iphi
245 <<"; Global: Cell Eta "<<cellEta<<"; Cell Phi "<<cellPhi*TMath::RadToDeg()<<endl;
246 if(i==0) sm = iSupMod;
247 else{
248 if(sm!=iSupMod) printf("******CLUSTER SHARED BY 2 SuperModules!!!!\n");
249 }
250 }// geometry on
251 }// cluster cell loop
252 }// print cell clusters
253 } //cluster loop
254 } // event loop
255
256
257 //Write histograms in a file
258 TFile * fhisto = (TFile*) new TFile("histos.root","recreate");
259 hEta->Write();
260 hPhi->Write();
261 hE ->Write();
262 hTime->Write();
263 if(kPrintKine){
264 hMCEta->Write();
265 hMCPhi->Write();
266 hMCE ->Write();
267 hDEta->Write();
268 hDPhi->Write();
269 hDE ->Write();
cd259404 270 }
c8fe2783 271 fhisto->Close();
cd259404 272}