8 #include "AliAnalysisTask.h"
9 #include "AliAnalysisManager.h"
11 #include "AliESDEvent.h"
12 #include "AliESDInputHandler.h"
13 #include "AliAODEvent.h"
14 #include "AliAODInputHandler.h"
16 #include "AliVEvent.h"
17 #include "AliVCaloCells.h"
18 #include "AliVCluster.h"
20 #include "AliEMCALGeometry.h"
21 #include "AliGeomManager.h"
23 #include "AliAnalysisTaskEMCALIsolation.h"
26 #include "AliMCEventHandler.h"
27 #include "AliMCEvent.h"
28 #include "AliMCParticle.h"
30 #include "TDatabasePDG.h"
31 #include "TParticlePDG.h"
33 // Task for isolating gammas with EMCAL
34 // Author: Marco Marquard
36 ClassImp(AliAnalysisTaskEMCALIsolation)
38 //________________________________________________________________________
40 AliAnalysisTaskEMCALIsolation::AliAnalysisTaskEMCALIsolation()
41 : AliAnalysisTaskSE(),
59 fGeoName("EMCAL_COMPLETE"),
66 DefineInput(0, TChain::Class());
67 DefineOutput(1, TList::Class());
70 //________________________________________________________________________
71 AliAnalysisTaskEMCALIsolation::AliAnalysisTaskEMCALIsolation(const char *name)
72 : AliAnalysisTaskSE(name),
90 fGeoName("EMCAL_COMPLETE"),
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());
106 //________________________________________________________________________
107 AliAnalysisTaskEMCALIsolation::~AliAnalysisTaskEMCALIsolation()
114 delete fHistGlobalHmap;
115 delete fHistGlobalHmap0;
120 //________________________________________________________________________
121 void AliAnalysisTaskEMCALIsolation::UserCreateOutputObjects()
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};
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};
129 fOutputList = new TList();
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);
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);
141 fHist = new TH1F("","",,,,);
142 fHist->GetXaxis()->SetTitle("");
143 fHist->GetYaxis()->SetTitle("");
144 fHist->SetMarkerStyle(kFullCircle);
145 fOutputList->Add(fHist);
148 fTreeEvent = new TTree("eventsT","eventsT");
149 fTreeEvent->Branch("emclus", &emclus);
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);
157 fOutputList->Add(fHistGlobalHmap);
158 fOutputList->Add(fHistGlobalHmap0);
159 fOutputList->Add(fTreeEvent);
160 fOutputList->Add(fTreeCluster);
166 //________________________________________________________________________
167 void AliAnalysisTaskEMCALIsolation::UserExec(Option_t *)
170 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
173 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
176 Printf("ERROR: Could not get ESDInputHandler");
179 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
182 Printf("ERROR: Could not get AODInputHandler");
186 fESD = esdH->GetEvent();
187 am->LoadBranch("AliESDRun.");
188 am->LoadBranch("AliESDHeader.");
191 fAOD = aodH->GetEvent();
194 AliFatal("Neither ESD nor AOD event found");
198 AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
200 Printf("ERROR: Could not retrieve MC event handler");
203 fMC = mcH->MCEvent();
205 Printf("ERROR: Could not retrieve MC event");
208 // stack to get truth infos of particles:
209 fStack = fMC->Stack();
211 Printf("ERROR: Could not retrieve stack from MC event");
214 //Printf(" MC particles: %d", fMC->GetNumberOfTracks());
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);
223 AliFatal("EMCAL(fGeom) geometry is not initialized");
227 AliVCaloCells *cells = 0;
230 am->LoadBranch("EMCALCells.");
231 fESDCells = fESD->GetEMCALCells();
235 am->LoadBranch("emcalCells");
236 fAODCells = fAOD->GetEMCALCells();
244 AliFatal("No cells loaded");
249 am->LoadBranch("CaloClusters");
250 TList *l = fESD->GetList();
251 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
254 am->LoadBranch("caloClusters");
255 TList *l = fAOD->GetList();
256 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
259 AliFatal("Impossible to not have either pointer to ESD or AOD event");
262 TObjArray *clusters = fEsdClusters;
264 clusters = fAodClusters;
266 Printf("ERROR: no clusters node in event!");
270 //========================================= start of analysis ============================================
273 //============> DATA LOOP =================>
275 //------------> cell loop ----------------->
276 Int_t ncells = cells->GetNumberOfCells();
277 //for (Int_t i = 0; i<15000; ++i){
279 for (Int_t i = 0; i<ncells; ++i){
280 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
282 //Check if this absId exists
283 if(!(fGeom->CheckAbsCellId(absID))){
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);
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);
294 printf("ID: %d ; supermodule: %d ; module %d ; phi %d ; eta %d \n",absID, nSupMod, nModule, nIphi, nIeta);
297 Int_t etainv = 47-ieta;
298 //fHSMCellE[nSupMod]->Fill(etainv,phi,i);
301 Int_t binx = etainv + 48*((nSupMod+1) % 2);
302 Int_t biny = iphi + 24*nSupMod/2;
306 Double_t cellE = cells->GetAmplitude(i);
307 //fHistGlobalHmap0->Fill(binx,biny,cellE);
309 Double_t nEtaGlobal, nPhiGlobal;
312 nSupMod2 = fGeom->GetSuperModuleNumber(absID);
313 fGeom->EtaPhiFromIndex(absID,nEtaGlobal,nPhiGlobal);
316 printf("SM number: %d ; global eta: \t%6.4f ; global phi: \t%6.4f \n",nSupMod2,nEtaGlobal,nPhiGlobal);
321 //if(nPhiGlobal < 0){
322 // nPhi = TMath::Pi()+TMath::Abs(TMath::Pi()+nPhiGlobal);
327 fHistGlobalHmap0->Fill(binx,biny,cellE);
328 fHistGlobalHmap->Fill(nEtaGlobal,nPhiGlobal,cellE);
333 printf("\nThere are %d cells in this event\n", ncells);
334 printf("There are %d tracks in this event\n", fESD->GetNumberOfTracks());
338 //<------------ cell loop <-----------------
340 //------------> track loop ----------------->
341 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
342 AliESDtrack* track = fESD->GetTrack(iTracks);
344 printf("ERROR: Could not receive track %d\n", iTracks);
350 //<------------ track loop <-----------------
352 //------------> cluster loop ----------------->
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));
360 //printf("Cluster empty!\n");
363 if (!clus->IsEMCAL()){
364 //printf("Cluster %d is not EMCAL!\n",i);
365 // cout << "at " << i << " cluster not EMCAL" << endl;
372 //<------------ cluster loop <-----------------
375 //<============ DATA LOOP <=================
378 //============> MC LOOP =================>
381 //------------> cell loop ----------------->
382 //<------------ cell loop <-----------------
384 //------------> track loop ----------------->
385 //<------------ track loop <-----------------
387 //------------> cluster loop ----------------->
388 for(int i = 0; i < nclus; i++){
389 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
391 //printf("Cluster empty!\n");
394 if (!clus->IsEMCAL()){
395 //printf("Cluster %d is not EMCAL!\n",i);
396 // cout << "at " << i << " cluster not EMCAL" << endl;
403 Int_t* mcarr = clus->GetLabels();
404 // number of contributors
405 Int_t nl = clus->GetNLabels();
410 for(int ip=0;ip<nl;ip++){
411 Int_t entry = mcarr[ip];
412 AliMCParticle *mcPart = static_cast<AliMCParticle*>(fMC->GetTrack(entry));
416 contPID = mcPart->PdgCode();
418 printf("%d \t %d \n",clslabel,mcPart->PdgCode());
420 prodrad = mcPart->Particle()->R();
421 Int_t indexMoth = mcPart->GetMother();
424 AliMCParticle* moth = static_cast<AliMCParticle*>(fMC->GetTrack(indexMoth));
425 // //AliMCParticle* moth = (AliMCParticle*)stack->At(indexMoth);
426 mothPID = moth->PdgCode();
433 //printf("%d \t %d \n",clslabel,nl);
436 fTreeCluster->Fill();
440 //<------------ cluster loop <-----------------
443 //<============ MC LOOP <=================
451 //cluster loop 2 for analysis
458 PostData(1, fOutputList);
461 //________________________________________________________________________
462 void AliAnalysisTaskEMCALIsolation::Terminate(Option_t *)
464 // Draw result to the screen
465 // Called once at the end of the query
467 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
469 printf("ERROR: Output list not available\n");
473 fHistGlobalHmap = dynamic_cast<TH2F*> (fOutputList->At(2));
474 if (!fHistGlobalHmap) {
475 printf("ERROR: fHistGlobalHmap not available\n");
479 fHistGlobalHmap->DrawCopy("colz");
484 const char * AliAnalysisTaskEMCALIsolation::GetParticleName(Int_t pdg)
486 TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg);
487 if(p1) return p1->GetName();
488 return Form("%d", pdg);