]>
Commit | Line | Data |
---|---|---|
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" | |
754a26d5 | 19 | #include "AliESDCaloTrigger.h" |
cd259404 | 20 | #include "AliPID.h" |
c8fe2783 | 21 | #include "AliEMCALGeometry.h" |
cd259404 | 22 | |
23 | #endif | |
24 | ||
9480984d | 25 | //Change the bool depending on what information you want to print |
26 | // when all FALSE, prints minimum cluster information. | |
c8fe2783 | 27 | Bool_t kPrintKine = kFALSE; //Do not use for raw data. |
28 | Bool_t kPrintCaloCells = kFALSE; | |
754a26d5 | 29 | Bool_t kPrintCaloTrigger = kFALSE; |
9480984d | 30 | Bool_t kPrintTrackMatches = kFALSE; |
31 | Bool_t kPrintClusterCells = kFALSE; | |
c8fe2783 | 32 | Bool_t kPrintClusterPID = kFALSE; |
9480984d | 33 | |
cd259404 | 34 | void TestESD() { |
c8fe2783 | 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 | |
cd259404 | 65 | TFile* f = new TFile("AliESDs.root"); |
66 | TTree* esdTree = (TTree*)f->Get("esdTree"); | |
cd259404 | 67 | AliESDEvent* esd = new AliESDEvent(); |
68 | esd->ReadFromTree(esdTree); | |
69 | ||
c8fe2783 | 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 | |
cd259404 | 76 | Int_t nEvt = esdTree->GetEntries(); |
cd259404 | 77 | for(Int_t iev = 0; iev < nEvt; iev++) { |
9480984d | 78 | cout << "<<<< Event: " << iev+1 << "/" << nEvt << " >>>>"<<endl; |
cd259404 | 79 | esdTree->GetEvent(iev); |
c8fe2783 | 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 | //AliESDCaloCells &cells= *(esd->GetEMCALCells()); | |
cd259404 | 109 | |
c8fe2783 | 110 | if(kPrintCaloCells){ |
111 | Int_t nTotalCells = cells.GetNumberOfCells() ; | |
112 | //Int_t type = cells.GetType(); | |
113 | for (Int_t icell= 0; icell < nTotalCells; icell++) { | |
114 | cout<<"Cell : "<<icell<<"/"<<nTotalCells<<" ID: "<<cells.GetCellNumber(icell)<<" Amplitude: "<<cells.GetAmplitude(icell)<<" Time: "<<cells.GetTime(icell)*1e9<<endl; | |
115 | }// cell loop | |
116 | } | |
117 | ||
754a26d5 | 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 | ||
c8fe2783 | 142 | //------------------------------------------------------ |
143 | // Calo Clusters | |
144 | //------------------------------------------------------ | |
145 | ||
146 | //Get CaloClusters Array | |
147 | caloClusters->Clear(); | |
148 | esd->GetEMCALClusters(caloClusters); | |
149 | ||
cd259404 | 150 | //loop over clusters |
151 | Int_t nclus = caloClusters->GetEntries(); | |
152 | for (Int_t icl = 0; icl < nclus; icl++) { | |
c8fe2783 | 153 | |
154 | AliVCluster* clus = (AliVCluster*)caloClusters->At(icl); | |
155 | //AliESDCluster* clus = (AliESDCluster*)caloClusters->At(icl); | |
cd259404 | 156 | Float_t energy = clus->E(); |
157 | clus->GetPosition(pos); | |
158 | TVector3 vpos(pos[0],pos[1],pos[2]); | |
c8fe2783 | 159 | |
160 | //We can get a momentum TLorentzVector per cluster, corrected by the vertex position | |
161 | //TLorentzVector p; | |
162 | //clus->GetMomentum(p,vertex_position); | |
163 | ||
cd259404 | 164 | Double_t cphi = vpos.Phi(); |
165 | Double_t ceta = vpos.Eta(); | |
c8fe2783 | 166 | |
167 | Int_t nMatched = clus->GetNTracksMatched(); | |
168 | Int_t trackIndex = clus->GetTrackMatchedIndex(); | |
169 | Int_t nLabels = clus->GetNLabels(); | |
cd259404 | 170 | Int_t labelIndex = clus->GetLabel(); |
c8fe2783 | 171 | Int_t nCells = clus->GetNCells(); |
9480984d | 172 | |
c8fe2783 | 173 | //Fill some histograms |
174 | hEta->Fill(ceta); | |
175 | hPhi->Fill(cphi*TMath::RadToDeg()); | |
176 | hE ->Fill(energy); | |
177 | hTime->Fill(clus->GetTOF()*1e9); | |
178 | ||
179 | //Print basic cluster information | |
180 | cout << "Cluster: " << icl+1 << "/" << nclus << " Energy: " << energy << "; Phi: " | |
181 | << cphi*TMath::RadToDeg() << "; Eta: " << ceta << "; NCells: " << nCells << " ;#Matches: " << nMatched | |
182 | << "; Index: " << trackIndex << "; #Labels: " << nLabels << " Index: " | |
183 | << labelIndex << "; Time "<<clus->GetTOF()*1e9<<" ns "<<endl; | |
184 | ||
185 | //Print primary info | |
186 | if(stack && kPrintKine) { | |
9480984d | 187 | if(labelIndex >= 0 && labelIndex < stack->GetNtrack()){ |
c8fe2783 | 188 | TParticle * particle = stack->Particle(labelIndex); |
189 | //Fill histograms with primary info | |
190 | hMCEta->Fill(particle->Eta()); | |
191 | hMCPhi->Fill(particle->Phi()*TMath::RadToDeg()); | |
192 | hMCE ->Fill(particle->Energy()); | |
193 | hDEta ->Fill(ceta-particle->Eta()); | |
194 | hDPhi ->Fill(cphi*TMath::RadToDeg()-particle->Phi()*TMath::RadToDeg()); | |
195 | hDE ->Fill(energy-particle->Energy()); | |
196 | //Print primary values | |
197 | cout<<" More contributing primary: "<<particle->GetName()<<"; with kinematics: "<<endl; | |
198 | cout<<" \t Energy: "<<particle->Energy()<<"; Phi: "<<particle->Phi()*TMath::RadToDeg()<<"; Eta: "<<particle->Eta()<<endl; | |
199 | for(Int_t i = 1; i < nLabels; i++){ | |
200 | //particle = stack->Particle((((AliESDCaloCluster*)clus)->GetLabelsArray())->At(i)); | |
201 | particle = stack->Particle((clus->GetLabels())[i]); | |
202 | //or Int_t *labels = clus->GetLabels(); | |
203 | //particle = stack->Particle(labels[i]); | |
204 | cout<<" Other contributing primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl; | |
205 | } | |
9480984d | 206 | } |
207 | else if( labelIndex >= stack->GetNtrack()) cout <<"PROBLEM, label is too large : "<<labelIndex<<" >= particles in stack "<< stack->GetNtrack() <<endl; | |
c8fe2783 | 208 | else cout<<"Negative label!!! : "<<labelIndex<<endl; |
209 | } // play with stack | |
210 | ||
211 | // Matching results | |
212 | if(kPrintTrackMatches && trackIndex >= 0) { | |
213 | AliESDtrack* track = esd->GetTrack(trackIndex); | |
214 | Double_t tphi = track->GetOuterParam()->Phi(); | |
215 | Double_t teta = track->GetOuterParam()->Eta(); | |
216 | Double_t tmom = track->GetOuterParam()->P(); | |
217 | cout << "\t Track Momentum: " << tmom << " phi: " << tphi << " eta: " << teta << endl; | |
218 | ||
219 | Double_t deta = teta - ceta; | |
220 | Double_t dphi = tphi - cphi; | |
221 | if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi(); | |
222 | if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi(); | |
223 | Double_t dR = sqrt(dphi*dphi + deta*deta); | |
224 | ||
225 | Double_t pOverE = tmom/energy; | |
226 | ||
227 | if(dR < 0.02 && pOverE < 1.8 && nCells > 1) { | |
228 | cout << "\n\t Excellent MATCH! dR = " << dR << " p/E = " << pOverE << " nCells = " << nCells << endl; | |
9480984d | 229 | } |
c8fe2783 | 230 | }// matching |
231 | ||
232 | //Get PID weights and print them | |
233 | if(kPrintClusterPID){ | |
234 | const Double_t *pid = clus->GetPID(); | |
235 | 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", | |
236 | pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0], | |
237 | pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon], | |
238 | pid[AliVCluster::kPion], pid[AliVCluster::kKaon], pid[AliVCluster::kProton], | |
239 | pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]); | |
240 | }//PID | |
241 | ||
242 | //Get CaloCells of cluster and print their info, position. | |
243 | if(kPrintClusterCells){ | |
244 | UShort_t * index = clus->GetCellsAbsId() ; | |
245 | Double_t * fraction = clus->GetCellsAmplitudeFraction() ; | |
246 | Int_t sm = -1; | |
247 | for(Int_t i = 0; i < nCells ; i++){ | |
248 | Int_t absId = index[i]; // or clus->GetCellNumber(i) ; | |
249 | Double_t ampFract = fraction[i]; | |
250 | Float_t amp = cells.GetCellAmplitude(absId) ; | |
251 | Double_t time = cells.GetCellTime(absId); | |
252 | cout<<"\t Cluster Cell: AbsID : "<< absId << " == "<<clus->GetCellAbsId(i) <<"; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time*1e9<<endl; | |
253 | //Geometry methods | |
254 | Int_t iSupMod = 0 ; | |
255 | Int_t iTower = 0 ; | |
256 | Int_t iIphi = 0 ; | |
257 | Int_t iIeta = 0 ; | |
258 | Int_t iphi = 0 ; | |
259 | Int_t ieta = 0 ; | |
260 | if(geom){ | |
261 | geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); | |
262 | //Gives SuperModule and Tower numbers | |
263 | geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, | |
264 | iIphi, iIeta,iphi,ieta); | |
265 | //Gives label of cell in eta-phi position per each supermodule | |
266 | Float_t cellPhi = 0; | |
267 | Float_t cellEta = 0; | |
268 | geom->EtaPhiFromIndex(absId,cellEta,cellPhi); | |
269 | cout<< " SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta | |
270 | <<"; Phi "<<iIphi<<"; Index: Cell Eta "<<ieta<<"; Cell Phi "<<iphi | |
271 | <<"; Global: Cell Eta "<<cellEta<<"; Cell Phi "<<cellPhi*TMath::RadToDeg()<<endl; | |
272 | if(i==0) sm = iSupMod; | |
273 | else{ | |
274 | if(sm!=iSupMod) printf("******CLUSTER SHARED BY 2 SuperModules!!!!\n"); | |
275 | } | |
276 | }// geometry on | |
277 | }// cluster cell loop | |
278 | }// print cell clusters | |
279 | } //cluster loop | |
280 | } // event loop | |
281 | ||
282 | ||
283 | //Write histograms in a file | |
284 | TFile * fhisto = (TFile*) new TFile("histos.root","recreate"); | |
285 | hEta->Write(); | |
286 | hPhi->Write(); | |
287 | hE ->Write(); | |
288 | hTime->Write(); | |
289 | if(kPrintKine){ | |
290 | hMCEta->Write(); | |
291 | hMCPhi->Write(); | |
292 | hMCE ->Write(); | |
293 | hDEta->Write(); | |
294 | hDPhi->Write(); | |
295 | hDE ->Write(); | |
cd259404 | 296 | } |
c8fe2783 | 297 | fhisto->Close(); |
cd259404 | 298 | } |