generatio of various plaes in on the fly
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALPhoton.cxx
1 // $Id$
2 //
3 //
4 //
5 //
6
7 #include "TChain.h"
8 #include "TTree.h"
9 #include "TH1F.h"
10 #include "TH2F.h"
11 #include "TCanvas.h"
12
13 #include "AliAnalysisTask.h"
14 #include "AliAnalysisManager.h"
15
16 #include "AliESDEvent.h"
17 #include "AliESDHeader.h"
18 #include "AliESDUtils.h"
19 #include "AliESDInputHandler.h"
20 #include "AliCentrality.h"
21 #include "AliESDpid.h"
22 #include "AliKFParticle.h"
23
24 #include "AliMCEventHandler.h"
25 #include "AliMCEvent.h"
26 #include "AliStack.h"
27 #include "TParticle.h"
28
29
30 #include "AliESDtrackCuts.h"
31 #include "AliESDv0.h"
32 #include "AliV0vertexer.h"
33 #include "AliESDCaloCluster.h"
34 #include "AliESDCaloCells.h"
35 #include "AliEMCALGeometry.h"
36 #include "AliEMCALRecoUtils.h"
37 #include "AliVCluster.h"
38 #include "AliAnalysisTaskEMCALClusterizeFast.h"
39 #include "TLorentzVector.h"
40
41
42 #include "AliAnalysisTaskEMCALPhoton.h"
43 #include "TFile.h"
44
45
46 using std::cout;
47 using std::endl;
48
49 ClassImp(AliAnalysisTaskEMCALPhoton)
50
51 //________________________________________________________________________
52 AliAnalysisTaskEMCALPhoton::AliAnalysisTaskEMCALPhoton() : 
53   AliAnalysisTaskSE(), 
54   fTrCuts(0),
55   fPrTrCuts(0),
56   fSelTracks(0),
57   fSelPrimTracks(0),
58   fPhotConvArray(0),
59   fMyClusts(0),
60   fMyAltClusts(0),
61   fMyCells(0),
62   fMyTracks(0),
63   fMyMcParts(0),
64   fHeader(0x0),
65   fCaloClusters(0),
66   fCaloClustersNew(0),
67   fEMCalCells(0),
68   fGeom(0x0),
69   fTimeResTOF(0),
70   fMipResponseTPC(0),
71   fGeoName("EMCAL_COMPLETEV1"),
72   fPeriod("LHC11d"),
73   fIsTrain(0),
74   fIsMC(0),
75   fDebug(0),
76   fIsGrid(0),
77   fClusThresh(2.0),
78   fClusterizer(0),
79   fCaloClustersName("EmcalClusterv2"),
80   fESD(0),
81   fMCEvent(0),
82   fStack(0x0),
83   fOutputList(0),
84   fTree(0),
85   fNV0sBefAndAftRerun(0),
86   fConversionVtxXY(0),
87   fInvMassV0(0),
88   fInvMassV0KF(0),
89   fInvMassV0SS(0),
90   fDedxPAll(0)
91 {
92   // Default constructor.
93 }
94
95 //________________________________________________________________________
96 AliAnalysisTaskEMCALPhoton::AliAnalysisTaskEMCALPhoton(const char *name) : 
97   AliAnalysisTaskSE(name), 
98   fTrCuts(0),
99   fPrTrCuts(0),
100   fSelTracks(0),
101   fSelPrimTracks(0),
102   fPhotConvArray(0),
103   fMyClusts(0),
104   fMyAltClusts(0),
105   fMyCells(0),
106   fMyTracks(0),
107   fMyMcParts(0),
108   fHeader(0),
109   fCaloClusters(0),
110   fCaloClustersNew(0),
111   fEMCalCells(0),
112   fGeom(0x0),
113   fTimeResTOF(0),
114   fMipResponseTPC(0),
115   fGeoName("EMCAL_COMPLETEV1"),
116   fPeriod("LHC11c"),
117   fIsTrain(0),
118   fIsMC(0),
119   fDebug(0),
120   fIsGrid(0),
121   fClusThresh(2.),
122   fClusterizer(0),
123   fCaloClustersName("EmcalClusterv2"),
124   fESD(0),
125   fMCEvent(0),
126   fStack(0x0),
127   fOutputList(0),
128   fTree(0),
129   fNV0sBefAndAftRerun(0),
130   fConversionVtxXY(0),
131   fInvMassV0(0),
132   fInvMassV0KF(0),
133   fInvMassV0SS(0),
134   fDedxPAll(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   DefineOutput(2, TTree::Class());
145 }
146
147 //________________________________________________________________________
148 void AliAnalysisTaskEMCALPhoton::UserCreateOutputObjects()
149 {
150   // Create histograms, called once.
151     
152   fSelTracks = new TObjArray();
153   
154   fSelPrimTracks = new TObjArray();
155   
156   if (TClass::GetClass("AliPhotonHeaderObj"))
157       TClass::GetClass("AliPhotonHeaderObj")->IgnoreTObjectStreamer();
158   fHeader = new AliPhotonHeaderObj;
159
160   if (TClass::GetClass("AliPhotonConvObj"))
161       TClass::GetClass("AliPhotonConvObj")->IgnoreTObjectStreamer();
162   fPhotConvArray = new TClonesArray("AliPhotonConvObj");
163   
164   if (TClass::GetClass("AliPhotonClusterObj"))
165       TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer();
166   fMyClusts = new TClonesArray("AliPhotonClusterObj");
167   
168   if (TClass::GetClass("AliPhotonCellObj"))
169       TClass::GetClass("AliPhotonCellObj")->IgnoreTObjectStreamer();
170   fMyCells = new TClonesArray("AliPhotonCellObj");
171   
172   if (TClass::GetClass("AliPhotonTrackObj"))
173       TClass::GetClass("AliPhotonTrackObj")->IgnoreTObjectStreamer();
174   fMyTracks = new TClonesArray("AliPhotonTrackObj");
175
176   if (fClusterizer || fIsTrain){
177     if(fClusterizer)
178       fCaloClustersName = fClusterizer->GetNewClusterArrayName();
179     else {
180       if(fPeriod.Contains("10h") || fPeriod.Contains("11h"))
181         fCaloClustersName = "EmcalClustersv1";
182       else
183         fCaloClustersName = "EmcalClustersv2";
184     }
185     if (TClass::GetClass("AliPhotonClusterObj"))
186         TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer();
187     fMyAltClusts = new TClonesArray("AliPhotonClusterObj");
188   }
189   cout<<fCaloClustersName<<endl;
190   if(fIsMC){
191     if (TClass::GetClass("AliPhotonMcPartObj"))
192         TClass::GetClass("AliPhotonMcPartObj")->IgnoreTObjectStreamer();
193     fMyMcParts = new TClonesArray("AliPhotonMcPartObj");
194   }
195  
196   fCaloClusters = new TRefArray();
197     
198   fOutputList = new TList();
199   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging) 
200
201   if( !fTree){
202     TFile *f = OpenFile(2); 
203     TDirectory::TContext context(f);
204     if (f) {
205       f->SetCompressionLevel(2);
206       fTree = new TTree("photon_ana_out", "StandaloneTree");
207       fTree->SetDirectory(f);
208       if (fIsTrain) {
209         fTree->SetAutoFlush(-2*1024*1024);
210         fTree->SetAutoSave(0);
211       } else {
212         fTree->SetAutoFlush(-32*1024*1024);
213         fTree->SetAutoSave(0);
214       }
215     
216     fTree->Branch("header", &fHeader, 16*1024, 99);    
217     fTree->Branch("conversions", &fPhotConvArray, 8*16*1024, 99);
218     fTree->Branch("clusters", &fMyClusts, 8*16*1024, 99);
219     if(fClusterizer || fIsTrain)
220       fTree->Branch(fCaloClustersName, &fMyAltClusts, 8*16*1024, 99);
221     fTree->Branch("cells", &fMyCells, 8*16*1024, 99);
222     fTree->Branch("IsoTracks", &fMyTracks, 8*16*1024, 99);
223     if(fIsMC)
224       fTree->Branch("mcparts", &fMyMcParts, 8*16*1024, 99);
225     }
226   }
227   //if(fIsGrid)fOutputList->Add(fTree);
228   fGeom = AliEMCALGeometry::GetInstance(fGeoName);
229   
230   
231   fNV0sBefAndAftRerun = new TH2F("hNV0sBefAndAftRerun","check if the number of v0s change with rerun;old v0 n;new v0 n",50,0.5,50.5,50,0.5,50.5);
232   fOutputList->Add(fNV0sBefAndAftRerun);
233   
234   fConversionVtxXY = new TH2F("hConversionVtxXY","x and y of conversion vertex candidates;x;y",1000,-100,100,1000,-100,100);
235   fOutputList->Add(fConversionVtxXY);
236   
237   fInvMassV0 = new TH1F("hInvMassV0","v0->GetEffMass();v0->GetEffMass();dN/dM",400,0,4);
238   fOutputList->Add(fInvMassV0);
239   
240   fInvMassV0KF = new TH1F("hInvMassV0KF","Inv. mass calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4);
241   fOutputList->Add(fInvMassV0KF);
242       
243   fInvMassV0SS = new TH1F("hInvMassV0SS","Inv. mass (same sign) calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4);
244   fOutputList->Add(fInvMassV0SS);
245       
246   fDedxPAll = new TH2F("hDedxPAll","dE/dx vs p (all selected tracks);p (GeV/c);dE/dx (a.u.)",400,0,40, 150,0,150);
247   fOutputList->Add(fDedxPAll);
248   
249   
250   PostData(1, fOutputList);
251   PostData(2, fTree);
252 }
253
254 //________________________________________________________________________
255 void AliAnalysisTaskEMCALPhoton::UserExec(Option_t *) 
256 {
257   // User exec, called once per event.
258
259   Bool_t isSelected = 0;
260   if(fPeriod.Contains("11")){
261     if(fPeriod.Contains("11a"))
262       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
263     if(fPeriod.Contains("11c") ||fPeriod.Contains("11d") )
264       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
265     if(fPeriod.Contains("11h") )
266       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA);
267
268   }
269   if(fIsMC){
270     isSelected = kTRUE;  
271   }
272
273
274   // Post output data.
275   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
276   if (!fESD) {
277     printf("ERROR: fESD not available, returning...\n");
278     return;
279   }
280   
281   AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
282   if(!pv) 
283     return;
284   if(TMath::Abs(pv->GetZ())>15)
285     return;
286
287   TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
288   TFile *inpfile = (TFile*)tree->GetCurrentFile();
289
290   // Track loop to fill a pT spectrum
291   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
292     AliESDtrack* track = fESD->GetTrack(iTracks);
293     if (!track)
294       continue;
295     
296     
297     if (fTrCuts && fTrCuts->IsSelected(track)){
298       fSelTracks->Add(track);
299       fDedxPAll->Fill(track->P(), track->GetTPCsignal());
300     }
301     if (fPrTrCuts && fPrTrCuts->IsSelected(track))
302       fSelPrimTracks->Add(track);
303   } //track loop 
304
305   fHeader->fInputFileName  = inpfile->GetName();
306   fHeader->fTrClassMask    = fESD->GetHeader()->GetTriggerMask();
307   fHeader->fTrCluster      = fESD->GetHeader()->GetTriggerCluster();
308   AliCentrality *cent = InputEvent()->GetCentrality();
309   fHeader->fV0Cent    = cent->GetCentralityPercentileUnchecked("V0M");
310   fHeader->fCl1Cent   = cent->GetCentralityPercentileUnchecked("CL1");
311   fHeader->fTrCent    = cent->GetCentralityPercentileUnchecked("TRK");
312   if(!fIsTrain){
313     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
314       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
315         break;
316       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
317     }
318   }
319   fESD->GetEMCALClusters(fCaloClusters);
320   fHeader->fNClus = fCaloClusters->GetEntries();
321   fEMCalCells = fESD->GetEMCALCells();
322   fHeader->fNCells = fEMCalCells->GetNumberOfCells();
323   AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
324   fHeader->fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); 
325
326   fMCEvent = MCEvent();
327   if(fMCEvent)
328     fStack = (AliStack*)fMCEvent->Stack();
329
330   
331   FindConversions();
332   FillMyCells();
333   FillMyClusters();
334   if(fCaloClustersNew)
335     FillMyAltClusters();
336   FillIsoTracks();
337   if(fIsMC)
338     GetMcParts();
339   
340   fTree->Fill();
341   fSelTracks->Clear();
342   fSelPrimTracks->Clear();
343   fPhotConvArray->Clear();
344   fMyClusts->Clear();
345   if(fMyAltClusts)
346     fMyAltClusts->Clear();
347   fMyCells->Clear();
348   fMyTracks->Clear();
349   if(fMyMcParts)
350     fMyMcParts->Clear();
351   fCaloClusters->Clear();
352   if(fCaloClustersNew)
353     fCaloClustersNew->Clear();
354   PostData(1, fOutputList);
355   PostData(2, fTree);
356 }      
357
358 //________________________________________________________________________
359 void AliAnalysisTaskEMCALPhoton::FindConversions() 
360 {
361   // Find conversion.
362
363   if(!fTrCuts)
364     return;
365   Int_t iconv = 0;
366   Int_t nV0Orig = fESD->GetNumberOfV0s();
367   Int_t ntracks = fESD->GetNumberOfTracks();
368   fESD->ResetV0s();
369   AliV0vertexer lV0vtxer;
370   lV0vtxer.Tracks2V0vertices(fESD);
371   Int_t nV0s = fESD->GetNumberOfV0s();
372   fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
373   for(Int_t iv0=0; iv0<nV0s; iv0++){
374     AliESDv0 * v0 = fESD->GetV0(iv0);
375     if(!v0)
376       continue;
377     Double_t vpos[3];
378     fInvMassV0->Fill(v0->GetEffMass());
379     v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
380     Int_t ipos = v0->GetPindex();
381     Int_t ineg = v0->GetNindex();
382     if(ipos<0 || ipos> ntracks)
383       continue;
384     if(ineg<0 || ineg> ntracks)
385       continue;
386     AliESDtrack *pos = static_cast<AliESDtrack*>(fESD->GetTrack(ipos));
387     AliESDtrack *neg = static_cast<AliESDtrack*>(fESD->GetTrack(ineg));
388     if(!pos || !neg)
389       continue;
390     if (!fTrCuts->IsSelected(pos) || !fTrCuts->IsSelected(neg))
391       continue;
392     /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
393       continue;*/
394     const AliExternalTrackParam * paramPos = v0->GetParamP() ;
395     const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
396     if(!paramPos || !paramNeg)
397       continue;
398     if(pos->GetSign() <0){//change tracks
399       pos=neg ;
400       neg=fESD->GetTrack(v0->GetPindex()) ;
401       paramPos=paramNeg ;
402       paramNeg=v0->GetParamP() ;
403     }
404     AliKFParticle negKF(*paramNeg,11);
405     AliKFParticle posKF(*paramPos,-11);
406     AliKFParticle photKF(negKF,posKF) ;
407    
408     if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0) 
409       continue ;
410     
411     if(pos->GetSign() == neg->GetSign()){ 
412       fInvMassV0SS->Fill(photKF.GetMass());
413       continue;
414     }
415     fInvMassV0KF->Fill(photKF.GetMass());
416     TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE()); 
417     if(photLV.M()>0.05 || photLV.M()<0)
418       continue;
419     fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
420     Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
421     if(convPhi<0)
422       convPhi+=TMath::TwoPi();
423     TVector3 vecpos(vpos);
424     Double_t v0Phi = 0;
425     if(vecpos.Perp()>0)
426       vecpos.Phi();
427     if(v0Phi<0)
428       v0Phi+=TMath::TwoPi();
429     AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
430     myconv->fPt = photLV.Pt();
431     myconv->fEta = photLV.Eta();
432     myconv->fPhi = convPhi;
433     myconv->fVR = vecpos.Perp();
434     if(vecpos.Perp()>0)
435       myconv->fVEta = vecpos.Eta();
436     myconv->fVPhi = v0Phi;
437     myconv->fMass = photLV.M();
438     myconv->fMcLabel = -3; //WARNING: include the correct labeling
439     //negative daughter
440    myconv->fNegPt      =  negKF.GetPt();
441    myconv->fNegEta     =  negKF.GetEta();
442    Double_t trackPhi   =  negKF.GetPhi();
443    if(trackPhi<0)
444      trackPhi+=TMath::TwoPi();
445    myconv->fNegPhi     =  trackPhi;
446    myconv->fNegDedx    =  neg->GetTPCsignal();
447    myconv->fNegMcLabel =  neg->GetLabel();
448     //negative daughter
449    myconv->fPosPt      =  posKF.GetPt();
450    myconv->fPosEta     =  posKF.GetEta();
451    trackPhi            =  posKF.GetPhi();
452    if(trackPhi<0)
453      trackPhi+=TMath::TwoPi();
454    myconv->fPosPhi     =  trackPhi;
455    myconv->fPosDedx    =  pos->GetTPCsignal();
456    myconv->fPosMcLabel =  pos->GetLabel();
457   }
458   return;
459 }
460
461 //________________________________________________________________________
462 void AliAnalysisTaskEMCALPhoton::FillMyCells() 
463 {
464   // Fill cells.
465
466   if (!fEMCalCells)
467     return;
468   Int_t ncells = fEMCalCells->GetNumberOfCells();
469   Int_t mcel = 0;
470   for(Int_t icell = 0; icell<ncells; icell++){
471     Int_t absID = TMath::Abs(fEMCalCells->GetCellNumber(icell));
472     AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
473     Float_t eta=-1, phi=-1;
474     if(!fGeom){
475       std::cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n";
476       return;
477     }
478     if(!fGeom)
479       return;
480     /*if(!fIsMC)*/fGeom->EtaPhiFromIndex(absID,eta,phi);
481     Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
482     mycell->fAbsID = absID;
483     mycell->fE = fEMCalCells->GetCellAmplitude(absID);
484     mycell->fEt = fEMCalCells->GetCellAmplitude(absID)*TMath::Sin(theta);
485     mycell->fEta = eta;
486     mycell->fPhi = phi;
487     mycell->fTime = fEMCalCells->GetCellTime(absID);
488   }
489 }
490
491 //________________________________________________________________________
492 void AliAnalysisTaskEMCALPhoton::FillMyClusters() 
493 {
494   // Fill clusters.
495
496   if (!fCaloClusters)
497     return;
498   Int_t nclus = fCaloClusters->GetEntries();
499   Int_t mcl = 0;
500   for(Int_t ic=0; ic < nclus; ic++){
501     AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
502     if(!clus)
503       continue;
504     if(!clus->IsEMCAL())
505       continue;
506     if(clus->E() < fClusThresh)
507       continue;
508     Float_t pos[3];
509     clus->GetPosition(pos);
510     TVector3 cpos(pos);
511     TString cellsAbsID;
512     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
513     myclus->fE       = clus->E();
514     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
515     myclus->fR       = cpos.Perp();
516     myclus->fEta     = cpos.Eta();
517     myclus->fPhi     = cpos.Phi();
518     if(cpos.Phi()<0){
519       myclus->fPhi+=TMath::TwoPi();
520     }
521     myclus->fN       = clus->GetNCells();
522     Short_t  id = -1;
523     myclus->fEmax    = GetMaxCellEnergy( clus, id); 
524     myclus->fIdmax   = id;
525     myclus->fTmax    = fEMCalCells->GetCellTime(id);
526     myclus->fEcross  = GetCrossEnergy( clus, id);
527     myclus->fDisp    = clus->GetDispersion();
528     myclus->fM20     = clus->GetM20();
529     myclus->fM02     = clus->GetM02();
530     myclus->fTrDEta  = clus->GetTrackDz();
531     myclus->fTrDPhi  = clus->GetTrackDx();
532     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
533     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
534     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
535     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
536     myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
537     myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
538     myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
539     myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
540     for(Int_t icell=0;icell<myclus->fN;icell++){
541       Int_t absID = clus->GetCellAbsId(icell);
542       cellsAbsID.Append(Form("%d",absID));
543       cellsAbsID.Append(";");
544     }
545     myclus->fCellsAbsId = cellsAbsID;
546     myclus->fMcLabel = clus->GetLabel(); 
547     Int_t matchIndex = clus->GetTrackMatchedIndex();
548     if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
549       myclus->fTrEp = -1;
550       continue;
551     }
552     AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
553     if(!track){
554       myclus->fTrEp = -1;
555       continue;
556     }
557     if(!fPrTrCuts){
558       myclus->fTrEp = -1;
559       continue;
560     }
561     if(!fPrTrCuts->IsSelected(track)){
562       myclus->fTrEp = -1;
563       continue;
564     }
565     myclus->fTrEp = clus->E()/track->P();
566   }
567   
568 }
569 //________________________________________________________________________
570 void AliAnalysisTaskEMCALPhoton::FillMyAltClusters() 
571 {
572   // Fill clusters.
573
574   if(!fCaloClustersNew)
575     return;
576   Int_t nclus = fCaloClustersNew->GetEntries();
577   Int_t mcl = 0;
578   for(Int_t ic=0; ic < nclus; ic++){
579     AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClustersNew->At(ic));
580     if(!clus)
581       continue;
582     if(!clus->IsEMCAL())
583       continue;
584     if(clus->E() < fClusThresh)
585       continue;
586     Float_t pos[3];
587     clus->GetPosition(pos);
588     TVector3 cpos(pos);
589     TString cellsAbsID;
590     AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
591     myclus->fE       = clus->E();
592     myclus->fEt      = clus->E()*TMath::Sin(cpos.Theta());
593     myclus->fR       = cpos.Perp();
594     myclus->fEta     = cpos.Eta();
595     myclus->fPhi     = cpos.Phi();
596     if(cpos.Phi()<0){
597       myclus->fPhi+=TMath::TwoPi();
598     }
599     myclus->fN       = clus->GetNCells();
600     myclus->fDisp    = clus->GetDispersion();
601     myclus->fM20     = clus->GetM20();
602     myclus->fM02     = clus->GetM02();
603     myclus->fTrDEta  = clus->GetTrackDz();
604     myclus->fTrDPhi  = clus->GetTrackDx();
605     myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
606     myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
607     myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
608     myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
609     for(Int_t icell=0;icell<myclus->fN;icell++){
610       Int_t absID = clus->GetCellAbsId(icell);
611       cellsAbsID.Append(Form("%d",absID));
612       cellsAbsID.Append(";");
613     }
614     myclus->fCellsAbsId = cellsAbsID;
615     myclus->fMcLabel = clus->GetLabel(); 
616     Int_t matchIndex = clus->GetTrackMatchedIndex();
617     if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
618       myclus->fTrEp = -1;
619       continue;
620     }
621     AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
622     if(!track){
623       myclus->fTrEp = -1;
624       continue;
625     }
626     if(!fPrTrCuts){
627       myclus->fTrEp = -1;
628       continue;
629     }
630     if(!fPrTrCuts->IsSelected(track)){
631       myclus->fTrEp = -1;
632       continue;
633     }
634     myclus->fTrEp = clus->E()/track->P();
635   }
636   
637 }
638 //________________________________________________________________________
639 void  AliAnalysisTaskEMCALPhoton::FillIsoTracks()
640 {
641   // Fill high pt tracks.
642
643   if(!fSelPrimTracks)
644     return;
645   Int_t ntracks = fSelPrimTracks->GetEntries();
646   Int_t imtr = 0;
647   for(Int_t it=0;it<ntracks; it++){
648     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(it));
649     if(!track)
650       continue;
651     /*if(track->Pt()<3)
652       continue;*/
653     AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
654     if(track->Phi()<1.0 || track->Phi()>3.55)
655       continue;
656     if(TMath::Abs(track->Eta())>1.1)
657       continue;
658     mtr->fPt = track->Pt();
659     mtr->fEta = track->Eta();
660     mtr->fPhi = track->Phi();
661     mtr->fCharge = track->Charge();
662     mtr->fDedx = track->GetTPCsignal();
663     mtr->fMcLabel = track->GetLabel();
664   }
665 }
666
667 //________________________________________________________________________
668 void  AliAnalysisTaskEMCALPhoton::GetMcParts()
669 {
670   // Get MC particles.
671
672   if (!fStack)
673     return;
674
675   const AliVVertex *evtVtx = fMCEvent->GetPrimaryVertex();
676   if (!evtVtx)
677     return;
678   Int_t ipart = 0;
679   Int_t nTracks = fStack->GetNtrack();
680   for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
681     TParticle *mcP = static_cast<TParticle*>(fStack->Particle(iTrack));
682     if (!mcP)
683       continue;
684     // primary particle
685     Double_t dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) + 
686                               (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
687                               (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
688     if(dR > 0.5)
689       continue;
690     
691     // kinematic cuts
692     Double_t pt = mcP->Pt() ;
693     if (pt<0.5)
694       continue;
695     Double_t eta = mcP->Eta();
696     if (eta<-0.7||eta>0.7)
697       continue;
698     Double_t phi  = mcP->Phi();
699     if (phi<1.0||phi>3.3)
700       continue;
701     // pion or eta meson or direct photon
702     if(mcP->GetPdgCode() == 111) {
703     } else if(mcP->GetPdgCode() == 221) {
704     } else if(mcP->GetPdgCode() == 22 ) {
705     } else
706       continue;
707     bool checkIfAlreadySaved = false;
708     for(Int_t imy=0;imy<fMyMcParts->GetEntries();imy++){
709       AliPhotonMcPartObj *mymc = static_cast<AliPhotonMcPartObj*>(fMyMcParts->At(imy));
710       if(!mymc)
711         continue;
712       if(mymc->fLabel == iTrack)
713         checkIfAlreadySaved = true;
714     }
715     if(!checkIfAlreadySaved)
716       FillMcPart(mcP, ipart++, iTrack);
717     for(Int_t id=mcP->GetFirstDaughter(); id <= mcP->GetLastDaughter(); id++){
718       if(id<=mcP->GetMother(0))
719         continue;
720       if(id<0 || id>nTracks)
721         continue;
722       TParticle *mcD = static_cast<TParticle*>(fStack->Particle(id));
723       if(!mcD)
724         continue;
725       FillMcPart(mcD, ipart++, id);
726     }
727   }
728 }
729
730 //________________________________________________________________________
731 void  AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t ipart, Int_t itrack)
732 {
733   // Fill MC particles.
734
735   if(!mcP)
736     return;
737   TVector3 vmcv(mcP->Vx(),mcP->Vy(), mcP->Vz());
738   AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(ipart));
739   mcp->fLabel = itrack ;
740   mcp->fPdg = mcP->GetPdgCode() ;
741   mcp->fPt = mcP->Pt() ;
742   mcp->fEta = mcP->Eta() ;
743   mcp->fPhi = mcP->Phi() ;
744   mcp->fVR = vmcv.Perp();
745   if(vmcv.Perp()>0){
746     mcp->fVEta = vmcv.Eta() ;
747     mcp->fVPhi = vmcv.Phi() ;
748   }
749   mcp->fMother = mcP->GetMother(0) ;
750 }
751
752 //________________________________________________________________________
753 Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
754 {
755   // Compute isolation based on tracks.
756   
757   Double_t trkIsolation = 0;
758   Double_t rad2 = radius*radius;
759   Int_t ntrks = fSelPrimTracks->GetEntries();
760   for(Int_t j = 0; j<ntrks; ++j) {
761     AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
762     if (!track)
763       continue;
764     if (track->Pt()<pt)
765       continue;
766     
767     Float_t eta = track->Eta();
768     Float_t phi = track->Phi();
769     Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
770     Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
771     if(dist>rad2)
772       continue;
773     trkIsolation += track->Pt();
774   } 
775   return trkIsolation;
776 }
777
778 //________________________________________________________________________
779 Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
780 {
781   // Get phi band.
782
783   if(!fSelPrimTracks)
784     return 0;
785   Int_t ntracks = fSelPrimTracks->GetEntries();
786   Double_t loweta = eta - R;
787   Double_t upeta = eta + R;
788   Double_t upphi = phi + R;
789   Double_t lowphi = phi - R;
790   Double_t et = 0;
791   for(Int_t itr=0; itr<ntracks; itr++){
792     AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(itr));
793     if(!track)
794       continue;
795     if(track->Pt()<minpt)
796       continue;
797     if((track->Eta() < upeta) && (track->Eta() > loweta))
798       continue;
799     if((track->Phi() > upphi) || (track->Phi() < lowphi))
800       continue;
801     et+=track->Pt();
802   }
803   return et;
804 }
805
806 //________________________________________________________________________
807 Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
808 {
809   // Calculate the energy of cross cells around the leading cell.
810
811   AliVCaloCells *cells = 0;
812   cells = fESD->GetEMCALCells();
813   if (!cells)
814     return 0;
815
816   if (!fGeom)
817     return 0;
818
819   Int_t iSupMod = -1;
820   Int_t iTower  = -1;
821   Int_t iIphi   = -1;
822   Int_t iIeta   = -1;
823   Int_t iphi    = -1;
824   Int_t ieta    = -1;
825   Int_t iphis   = -1;
826   Int_t ietas   = -1;
827
828   Double_t crossEnergy = 0;
829
830   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
831   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
832
833   Int_t ncells = cluster->GetNCells();
834   for (Int_t i=0; i<ncells; i++) {
835     Int_t cellAbsId = cluster->GetCellAbsId(i);
836     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
837     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
838     Int_t aphidiff = TMath::Abs(iphi-iphis);
839     if (aphidiff>1)
840       continue;
841     Int_t aetadiff = TMath::Abs(ieta-ietas);
842     if (aetadiff>1)
843       continue;
844     if ( (aphidiff==1 && aetadiff==0) ||
845         (aphidiff==0 && aetadiff==1) ) {
846       crossEnergy += cells->GetCellAmplitude(cellAbsId);
847     }
848   }
849
850   return crossEnergy;
851 }
852
853 //________________________________________________________________________
854 Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
855 {
856   // Get maximum energy of attached cell.
857
858   id = -1;
859
860   AliVCaloCells *cells = 0;
861   cells = fESD->GetEMCALCells();
862   if (!cells)
863     return 0;
864
865   Double_t maxe = 0;
866   Int_t ncells = cluster->GetNCells();
867   for (Int_t i=0; i<ncells; i++) {
868     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
869     if (e>maxe) {
870       maxe = e;
871       id   = cluster->GetCellAbsId(i);
872     }
873   }
874   return maxe;
875 }
876
877 //________________________________________________________________________
878 void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *) 
879 {
880   // Called once at the end of the query
881 /*  if(fIsGrid)
882     return;*/
883   if (fTree) {
884     printf("***tree %s being saved***\n",fTree->GetName());
885     TFile *f = OpenFile(2);
886     TDirectory::TContext context(f);
887     if (f) 
888       fTree->Write();
889   }
890 }