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