e026a16039ec21c649bd60e616f9209f0a022ad3
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALIsoPhoton.cxx
1 // $Id$
2
3 #include "AliAnalysisTaskEMCALIsoPhoton.h"
4
5 #include <TCanvas.h>
6 #include <TChain.h>
7 #include <TFile.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <THnSparse.h>
12 #include <TLorentzVector.h>
13 #include <TList.h>
14
15 #include "AliAnalysisManager.h"
16 #include "AliAnalysisTask.h"
17 #include "AliEMCALGeometry.h"
18 #include "AliEMCALRecoUtils.h"
19 #include "AliESDCaloCells.h"
20 #include "AliESDCaloCluster.h"
21 #include "AliESDEvent.h"
22 #include "AliESDHeader.h"
23 #include "AliESDInputHandler.h"
24 #include "AliESDUtils.h"
25 #include "AliESDtrack.h"
26 #include "AliESDtrackCuts.h"
27 #include "AliMCEvent.h"
28 #include "AliMCEventHandler.h"
29 #include "AliStack.h"
30 #include "AliV0vertexer.h"
31 #include "AliVCluster.h"
32
33 ClassImp(AliAnalysisTaskEMCALIsoPhoton)
34
35 //________________________________________________________________________
36 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : 
37   AliAnalysisTaskSE(), 
38   fCaloClusters(0),
39   fSelPrimTracks(0),
40   fTracks(0),
41   fEMCalCells(0),
42   fPrTrCuts(0),
43   fGeom(0x0),
44   fGeoName("EMCAL_COMPLETEV1"),
45   fPeriod("LHC11c"),
46   fTrigBit("kEMC7"),
47   fIsTrain(0),
48   fIsMc(0),
49   fDebug(0),
50   fPathStrOpt("/"),
51   fExoticCut(0.97),
52   fIsoConeR(0.4),
53   fNDimensions(7),
54   fECut(3.),
55   fTrackMult(0),        
56   fMcIdFamily(""),
57   fNClusForDirPho(0),
58   fDirPhoPt(0),
59   fHigherPtCone(0),
60   fESD(0),
61   fMCEvent(0),
62   fStack(0),
63   fOutputList(0),
64   fEvtSel(0),
65   fNClusEt10(0),
66   fRecoPV(0),
67   fPVtxZ(0),
68   fTrMultDist(0),
69   fMCDirPhotonPtEtaPhi(0),
70   fDecayPhotonPtMC(0),
71   fCellAbsIdVsAmpl(0),       
72   fNClusHighClusE(0),
73   fHigherPtConeM02(0),
74   fClusEtMcPt(0),
75   fClusMcDetaDphi(0),
76   fNClusPerPho(0),
77   fMcPtInConeBG(0),
78   fMcPtInConeSBG(0),
79   fMcPtInConeBGnoUE(0),
80   fMcPtInConeSBGnoUE(0),
81   fMCDirPhotonPtEtaPhiNoClus(0),
82   fHnOutput(0)
83 {
84   // Default constructor.
85 }
86
87 //________________________________________________________________________
88 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : 
89   AliAnalysisTaskSE(name), 
90   fCaloClusters(0),
91   fSelPrimTracks(0),
92   fTracks(0),
93   fEMCalCells(0),
94   fPrTrCuts(0),
95   fGeom(0x0),
96   fGeoName("EMCAL_COMPLETEV1"),
97   fPeriod("LHC11c"),
98   fTrigBit("kEMC7"),
99   fIsTrain(0),
100   fIsMc(0),
101   fDebug(0),
102   fPathStrOpt("/"),
103   fExoticCut(0.97),
104   fIsoConeR(0.4),
105   fNDimensions(7),
106   fECut(3.),
107   fTrackMult(0),        
108   fMcIdFamily(""),
109   fNClusForDirPho(0),
110   fDirPhoPt(0),
111   fHigherPtCone(0),
112   fESD(0),
113   fMCEvent(0),
114   fStack(0),
115   fOutputList(0),
116   fEvtSel(0),
117   fNClusEt10(0),
118   fRecoPV(0),
119   fPVtxZ(0),            
120   fTrMultDist(0),
121   fMCDirPhotonPtEtaPhi(0),
122   fDecayPhotonPtMC(0),
123   fCellAbsIdVsAmpl(0),       
124   fNClusHighClusE(0),   
125   fHigherPtConeM02(0),
126   fClusEtMcPt(0),
127   fClusMcDetaDphi(0),
128   fNClusPerPho(0),
129   fMcPtInConeBG(0),
130   fMcPtInConeSBG(0),
131   fMcPtInConeBGnoUE(0),
132   fMcPtInConeSBGnoUE(0),
133   fMCDirPhotonPtEtaPhiNoClus(0),
134   fHnOutput(0)
135 {
136   // Constructor
137
138   // Define input and output slots here
139   // Input slot #0 works with a TChain
140   DefineInput(0, TChain::Class());
141   // Output slot #0 id reserved by the base class for AOD
142   // Output slot #1 writes into a TH1 container
143   DefineOutput(1, TList::Class());
144 }
145
146 //________________________________________________________________________
147 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
148 {
149   // Create histograms, called once.
150     
151   fCaloClusters = new TRefArray();
152   fSelPrimTracks = new TObjArray();
153
154   
155   fOutputList = new TList();
156   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging) 
157   
158   fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
159   
160   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
161   fOutputList->Add(fEvtSel);
162   
163   fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
164   fOutputList->Add(fNClusEt10);
165   
166   fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
167   fOutputList->Add(fRecoPV);
168
169   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
170   fOutputList->Add(fPVtxZ);
171
172   fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5);
173   fOutputList->Add(fTrMultDist);
174
175   fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20);
176   fMCDirPhotonPtEtaPhi->Sumw2();
177   fOutputList->Add(fMCDirPhotonPtEtaPhi);
178
179   fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100);
180   fDecayPhotonPtMC->Sumw2();
181   fOutputList->Add(fDecayPhotonPtMC);
182
183   fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,24*48*10,-0.5,24*48*10-0.5);
184   fOutputList->Add(fCellAbsIdVsAmpl);
185
186   fNClusHighClusE = new TH2F("hNClusHighClusE","total number of clusters vs. highest clus energy in the event;E (GeV);NClus",200,0,100,301,-0.5,300.5);
187   fOutputList->Add(fNClusHighClusE);
188
189   fHigherPtConeM02 = new TH2F("hHigherPtConeM02","#lambda_{0}^{2} vs. in-cone-p_{T}^{max};p_{T}^{max} (GeV/c, in the cone);#lambda_{0}^{2}",1000,0,100,400,0,4);
190   fOutputList->Add(fHigherPtConeM02);
191
192   fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100);
193   fOutputList->Add(fClusEtMcPt);
194
195   fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7);
196   fOutputList->Add(fClusMcDetaDphi);
197
198   fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5);
199   fOutputList->Add(fNClusPerPho);
200
201   fMcPtInConeBG = new TH2F("hMcPtInConeBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,500,0,50);
202   fOutputList->Add(fMcPtInConeBG);
203
204   fMcPtInConeSBG  = new TH2F("hMcPtInConeSBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,500,0,50);
205   fOutputList->Add(fMcPtInConeSBG);
206
207   fMcPtInConeBGnoUE = new TH2F("hMcPtInConeBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,500,0,50);
208   fOutputList->Add(fMcPtInConeBGnoUE);
209
210   fMcPtInConeSBGnoUE  = new TH2F("hMcPtInConeSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,500,0,50);
211   fOutputList->Add(fMcPtInConeSBGnoUE);
212
213
214   fMCDirPhotonPtEtaPhiNoClus = new TH3F("hMCDirPhotonPhiEtaNoClus","p_{T}, #eta and  #phi of prompt photons with no reco clusters;p_{T};#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20);
215   fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus);
216
217   Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000,  nAllIso=1000,  nCeIsoNoUE=1000,  nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=100, nPhoMcPt=101;
218   Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt};
219   fNDimensions = sizeof(bins)/sizeof(Int_t);
220   const Int_t ndims =   fNDimensions;
221   Double_t xmin[] = { 0.,   0.,  -10.,   -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-1.5};
222   Double_t xmax[] = { 100., 4., 190., 190., 190.,  190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5};
223   fHnOutput =  new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult,mc-p_{T}^{#gamma}", ndims, bins, xmin, xmax);
224   fOutputList->Add(fHnOutput);
225
226
227
228   PostData(1, fOutputList);
229 }
230
231 //________________________________________________________________________
232 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) 
233 {
234   // Main loop, called for each event.
235
236   // event trigger selection
237   Bool_t isSelected = 0;
238   if(fPeriod.Contains("11a")){
239     if(fTrigBit.Contains("kEMC"))
240       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
241     else
242       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
243   }
244   else{
245     if(fTrigBit.Contains("kEMC"))
246       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
247     else
248       if(fTrigBit.Contains("kMB"))
249         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
250       else
251         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
252   }
253   if(fIsMc)
254     isSelected = kTRUE;
255   if(fDebug)
256     printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
257   if(!isSelected )
258         return; 
259   if(fIsMc){
260     TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
261     TFile *file = (TFile*)tree->GetCurrentFile();
262     TString filename = file->GetName();
263     if(!filename.Contains(fPathStrOpt.Data()))
264       return;
265   }
266   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
267   if (!fESD) {
268     printf("ERROR: fESD not available\n");
269     return;
270   }
271   
272   fEvtSel->Fill(0);
273   if(fDebug)
274     printf("fESD is ok\n");
275   
276   AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
277   if(!pv)
278     return;
279   if(!pv->GetStatus())
280     fRecoPV->Fill(0);
281   else
282     fRecoPV->Fill(1);
283   fPVtxZ->Fill(pv->GetZ());
284   if(TMath::Abs(pv->GetZ())>15)
285     return;
286   if(fDebug)
287     printf("passed vertex cut\n");
288
289   fEvtSel->Fill(1);
290
291   fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
292   if(!fTracks){
293     AliError("Track array in event is NULL!");
294     return;
295   }
296   // Track loop to fill a pT spectrum
297   const Int_t Ntracks = fTracks->GetEntriesFast();
298   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
299     //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
300     //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
301     AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(iTracks));
302     if (!track)
303       continue;
304     if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
305       fSelPrimTracks->Add(track);
306       //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
307     }
308   }
309
310   if(!fIsTrain){
311     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
312       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
313         break;
314       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
315     }
316   }
317   AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
318   fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); 
319   fTrMultDist->Fill(fTrackMult);
320
321   fESD->GetEMCALClusters(fCaloClusters);
322   fEMCalCells = fESD->GetEMCALCells();
323   
324     
325
326   fMCEvent = MCEvent();
327   if(fMCEvent){
328     if(fDebug)
329       std::cout<<"MCevent exists\n";
330     fStack = (AliStack*)fMCEvent->Stack();
331   }
332   else{
333     if(fDebug)
334       std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
335   }
336   FollowGamma();
337   if(fDebug)
338     printf("passed calling of FollowGamma\n");
339   FillClusHists(); 
340   if(fDebug)
341     printf("passed calling of FillClusHists\n");
342   FillMcHists();
343   if(fDebug)
344     printf("passed calling of FillMcHists\n");
345
346   fCaloClusters->Clear();
347   fSelPrimTracks->Clear();
348   fNClusForDirPho = 0;
349
350   PostData(1, fOutputList);
351 }      
352
353 //________________________________________________________________________
354 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
355 {
356   // Fill cluster histograms.
357
358   if(!fCaloClusters)
359     return;
360   const Int_t nclus = fCaloClusters->GetEntries();
361   if(nclus==0)
362     return;
363   if(fDebug)
364     printf("Inside FillClusHists and there are %d clusters\n",nclus);
365   Double_t maxE = 0;
366   Int_t nclus10 = 0;
367   Double_t ptmc=-1;
368   for(Int_t ic=0;ic<nclus;ic++){
369     maxE=0;
370     AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
371     if(!c)
372       continue;
373     if(!c->IsEMCAL())
374       continue;
375     if(c->E()<fECut)
376       continue;
377     Short_t id;
378     Double_t Emax = GetMaxCellEnergy( c, id);
379     Double_t Ecross = GetCrossEnergy( c, id);
380     if((1-Ecross/Emax)>fExoticCut)
381       continue;
382     Float_t clsPos[3] = {0,0,0};
383     c->GetPosition(clsPos);
384     TVector3 clsVec(clsPos);
385     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
386     if(fDebug)
387       printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E());
388     if(Et>10)
389       nclus10++;
390     Float_t ceiso, cephiband, cecore;
391     Float_t triso, trphiband, trcore;
392     Float_t alliso, allphiband, allcore;
393     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
394     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
395     GetCeIso(clsVec, ceiso, cephiband, cecore);
396     GetTrIso(clsVec, triso, trphiband, trcore);
397     Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
398     if(Et>10 && Et<15 && dr>0.025){
399       fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
400       if(fDebug)
401         printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
402     }
403     alliso = ceiso + triso;
404     allphiband = cephiband + trphiband;
405     allcore = cecore + trcore;
406     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
407     Float_t trisoue =  trphiband/phibandArea*netConeArea;
408     Float_t allisoue =  allphiband/phibandArea*netConeArea;
409     Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR); 
410     if(fDebug && Et>10)
411       printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et);
412     if(c->GetM02()>0.5 && c->GetM02()<2.0){
413       fMcPtInConeBG->Fill(alliso-Et-allisoue, mcptsum);
414       fMcPtInConeBGnoUE->Fill(alliso-Et, mcptsum);
415     }
416     if(c->GetM02()>0.1 && c->GetM02()<0.3){
417       fMcPtInConeSBG->Fill(alliso-Et-allisoue, mcptsum);
418       fMcPtInConeSBGnoUE->Fill(alliso-Et, mcptsum);
419     }
420     const Int_t ndims =   fNDimensions;
421     Double_t outputValues[ndims];
422     ptmc = GetClusSource(c);
423     outputValues[0] = Et;
424     outputValues[1] = c->GetM02();
425     outputValues[2] = ceiso-Et/*cecore*/-ceisoue;
426     outputValues[3] = triso-trisoue;
427     outputValues[4] = alliso-Et/*cecore*/-allisoue;
428     outputValues[5] = ceiso-Et;
429     outputValues[6] = alliso-Et;
430     outputValues[7] = c->GetTrackDx();
431     outputValues[8] = c->GetTrackDz();
432     outputValues[9] = clsVec.Eta();
433     outputValues[10] = clsVec.Phi();
434     outputValues[11] = fEMCalCells->GetCellTime(id);
435     outputValues[12] = fTrackMult;
436     outputValues[13] = ptmc;
437     fHnOutput->Fill(outputValues);
438     if(c->E()>maxE)
439       maxE = c->E();
440   }
441   fNClusHighClusE->Fill(maxE,nclus);
442   fNClusEt10->Fill(nclus10);
443   fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
444
445
446 //________________________________________________________________________
447 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
448 {
449   // Get cell isolation.
450
451   if(!fEMCalCells)
452     return;
453   const Int_t ncells = fEMCalCells->GetNumberOfCells();
454   Float_t totiso=0;
455   Float_t totband=0;
456   Float_t totcore=0;
457   Float_t etacl = vec.Eta();
458   Float_t phicl = vec.Phi();
459   Float_t thetacl = vec.Theta();
460   if(phicl<0)
461     phicl+=TMath::TwoPi();
462   Int_t absid = -1;
463   Float_t eta=-1, phi=-1;  
464   for(int icell=0;icell<ncells;icell++){
465     absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
466     Float_t celltime = fEMCalCells->GetCellTime(absid);
467     if(celltime>2e-8 || celltime<-2e-8)
468       continue;
469     if(!fGeom)
470       return;
471     fGeom->EtaPhiFromIndex(absid,eta,phi);
472     Float_t dphi = TMath::Abs(phi-phicl);
473     Float_t deta = TMath::Abs(eta-etacl);
474     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
475     Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
476     if(R<fIsoConeR){
477       totiso += etcell;
478       if(R<0.04)
479         totcore += etcell;
480     }
481     else{
482       if(dphi>fIsoConeR)
483         continue;
484       if(deta<fIsoConeR)
485         continue;
486       totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
487     }
488   }
489   iso = totiso;
490   phiband = totband;
491   core = totcore;
492
493
494 //________________________________________________________________________
495 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
496 {
497   // Get track isolation.
498
499   if(!fSelPrimTracks)
500     return;
501   fHigherPtCone = 0;
502   const Int_t ntracks = fSelPrimTracks->GetEntries();
503   Double_t totiso=0;
504   Double_t totband=0;
505   Double_t totcore=0;
506   Double_t etacl = vec.Eta();
507   Double_t phicl = vec.Phi();
508   if(phicl<0)
509     phicl+=TMath::TwoPi();
510   for(int itrack=0;itrack<ntracks;itrack++){
511     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
512     if(!track)
513       continue;
514     Double_t dphi = TMath::Abs(track->Phi()-phicl);
515     Double_t deta = TMath::Abs(track->Eta()-etacl);
516     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
517     Double_t pt = track->Pt();
518     if(pt>fHigherPtCone)
519       fHigherPtCone = pt;
520     if(R<fIsoConeR){
521       totiso += track->Pt();
522       if(R<0.04)
523         totcore += pt;
524     }
525     else{
526       if(dphi>fIsoConeR)
527         continue;
528       if(deta<fIsoConeR)
529         continue;
530       totband += track->Pt();
531     }
532   }
533   iso = totiso;
534   phiband = totband;
535   core = totcore;
536
537
538 //________________________________________________________________________
539 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
540 {
541   // Calculate the energy of cross cells around the leading cell.
542
543   AliVCaloCells *cells = 0;
544   cells = fESD->GetEMCALCells();
545   if (!cells)
546     return 0;
547
548   if (!fGeom)
549     return 0;
550
551   Int_t iSupMod = -1;
552   Int_t iTower  = -1;
553   Int_t iIphi   = -1;
554   Int_t iIeta   = -1;
555   Int_t iphi    = -1;
556   Int_t ieta    = -1;
557   Int_t iphis   = -1;
558   Int_t ietas   = -1;
559
560   Double_t crossEnergy = 0;
561
562   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
563   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
564
565   Int_t ncells = cluster->GetNCells();
566   for (Int_t i=0; i<ncells; i++) {
567     Int_t cellAbsId = cluster->GetCellAbsId(i);
568     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
569     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
570     Int_t aphidiff = TMath::Abs(iphi-iphis);
571     if (aphidiff>1)
572       continue;
573     Int_t aetadiff = TMath::Abs(ieta-ietas);
574     if (aetadiff>1)
575       continue;
576     if ( (aphidiff==1 && aetadiff==0) ||
577         (aphidiff==0 && aetadiff==1) ) {
578       crossEnergy += cells->GetCellAmplitude(cellAbsId);
579     }
580   }
581
582   return crossEnergy;
583 }
584
585 //________________________________________________________________________
586 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
587 {
588   // Get maximum energy of attached cell.
589
590   id = -1;
591
592   AliVCaloCells *cells = 0;
593   cells = fESD->GetEMCALCells();
594   if (!cells)
595     return 0;
596
597   Double_t maxe = 0;
598   Int_t ncells = cluster->GetNCells();
599   for (Int_t i=0; i<ncells; i++) {
600     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
601     if (e>maxe) {
602       maxe = e;
603       id   = cluster->GetCellAbsId(i);
604     }
605   }
606   return maxe;
607 }
608
609 //________________________________________________________________________
610 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
611 {
612   if(!fStack)
613     return;
614   Int_t nTracks = fStack->GetNtrack();
615   if(fDebug)
616     printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
617   for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
618     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));  
619     if(!mcp)
620       continue;  
621     Int_t pdg = mcp->GetPdgCode();
622     if(pdg!=22)
623       continue;
624     if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
625       continue;
626     Int_t imom = mcp->GetMother(0);
627     if(imom<0 || imom>nTracks)
628       continue;
629     TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));  
630     if(!mcmom)
631       continue;
632     Int_t pdgMom = mcmom->GetPdgCode();
633     if((imom==6 || imom==7) && pdgMom==22) {
634       fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
635       if(fNClusForDirPho==0)
636         fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
637       if(fDebug){
638         printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatusCode());
639         printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatusCode());
640       }
641     }
642     else{
643       if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
644         fDecayPhotonPtMC->Fill(mcp->Pt());
645     }
646   }
647 }
648 //________________________________________________________________________
649 Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
650 {
651   if(!c)
652     return -1;
653   if(!fStack)
654     return -1;
655   Int_t nmcp = fStack->GetNtrack();
656   Int_t clabel = c->GetLabel();
657   if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
658     printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
659   if(!fMcIdFamily.Contains(Form("%d",clabel)))
660     return -1;
661   fNClusForDirPho++;
662   TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
663   Int_t partonpos = partonposstr.Atoi();
664   if(fDebug)
665     printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
666   if(clabel<0 || clabel>nmcp)
667     return 0;
668   Float_t clsPos[3] = {0,0,0};
669   c->GetPosition(clsPos);
670   TVector3 clsVec(clsPos);
671   Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
672   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
673   if(!mcp)
674     return -1;
675   if(fDebug){
676     printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),clabel);
677   }
678   Int_t lab1 =  mcp->GetFirstDaughter();
679   if(lab1<0 || lab1>nmcp)
680     return -1;
681   TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
682   if(!mcd)
683     return -1;
684   if(fDebug)
685     printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->Energy(),mcd->GetPdgCode(),lab1);
686   if(mcd->GetPdgCode()==22){
687     fClusEtMcPt->Fill(mcd->Pt(), Et);
688     fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
689   }
690   else{
691     printf("Warning: daughter of photon parton is not a photon\n");
692     return -1;
693   }
694   return fDirPhoPt;
695 }
696 //________________________________________________________________________
697 void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
698 {
699   if(!fStack)
700     return;
701   Int_t selfid = 6;
702   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
703   if(!mcp)
704     return;
705   if(mcp->GetPdgCode()!=22){
706     mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
707     if(!mcp)
708       return;
709   }  
710   Int_t daug0f =  mcp->GetFirstDaughter();
711   Int_t daug0l =  mcp->GetLastDaughter();
712   Int_t nd0 = daug0l - daug0f;
713   if(fDebug)
714     printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),nd0+1);
715   fMcIdFamily = Form("%d,",selfid);
716   GetDaughtersInfo(daug0f,daug0l, selfid,"");
717   if(fDebug){
718     printf("\t---- end of the prompt  gamma's family tree ----\n\n");
719     printf("Family id string = %s,\n\n",fMcIdFamily.Data());
720   }
721   TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
722   if(!md)
723     return;
724   fDirPhoPt = md->Pt();
725 }
726 //________________________________________________________________________
727 void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
728 {
729   int nmcp = fStack->GetNtrack();
730   if(firstd<0 || firstd>nmcp)
731     return;
732   if(lastd<0 || lastd>nmcp)
733     return;
734   if(firstd>lastd){
735     printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
736     return;
737   }
738   TString indenter = Form("\t%s",inputind);
739   TParticle *dp = 0x0;
740   if(fDebug)
741     printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
742   for(int id=firstd; id<lastd+1; id++){
743     dp =  static_cast<TParticle*>(fStack->Particle(id));
744     if(!dp)
745       continue;
746     Int_t dfd = dp->GetFirstDaughter(); 
747     Int_t dld = dp->GetLastDaughter();
748     Int_t dnd =  dld - dfd + 1;
749     Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
750     if(fDebug)
751       printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->Energy(), vr, dp->GetPdgCode(), dnd, dfd, dld);
752     fMcIdFamily += Form("%d,",id);
753     GetDaughtersInfo(dfd,dld,id,indenter.Data());
754   }
755 }
756
757 //________________________________________________________________________
758 Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
759 {
760   if(!fStack)
761     return 0;
762   if(fDebug)
763     printf("Inside GetMcPtSumInCone!!\n");
764   Int_t nTracks = fStack->GetNtrack();
765   Float_t ptsum = 0;
766   for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
767     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));  
768     if(!mcp)
769       continue;  
770     if(mcp->Rho()>2.5)
771       continue;
772     else {
773       if(fDebug)
774         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
775     }
776     Float_t dphi = mcp->Phi() - phiclus;
777     Float_t deta = mcp->Eta() - etaclus;
778     if(fDebug)
779       printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
780     if(deta>10)
781       continue;
782     Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
783     if(dR>R)
784       continue;
785     ptsum += mcp->Pt();
786   }
787   return ptsum;
788 }
789 //________________________________________________________________________
790 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *) 
791 {
792   // Called once at the end of the query.
793 }