updated for Non-HFE systemaitc study
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFECal.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *               
3  *                                                                        *               
4  * Author: The ALICE Off-line Project.                                    *               
5  * Contributors are mentioned in the code where appropriate.              *               
6  *                                                                        *               
7  * Permission to use, copy, modify and distribute this software and its   *               
8  * documentation strictly for non-commercial purposes is hereby granted   *               
9  * without fee, provided that the above copyright notice appears in all   *               
10  * copies and that both the copyright notice and this permission notice   *               
11  * appear in the supporting documentation. The authors make no claims     *               
12  * about the suitability of this software for any purpose. It is          *               
13  * provided "as is" without express or implied warranty.                  *               
14  **************************************************************************/
15
16 // Class for heavy-flavour electron  with EMCal triggered events
17 // Author: Shingo Sakai
18
19
20 #include "TChain.h"
21 #include "TTree.h"
22 #include "TH2F.h"
23 #include "TF1.h"
24 #include "TMath.h"
25 #include "TCanvas.h"
26 #include "THnSparse.h"
27 #include "TLorentzVector.h"
28 #include "TString.h"
29 #include "TFile.h"
30 #include "TGraphErrors.h"
31
32 #include "AliAnalysisTask.h"
33 #include "AliAnalysisManager.h"
34
35 #include "AliESDEvent.h"
36 #include "AliESDHandler.h"
37 #include "AliAODEvent.h"
38 #include "AliAODHandler.h"
39
40 #include "AliMCEventHandler.h"
41 #include "AliMCEvent.h"
42 #include "AliMCParticle.h"
43
44 #include "AliAnalysisTaskHFECal.h"
45 #include "TGeoGlobalMagField.h"
46 #include "AliLog.h"
47 #include "AliAnalysisTaskSE.h"
48 #include "TRefArray.h"
49 #include "TVector.h"
50 #include "AliESDInputHandler.h"
51 #include "AliESDpid.h"
52 #include "AliESDtrackCuts.h"
53 #include "AliPhysicsSelection.h"
54 #include "AliESDCaloCluster.h"
55 #include "AliAODCaloCluster.h"
56 #include "AliEMCALRecoUtils.h"
57 #include "AliEMCALGeometry.h"
58 #include "AliGeomManager.h"
59 #include "stdio.h"
60 #include "TGeoManager.h"
61 #include "iostream"
62 #include "fstream"
63
64 #include "AliEMCALTrack.h"
65 #include "AliMagF.h"
66
67 #include "AliKFParticle.h"
68 #include "AliKFVertex.h"
69
70 #include "AliPID.h"
71 #include "AliPIDResponse.h"
72 #include "AliHFEcontainer.h"
73 #include "AliHFEcuts.h"
74 #include "AliHFEpid.h"
75 #include "AliHFEpidBase.h"
76 #include "AliHFEpidQAmanager.h"
77 #include "AliHFEtools.h"
78 #include "AliCFContainer.h"
79 #include "AliCFManager.h"
80
81 #include "AliStack.h"
82
83 #include "AliCentrality.h"
84
85 using namespace std;
86
87 ClassImp(AliAnalysisTaskHFECal)
88 //________________________________________________________________________
89 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name) 
90   : AliAnalysisTaskSE(name)
91   ,fESD(0)
92   ,fMC(0)
93   ,stack(0)
94   ,fGeom(0)
95   ,fOutputList(0)
96   ,fqahist(1) 
97   ,fTrackCuts(0)
98   ,fCuts(0)
99   ,fIdentifiedAsOutInz(kFALSE)
100   ,fPassTheEventCut(kFALSE)
101   ,fRejectKinkMother(kFALSE)
102   ,fmcData(kFALSE)
103   ,fVz(0.0)
104   ,fCFM(0)      
105   ,fPID(0)
106   ,fPIDqa(0)           
107   ,fOpeningAngleCut(0.1)
108   ,fInvmassCut(0)        // no mass
109   ,fSetMassConstraint(kTRUE)
110   ,fSetMassWidthCut(kFALSE)
111   ,fNoEvents(0)
112   ,fEMCAccE(0)
113   ,hEMCAccE(0)
114   ,fTrkpt(0)
115   ,fTrkEovPBef(0)        
116   ,fTrkEovPAft(0)       
117   ,fdEdxBef(0)   
118   ,fdEdxAft(0)   
119   ,fIncpT(0)    
120   ,fIncpTM20(0) 
121   ,fInvmassLS(0)                
122   ,fInvmassULS(0)               
123   ,fInvmassLSmc(0)              
124   ,fInvmassULSmc(0)             
125   ,fInvmassLSreco(0)            
126   ,fInvmassULSreco(0)           
127   ,fInvmassLSmc0(0)             
128   ,fInvmassLSmc1(0)             
129   ,fInvmassLSmc2(0)             
130   ,fInvmassLSmc3(0)             
131   ,fInvmassULSmc0(0)            
132   ,fInvmassULSmc1(0)            
133   ,fInvmassULSmc2(0)            
134   ,fInvmassULSmc3(0)            
135   ,fOpeningAngleLS(0)   
136   ,fOpeningAngleULS(0)  
137   ,fPhotoElecPt(0)
138   ,fPhoElecPt(0)
139   ,fPhoElecPtM20(0)
140   ,fSameElecPt(0)
141   ,fSameElecPtM20(0)
142   ,fTrackPtBefTrkCuts(0)         
143   ,fTrackPtAftTrkCuts(0)
144   ,fTPCnsigma(0)
145   ,fCent(0)
146   ,fEleInfo(0)
147   ,fElenSigma(0)
148   /*
149   ,fClsEBftTrigCut(0)    
150   ,fClsEAftTrigCut(0)    
151   ,fClsEAftTrigCut1(0)   
152   ,fClsEAftTrigCut2(0)   
153   ,fClsEAftTrigCut3(0)   
154   ,fClsEAftTrigCut4(0)   
155   ,fClsETime(0)       
156   ,fClsETime1(0)       
157   ,fTrigTimes(0)
158   ,fCellCheck(0)
159   */
160   ,fInputHFEMC(0)
161   ,fInputAlle(0)
162   ,fIncpTMChfe(0)       
163   ,fIncpTMChfeAll(0)    
164   ,fIncpTMCM20hfe(0)    
165   ,fIncpTMCM20hfeAll(0) 
166   ,fIncpTMCM20hfeCheck(0)       
167   ,fInputHFEMC_weight(0)
168   ,fIncpTMCM20hfeCheck_weight(0)
169   ,fIncpTMCpho(0)       
170   ,fIncpTMCM20pho(0)    
171   ,fPhoElecPtMC(0)
172   ,fPhoElecPtMCM20(0)
173   ,fSameElecPtMC(0)
174   ,fSameElecPtMCM20(0)
175   ,fIncpTMCM20pho_pi0e(0)       
176   ,fPhoElecPtMCM20_pi0e(0)
177   ,fSameElecPtMCM20_pi0e(0)
178   ,fIncpTMCM20pho_eta(0)        
179   ,fPhoElecPtMCM20_eta(0)
180   ,fSameElecPtMCM20_eta(0)
181   ,fIncpTMCpho_pi0e_TPC(0)      
182   ,fPhoElecPtMC_pi0e_TPC(0)
183   ,fSameElecPtMC_pi0e_TPC(0)
184   ,CheckNclust(0)
185   ,CheckNits(0)
186   ,Hpi0pTcheck(0)
187   ,HETApTcheck(0)
188   ,HphopTcheck(0)
189   ,HDpTcheck(0)
190   ,HBpTcheck(0)
191   ,fpTCheck(0)
192   ,fMomDtoE(0) 
193   ,fLabelCheck(0)
194   ,fgeoFake(0)
195   ,fFakeTrk0(0)
196   ,fFakeTrk1(0)
197   ,ftimingEle(0)
198   ,fIncMaxE(0)
199   ,fIncReco(0)
200   ,fPhoReco(0)
201   ,fSamReco(0) 
202   ,fIncRecoMaxE(0)
203   ,fPhoRecoMaxE(0)
204   ,fSamRecoMaxE(0) 
205   //,fnSigEtaCorr(NULL)
206 {
207   //Named constructor
208   
209   fPID = new AliHFEpid("hfePid");
210   fTrackCuts = new AliESDtrackCuts();
211   
212   for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
213
214   // Define input and output slots here
215   // Input slot #0 works with a TChain
216   DefineInput(0, TChain::Class());
217   // Output slot #0 id reserved by the base class for AOD
218   // Output slot #1 writes into a TH1 container
219   // DefineOutput(1, TH1I::Class());
220   DefineOutput(1, TList::Class());
221   //  DefineOutput(3, TTree::Class());
222 }
223
224 //________________________________________________________________________
225 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() 
226   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
227   ,fESD(0)
228   ,fMC(0)
229   ,stack(0)
230   ,fGeom(0)
231   ,fOutputList(0)
232   ,fqahist(1)
233   ,fTrackCuts(0)
234   ,fCuts(0)
235   ,fIdentifiedAsOutInz(kFALSE)
236   ,fPassTheEventCut(kFALSE)
237   ,fRejectKinkMother(kFALSE)
238   ,fmcData(kFALSE)
239   ,fVz(0.0)
240   ,fCFM(0)      
241   ,fPID(0)       
242   ,fPIDqa(0)           
243   ,fOpeningAngleCut(0.1)
244   ,fInvmassCut(0)        // no mass
245   ,fSetMassConstraint(kTRUE)
246   ,fSetMassWidthCut(kFALSE)
247   ,fNoEvents(0)
248   ,fEMCAccE(0)
249   ,hEMCAccE(0)
250   ,fTrkpt(0)
251   ,fTrkEovPBef(0)        
252   ,fTrkEovPAft(0)        
253   ,fdEdxBef(0)   
254   ,fdEdxAft(0)   
255   ,fIncpT(0)     
256   ,fIncpTM20(0)  
257   ,fInvmassLS(0)                
258   ,fInvmassULS(0)               
259   ,fInvmassLSmc(0)              
260   ,fInvmassULSmc(0)             
261   ,fInvmassLSreco(0)            
262   ,fInvmassULSreco(0)           
263   ,fInvmassLSmc0(0)             
264   ,fInvmassLSmc1(0)             
265   ,fInvmassLSmc2(0)             
266   ,fInvmassLSmc3(0)             
267   ,fInvmassULSmc0(0)            
268   ,fInvmassULSmc1(0)            
269   ,fInvmassULSmc2(0)            
270   ,fInvmassULSmc3(0)            
271   ,fOpeningAngleLS(0)   
272   ,fOpeningAngleULS(0)  
273   ,fPhotoElecPt(0)
274   ,fPhoElecPt(0)
275   ,fPhoElecPtM20(0)
276   ,fSameElecPt(0)
277   ,fSameElecPtM20(0)
278   ,fTrackPtBefTrkCuts(0)         
279   ,fTrackPtAftTrkCuts(0)                  
280   ,fTPCnsigma(0)
281   ,fCent(0)
282   ,fEleInfo(0)
283   ,fElenSigma(0)
284   /*
285   ,fClsEBftTrigCut(0)    
286   ,fClsEAftTrigCut(0)    
287   ,fClsEAftTrigCut1(0)   
288   ,fClsEAftTrigCut2(0)   
289   ,fClsEAftTrigCut3(0)   
290   ,fClsEAftTrigCut4(0)   
291   ,fClsETime(0)       
292   ,fClsETime1(0)       
293   ,fTrigTimes(0)
294   ,fCellCheck(0)
295   */
296   ,fInputHFEMC(0)
297   ,fInputAlle(0)
298   ,fIncpTMChfe(0)       
299   ,fIncpTMChfeAll(0)    
300   ,fIncpTMCM20hfe(0)    
301   ,fIncpTMCM20hfeAll(0) 
302   ,fIncpTMCM20hfeCheck(0)       
303   ,fInputHFEMC_weight(0)
304   ,fIncpTMCM20hfeCheck_weight(0)
305   ,fIncpTMCpho(0)       
306   ,fIncpTMCM20pho(0)    
307   ,fPhoElecPtMC(0)
308   ,fPhoElecPtMCM20(0)
309   ,fSameElecPtMC(0)
310   ,fSameElecPtMCM20(0)
311   ,fIncpTMCM20pho_pi0e(0)       
312   ,fPhoElecPtMCM20_pi0e(0)
313   ,fSameElecPtMCM20_pi0e(0)
314   ,fIncpTMCM20pho_eta(0)        
315   ,fPhoElecPtMCM20_eta(0)
316   ,fSameElecPtMCM20_eta(0)
317   ,fIncpTMCpho_pi0e_TPC(0)      
318   ,fPhoElecPtMC_pi0e_TPC(0)
319   ,fSameElecPtMC_pi0e_TPC(0)
320   ,CheckNclust(0)
321   ,CheckNits(0)
322   ,Hpi0pTcheck(0)
323   ,HETApTcheck(0)
324   ,HphopTcheck(0)
325   ,HDpTcheck(0)
326   ,HBpTcheck(0)
327   ,fpTCheck(0)
328   ,fMomDtoE(0)
329   ,fLabelCheck(0)
330   ,fgeoFake(0)
331   ,fFakeTrk0(0)
332   ,fFakeTrk1(0)
333   ,ftimingEle(0)
334   ,fIncMaxE(0)
335   ,fIncReco(0)
336   ,fPhoReco(0)
337   ,fSamReco(0)
338   ,fIncRecoMaxE(0)
339   ,fPhoRecoMaxE(0)
340   ,fSamRecoMaxE(0)
341   //,fnSigEtaCorr(NULL)
342 {
343         //Default constructor
344         fPID = new AliHFEpid("hfePid");
345
346         fTrackCuts = new AliESDtrackCuts();
347         
348        for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
349
350         // Constructor
351         // Define input and output slots here
352         // Input slot #0 works with a TChain
353         DefineInput(0, TChain::Class());
354         // Output slot #0 id reserved by the base class for AOD
355         // Output slot #1 writes into a TH1 container
356         // DefineOutput(1, TH1I::Class());
357         DefineOutput(1, TList::Class());
358         //DefineOutput(3, TTree::Class());
359 }
360 //_________________________________________
361
362 AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
363 {
364   //Destructor 
365   
366   delete fOutputList;
367   delete fGeom;
368   delete fPID;
369   delete fCuts;
370   delete fCFM;
371   delete fPIDqa;
372   delete fTrackCuts;
373 }
374 //_________________________________________
375
376 void AliAnalysisTaskHFECal::UserExec(Option_t*)
377 {
378   //Main loop
379   //Called for each event
380   
381   // create pointer to event
382   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
383   if (!fESD) {
384     printf("ERROR: fESD not available\n");
385     return;
386   }
387   
388   if(!fCuts){
389     AliError("HFE cuts not available");
390     return;
391   }
392   
393   if(!fPID->IsInitialized()){ 
394     // Initialize PID with the given run number
395     AliWarning("PID not initialised, get from Run no");
396     fPID->InitializePID(fESD->GetRunNumber());
397   }
398  
399   if(fmcData)fMC = MCEvent();
400   //AliStack* stack = NULL;
401   if(fmcData && fMC)stack = fMC->Stack();
402
403   Float_t cent = -1.;
404   AliCentrality *centrality = fESD->GetCentrality(); 
405   cent = centrality->GetCentralityPercentile("V0M");
406
407   //---- fill MC track info
408   if(fmcData && fMC)
409     {
410     Int_t nParticles = stack->GetNtrack();
411     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
412       TParticle* particle = stack->Particle(iParticle);
413       int fPDG = particle->GetPdgCode(); 
414       double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
415       double pTMC = particle->Pt();
416       double proR = particle->R();
417       double etaMC = particle->Eta();
418       if(fabs(etaMC)>0.6)continue;
419       Bool_t mcInDtoE= kFALSE;
420       Bool_t mcInBtoE= kFALSE;
421
422       Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
423       //if(!MChijing)printf("not MC hijing");
424       int iHijing = 1;
425       if(!MChijing)iHijing = 0;
426       double mcphoinfo[5];
427       mcphoinfo[0] = cent;
428       mcphoinfo[1] = pTMC;
429       mcphoinfo[2] = iHijing;
430       //if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
431       //if(fPDG==221)HETApTcheck->Fill(pTMC,iHijing);
432       if(fPDG==111)Hpi0pTcheck->Fill(mcphoinfo);
433       if(fPDG==221)HETApTcheck->Fill(mcphoinfo);
434       if(fabs(fPDG)==411 || fabs(fPDG)==413 || fabs(fPDG)==421 || fabs(fPDG)==423 || fabs(fPDG)==431)HDpTcheck->Fill(pTMC,iHijing);
435       if(fabs(fPDG)==511 || fabs(fPDG)==513 || fabs(fPDG)==521 || fabs(fPDG)==523 || fabs(fPDG)==531)HBpTcheck->Fill(pTMC,iHijing);
436
437       if(particle->GetFirstMother()>-1 && fPDG==22)
438         {
439          int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();  
440          if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
441         } 
442
443       if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
444         {
445             int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();  
446             double pTMCparent = stack->Particle(particle->GetFirstMother())->Pt();  
447             if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
448             if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
449             if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)
450                {
451                 fInputHFEMC->Fill(cent,pTMC);
452                 double mcinfo[5];
453                 mcinfo[0] = cent;
454                 mcinfo[1] = pTMC;
455                 mcinfo[2] = 0.0;
456                 mcinfo[3] = iHijing;
457                 mcinfo[4] = pTMCparent;
458                 fInputHFEMC_weight->Fill(mcinfo);
459                }
460          }
461
462
463          if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
464
465       }
466     } 
467
468   fNoEvents->Fill(0);
469
470   Int_t fNOtrks =  fESD->GetNumberOfTracks();
471   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
472   
473   Double_t pVtxZ = -999;
474   pVtxZ = pVtx->GetZ();
475   
476   if(TMath::Abs(pVtxZ)>10) return;
477   fNoEvents->Fill(1);
478   
479   if(fNOtrks<1) return;
480   fNoEvents->Fill(2);
481   
482   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
483   if(!pidResponse){
484     AliDebug(1, "Using default PID Response");
485     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); 
486   }
487   
488   fPID->SetPIDResponse(pidResponse);
489   
490   fCFM->SetRecEventInfo(fESD);
491   
492   //Float_t cent = -1.;
493   //AliCentrality *centrality = fESD->GetCentrality(); 
494   //cent = centrality->GetCentralityPercentile("V0M");
495   fCent->Fill(cent);
496   
497   //if(cent>90.) return;
498         
499  // Calorimeter info.
500  
501    FindTriggerClusters();
502
503   // make EMCAL array 
504   double maxE = 0.0;
505   for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
506      {
507       AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
508       if(clust && clust->IsEMCAL())
509         {
510          double clustE = clust->E();
511          float  emcx[3]; // cluster pos
512          clust->GetPosition(emcx);
513          TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
514          double emcphi = clustpos.Phi(); 
515          double emceta = clustpos.Eta();
516          double calInfo[5];
517          calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2(); 
518          //fEMCAccE->Fill(calInfo); 
519          //if(clustE>3.0)fEMCAccE->Fill(calInfo); 
520          //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo); 
521          hEMCAccE->Fill(cent,clustE); 
522          if(clustE>maxE)maxE = clustE; 
523         }
524    }
525
526   // Track loop 
527   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
528     AliESDtrack* track = fESD->GetTrack(iTracks);
529     if (!track) {
530       printf("ERROR: Could not receive track %d\n", iTracks);
531       continue;
532     }
533    
534     int parentlabel = 99999;
535     int parentPID = 99999;
536     int grand_parentlabel = 99999;
537     int grand_parentPID = 99999;
538     Bool_t mcPho = kFALSE;
539     Bool_t mcDtoE= kFALSE;
540     Bool_t mcBtoE= kFALSE;
541     Bool_t mcOrgPi0 = kFALSE;
542     Bool_t mcOrgEta = kFALSE;
543     double mcele = -1.;
544     double mcpT = 0.0;
545     double mcMompT = 0.0;
546     //double mcGrandMompT = 0.0;
547     double mcWeight = -10.0;
548
549     int iHijing = 1;
550     int mcLabel = -1;
551
552     if(fmcData && fMC && stack)
553       {
554        Int_t label = TMath::Abs(track->GetLabel());
555        mcLabel = track->GetLabel();
556        
557        if(mcLabel>-1)
558        {
559       
560                Bool_t MChijing = fMC->IsFromBGEvent(label);
561                if(!MChijing)iHijing = 0;
562
563                TParticle* particle = stack->Particle(label);
564                int mcpid = particle->GetPdgCode();
565                mcpT = particle->Pt();
566                //printf("MCpid = %d",mcpid);
567                if(particle->GetFirstMother()>-1)
568                {
569                        //int parentlabel = particle->GetFirstMother();
570                        //int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
571                        //mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
572                        FindMother(particle, parentlabel, parentPID);
573                        mcMompT = stack->Particle(parentlabel)->Pt();
574                        if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
575                        if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
576                        if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
577
578                        // make D->e pT correlation
579                        if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT); 
580
581                        //cout << "check PID = " << parentPID << endl;
582                        //cout << "check pho = " << mcPho << endl;
583                        //cout << "check D or B = " << mcDtoE << endl;
584                        // pi->e (Dalitz)
585                        if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
586                        {
587                                //cout << "find pi0->e " <<  endl;
588                                mcOrgPi0 = kTRUE;
589                                mcWeight = GetMCweight(mcMompT); 
590                        }
591                        // eta->e (Dalitz)
592                        if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
593                        {
594                                //cout << "find Eta->e " <<  endl;
595                                mcOrgEta = kTRUE;
596                                mcWeight = GetMCweightEta(mcMompT); 
597                        }
598
599                        // access grand parent 
600                        TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
601                        //if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
602                        if(particle_parent->GetFirstMother()>-1 && (parentPID==22 || parentPID==111) && fabs(mcpid)==11) // get grand parent g->e
603                        {
604                                //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
605                                //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
606                                FindMother(particle_parent, grand_parentlabel, grand_parentPID);
607                                double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
608                                if(grand_parentPID==111 && mcGrandpT>0.0)
609                                {
610                                        // check eta->pi0 decay !
611                                        int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
612                                        TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
613                                        FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
614                                        if(grand2_parentPID==221)
615                                        {
616                                                //cout << "find Eta->e " <<  endl;
617                                                double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
618                                                mcOrgEta = kTRUE;
619                                                mcWeight = GetMCweight(mcGrandpT2);  
620                                                mcMompT = mcGrandpT2; 
621                                        }
622                                        else
623                                        {
624                                                //cout << "find pi0->e " <<  endl;
625                                                mcOrgPi0 = kTRUE;
626                                                mcWeight = GetMCweight(mcGrandpT);  
627                                                mcMompT = mcGrandpT; 
628                                        }
629                                }
630
631                                if(grand_parentPID==221 && mcGrandpT>0.0)
632                                {
633                                        //cout << "find Eta->e " <<  endl;
634                                        mcOrgEta = kTRUE;
635                                        mcOrgPi0 = kFALSE;
636                                        mcWeight = GetMCweightEta(mcGrandpT); 
637                                        mcMompT = mcGrandpT; 
638                                }
639                        }
640                }
641
642                //cout << "===================="<<endl;
643                //cout << "mcDtoE : " << mcDtoE << endl; 
644                //cout << "mcBtoE : " << mcBtoE << endl; 
645                //cout << "mcPho : " << mcPho << endl; 
646
647                if(fabs(mcpid)==11)mcele= 0.; 
648                //cout << "check e: " << mcele << endl; 
649                if(fabs(mcpid)==11 && mcDtoE)mcele= 1.; 
650                //cout << "check D->e: " << mcele << endl; 
651                if(fabs(mcpid)==11 && mcBtoE)mcele= 2.; 
652                //cout << "check B->e: " << mcele << endl; 
653                if(fabs(mcpid)==11 && mcPho)mcele= 3.; 
654                //cout << "check Pho->e: " << mcele << endl; 
655
656                //cout << "check PID " << endl;
657                if(fabs(mcpid)!=11)
658                  {
659                   //cout << "!= 11" << endl;
660                   //cout << mcpid << endl;
661                  }
662                if(mcele==-1)
663                  {
664                   //cout << "mcele==-1" << endl;
665                   //cout << mcele << endl;
666                   //cout << mcpid << endl;
667                  }
668  
669        } // end of mcLabel>-1
670
671       } // end of MC info.
672  
673     //cout << "Pi0 = " << mcOrgPi0 << " ; Eta = " << mcOrgEta << endl; 
674     //printf("weight = %f\n",mcWeight);
675
676     if(TMath::Abs(track->Eta())>0.6) continue;
677     if(TMath::Abs(track->Pt()<2.0)) continue;
678     
679     fTrackPtBefTrkCuts->Fill(track->Pt());              
680     // RecKine: ITSTPC cuts  
681     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
682     
683     //RecKink
684     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
685       if(track->GetKinkIndex(0) != 0) continue;
686     } 
687     
688     // RecPrim
689     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
690     
691     // HFEcuts: ITS layers cuts
692     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
693     
694     // HFE cuts: TPC PID cleanup
695     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
696
697     int nTPCcl = track->GetTPCNcls();
698     //int nTPCclF = track->GetTPCNclsF(); // warnings
699     int nITS = track->GetNcls(0);
700     
701     fTrackPtAftTrkCuts->Fill(track->Pt());              
702     
703     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
704     pt = track->Pt();
705     if(pt<2.0)continue;
706     
707     // Track extrapolation
708     
709     //Int_t charge = track->Charge();
710     fTrkpt->Fill(pt);
711     mom = track->P();
712     phi = track->Phi();
713     eta = track->Eta();
714     dEdx = track->GetTPCsignal();
715     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
716
717     //cout << "nSigma correctoon-----" << endl;
718     //cout << "org = " << fTPCnSigma << endl; 
719     if(!fmcData) // nsigma eta correction
720        {
721         double nSigexpCorr = NsigmaCorrection(eta,cent);
722         fTPCnSigma -= nSigexpCorr;
723        }
724
725     //cout << "correction = " << fTPCnSigma << endl; 
726
727         double ncells = -1.0;
728         double m20 = -1.0;
729         double m02 = -1.0;
730         double disp = -1.0;
731         double rmatch = -1.0;  
732         double nmatch = -1.0;
733         double oppstatus = 0.0;
734         double emctof = 0.0;
735         Bool_t MaxEmatch = kFALSE;
736
737     Bool_t fFlagPhotonicElec = kFALSE;
738     Bool_t fFlagConvinatElec = kFALSE;
739
740     Int_t clsId = track->GetEMCALcluster();
741     if (clsId>0){
742       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
743       if(clust && clust->IsEMCAL()){
744
745                 double clustE = clust->E();
746                 if(clustE==maxE)MaxEmatch = kTRUE;
747                 eop = clustE/fabs(mom);
748                  //cout << "eop org = "<< eop << endl;
749                 if(mcLabel>-1.0)
750                   {
751                    double mceopcorr = MCEopMeanCorrection(pt,cent);
752                    eop += mceopcorr;
753                   }
754                 //cout << "eop corr = " << eop << endl;
755
756                 //double clustT = clust->GetTOF();
757                 ncells = clust->GetNCells();
758                 m02 = clust->GetM02();
759                 m20 = clust->GetM20();
760                 disp = clust->GetDispersion();
761                 double delphi = clust->GetTrackDx(); 
762                 double deleta = clust->GetTrackDz(); 
763                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
764                 nmatch = clust->GetNTracksMatched();
765                 emctof = clust->GetTOF();
766                 //cout << "emctof = " << emctof << endl;
767
768                 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
769                 {
770                   SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
771                 }
772                 if(fFlagPhotonicElec)oppstatus = 1.0;
773                 if(fFlagConvinatElec)oppstatus = 2.0;
774                 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
775
776                   double valdedx[16];
777                   valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
778                   //valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
779                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = nmatch; valdedx[9] = m20; valdedx[10] = mcpT;
780                   valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
781                   valdedx[15] = mcele;
782                   //if(fqahist==1)fEleInfo->Fill(valdedx);
783                  
784
785       }
786     }
787
788     //Get Cal info PID response
789     double eop2;
790     double ss[4];
791     Double_t nSigmaEop = fPID->GetPIDResponse()->NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
792     if(fTPCnSigma>-1.5 && fTPCnSigma<3.0 && nITS>2.5 && nTPCcl>100)
793       {
794        double valEop[3];
795        valEop[0] = cent;
796        valEop[1] = pt;
797        valEop[2] = nSigmaEop;
798        fElenSigma->Fill(valEop);
799       }
800
801    // ============ PID
802
803     if(nITS<2.5)continue;
804     if(nTPCcl<100)continue;
805    
806     CheckNclust->Fill(nTPCcl); 
807     CheckNits->Fill(nITS); 
808
809     fdEdxBef->Fill(mom,fTPCnSigma);
810     fTPCnsigma->Fill(mom,fTPCnSigma);
811     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
812
813     Int_t pidpassed = 1;
814  
815     // check reco eff. with TPC
816
817     double phoval[5];
818     phoval[0] = cent;
819     phoval[1] = pt;
820     phoval[2] = fTPCnSigma;
821     phoval[3] = iHijing;
822     phoval[4] = mcMompT;
823    
824     if((fTPCnSigma >= -1.0 && fTPCnSigma <= 3) && mcele>-1 && mcPho && mcOrgPi0)
825       {
826         if(iHijing==1)mcWeight = 1.0; 
827         fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight);    
828         if(fFlagPhotonicElec) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
829         if(fFlagConvinatElec) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
830       }
831
832     //--- track accepted
833     AliHFEpidObject hfetrack;
834     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
835     hfetrack.SetRecTrack(track);
836     double binct = 10.5;
837     if((0.0< cent) && (cent<5.0)) binct = 0.5;
838     if((5.0< cent) && (cent<10.0)) binct = 1.5;
839     if((10.0< cent) && (cent<20.0)) binct = 2.5;
840     if((20.0< cent) && (cent<30.0)) binct = 3.5;
841     if((30.0< cent) && (cent<40.0)) binct = 4.5;
842     if((40.0< cent) && (cent<50.0)) binct = 5.5;
843     if((50.0< cent) && (cent<60.0)) binct = 6.5;
844     if((60.0< cent) && (cent<70.0)) binct = 7.5;
845     if((70.0< cent) && (cent<80.0)) binct = 8.5;
846     if((80.0< cent) && (cent<90.0)) binct = 9.5;
847     if((90.0< cent) && (cent<100.0)) binct = 10.5;
848
849     hfetrack.SetCentrality((int)binct); //added
850     hfetrack.SetPbPb();
851     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
852
853     if(pidpassed==0) continue;
854     
855     fTrkEovPAft->Fill(pt,eop);
856     fdEdxAft->Fill(mom,fTPCnSigma);
857
858     // Fill real data
859     fIncpT->Fill(cent,pt);    
860     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
861     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
862
863     if(m20>0.0 && m20<0.3)
864       { 
865        fIncpTM20->Fill(cent,pt);   
866        ftimingEle->Fill(pt,emctof); 
867        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
868        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
869      }
870     
871  
872     //--------
873      
874     double recopT =  SumpT(iTracks,track);
875
876     if(m20>0.0 && m20<0.3)
877       {
878        if(MaxEmatch)fIncMaxE->Fill(cent,pt);
879        if(pt>5.0)
880          {
881           fIncReco->Fill(cent,recopT);
882           if(fFlagPhotonicElec) fPhoReco->Fill(cent,recopT);
883           if(fFlagConvinatElec) fSamReco->Fill(cent,recopT);
884           if(MaxEmatch)
885             {
886              fIncRecoMaxE->Fill(cent,recopT);
887              if(fFlagPhotonicElec) fPhoRecoMaxE->Fill(cent,recopT);
888              if(fFlagConvinatElec) fSamRecoMaxE->Fill(cent,recopT);
889             }
890          }
891      }
892
893     // MC
894     // check label for electron candidiates
895
896     int idlabel = 1;
897     if(mcLabel==0)idlabel = 0;
898     fLabelCheck->Fill(pt,idlabel);
899     if(mcLabel==0)fgeoFake->Fill(phi,eta);
900
901     if(mcLabel<0 && m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)
902        {
903         fFakeTrk0->Fill(cent,pt);
904        }
905
906     if(mcele>-1) // select MC electrons
907       {
908
909           fIncpTMChfeAll->Fill(cent,pt);    
910           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);    
911           if(m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)fFakeTrk1->Fill(cent,pt);
912
913        if(mcBtoE || mcDtoE) // select B->e & D->e
914          {
915           fIncpTMChfe->Fill(cent,pt);    
916           //if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);    
917           //if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT);    
918           if(m20>0.0 && m20<0.3)
919             {
920                 //cout << "MC label = " << mcLabel << endl;
921                 fIncpTMCM20hfe->Fill(cent,pt);    
922                 fIncpTMCM20hfeCheck->Fill(cent,mcpT);    
923                 fIncpTMCM20hfeCheck_weight->Fill(phoval);    
924             }
925          }
926       
927        if(mcPho) // select photonic electrons
928         {
929
930          fIncpTMCpho->Fill(phoval);    
931          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
932          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
933
934          if(m20>0.0 && m20<0.3) 
935            {
936             fIncpTMCM20pho->Fill(phoval);    
937             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
938             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
939             // pi0->g->e
940             if(mcWeight>-1)
941               {
942                if(iHijing==1)mcWeight = 1.0; 
943                if(mcOrgPi0)
944                  {
945                   fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);    
946                   if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
947                   if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
948                   //fIncpTMCM20pho_pi0e->Fill(phoval);   // v5-04-02-AN & v5-04-06-AN 
949                   //if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval);
950                   //if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval);
951                  }
952                if(mcOrgEta)
953                  {
954                   fIncpTMCM20pho_eta->Fill(phoval,mcWeight);    
955                   if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
956                   if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
957                   //fIncpTMCM20pho_eta->Fill(phoval);    
958                   //if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval);
959                   //if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval);
960                  }
961                // --- eta
962               }
963            }
964         } 
965       } 
966  }
967  PostData(1, fOutputList);
968 }
969 //_________________________________________
970 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
971 {
972   //--- Check MC
973  
974   //Bool_t mcData = kFALSE;
975   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
976     {
977      fmcData = kTRUE;
978      printf("+++++ MC Data available");
979     }
980   if(fmcData)
981     {
982      printf("++++++++= MC analysis \n");
983     }
984   else
985    {
986      printf("++++++++ real data analysis \n");
987    }
988
989    printf("+++++++ QA hist %d",fqahist);
990
991   //---- Geometry
992   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
993
994   //--------Initialize PID
995   //fPID->SetHasMCData(kFALSE);
996   fPID->SetHasMCData(fmcData);
997   if(!fPID->GetNumberOfPIDdetectors()) 
998     {
999       fPID->AddDetector("TPC", 0);
1000       fPID->AddDetector("EMCAL", 1);
1001     }
1002   
1003   Double_t params[4];
1004   const char *cutmodel;
1005   cutmodel = "pol0";
1006   params[0] = -1.0; //sigma min
1007   double maxnSig = 3.0;
1008   if(fmcData)
1009     {
1010      params[0] = -5.0; //sigma min
1011      maxnSig = 5.0; 
1012     } 
1013
1014   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
1015
1016   fPID->SortDetectors(); 
1017   fPIDqa = new AliHFEpidQAmanager();
1018   fPIDqa->Initialize(fPID);
1019  
1020   //------- fcut --------------  
1021   fCuts = new AliHFEcuts();
1022   fCuts->CreateStandardCuts();
1023   //fCuts->SetMinNClustersTPC(100);
1024   fCuts->SetMinNClustersTPC(90);
1025   fCuts->SetMinRatioTPCclusters(0.6);
1026   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
1027   //fCuts->SetMinNClustersITS(3);
1028   fCuts->SetMinNClustersITS(2);
1029   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
1030   fCuts->SetCheckITSLayerStatus(kFALSE);
1031   fCuts->SetVertexRange(10.);
1032   fCuts->SetTOFPIDStep(kFALSE);
1033   fCuts->SetPtRange(2, 50);
1034   fCuts->SetMaxImpactParam(3.,3.);
1035
1036   //--------Initialize correction Framework and Cuts
1037   fCFM = new AliCFManager;
1038   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
1039   fCFM->SetNStepParticle(kNcutSteps);
1040   for(Int_t istep = 0; istep < kNcutSteps; istep++)
1041     fCFM->SetParticleCutsList(istep, NULL);
1042   
1043   if(!fCuts){
1044     AliWarning("Cuts not available. Default cuts will be used");
1045     fCuts = new AliHFEcuts;
1046     fCuts->CreateStandardCuts();
1047   }
1048   fCuts->Initialize(fCFM);
1049   
1050   //---------Output Tlist
1051   fOutputList = new TList();
1052   fOutputList->SetOwner();
1053   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
1054   
1055   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
1056   fOutputList->Add(fNoEvents);
1057   
1058   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
1059   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5}; 
1060   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5}; 
1061   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
1062   //if(fqahist==1)fOutputList->Add(fEMCAccE);
1063
1064   hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
1065   fOutputList->Add(hEMCAccE);
1066
1067   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
1068   fOutputList->Add(fTrkpt);
1069   
1070   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
1071   fOutputList->Add(fTrackPtBefTrkCuts);
1072   
1073   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
1074   fOutputList->Add(fTrackPtAftTrkCuts);
1075   
1076   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
1077   fOutputList->Add(fTPCnsigma);
1078   
1079   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
1080   fOutputList->Add(fTrkEovPBef);
1081   
1082   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
1083   fOutputList->Add(fTrkEovPAft);
1084   
1085   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
1086   fOutputList->Add(fdEdxBef);
1087   
1088   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
1089   fOutputList->Add(fdEdxAft);
1090   
1091   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
1092   fOutputList->Add(fIncpT);
1093
1094   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1095   fOutputList->Add(fIncpTM20);
1096   
1097   Int_t nBinspho[9] =  { 10,  30,   600, 60,   50,    4,  40,   8,  30};
1098   Double_t minpho[9] = {  0.,  0., -0.1, 40,  0, -0.5,   0,-1.5,   0};   
1099   Double_t maxpho[9] = {100., 30.,  0.5, 100,   1,  3.5,   2, 6.5,  30};   
1100
1101   fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
1102   if(fqahist==1)fOutputList->Add(fInvmassLS);
1103   
1104   fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
1105   if(fqahist==1)fOutputList->Add(fInvmassULS);
1106   
1107   fInvmassLSmc = new THnSparseD("fInvmassLSmc", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
1108   if(fqahist==1)fOutputList->Add(fInvmassLSmc);
1109   
1110   fInvmassULSmc = new THnSparseD("fInvmassULSmc", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
1111   if(fqahist==1)fOutputList->Add(fInvmassULSmc);
1112
1113   fInvmassLSreco = new TH2D("fInvmassLSreco", "Inv mass of LS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1114   fInvmassLSreco->Sumw2();
1115   fOutputList->Add(fInvmassLSreco);
1116
1117   fInvmassULSreco = new TH2D("fInvmassULSreco", "Inv mass of ULS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1118   fInvmassULSreco->Sumw2();
1119   fOutputList->Add(fInvmassULSreco);
1120
1121   fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1122   fInvmassLSmc0->Sumw2();
1123   fOutputList->Add(fInvmassLSmc0);
1124   
1125   fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1126   fInvmassLSmc1->Sumw2();
1127   fOutputList->Add(fInvmassLSmc1);
1128
1129   fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1130   fInvmassLSmc2->Sumw2();
1131   fOutputList->Add(fInvmassLSmc2);
1132
1133   fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1134   fInvmassLSmc3->Sumw2();
1135   fOutputList->Add(fInvmassLSmc3);
1136
1137   fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1138   fInvmassULSmc0->Sumw2();
1139   fOutputList->Add(fInvmassULSmc0);
1140
1141   fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1142   fInvmassULSmc1->Sumw2();
1143   fOutputList->Add(fInvmassULSmc1);
1144
1145   fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1146   fInvmassULSmc2->Sumw2();
1147   fOutputList->Add(fInvmassULSmc2);
1148
1149   fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1150   fInvmassULSmc3->Sumw2();
1151   fOutputList->Add(fInvmassULSmc3);
1152
1153   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1154   fOutputList->Add(fOpeningAngleLS);
1155   
1156   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1157   fOutputList->Add(fOpeningAngleULS);
1158   
1159   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1160   fOutputList->Add(fPhotoElecPt);
1161   
1162   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
1163   fOutputList->Add(fPhoElecPt);
1164   
1165   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
1166   fOutputList->Add(fPhoElecPtM20);
1167
1168   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
1169   fOutputList->Add(fSameElecPt);
1170
1171   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
1172   fOutputList->Add(fSameElecPtM20);
1173
1174   fCent = new TH1F("fCent","Centrality",200,0,100) ;
1175   fOutputList->Add(fCent);
1176  
1177   // Make common binning
1178   const Double_t kMinP = 0.;
1179   const Double_t kMaxP = 20.;
1180
1181   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
1182   Int_t nBins[16] =  {  100,     7,  60,    20,    90,  250,   25,   40,   10, 200,  100, 100,  500,    5, 100,    8};
1183   Double_t min[16] = {kMinP,  -0.5, 1.0,  -1.0,  -5.0,    0,    0,    0,  0.0, 0.0,  0.0,   0,    0, -0.5,  80, -1.5};
1184   Double_t max[16] = {kMaxP,   6.5, 4.0,   1.0,   4.0,  2.5, 0.05,   40,   10, 2.0, 20.0, 100,  100,  4.5, 180,  6.5};
1185   fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;clsF;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
1186   if(fqahist==1)fOutputList->Add(fEleInfo);
1187
1188   // Make common binning
1189   Int_t nBinsEop[3] =  { 10, 50, 100};
1190   Double_t minEop[3] = {  0,  0,  -5};
1191   Double_t maxEop[3] = {100, 50,   5};
1192   fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop);
1193   fOutputList->Add(fElenSigma);
1194
1195
1196   //<---  Trigger info
1197   /*
1198   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1199   fOutputList->Add(fClsEBftTrigCut);
1200
1201   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1202   fOutputList->Add(fClsEAftTrigCut);
1203
1204   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1205   fOutputList->Add(fClsEAftTrigCut1);
1206
1207   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1208   fOutputList->Add(fClsEAftTrigCut2);
1209
1210   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1211   fOutputList->Add(fClsEAftTrigCut3);
1212
1213   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1214   fOutputList->Add(fClsEAftTrigCut4);
1215
1216   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1217   fOutputList->Add(fClsETime);
1218
1219   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1220   fOutputList->Add(fClsETime1);
1221
1222   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1223   fOutputList->Add(fTrigTimes);
1224
1225   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1226   fOutputList->Add(fCellCheck);
1227   */
1228   //<---------- MC
1229
1230   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1231   fOutputList->Add(fInputHFEMC);
1232
1233   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
1234   fOutputList->Add(fInputAlle);
1235
1236   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1237   fOutputList->Add(fIncpTMChfe);
1238
1239   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
1240   fOutputList->Add(fIncpTMChfeAll);
1241
1242   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1243   fOutputList->Add(fIncpTMCM20hfe);
1244
1245   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
1246   fOutputList->Add(fIncpTMCM20hfeAll);
1247
1248   fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1249   fOutputList->Add(fIncpTMCM20hfeCheck);
1250
1251   Int_t nBinspho2[5] =  { 200, 100,    7,    3, 700};
1252   Double_t minpho2[5] = {  0.,  0., -2.5, -0.5, 0.};   
1253   Double_t maxpho2[5] = {100., 50.,  4.5,  2.5, 70.};   
1254   
1255   fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1256   fOutputList->Add(fInputHFEMC_weight);
1257   
1258   fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1259   fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1260   
1261   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
1262   fOutputList->Add(fIncpTMCpho);
1263
1264   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1265   fOutputList->Add(fIncpTMCM20pho);
1266
1267   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1268   fOutputList->Add(fPhoElecPtMC);
1269   
1270   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1271   fOutputList->Add(fPhoElecPtMCM20);
1272
1273   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1274   fOutputList->Add(fSameElecPtMC);
1275
1276   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1277   fOutputList->Add(fSameElecPtMCM20);
1278
1279   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1280   fIncpTMCM20pho_pi0e->Sumw2();
1281   fOutputList->Add(fIncpTMCM20pho_pi0e);
1282
1283   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1284   fPhoElecPtMCM20_pi0e->Sumw2();
1285   fOutputList->Add(fPhoElecPtMCM20_pi0e);
1286
1287   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1288   fSameElecPtMCM20_pi0e->Sumw2();
1289   fOutputList->Add(fSameElecPtMCM20_pi0e);
1290  // 
1291   fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1292   fIncpTMCM20pho_eta->Sumw2();
1293   fOutputList->Add(fIncpTMCM20pho_eta);
1294
1295   fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1296   fPhoElecPtMCM20_eta->Sumw2();
1297   fOutputList->Add(fPhoElecPtMCM20_eta);
1298
1299   fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1300   fSameElecPtMCM20_eta->Sumw2();
1301   fOutputList->Add(fSameElecPtMCM20_eta);
1302   // ------------
1303   fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1304   fIncpTMCpho_pi0e_TPC->Sumw2();
1305   fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1306
1307   fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1308   fPhoElecPtMC_pi0e_TPC->Sumw2();
1309   fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1310
1311   fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1312   fSameElecPtMC_pi0e_TPC->Sumw2();
1313   fOutputList->Add(fSameElecPtMC_pi0e_TPC);
1314   //-------------
1315
1316   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1317   fOutputList->Add(CheckNclust);
1318
1319   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1320   fOutputList->Add(CheckNits);
1321   /*
1322   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1323   fOutputList->Add(Hpi0pTcheck);
1324
1325   HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1326   fOutputList->Add(HETApTcheck);
1327   */
1328
1329   Int_t nBinspho3[3] =  { 200, 100, 3};
1330   Double_t minpho3[3] = {  0.,  0., -0.5};   
1331   Double_t maxpho3[3] = {100., 50., 2.5};   
1332
1333   Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1334   fOutputList->Add(Hpi0pTcheck);
1335
1336   HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1337   fOutputList->Add(HETApTcheck);
1338   //--
1339   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1340   fOutputList->Add(HphopTcheck);
1341   //
1342   HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1343   fOutputList->Add(HDpTcheck);
1344
1345   HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1346   fOutputList->Add(HBpTcheck);
1347   //
1348
1349   fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1350   fOutputList->Add(fpTCheck);
1351
1352   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1353   fOutputList->Add(fMomDtoE);
1354
1355   fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1356   fOutputList->Add(fLabelCheck);
1357
1358   fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1359   fOutputList->Add(fgeoFake);
1360
1361   fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20);
1362   fOutputList->Add(fFakeTrk0);
1363
1364   fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20);
1365   fOutputList->Add(fFakeTrk1);
1366
1367   ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1368   fOutputList->Add(ftimingEle);
1369
1370   // eta correction
1371   // note: parameters 01/31new.TPCnSigmaEtaDep
1372   // 70-90 delta_eta = 0.2
1373
1374   double etaval[12] = {-0.55,-0.45,-0.35,-0.25,-0.15,-0.05,0.05,0.15,0.25,0.35,0.45,0.55};
1375   double corr0[12]= {-0.569177,-0.528844,-0.391979,-0.165494,0.0283495,0.156171,0.266353,0.13103,-0.0250842,-0.274089,-0.45481,-0.536291}; // 0-10 (done)
1376   double corr1[12]= {-0.404742,-0.278953,-0.218069,0.00139927,0.191412,0.354403,0.524594,0.341778,0.244199,-0.112146,-0.160692,-0.352832}; // 10-20 (done)
1377   double corr2[12] = {-0.306007,-0.16821,-0.0248635,0.202233,0.447051,0.497197,0.712251,0.433482,0.337907,0.168426,-0.0693229,-0.0728351}; // 20-30 (done)
1378   double corr3[12] = {-0.13884,-0.0503553,0.104403,0.389773,0.50697,0.539048,0.751642,0.655636,0.518563,0.308156,0.0361159,-0.0491439}; // 30-40 (done)
1379   double corr4[12] = {-0.0319431,0.0808711,0.208774,0.443217,0.557762,0.61453,0.889519,0.808282,0.620394,0.267092,0.15241,-0.0458664}; // 40-50 (done)
1380   double corr5[12] = {-0.130625,0.0189124,0.190344,0.467431,0.546353,0.672251,0.731541,0.802101,0.437108,0.294081,0.193682,0.159074}; // 50-70(done)
1381   double corr6[12] = {0.0600197,0.0600197,0.358366,0.358366,0.973734,0.973734,0.759812,0.759812,0.667861,0.667861,0.415635,0.415635}; // 70-90(done)
1382  
1383   fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10
1384   fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20
1385   fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30
1386   fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40 
1387   fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50
1388   fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70
1389   fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90
1390
1391   fIncMaxE = new TH2D("fIncMaxE","Inc",10,0,100,10,0,100);
1392   fOutputList->Add(fIncMaxE);
1393
1394   fIncReco = new TH2D("fIncReco","Inc",10,0,100,100,0,500);
1395   fOutputList->Add(fIncReco);
1396
1397   fPhoReco = new TH2D("fPhoReco","Pho",10,0,100,100,0,500);
1398   fOutputList->Add(fPhoReco);
1399
1400   fSamReco = new TH2D("fSamReco","Same",10,0,100,100,0,500);
1401   fOutputList->Add(fSamReco);
1402
1403   fIncRecoMaxE = new TH2D("fIncRecoMaxE","Inc",10,0,100,100,0,500);
1404   fOutputList->Add(fIncRecoMaxE);
1405
1406   fPhoRecoMaxE = new TH2D("fPhoRecoMaxE","Pho",10,0,100,100,0,500);
1407   fOutputList->Add(fPhoRecoMaxE);
1408
1409   fSamRecoMaxE = new TH2D("fSamRecoMaxE","Same",10,0,100,100,0,500);
1410   fOutputList->Add(fSamRecoMaxE);
1411
1412   PostData(1,fOutputList);
1413 }
1414
1415 //________________________________________________________________________
1416 void AliAnalysisTaskHFECal::Terminate(Option_t *)
1417 {
1418   // Info("Terminate");
1419         AliAnalysisTaskSE::Terminate();
1420 }
1421
1422 //________________________________________________________________________
1423 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1424 {
1425   // Check single track cuts for a given cut step
1426   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1427   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1428   return kTRUE;
1429 }
1430 //_________________________________________
1431 void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent, Bool_t tagpi0, Bool_t tageta)
1432 {
1433   //Identify non-heavy flavour electrons using Invariant mass method
1434   
1435   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1436   fTrackCuts->SetRequireTPCRefit(kTRUE);
1437   fTrackCuts->SetRequireITSRefit(kTRUE);
1438   fTrackCuts->SetEtaRange(-0.9,0.9);
1439   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1440   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1441   fTrackCuts->SetMinNClustersTPC(90);
1442   
1443   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1444   
1445   double ptEle = track->Pt();  //add
1446   if(ibgevent==0 && w > 0.0)
1447      {
1448       fpTCheck->Fill(ptEle,w);
1449      }
1450
1451   Bool_t flagPhotonicElec = kFALSE;
1452   Bool_t flagConvinatElec = kFALSE;
1453   
1454   int p1 = 0;
1455   if(mce==3)
1456      {
1457        Int_t label = TMath::Abs(track->GetLabel());
1458        TParticle* particle = stack->Particle(label);
1459        p1 = particle->GetFirstMother();
1460      }
1461
1462   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1463   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1464     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1465     if (!trackAsso) {
1466       printf("ERROR: Could not receive track %d\n", jTracks);
1467       continue;
1468     }
1469     if(itrack==jTracks)continue;    
1470     int jbgevent = 0;    
1471
1472     int p2 = 0;
1473     if(mce==3)
1474     {
1475       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1476       TParticle* particle2 = stack->Particle(label2);
1477       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1478       if(MChijing_ass)jbgevent =1;
1479       if(particle2->GetFirstMother()>-1)
1480          p2 = particle2->GetFirstMother();
1481     }
1482
1483     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1484     Double_t mass=999., width = -999;
1485     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1486     
1487     //ptPrim = track->Pt();
1488     ptPrim = ptEle;
1489
1490     dEdxAsso = trackAsso->GetTPCsignal();
1491     ptAsso = trackAsso->Pt();
1492     Int_t chargeAsso = trackAsso->Charge();
1493     Int_t charge = track->Charge();
1494     
1495
1496     if(ptAsso <0.5) continue;
1497     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1498     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1499     
1500     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1501     if(charge>0) fPDGe1 = -11;
1502     if(chargeAsso>0) fPDGe2 = -11;
1503  
1504     //printf("chargeAsso = %d\n",chargeAsso);
1505     //printf("charge = %d\n",charge);
1506     if(charge == chargeAsso) fFlagLS = kTRUE;
1507     if(charge != chargeAsso) fFlagULS = kTRUE;
1508     
1509     //printf("fFlagLS = %d\n",fFlagLS);
1510     //printf("fFlagULS = %d\n",fFlagULS);
1511     printf("\n");
1512
1513     AliKFParticle ge1(*track, fPDGe1);
1514     AliKFParticle ge2(*trackAsso, fPDGe2);
1515     AliKFParticle recg(ge1, ge2);
1516     
1517     // vertex 
1518     AliKFVertex primV(*pVtx);
1519     primV += recg;
1520     recg.SetProductionVertex(primV);
1521     
1522     // check chi2
1523     if(recg.GetNDF()<1) continue;
1524     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1525     Double_t chi2cut = 3.0;
1526
1527     // mass.
1528     if(fSetMassConstraint)
1529       {
1530        recg.SetMassConstraint(0,0.0001);
1531        chi2cut = 30.0;
1532       }
1533     recg.GetMass(mass,width);
1534
1535     if(fSetMassWidthCut && width>1e10)continue;
1536
1537         // angle   
1538     openingAngle = ge1.GetAngle(ge2);
1539     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1540     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1541     
1542     double ishower = 0;
1543     if(shower>0.0 && shower<0.3)ishower = 1;
1544
1545     double phoinfo[9];
1546     phoinfo[0] = cent;
1547     phoinfo[1] = ptPrim;
1548     phoinfo[2] = mass;
1549     phoinfo[3] = nSig;
1550     //phoinfo[3] = dEdxAsso;
1551     phoinfo[4] = openingAngle;
1552     phoinfo[5] = ishower;
1553     phoinfo[6] = ep;
1554     phoinfo[7] = mce;
1555     phoinfo[8] = ptAsso;
1556
1557     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1558     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1559     if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1560     if(fFlagULS && ibgevent==0 && jbgevent==0)
1561        {
1562          fInvmassULSmc->Fill(phoinfo,w);
1563        }
1564     //printf("fInvmassCut %f\n",fInvmassCut);
1565     //printf("openingAngle %f\n",fOpeningAngleCut);
1566
1567     // angle cut
1568     if(openingAngle > fOpeningAngleCut) continue;
1569     // chi2 cut
1570     //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1571     if(chi2recg>chi2cut) continue;
1572  
1573     if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1574     if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1575   
1576     // for real data  
1577     //printf("mce =%f\n",mce);
1578     if(mce<-0.5) // mce==-1. is real
1579        {
1580          //printf("Real data\n");
1581          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1582                flagPhotonicElec = kTRUE;
1583               }
1584          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1585                flagConvinatElec = kTRUE;
1586               }
1587         }
1588     // for MC data  
1589     else
1590        {
1591          //printf("MC data\n");
1592
1593          if(w>0.0)
1594            {
1595            //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1596            if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1597            if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1598            if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1599            if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1600            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1601            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1602            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1603            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1604           }
1605
1606          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1607          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1608                flagPhotonicElec = kTRUE;
1609               }
1610          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1611                flagConvinatElec = kTRUE;
1612               }
1613         }
1614
1615   }
1616   fFlagPhotonicElec = flagPhotonicElec;
1617   fFlagConvinatElec = flagConvinatElec;
1618   
1619 }
1620 //-------------------------------------------
1621
1622 void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1623 {
1624  //int label = 99999;
1625  //int pid = 99999;
1626
1627  if(part->GetFirstMother()>-1)
1628    {
1629     label = part->GetFirstMother();
1630     pid = stack->Particle(label)->GetPdgCode();
1631    }
1632    //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
1633 }
1634
1635 double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1636 {
1637         double weight = 1.0;
1638
1639         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1640         {
1641                 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1642         }
1643         else
1644         {
1645                 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1646         }
1647   return weight;
1648 }
1649
1650 double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1651 {
1652   double weight = 1.0;
1653
1654   weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1655   return weight;
1656 }
1657
1658
1659 //_________________________________________
1660 void AliAnalysisTaskHFECal::FindTriggerClusters()
1661 {
1662   //cout << "finding trigger patch" << endl; 
1663   // constants
1664   const int nModuleCols = 2;
1665   const int nModuleRows = 5;
1666   const int nColsFeeModule = 48;
1667   const int nRowsFeeModule = 24;
1668   const int nColsFaltroModule = 24;
1669   const int nRowsFaltroModule = 12;
1670   //const int faltroWidthMax = 20;
1671
1672   // part 1, trigger extraction -------------------------------------
1673   Int_t globCol, globRow;
1674   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1675   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1676
1677   //Int_t trigtimes[faltroWidthMax];
1678   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1679   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1680   //Double_t fTrigCutLow = 6;
1681   //Double_t fTrigCutHigh = 10;
1682   Double_t fTimeCutLow =  469e-09;
1683   Double_t fTimeCutHigh = 715e-09;
1684
1685   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1686
1687   // erase trigger maps
1688   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1689   {
1690     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1691     {
1692       ftriggersCut[i][j] = 0;
1693       ftriggers[i][j] = 0;
1694       ftriggersTime[i][j] = 0;
1695     }
1696   }
1697
1698   Int_t iglobCol=0, iglobRow=0;
1699   // go through triggers
1700   if( fCaloTrigger->GetEntries() > 0 )
1701   {
1702     // needs reset
1703     fCaloTrigger->Reset();
1704     while( fCaloTrigger->Next() )
1705     {
1706       fCaloTrigger->GetPosition( globCol, globRow );
1707       fCaloTrigger->GetNL0Times( ntimes );
1708       /*
1709       // no L0s
1710       if( ntimes < 1 )   continue;
1711       // get precise timings
1712       fCaloTrigger->GetL0Times( trigtimes );
1713       trigInCut = 0;
1714       for(Int_t i = 0; i < ntimes; i++ )
1715       {
1716         // save the first trigger time in channel
1717         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1718           triggersTime[globCol][globRow] = trigtimes[i];
1719         //printf("trigger times: %d\n",trigtimes[i]);
1720         // check if in cut
1721         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1722           trigInCut = 1;
1723
1724         fTrigTimes->Fill(trigtimes[i]);
1725       }
1726       */
1727    
1728       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1729       Int_t bit = 0;
1730       fCaloTrigger->GetTriggerBits(bit);
1731       //cout << "bit = " << bit << endl;
1732
1733       Int_t ts = 0;
1734       fCaloTrigger->GetL1TimeSum(ts);
1735       //cout << "ts = " << ts << endl;
1736       if (ts > 0)ftriggers[globCol][globRow] = 1;
1737       // number of triggered channels in event
1738       nTrigChannel++;
1739       // ... inside cut
1740       if(ts>0 && (bit >> 6 & 0x1))
1741       {
1742         iglobCol = globCol;
1743         iglobRow = globRow;
1744         nTrigChannelCut++;
1745         //cout << "ts cut = " << ts << endl;
1746         //cout << "globCol = " << globCol << endl;
1747         //cout << "globRow = " << globRow << endl;
1748         ftriggersCut[globCol][globRow] = 1;
1749       }
1750
1751     } // calo trigger entries
1752   } // has calo trigger entries
1753
1754   // part 2 go through the clusters here -----------------------------------
1755   //cout << " part 2 go through the clusters here ----------------------------------- " << endl; 
1756   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1757   Short_t cellAddr, nSACell;
1758   Int_t mclabel;
1759   //Int_t nSACell, iSACell, mclabel;
1760   Int_t iSACell;
1761   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
1762   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1763
1764   TRefArray *fCaloClusters = new TRefArray();
1765   fESD->GetEMCALClusters( fCaloClusters ); 
1766   nCluster = fCaloClusters->GetEntries();
1767
1768
1769   // save all cells times if there are clusters  
1770   if( nCluster > 0 ){
1771     // erase time map
1772     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){ 
1773       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1774         cellTime[i][j] = 0.;
1775         cellEnergy[i][j] = 0.;
1776       }
1777     }
1778
1779     // get cells
1780     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1781     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1782     nSACell = fCaloCells->GetNumberOfCells();
1783     for(iSACell = 0; iSACell < nSACell; iSACell++ ){ 
1784       // get the cell info *fCal
1785       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1786
1787       // get cell position 
1788       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta ); 
1789       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1790
1791       // convert co global phi eta  
1792       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1793       geta = ieta + nColsFeeModule*(nSupMod%2);
1794
1795       // save cell time and energy
1796       cellTime[geta][gphi] = cellTimeT;
1797       cellEnergy[geta][gphi] = cellAmp;
1798
1799     }
1800   }
1801
1802   Int_t nClusterTrig, nClusterTrigCut;
1803   UShort_t *cellAddrs;
1804   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1805   Float_t clsPos[3] = {0.,0.,0.};
1806
1807   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1808   {
1809     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1810     if(!cluster || !cluster->IsEMCAL()) continue;
1811
1812     // get cluster cells
1813     nCell = cluster->GetNCells();
1814
1815     // get cluster energy
1816     clsE = cluster->E();
1817
1818     // get cluster position
1819     cluster->GetPosition(clsPos);
1820     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1821     clsEta = clsPosVec.Eta();
1822     clsPhi = clsPosVec.Phi();
1823
1824     // get the cell addresses
1825     cellAddrs = cluster->GetCellsAbsId();
1826
1827     // check if the cluster contains cell, that was marked as triggered
1828     nClusterTrig = 0;
1829     nClusterTrigCut = 0;
1830
1831     // loop the cells to check, if cluser in acceptance
1832     // any cluster with a cell outside acceptance is not considered
1833     for( iCell = 0; iCell < nCell; iCell++ )
1834     {
1835      // check hot cell
1836      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]); 
1837
1838       // get cell position
1839       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1840       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1841
1842       // convert co global phi eta
1843       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1844       geta = ieta + nColsFeeModule*(nSupMod%2);
1845
1846       if( cellTime[geta][gphi] > 0. ){ 
1847         clusterTime += cellTime[geta][gphi];
1848         gCell++;
1849       }
1850
1851       // get corresponding FALTRO
1852       fphi = gphi / 2;
1853       feta = geta / 2;
1854
1855         //cout << "fphi = " << fphi << endl;
1856         //cout << "feta = " << feta << endl;
1857
1858       // try to match with a triggered
1859       if( ftriggers[feta][fphi]==1)
1860       {  nClusterTrig++;
1861       }
1862       if( ftriggersCut[feta][fphi]==1)
1863       { nClusterTrigCut++;
1864       }
1865
1866       //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
1867      
1868     } // cells
1869
1870
1871     if( gCell > 0 ) 
1872       clusterTime = clusterTime / (Double_t)gCell;
1873     // fix the reconstriction code time 100ns jumps
1874     if( fESD->GetBunchCrossNumber() % 4 < 2 )
1875       clusterTime -= 0.0000001;
1876
1877     //fClsETime->Fill(clsE,clusterTime);
1878     //fClsEBftTrigCut->Fill(clsE);
1879
1880     if(nClusterTrig>0){
1881       //fClsETime1->Fill(clsE,clusterTime);
1882     }
1883
1884     if(nClusterTrig>0){
1885       cluster->SetChi2(1);
1886       //fClsEAftTrigCut1->Fill(clsE);                                               
1887     }
1888
1889     if(nClusterTrigCut>0){
1890       cluster->SetChi2(2);
1891       //fClsEAftTrigCut2->Fill(clsE);
1892     }
1893
1894     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1895     {
1896       cluster->SetChi2(3);
1897       //fClsEAftTrigCut3->Fill(clsE);
1898     }
1899
1900     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1901     {
1902       // cluster->SetChi2(4);
1903       //fClsEAftTrigCut4->Fill(clsE);
1904     }
1905     if(nClusterTrigCut<1)
1906     {
1907       cluster->SetChi2(0);
1908
1909       //fClsEAftTrigCut->Fill(clsE);
1910     }
1911
1912   } // clusters
1913 }
1914
1915 // <-------- only MC correction
1916 double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
1917 {
1918   TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)"); 
1919   TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)"); 
1920
1921  double shift = 0.0;
1922   
1923  if(central>0 && central<=10)
1924    {
1925     fcorr0->SetParameters(1.045,1.288,3.18e-01); //
1926     fcorr1->SetParameters(9.91e-01,3.466,2.344);
1927    }
1928  else if(central>10 && central<=20)
1929    {
1930     fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
1931     fcorr1->SetParameters(0.975,2.276,1.501e-01);
1932    }
1933  else if(central>20 && central<=30)
1934    {
1935     fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
1936     fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
1937    }
1938  else if(central>30 && central<=40)
1939    {
1940     fcorr0->SetParameters(1.00,1.466,2.305e-1);
1941     fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
1942    }
1943  else if(central>40 && central<=50)
1944    {
1945     fcorr0->SetParameters(1.00,1.422,1.518e-01);
1946     fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
1947    }
1948  
1949  else if(central>50 && central<=70)
1950    {
1951     fcorr0->SetParameters(0.989,2.495,2.167);
1952     fcorr1->SetParameters(0.961,1.734,1.438e-01);
1953    }
1954  else if(central>70 && central<=100)
1955    {
1956     fcorr0->SetParameters(0.981,-3.373,3.93327);
1957     fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
1958    }
1959  
1960
1961  shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
1962
1963  return shift;
1964 }
1965
1966 // <-------- only Data correction
1967 double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central)
1968 {
1969  int icent = 0;
1970
1971  if(central>=0 && central<10)
1972    {
1973     icent = 0;
1974    }
1975  else if(central>=10 && central<20)
1976   {
1977    icent = 1;
1978   }
1979  else if(central>=20 && central<30)
1980   {
1981    icent = 2;
1982   }
1983  else if(central>=30 && central<40)
1984   {
1985    icent = 3;
1986   }
1987  else if(central>=40 && central<50)
1988   {
1989    icent = 4;
1990   }
1991  else if(central>=50 && central<70)
1992   {
1993    icent = 5;
1994   }
1995  else
1996   {
1997    icent = 6;
1998   }
1999
2000  double shift = fnSigEtaCorr[icent]->Eval(tmpeta);
2001  
2002  //cout << "eta correction"<< endl;
2003  //cout << "cent = "<< central<< endl;
2004  //cout << "icent = "<< icent << endl;
2005  //cout << "shift = "<< shift << endl;
2006
2007  return shift;
2008
2009 }
2010
2011
2012 double AliAnalysisTaskHFECal::SumpT(Int_t itrack, AliESDtrack* track)
2013 {
2014  
2015   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
2016   fTrackCuts->SetRequireTPCRefit(kTRUE);
2017   fTrackCuts->SetRequireITSRefit(kTRUE);
2018   fTrackCuts->SetEtaRange(-0.9,0.9);
2019   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
2020   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
2021   fTrackCuts->SetMinNClustersTPC(90);
2022  
2023   double pTrecp = track->Pt();
2024   double phiorg = track->Phi();
2025   double etaorg = track->Eta();
2026
2027   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
2028     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
2029     if (!trackAsso) {
2030       printf("ERROR: Could not receive track %d\n", jTracks);
2031       continue;
2032     }
2033     if(itrack==jTracks)continue;
2034     double pTAss = trackAsso->Pt();
2035     double etaAss = trackAsso->Eta();
2036     double phiAss = trackAsso->Phi();
2037
2038     double delphi = phiorg - phiAss;
2039     double deleta = etaorg - etaAss;
2040
2041     double R = sqrt(pow(deleta,2)+pow(delphi,2));
2042     if(pTAss<0.5)continue;
2043     if(R<0.4)pTrecp+=pTAss;
2044
2045     }
2046  
2047    return pTrecp;
2048 }
2049