]>
Commit | Line | Data |
---|---|---|
62246135 | 1 | #include "TChain.h" |
2 | #include "TTree.h" | |
3 | #include "TH1F.h" | |
4 | #include "TCanvas.h" | |
5 | #include "TFile.h" | |
6 | #include "TObjArray.h" | |
7 | ||
8 | #include "AliAnalysisTask.h" | |
9 | #include "AliAnalysisManager.h" | |
10 | ||
11 | #include "AliESDEvent.h" | |
12 | #include "AliESDInputHandler.h" | |
13 | #include "AliAODEvent.h" | |
14 | #include "AliAODInputHandler.h" | |
15 | ||
16 | #include "AliVEvent.h" | |
17 | #include "AliVCaloCells.h" | |
18 | #include "AliVCluster.h" | |
19 | ||
20 | #include "AliEMCALGeometry.h" | |
21 | #include "AliGeomManager.h" | |
22 | ||
23 | #include "AliAnalysisTaskEMCALIsolation.h" | |
24 | ||
25 | #include "AliStack.h" | |
26 | #include "AliMCEventHandler.h" | |
27 | #include "AliMCEvent.h" | |
28 | #include "AliMCParticle.h" | |
29 | ||
30 | #include "TDatabasePDG.h" | |
31 | #include "TParticlePDG.h" | |
32 | ||
33 | // Task for isolating gammas with EMCAL | |
34 | // Author: Marco Marquard | |
35 | ||
36 | ClassImp(AliAnalysisTaskEMCALIsolation) | |
37 | ||
38 | //________________________________________________________________________ | |
39 | ||
40 | AliAnalysisTaskEMCALIsolation::AliAnalysisTaskEMCALIsolation() | |
41 | : AliAnalysisTaskSE(), | |
42 | bVerbose(0), | |
43 | bMC(0), | |
44 | fESD(0), | |
45 | fAOD(0), | |
46 | fMC(0), | |
47 | fStack(0), | |
48 | fOutputList(0), | |
49 | fHistGlobalHmap(0), | |
50 | fHistGlobalHmap0(0), | |
51 | fTreeEvent(0), | |
52 | fTreeCluster(0), | |
53 | emclus(0), | |
54 | prodrad(0), | |
55 | contPID(-999), | |
56 | mothPID(-999), | |
57 | trackmatch(0), | |
58 | fGeom(0), | |
59 | fGeoName("EMCAL_COMPLETE"), | |
60 | fESDCells(0), | |
61 | fAODCells(0), | |
62 | fEsdClusters(0), | |
63 | fAodClusters(0) | |
64 | { | |
65 | //default constructor | |
66 | DefineInput(0, TChain::Class()); | |
67 | DefineOutput(1, TList::Class()); | |
68 | } | |
69 | ||
70 | //________________________________________________________________________ | |
71 | AliAnalysisTaskEMCALIsolation::AliAnalysisTaskEMCALIsolation(const char *name) | |
72 | : AliAnalysisTaskSE(name), | |
73 | bVerbose(0), | |
74 | bMC(0), | |
75 | fESD(0), | |
76 | fAOD(0), | |
77 | fMC(0), | |
78 | fStack(0), | |
79 | fOutputList(0), | |
80 | fHistGlobalHmap(0), | |
81 | fHistGlobalHmap0(0), | |
82 | fTreeEvent(0), | |
83 | fTreeCluster(0), | |
84 | emclus(0), | |
85 | prodrad(0), | |
86 | contPID(-999), | |
87 | mothPID(-999), | |
88 | trackmatch(0), | |
89 | fGeom(0), | |
90 | fGeoName("EMCAL_COMPLETE"), | |
91 | fESDCells(0), | |
92 | fAODCells(0), | |
93 | fEsdClusters(0), | |
94 | fAodClusters(0) | |
95 | { | |
96 | // Constructor | |
97 | ||
98 | // Define input and output slots here | |
99 | // Input slot #0 works with a TChain | |
100 | DefineInput(0, TChain::Class()); | |
101 | // Output slot #0 id reserved by the base class for AOD | |
102 | // Output slot #1 writes into a TH1 container | |
103 | DefineOutput(1, TList::Class()); | |
104 | } | |
105 | ||
106 | //________________________________________________________________________ | |
107 | AliAnalysisTaskEMCALIsolation::~AliAnalysisTaskEMCALIsolation() | |
108 | { | |
109 | //Destructor | |
110 | ||
111 | delete fESD; | |
112 | delete fAOD; | |
113 | delete fOutputList; | |
114 | delete fHistGlobalHmap; | |
115 | delete fHistGlobalHmap0; | |
116 | delete fTreeCluster; | |
117 | delete fTreeEvent; | |
118 | } | |
119 | ||
120 | //________________________________________________________________________ | |
121 | void AliAnalysisTaskEMCALIsolation::UserCreateOutputObjects() | |
122 | { | |
123 | // Create histograms | |
124 | // Called once | |
125 | Double_t etagaps[97]= {-0.66687,-0.653,-0.63913,-0.62526,-0.61139,-0.59752,-0.58365,-0.56978,-0.55591,-0.54204,-0.52817,-0.5143,-0.50043,-0.48656,-0.47269,-0.45882,-0.44495,-0.43108,-0.41721,-0.40334,-0.38947,-0.3756,-0.36173,-0.34786,-0.33399,-0.32012,-0.30625,-0.29238,-0.27851,-0.26464,-0.25077,-0.2369,-0.22303,-0.20916,-0.19529,-0.18142,-0.16755,-0.15368,-0.13981,-0.12594,-0.11207,-0.0982,-0.08433,-0.07046,-0.05659,-0.04272,-0.02885,-0.01498,-0.00111,0.01276,0.02663,0.0405,0.05437,0.06824,0.08211,0.09598,0.10985,0.12372,0.13759,0.15146,0.16533,0.1792,0.19307,0.20694,0.22081,0.23468,0.24855,0.26242,0.27629,0.29016,0.30403,0.3179,0.33177,0.34564,0.35951,0.37338,0.38725,0.40112,0.41499,0.42886,0.44273,0.4566,0.47047,0.48434,0.49821,0.51208,0.52595,0.53982,0.55369,0.56756,0.58143,0.5953,0.60917,0.62304,0.63691,0.65078,0.66465}; | |
126 | ||
127 | Double_t phigaps[125]= {1.408,1.4215,1.435,1.4485,1.462,1.4755,1.489,1.5025,1.516,1.5295,1.543,1.5565,1.57,1.5835,1.597,1.6105,1.624,1.6375,1.651,1.6645,1.678,1.6915,1.705,1.7185,1.732,1.758,1.7715,1.785,1.7985,1.812,1.8255,1.839,1.8525,1.866,1.8795,1.893,1.9065,1.92,1.9335,1.947,1.9605,1.974,1.9875,2.001,2.0145,2.028,2.0415,2.055,2.0685,2.082,2.108,2.1215,2.135,2.1485,2.162,2.1755,2.189,2.2025,2.216,2.2295,2.243,2.2565,2.27,2.2835,2.297,2.3105,2.324,2.3375,2.351,2.3645,2.378,2.3915,2.405,2.4185,2.432,2.456,2.4695,2.483,2.4965,2.51,2.5235,2.537,2.5505,2.564,2.5775,2.591,2.6045,2.618,2.6315,2.645,2.6585,2.672,2.6855,2.699,2.7125,2.726,2.7395,2.753,2.7665,2.78,2.804,2.8175,2.831,2.8445,2.858,2.8715,2.885,2.8985,2.912,2.9255,2.939,2.9525,2.966,2.9795,2.993,3.0065,3.02,3.0335,3.047,3.0605,3.074,3.0875,3.101,3.1145,3.128}; | |
128 | ||
129 | fOutputList = new TList(); | |
130 | ||
131 | fHistGlobalHmap = new TH2F("fHistGlobalHmap","Hit map for EMCAL",96,etagaps,124,phigaps); | |
132 | fHistGlobalHmap->GetXaxis()->SetTitle("#eta"); | |
133 | fHistGlobalHmap->GetYaxis()->SetTitle("#phi"); | |
134 | fOutputList->Add(fHistGlobalHmap); | |
135 | ||
136 | fHistGlobalHmap0 = new TH2F("fHistGlobalHmap0","Hit map for EMCAL",96,0,96,120,0,120); | |
137 | fHistGlobalHmap0->GetXaxis()->SetTitle("#eta"); | |
138 | fHistGlobalHmap0->GetYaxis()->SetTitle("#phi"); | |
139 | fOutputList->Add(fHistGlobalHmap0); | |
140 | /* | |
141 | fHist = new TH1F("","",,,,); | |
142 | fHist->GetXaxis()->SetTitle(""); | |
143 | fHist->GetYaxis()->SetTitle(""); | |
144 | fHist->SetMarkerStyle(kFullCircle); | |
145 | fOutputList->Add(fHist); | |
146 | */ | |
147 | ||
148 | fTreeEvent = new TTree("eventsT","eventsT"); | |
149 | fTreeEvent->Branch("emclus", &emclus); | |
150 | ||
151 | fTreeCluster = new TTree("clusterT","clusterT"); | |
152 | fTreeCluster->Branch("emclus", &emclus); | |
153 | fTreeCluster->Branch("prodrad", &prodrad); | |
154 | fTreeCluster->Branch("contPID", &contPID); | |
155 | fTreeCluster->Branch("mothPID", &mothPID); | |
156 | ||
157 | fOutputList->Add(fHistGlobalHmap); | |
158 | fOutputList->Add(fHistGlobalHmap0); | |
159 | fOutputList->Add(fTreeEvent); | |
160 | fOutputList->Add(fTreeCluster); | |
161 | ||
162 | ||
163 | ||
164 | } | |
165 | ||
166 | //________________________________________________________________________ | |
167 | void AliAnalysisTaskEMCALIsolation::UserExec(Option_t *) | |
168 | { | |
169 | ||
170 | AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); | |
171 | ||
172 | ||
173 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler()); | |
174 | if (!esdH) | |
175 | { | |
176 | Printf("ERROR: Could not get ESDInputHandler"); | |
177 | } | |
178 | ||
179 | AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler()); | |
180 | if (!aodH && !esdH) | |
181 | { | |
182 | Printf("ERROR: Could not get AODInputHandler"); | |
183 | } | |
184 | ||
185 | if(esdH){ | |
186 | fESD = esdH->GetEvent(); | |
187 | am->LoadBranch("AliESDRun."); | |
188 | am->LoadBranch("AliESDHeader."); | |
189 | } | |
190 | else if(aodH){ | |
191 | fAOD = aodH->GetEvent(); | |
192 | } | |
193 | else{ | |
194 | AliFatal("Neither ESD nor AOD event found"); | |
195 | return; | |
196 | } | |
197 | ||
198 | AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); | |
199 | if (!mcH) { | |
200 | Printf("ERROR: Could not retrieve MC event handler"); | |
201 | return; | |
202 | } | |
203 | fMC = mcH->MCEvent(); | |
204 | if (!fMC) { | |
205 | Printf("ERROR: Could not retrieve MC event"); | |
206 | return; | |
207 | } | |
208 | // stack to get truth infos of particles: | |
209 | fStack = fMC->Stack(); | |
210 | if (!fStack) { | |
211 | Printf("ERROR: Could not retrieve stack from MC event"); | |
212 | return; | |
213 | } | |
214 | //Printf(" MC particles: %d", fMC->GetNumberOfTracks()); | |
215 | ||
216 | ||
217 | fGeom = AliEMCALGeometry::GetInstance(fGeoName); | |
218 | //AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry(); | |
219 | //Double_t phimin = emc->GetArm1PhiMin(); | |
220 | //Double_t phimax = emc->GetArm1PhiMax(); | |
221 | //printf("EMCAL coverage in phi from: %f to: %f \n",phimin,phimax); | |
222 | if (!fGeom){ | |
223 | AliFatal("EMCAL(fGeom) geometry is not initialized"); | |
224 | return; | |
225 | } | |
226 | ||
227 | AliVCaloCells *cells = 0; | |
228 | ||
229 | if (fESD){ | |
230 | am->LoadBranch("EMCALCells."); | |
231 | fESDCells = fESD->GetEMCALCells(); | |
232 | cells = fESDCells; | |
233 | } | |
234 | else{ | |
235 | am->LoadBranch("emcalCells"); | |
236 | fAODCells = fAOD->GetEMCALCells(); | |
237 | cells = fAODCells; | |
238 | } | |
239 | ||
240 | ||
241 | ||
242 | ||
243 | if (!cells){ | |
244 | AliFatal("No cells loaded"); | |
245 | return; | |
246 | } | |
247 | ||
248 | if(fESD){ | |
249 | am->LoadBranch("CaloClusters"); | |
250 | TList *l = fESD->GetList(); | |
251 | fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters")); | |
252 | } | |
253 | else if(fAOD){ | |
254 | am->LoadBranch("caloClusters"); | |
255 | TList *l = fAOD->GetList(); | |
256 | fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters")); | |
257 | } | |
258 | else{ | |
259 | AliFatal("Impossible to not have either pointer to ESD or AOD event"); | |
260 | } | |
261 | ||
262 | TObjArray *clusters = fEsdClusters; | |
263 | if (!clusters) | |
264 | clusters = fAodClusters; | |
265 | if (!clusters){ | |
266 | Printf("ERROR: no clusters node in event!"); | |
267 | //return; | |
268 | } | |
269 | ||
270 | //========================================= start of analysis ============================================ | |
271 | ||
272 | ||
273 | //============> DATA LOOP =================> | |
274 | ||
275 | //------------> cell loop -----------------> | |
276 | Int_t ncells = cells->GetNumberOfCells(); | |
277 | //for (Int_t i = 0; i<15000; ++i){ | |
278 | // Int_t absID = i; | |
279 | for (Int_t i = 0; i<ncells; ++i){ | |
280 | Int_t absID = TMath::Abs(cells->GetCellNumber(i)); | |
281 | ||
282 | //Check if this absId exists | |
283 | if(!(fGeom->CheckAbsCellId(absID))){ | |
284 | continue; | |
285 | } | |
286 | Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta; | |
287 | // Get from the absId the super module number, the module number and the eta-phi index (0 or 1) in the module | |
288 | fGeom->GetCellIndex(absID, nSupMod, nModule, nIphi, nIeta); | |
289 | ||
290 | // Get from the the super module number, the module number and the eta-phi index (0 or 1) in the module the tower row (iphi) and column (ieta) | |
291 | fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta); | |
292 | ||
293 | if(bVerbose) | |
294 | printf("ID: %d ; supermodule: %d ; module %d ; phi %d ; eta %d \n",absID, nSupMod, nModule, nIphi, nIeta); | |
295 | ||
296 | ||
297 | Int_t etainv = 47-ieta; | |
298 | //fHSMCellE[nSupMod]->Fill(etainv,phi,i); | |
299 | ||
300 | ||
301 | Int_t binx = etainv + 48*((nSupMod+1) % 2); | |
302 | Int_t biny = iphi + 24*nSupMod/2; | |
303 | if(nSupMod % 2 != 0) | |
304 | biny -= 12; | |
305 | ||
306 | Double_t cellE = cells->GetAmplitude(i); | |
307 | //fHistGlobalHmap0->Fill(binx,biny,cellE); | |
308 | ||
309 | Double_t nEtaGlobal, nPhiGlobal; | |
310 | Int_t nSupMod2; | |
311 | ||
312 | nSupMod2 = fGeom->GetSuperModuleNumber(absID); | |
313 | fGeom->EtaPhiFromIndex(absID,nEtaGlobal,nPhiGlobal); | |
314 | ||
315 | if(bVerbose){ | |
316 | printf("SM number: %d ; global eta: \t%6.4f ; global phi: \t%6.4f \n",nSupMod2,nEtaGlobal,nPhiGlobal); | |
317 | } | |
318 | ||
319 | ||
320 | //Double_t nPhi; | |
321 | //if(nPhiGlobal < 0){ | |
322 | // nPhi = TMath::Pi()+TMath::Abs(TMath::Pi()+nPhiGlobal); | |
323 | //} | |
324 | //else{ | |
325 | //nPhi = nPhiGlobal; | |
326 | //} | |
327 | fHistGlobalHmap0->Fill(binx,biny,cellE); | |
328 | fHistGlobalHmap->Fill(nEtaGlobal,nPhiGlobal,cellE); | |
329 | ||
330 | } | |
331 | ||
332 | if(bVerbose){ | |
333 | printf("\nThere are %d cells in this event\n", ncells); | |
334 | printf("There are %d tracks in this event\n", fESD->GetNumberOfTracks()); | |
335 | } | |
336 | ||
337 | ||
338 | //<------------ cell loop <----------------- | |
339 | ||
340 | //------------> track loop -----------------> | |
341 | for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) { | |
342 | AliESDtrack* track = fESD->GetTrack(iTracks); | |
343 | if (!track) { | |
344 | printf("ERROR: Could not receive track %d\n", iTracks); | |
345 | continue; | |
346 | } | |
347 | ||
348 | ||
349 | } | |
350 | //<------------ track loop <----------------- | |
351 | ||
352 | //------------> cluster loop -----------------> | |
353 | emclus = 0; | |
354 | Int_t nclus = clusters->GetEntries(); | |
355 | //printf("%d clusters\n",nclus); | |
356 | //cluster loop 1 for counting EMCal cluster | |
357 | for(int i = 0; i < nclus; i++){ | |
358 | AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i)); | |
359 | if (!clus){ | |
360 | //printf("Cluster empty!\n"); | |
361 | continue; | |
362 | } | |
363 | if (!clus->IsEMCAL()){ | |
364 | //printf("Cluster %d is not EMCAL!\n",i); | |
365 | // cout << "at " << i << " cluster not EMCAL" << endl; | |
366 | continue; | |
367 | } | |
368 | ||
369 | emclus++; | |
370 | } | |
371 | ||
372 | //<------------ cluster loop <----------------- | |
373 | ||
374 | ||
375 | //<============ DATA LOOP <================= | |
376 | ||
377 | ||
378 | //============> MC LOOP =================> | |
379 | if(bMC){ | |
380 | ||
381 | //------------> cell loop -----------------> | |
382 | //<------------ cell loop <----------------- | |
383 | ||
384 | //------------> track loop -----------------> | |
385 | //<------------ track loop <----------------- | |
386 | ||
387 | //------------> cluster loop -----------------> | |
388 | for(int i = 0; i < nclus; i++){ | |
389 | AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i)); | |
390 | if (!clus){ | |
391 | //printf("Cluster empty!\n"); | |
392 | continue; | |
393 | } | |
394 | if (!clus->IsEMCAL()){ | |
395 | //printf("Cluster %d is not EMCAL!\n",i); | |
396 | // cout << "at " << i << " cluster not EMCAL" << endl; | |
397 | continue; | |
398 | } | |
399 | //MCTruth | |
400 | int clslabel = -1; | |
401 | ||
402 | // all contributors | |
403 | Int_t* mcarr = clus->GetLabels(); | |
404 | // number of contributors | |
405 | Int_t nl = clus->GetNLabels(); | |
406 | ||
407 | ||
408 | ||
409 | ||
410 | for(int ip=0;ip<nl;ip++){ | |
411 | Int_t entry = mcarr[ip]; | |
412 | AliMCParticle *mcPart = static_cast<AliMCParticle*>(fMC->GetTrack(entry)); | |
413 | if(!mcPart){ | |
414 | continue; | |
415 | } | |
416 | contPID = mcPart->PdgCode(); | |
417 | if(bVerbose){ | |
418 | printf("%d \t %d \n",clslabel,mcPart->PdgCode()); | |
419 | } | |
420 | prodrad = mcPart->Particle()->R(); | |
421 | Int_t indexMoth = mcPart->GetMother(); | |
422 | if(indexMoth >= 0) | |
423 | { | |
424 | AliMCParticle* moth = static_cast<AliMCParticle*>(fMC->GetTrack(indexMoth)); | |
425 | // //AliMCParticle* moth = (AliMCParticle*)stack->At(indexMoth); | |
426 | mothPID = moth->PdgCode(); | |
427 | } | |
428 | ||
429 | ||
430 | ||
431 | ||
432 | ||
433 | //printf("%d \t %d \n",clslabel,nl); | |
434 | ||
435 | } | |
436 | fTreeCluster->Fill(); | |
437 | ||
438 | } | |
439 | ||
440 | //<------------ cluster loop <----------------- | |
441 | ||
442 | } | |
443 | //<============ MC LOOP <================= | |
444 | ||
445 | ||
446 | ||
447 | ||
448 | ||
449 | ||
450 | ||
451 | //cluster loop 2 for analysis | |
452 | //\cluster loop | |
453 | ||
454 | ||
455 | fTreeEvent->Fill(); | |
456 | ||
457 | ||
458 | PostData(1, fOutputList); | |
459 | } | |
460 | ||
461 | //________________________________________________________________________ | |
462 | void AliAnalysisTaskEMCALIsolation::Terminate(Option_t *) | |
463 | { | |
464 | // Draw result to the screen | |
465 | // Called once at the end of the query | |
466 | ||
467 | fOutputList = dynamic_cast<TList*> (GetOutputData(1)); | |
468 | if (!fOutputList) { | |
469 | printf("ERROR: Output list not available\n"); | |
470 | return; | |
471 | } | |
472 | ||
473 | fHistGlobalHmap = dynamic_cast<TH2F*> (fOutputList->At(2)); | |
474 | if (!fHistGlobalHmap) { | |
475 | printf("ERROR: fHistGlobalHmap not available\n"); | |
476 | return; | |
477 | } | |
478 | ||
479 | fHistGlobalHmap->DrawCopy("colz"); | |
480 | ||
481 | } | |
482 | ||
483 | ||
484 | const char * AliAnalysisTaskEMCALIsolation::GetParticleName(Int_t pdg) | |
485 | { | |
486 | TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg); | |
487 | if(p1) return p1->GetName(); | |
488 | return Form("%d", pdg); | |
489 | } |