update
[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 "TMath.h"
24 #include "TCanvas.h"
25 #include "THnSparse.h"
26 #include "TLorentzVector.h"
27 #include "TString.h"
28 #include "TFile.h"
29
30 #include "AliAnalysisTask.h"
31 #include "AliAnalysisManager.h"
32
33 #include "AliESDEvent.h"
34 #include "AliESDHandler.h"
35 #include "AliAODEvent.h"
36 #include "AliAODHandler.h"
37
38 #include "AliMCEventHandler.h"
39 #include "AliMCEvent.h"
40 #include "AliMCParticle.h"
41
42 #include "AliAnalysisTaskHFECal.h"
43 #include "TGeoGlobalMagField.h"
44 #include "AliLog.h"
45 #include "AliAnalysisTaskSE.h"
46 #include "TRefArray.h"
47 #include "TVector.h"
48 #include "AliESDInputHandler.h"
49 #include "AliESDpid.h"
50 #include "AliESDtrackCuts.h"
51 #include "AliPhysicsSelection.h"
52 #include "AliESDCaloCluster.h"
53 #include "AliAODCaloCluster.h"
54 #include "AliEMCALRecoUtils.h"
55 #include "AliEMCALGeometry.h"
56 #include "AliGeomManager.h"
57 #include "stdio.h"
58 #include "TGeoManager.h"
59 #include "iostream"
60 #include "fstream"
61
62 #include "AliEMCALTrack.h"
63 #include "AliMagF.h"
64
65 #include "AliKFParticle.h"
66 #include "AliKFVertex.h"
67
68 #include "AliPID.h"
69 #include "AliPIDResponse.h"
70 #include "AliHFEcontainer.h"
71 #include "AliHFEcuts.h"
72 #include "AliHFEpid.h"
73 #include "AliHFEpidBase.h"
74 #include "AliHFEpidQAmanager.h"
75 #include "AliHFEtools.h"
76 #include "AliCFContainer.h"
77 #include "AliCFManager.h"
78
79 #include "AliStack.h"
80
81 #include "AliCentrality.h"
82
83 ClassImp(AliAnalysisTaskHFECal)
84 //________________________________________________________________________
85 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name) 
86   : AliAnalysisTaskSE(name)
87   ,fESD(0)
88   ,fMC(0)
89   ,stack(0)
90   ,fGeom(0)
91   ,fOutputList(0)
92   ,fqahist(0) 
93   ,fTrackCuts(0)
94   ,fCuts(0)
95   ,fIdentifiedAsOutInz(kFALSE)
96   ,fPassTheEventCut(kFALSE)
97   ,fRejectKinkMother(kFALSE)
98   ,fmcData(kFALSE)
99   ,fVz(0.0)
100   ,fCFM(0)      
101   ,fPID(0)
102   ,fPIDqa(0)           
103   ,fOpeningAngleCut(0.1)
104   ,fInvmassCut(0.01)    
105   ,fNoEvents(0)
106   ,fEMCAccE(0)
107   ,fTrkpt(0)
108   ,fTrkEovPBef(0)        
109   ,fTrkEovPAft(0)       
110   ,fdEdxBef(0)   
111   ,fdEdxAft(0)   
112   ,fIncpT(0)    
113   ,fIncpTM20(0) 
114   ,fInvmassLS(0)                
115   ,fInvmassULS(0)               
116   ,fInvmassLSmc(0)              
117   ,fInvmassULSmc(0)             
118   ,fInvmassLSmc0(0)             
119   ,fInvmassLSmc1(0)             
120   ,fInvmassLSmc2(0)             
121   ,fInvmassLSmc3(0)             
122   ,fInvmassULSmc0(0)            
123   ,fInvmassULSmc1(0)            
124   ,fInvmassULSmc2(0)            
125   ,fInvmassULSmc3(0)            
126   ,fOpeningAngleLS(0)   
127   ,fOpeningAngleULS(0)  
128   ,fPhotoElecPt(0)
129   ,fPhoElecPt(0)
130   ,fPhoElecPtM20(0)
131   ,fSameElecPt(0)
132   ,fSameElecPtM20(0)
133   ,fTrackPtBefTrkCuts(0)         
134   ,fTrackPtAftTrkCuts(0)
135   ,fTPCnsigma(0)
136   ,fCent(0)
137   ,fEleInfo(0)
138   /*
139   ,fClsEBftTrigCut(0)    
140   ,fClsEAftTrigCut(0)    
141   ,fClsEAftTrigCut1(0)   
142   ,fClsEAftTrigCut2(0)   
143   ,fClsEAftTrigCut3(0)   
144   ,fClsEAftTrigCut4(0)   
145   ,fClsETime(0)       
146   ,fClsETime1(0)       
147   ,fTrigTimes(0)
148   ,fCellCheck(0)
149   */
150   ,fInputHFEMC(0)
151   ,fInputAlle(0)
152   ,fIncpTMChfe(0)       
153   ,fIncpTMChfeAll(0)    
154   ,fIncpTMCM20hfe(0)    
155   ,fIncpTMCM20hfeAll(0) 
156   ,fIncpTMCM20hfeCheck(0)       
157   ,fIncpTMCpho(0)       
158   ,fIncpTMCM20pho(0)    
159   ,fPhoElecPtMC(0)
160   ,fPhoElecPtMCM20(0)
161   ,fSameElecPtMC(0)
162   ,fSameElecPtMCM20(0)
163   ,fIncpTMCM20pho_pi0e(0)       
164   ,fPhoElecPtMCM20_pi0e(0)
165   ,fSameElecPtMCM20_pi0e(0)
166   ,fIncpTMCM20pho_eta(0)        
167   ,fPhoElecPtMCM20_eta(0)
168   ,fSameElecPtMCM20_eta(0)
169   ,CheckNclust(0)
170   ,CheckNits(0)
171   ,Hpi0pTcheck(0)
172   ,HETApTcheck(0)
173   ,HphopTcheck(0)
174   ,fpTCheck(0)
175   ,fMomDtoE(0) 
176 {
177   //Named constructor
178   
179   fPID = new AliHFEpid("hfePid");
180   fTrackCuts = new AliESDtrackCuts();
181   
182   // Define input and output slots here
183   // Input slot #0 works with a TChain
184   DefineInput(0, TChain::Class());
185   // Output slot #0 id reserved by the base class for AOD
186   // Output slot #1 writes into a TH1 container
187   // DefineOutput(1, TH1I::Class());
188   DefineOutput(1, TList::Class());
189   //  DefineOutput(3, TTree::Class());
190 }
191
192 //________________________________________________________________________
193 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() 
194   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
195   ,fESD(0)
196   ,fMC(0)
197   ,stack(0)
198   ,fGeom(0)
199   ,fOutputList(0)
200   ,fqahist(0)
201   ,fTrackCuts(0)
202   ,fCuts(0)
203   ,fIdentifiedAsOutInz(kFALSE)
204   ,fPassTheEventCut(kFALSE)
205   ,fRejectKinkMother(kFALSE)
206   ,fmcData(kFALSE)
207   ,fVz(0.0)
208   ,fCFM(0)      
209   ,fPID(0)       
210   ,fPIDqa(0)           
211   ,fOpeningAngleCut(0.1)
212   ,fInvmassCut(0.01)    
213   ,fNoEvents(0)
214   ,fEMCAccE(0)
215   ,fTrkpt(0)
216   ,fTrkEovPBef(0)        
217   ,fTrkEovPAft(0)        
218   ,fdEdxBef(0)   
219   ,fdEdxAft(0)   
220   ,fIncpT(0)     
221   ,fIncpTM20(0)  
222   ,fInvmassLS(0)                
223   ,fInvmassULS(0)               
224   ,fInvmassLSmc(0)              
225   ,fInvmassULSmc(0)             
226   ,fInvmassLSmc0(0)             
227   ,fInvmassLSmc1(0)             
228   ,fInvmassLSmc2(0)             
229   ,fInvmassLSmc3(0)             
230   ,fInvmassULSmc0(0)            
231   ,fInvmassULSmc1(0)            
232   ,fInvmassULSmc2(0)            
233   ,fInvmassULSmc3(0)            
234   ,fOpeningAngleLS(0)   
235   ,fOpeningAngleULS(0)  
236   ,fPhotoElecPt(0)
237   ,fPhoElecPt(0)
238   ,fPhoElecPtM20(0)
239   ,fSameElecPt(0)
240   ,fSameElecPtM20(0)
241   ,fTrackPtBefTrkCuts(0)         
242   ,fTrackPtAftTrkCuts(0)                  
243   ,fTPCnsigma(0)
244   ,fCent(0)
245   ,fEleInfo(0)
246   /*
247   ,fClsEBftTrigCut(0)    
248   ,fClsEAftTrigCut(0)    
249   ,fClsEAftTrigCut1(0)   
250   ,fClsEAftTrigCut2(0)   
251   ,fClsEAftTrigCut3(0)   
252   ,fClsEAftTrigCut4(0)   
253   ,fClsETime(0)       
254   ,fClsETime1(0)       
255   ,fTrigTimes(0)
256   ,fCellCheck(0)
257   */
258   ,fInputHFEMC(0)
259   ,fInputAlle(0)
260   ,fIncpTMChfe(0)       
261   ,fIncpTMChfeAll(0)    
262   ,fIncpTMCM20hfe(0)    
263   ,fIncpTMCM20hfeAll(0) 
264   ,fIncpTMCM20hfeCheck(0)       
265   ,fIncpTMCpho(0)       
266   ,fIncpTMCM20pho(0)    
267   ,fPhoElecPtMC(0)
268   ,fPhoElecPtMCM20(0)
269   ,fSameElecPtMC(0)
270   ,fSameElecPtMCM20(0)
271   ,fIncpTMCM20pho_pi0e(0)       
272   ,fPhoElecPtMCM20_pi0e(0)
273   ,fSameElecPtMCM20_pi0e(0)
274   ,fIncpTMCM20pho_eta(0)        
275   ,fPhoElecPtMCM20_eta(0)
276   ,fSameElecPtMCM20_eta(0)
277   ,CheckNclust(0)
278   ,CheckNits(0)
279   ,Hpi0pTcheck(0)
280   ,HETApTcheck(0)
281   ,HphopTcheck(0)
282   ,fpTCheck(0)
283   ,fMomDtoE(0)
284 {
285         //Default constructor
286         fPID = new AliHFEpid("hfePid");
287
288         fTrackCuts = new AliESDtrackCuts();
289         
290         // Constructor
291         // Define input and output slots here
292         // Input slot #0 works with a TChain
293         DefineInput(0, TChain::Class());
294         // Output slot #0 id reserved by the base class for AOD
295         // Output slot #1 writes into a TH1 container
296         // DefineOutput(1, TH1I::Class());
297         DefineOutput(1, TList::Class());
298         //DefineOutput(3, TTree::Class());
299 }
300 //_________________________________________
301
302 AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
303 {
304   //Destructor 
305   
306   delete fOutputList;
307   delete fGeom;
308   delete fPID;
309   delete fCuts;
310   delete fCFM;
311   delete fPIDqa;
312   delete fTrackCuts;
313 }
314 //_________________________________________
315
316 void AliAnalysisTaskHFECal::UserExec(Option_t*)
317 {
318   //Main loop
319   //Called for each event
320   
321   // create pointer to event
322   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
323   if (!fESD) {
324     printf("ERROR: fESD not available\n");
325     return;
326   }
327   
328   if(!fCuts){
329     AliError("HFE cuts not available");
330     return;
331   }
332   
333   if(!fPID->IsInitialized()){ 
334     // Initialize PID with the given run number
335     AliWarning("PID not initialised, get from Run no");
336     fPID->InitializePID(fESD->GetRunNumber());
337   }
338  
339   if(fmcData)fMC = MCEvent();
340   //AliStack* stack = NULL;
341   if(fmcData && fMC)stack = fMC->Stack();
342
343   Float_t cent = -1.;
344   AliCentrality *centrality = fESD->GetCentrality(); 
345   cent = centrality->GetCentralityPercentile("V0M");
346
347   //---- fill MC track info
348   if(fmcData && fMC)
349     {
350     Int_t nParticles = stack->GetNtrack();
351     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
352       TParticle* particle = stack->Particle(iParticle);
353       int fPDG = particle->GetPdgCode(); 
354       double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
355       double pTMC = particle->Pt();
356       double proR = particle->R();
357       double etaMC = particle->Eta();
358       if(fabs(etaMC)>0.6)continue;
359       Bool_t mcInDtoE= kFALSE;
360       Bool_t mcInBtoE= kFALSE;
361
362       Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
363       //if(!MChijing)printf("not MC hijing");
364       int iHijing = 1;
365       if(!MChijing)iHijing = 0;
366       if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
367       if(fPDG==221)HETApTcheck->Fill(pTMC,iHijing);
368
369       if(particle->GetFirstMother()>-1 && fPDG==22)
370         {
371          int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();  
372          if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
373         } 
374
375       if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
376         {
377             int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();  
378             if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
379             if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
380             if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)fInputHFEMC->Fill(cent,pTMC);
381          }
382
383          if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
384
385       }
386     } 
387
388   fNoEvents->Fill(0);
389
390   Int_t fNOtrks =  fESD->GetNumberOfTracks();
391   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
392   
393   Double_t pVtxZ = -999;
394   pVtxZ = pVtx->GetZ();
395   
396   if(TMath::Abs(pVtxZ)>10) return;
397   fNoEvents->Fill(1);
398   
399   if(fNOtrks<2) return;
400   fNoEvents->Fill(2);
401   
402   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
403   if(!pidResponse){
404     AliDebug(1, "Using default PID Response");
405     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); 
406   }
407   
408   fPID->SetPIDResponse(pidResponse);
409   
410   fCFM->SetRecEventInfo(fESD);
411   
412   //Float_t cent = -1.;
413   //AliCentrality *centrality = fESD->GetCentrality(); 
414   //cent = centrality->GetCentralityPercentile("V0M");
415   fCent->Fill(cent);
416   
417   //if(cent>90.) return;
418         
419  // Calorimeter info.
420  
421    FindTriggerClusters();
422
423   // make EMCAL array 
424   for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
425      {
426       AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
427       if(clust->IsEMCAL())
428         {
429          double clustE = clust->E();
430          float  emcx[3]; // cluster pos
431          clust->GetPosition(emcx);
432          TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
433          double emcphi = clustpos.Phi(); 
434          double emceta = clustpos.Eta();
435          double calInfo[5];
436          calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2(); 
437          //if(clustE>1.5)fEMCAccE->Fill(calInfo); 
438          //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo); 
439         }
440    }
441
442   // Track loop 
443   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
444     AliESDtrack* track = fESD->GetTrack(iTracks);
445     if (!track) {
446       printf("ERROR: Could not receive track %d\n", iTracks);
447       continue;
448     }
449    
450     int parentlabel = 99999;
451     int parentPID = 99999;
452     int grand_parentlabel = 99999;
453     int grand_parentPID = 99999;
454     Bool_t mcPho = kFALSE;
455     Bool_t mcDtoE= kFALSE;
456     Bool_t mcBtoE= kFALSE;
457     Bool_t mcOrgPi0 = kFALSE;
458     Bool_t mcOrgEta = kFALSE;
459     double mcele = -1.;
460     double mcpT = 0.0;
461     double mcMompT = 0.0;
462     //double mcGrandMompT = 0.0;
463     double mcWeight = -10.0;
464
465     int iHijing = 1;
466
467     if(fmcData && fMC && stack)
468       {
469        Int_t label = TMath::Abs(track->GetLabel());
470
471        Bool_t MChijing = fMC->IsFromBGEvent(label);
472        if(!MChijing)iHijing = 0;
473
474        TParticle* particle = stack->Particle(label);
475        int mcpid = particle->GetPdgCode();
476        mcpT = particle->Pt();
477        //printf("MCpid = %d",mcpid);
478        if(particle->GetFirstMother()>-1)
479          {
480           //int parentlabel = particle->GetFirstMother();
481           //int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
482           //mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
483           FindMother(particle, parentlabel, parentPID);
484           mcMompT = stack->Particle(parentlabel)->Pt();
485           if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
486           if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
487           if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
488
489           // make D->e pT correlation
490           if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT); 
491
492            //cout << "check PID = " << parentPID << endl;
493            //cout << "check pho = " << mcPho << endl;
494            //cout << "check D or B = " << mcDtoE << endl;
495           // pi->e (Dalitz)
496           if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
497               {
498                //cout << "find pi0->e " <<  endl;
499                mcOrgPi0 = kTRUE;
500                mcWeight = GetMCweight(mcMompT); 
501               }
502           // eta->e (Dalitz)
503           if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
504               {
505                //cout << "find Eta->e " <<  endl;
506                mcOrgEta = kTRUE;
507                mcWeight = GetMCweightEta(mcMompT); 
508               }
509
510           // access grand parent 
511           TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
512           //if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
513           if(particle_parent->GetFirstMother()>-1 && (parentPID==22 || parentPID==111) && fabs(mcpid)==11) // get grand parent g->e
514             {
515              //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
516              //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
517              FindMother(particle_parent, grand_parentlabel, grand_parentPID);
518              double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
519              if(grand_parentPID==111 && mcGrandpT>0.0)
520                 {
521                   // check eta->pi0 decay !
522                    int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
523                    TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
524                    FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
525                    if(grand2_parentPID==221)
526                      {
527                       //cout << "find Eta->e " <<  endl;
528                       double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
529                       mcOrgEta = kTRUE;
530                       mcWeight = GetMCweight(mcGrandpT2);  
531                       mcMompT = mcGrandpT2; 
532                      }
533                    else
534                      {
535                       //cout << "find pi0->e " <<  endl;
536                       mcOrgPi0 = kTRUE;
537                       mcWeight = GetMCweight(mcGrandpT);  
538                       mcMompT = mcGrandpT; 
539                      }
540                 }
541
542              if(grand_parentPID==221 && mcGrandpT>0.0)
543                 {
544                    //cout << "find Eta->e " <<  endl;
545                    mcOrgEta = kTRUE;
546                    mcOrgPi0 = kFALSE;
547                    mcWeight = GetMCweightEta(mcGrandpT); 
548                    mcMompT = mcGrandpT; 
549                 }
550             }
551          } 
552        if(fabs(mcpid)==11 && mcDtoE)mcele= 1.; 
553        if(fabs(mcpid)==11 && mcBtoE)mcele= 2.; 
554        if(fabs(mcpid)==11 && mcPho)mcele= 3.; 
555       }
556  
557     //cout << "Pi0 = " << mcOrgPi0 << " ; Eta = " << mcOrgEta << endl; 
558     //printf("weight = %f\n",mcWeight);
559
560     if(TMath::Abs(track->Eta())>0.6) continue;
561     if(TMath::Abs(track->Pt()<2.0)) continue;
562     
563     fTrackPtBefTrkCuts->Fill(track->Pt());              
564     // RecKine: ITSTPC cuts  
565     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
566     
567     //RecKink
568     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
569       if(track->GetKinkIndex(0) != 0) continue;
570     } 
571     
572     // RecPrim
573     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
574     
575     // HFEcuts: ITS layers cuts
576     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
577     
578     // HFE cuts: TPC PID cleanup
579     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
580
581     int nTPCcl = track->GetTPCNcls();
582     int nTPCclF = track->GetTPCNclsF();
583     int nITS = track->GetNcls(0);
584     
585     fTrackPtAftTrkCuts->Fill(track->Pt());              
586     
587     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
588     pt = track->Pt();
589     if(pt<2.0)continue;
590     
591     // Track extrapolation
592     
593     Int_t charge = track->Charge();
594     fTrkpt->Fill(pt);
595     mom = track->P();
596     phi = track->Phi();
597     eta = track->Eta();
598     dEdx = track->GetTPCsignal();
599     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
600     
601         double ncells = -1.0;
602         double m20 = -1.0;
603         double m02 = -1.0;
604         double disp = -1.0;
605         double rmatch = -1.0;  
606         double nmatch = -1.0;
607         double oppstatus = 0.0;
608
609     Bool_t fFlagPhotonicElec = kFALSE;
610     Bool_t fFlagConvinatElec = kFALSE;
611
612     Int_t clsId = track->GetEMCALcluster();
613     if (clsId>0){
614       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
615       if(clust && clust->IsEMCAL()){
616
617                 double clustE = clust->E();
618                 eop = clustE/fabs(mom);
619                 //double clustT = clust->GetTOF();
620                 ncells = clust->GetNCells();
621                 m02 = clust->GetM02();
622                 m20 = clust->GetM20();
623                 disp = clust->GetDispersion();
624                 double delphi = clust->GetTrackDx(); 
625                 double deleta = clust->GetTrackDz(); 
626                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
627                 nmatch = clust->GetNTracksMatched();
628
629                 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
630                 {
631                   SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
632                 }
633                 if(fFlagPhotonicElec)oppstatus = 1.0;
634                 if(fFlagConvinatElec)oppstatus = 2.0;
635                 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
636
637                   double valdedx[16];
638                   valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
639                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
640                   valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
641                   valdedx[15] = mcele;
642                   if(fqahist==1)fEleInfo->Fill(valdedx);
643                  
644
645       }
646     }
647      
648     if(nITS<2.5)continue;
649     if(nTPCcl<100)continue;
650    
651     CheckNclust->Fill(nTPCcl); 
652     CheckNits->Fill(nITS); 
653
654     fdEdxBef->Fill(mom,fTPCnSigma);
655     fTPCnsigma->Fill(mom,fTPCnSigma);
656     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
657
658     Int_t pidpassed = 1;
659     
660
661     //--- track accepted
662     AliHFEpidObject hfetrack;
663     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
664     hfetrack.SetRecTrack(track);
665     double binct = 10.5;
666     if((0.0< cent) && (cent<5.0)) binct = 0.5;
667     if((5.0< cent) && (cent<10.0)) binct = 1.5;
668     if((10.0< cent) && (cent<20.0)) binct = 2.5;
669     if((20.0< cent) && (cent<30.0)) binct = 3.5;
670     if((30.0< cent) && (cent<40.0)) binct = 4.5;
671     if((40.0< cent) && (cent<50.0)) binct = 5.5;
672     if((50.0< cent) && (cent<60.0)) binct = 6.5;
673     if((60.0< cent) && (cent<70.0)) binct = 7.5;
674     if((70.0< cent) && (cent<80.0)) binct = 8.5;
675     if((80.0< cent) && (cent<90.0)) binct = 9.5;
676     if((90.0< cent) && (cent<100.0)) binct = 10.5;
677
678     hfetrack.SetCentrality((int)binct); //added
679     hfetrack.SetPbPb();
680     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
681
682     if(pidpassed==0) continue;
683     
684     fTrkEovPAft->Fill(pt,eop);
685     fdEdxAft->Fill(mom,fTPCnSigma);
686
687     // Fill real data
688     fIncpT->Fill(cent,pt);    
689     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
690     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
691
692     if(m20>0.0 && m20<0.3)
693       { 
694        fIncpTM20->Fill(cent,pt);    
695        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
696        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
697      }
698     
699     // MC
700     if(mcele>0) // select MC electrons
701       {
702
703           fIncpTMChfeAll->Fill(cent,pt);    
704           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);    
705
706        if(mcBtoE || mcDtoE) // select B->e & D->e
707          {
708           fIncpTMChfe->Fill(cent,pt);    
709           if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);    
710           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT);    
711          }
712       
713        if(mcPho) // select photonic electrons
714         {
715          double phoval[5];
716          phoval[0] = cent;
717          phoval[1] = pt;
718          phoval[2] = fTPCnSigma;
719          phoval[3] = iHijing;
720          phoval[4] = mcMompT;
721
722          fIncpTMCpho->Fill(phoval);    
723          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
724          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
725
726          if(m20>0.0 && m20<0.3) 
727            {
728             fIncpTMCM20pho->Fill(phoval);    
729             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
730             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
731             // pi0->g->e
732             if(mcWeight>-1)
733               {
734                if(iHijing==1)mcWeight = 1.0; 
735                if(mcOrgPi0)
736                  {
737                   //fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);    
738                   //if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
739                   //if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
740                   fIncpTMCM20pho_pi0e->Fill(phoval);    
741                   if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval);
742                   if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval);
743                  }
744                if(mcOrgEta)
745                  {
746                   //fIncpTMCM20pho_eta->Fill(phoval,mcWeight);    
747                   //if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
748                   //if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
749                   fIncpTMCM20pho_eta->Fill(phoval);    
750                   if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval);
751                   if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval);
752                  }
753                // --- eta
754               }
755            }
756         } 
757       } 
758  }
759  PostData(1, fOutputList);
760 }
761 //_________________________________________
762 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
763 {
764   //--- Check MC
765  
766   //Bool_t mcData = kFALSE;
767   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
768     {
769      fmcData = kTRUE;
770      printf("+++++ MC Data available");
771     }
772   if(fmcData)
773     {
774      printf("++++++++= MC analysis \n");
775     }
776   else
777    {
778      printf("++++++++ real data analysis \n");
779    }
780
781    printf("+++++++ QA hist %d",fqahist);
782
783   //---- Geometry
784   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
785
786   //--------Initialize PID
787   //fPID->SetHasMCData(kFALSE);
788   fPID->SetHasMCData(fmcData);
789   if(!fPID->GetNumberOfPIDdetectors()) 
790     {
791       fPID->AddDetector("TPC", 0);
792       fPID->AddDetector("EMCAL", 1);
793     }
794   
795   Double_t params[4];
796   const char *cutmodel;
797   cutmodel = "pol0";
798   params[0] = -1.0; //sigma min
799   double maxnSig = 3.0;
800   if(fmcData)
801     {
802      params[0] = -5.0; //sigma min
803      maxnSig = 5.0; 
804     } 
805
806   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
807
808   fPID->SortDetectors(); 
809   fPIDqa = new AliHFEpidQAmanager();
810   fPIDqa->Initialize(fPID);
811  
812   //------- fcut --------------  
813   fCuts = new AliHFEcuts();
814   fCuts->CreateStandardCuts();
815   //fCuts->SetMinNClustersTPC(100);
816   fCuts->SetMinNClustersTPC(90);
817   fCuts->SetMinRatioTPCclusters(0.6);
818   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
819   //fCuts->SetMinNClustersITS(3);
820   fCuts->SetMinNClustersITS(2);
821   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
822   fCuts->SetCheckITSLayerStatus(kFALSE);
823   fCuts->SetVertexRange(10.);
824   fCuts->SetTOFPIDStep(kFALSE);
825   fCuts->SetPtRange(2, 50);
826   fCuts->SetMaxImpactParam(3.,3.);
827
828   //--------Initialize correction Framework and Cuts
829   fCFM = new AliCFManager;
830   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
831   fCFM->SetNStepParticle(kNcutSteps);
832   for(Int_t istep = 0; istep < kNcutSteps; istep++)
833     fCFM->SetParticleCutsList(istep, NULL);
834   
835   if(!fCuts){
836     AliWarning("Cuts not available. Default cuts will be used");
837     fCuts = new AliHFEcuts;
838     fCuts->CreateStandardCuts();
839   }
840   fCuts->Initialize(fCFM);
841   
842   //---------Output Tlist
843   fOutputList = new TList();
844   fOutputList->SetOwner();
845   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
846   
847   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
848   fOutputList->Add(fNoEvents);
849   
850   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
851   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5}; 
852   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5}; 
853   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
854   if(fqahist==1)fOutputList->Add(fEMCAccE);
855
856   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
857   fOutputList->Add(fTrkpt);
858   
859   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
860   fOutputList->Add(fTrackPtBefTrkCuts);
861   
862   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
863   fOutputList->Add(fTrackPtAftTrkCuts);
864   
865   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
866   fOutputList->Add(fTPCnsigma);
867   
868   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
869   fOutputList->Add(fTrkEovPBef);
870   
871   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
872   fOutputList->Add(fTrkEovPAft);
873   
874   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
875   fOutputList->Add(fdEdxBef);
876   
877   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
878   fOutputList->Add(fdEdxAft);
879   
880   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
881   fOutputList->Add(fIncpT);
882
883   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
884   fOutputList->Add(fIncpTM20);
885   
886   Int_t nBinspho[9] =  { 200, 100, 500, 12,   50,    4, 200,   8, 100};
887   Double_t minpho[9] = {  0.,  0.,  0., -2.5,  0, -0.5,   0,-1.5,   0};   
888   Double_t maxpho[9] = {100., 50., 0.5, 3.5,   1,  3.5,   2, 6.5,  50};   
889
890   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);
891   if(fqahist==1)fOutputList->Add(fInvmassLS);
892   
893   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);
894   if(fqahist==1)fOutputList->Add(fInvmassULS);
895   
896   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);
897   if(fqahist==1)fOutputList->Add(fInvmassLSmc);
898   
899   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);
900   if(fqahist==1)fOutputList->Add(fInvmassULSmc);
901
902   fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
903   fInvmassLSmc0->Sumw2();
904   fOutputList->Add(fInvmassLSmc0);
905   
906   fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
907   fInvmassLSmc1->Sumw2();
908   fOutputList->Add(fInvmassLSmc1);
909
910   fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
911   fInvmassLSmc2->Sumw2();
912   fOutputList->Add(fInvmassLSmc2);
913
914   fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
915   fInvmassLSmc3->Sumw2();
916   fOutputList->Add(fInvmassLSmc3);
917
918   fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
919   fInvmassULSmc0->Sumw2();
920   fOutputList->Add(fInvmassULSmc0);
921
922   fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
923   fInvmassULSmc1->Sumw2();
924   fOutputList->Add(fInvmassULSmc1);
925
926   fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
927   fInvmassULSmc2->Sumw2();
928   fOutputList->Add(fInvmassULSmc2);
929
930   fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
931   fInvmassULSmc3->Sumw2();
932   fOutputList->Add(fInvmassULSmc3);
933
934   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
935   fOutputList->Add(fOpeningAngleLS);
936   
937   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
938   fOutputList->Add(fOpeningAngleULS);
939   
940   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
941   fOutputList->Add(fPhotoElecPt);
942   
943   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
944   fOutputList->Add(fPhoElecPt);
945   
946   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
947   fOutputList->Add(fPhoElecPtM20);
948
949   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
950   fOutputList->Add(fSameElecPt);
951
952   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
953   fOutputList->Add(fSameElecPtM20);
954
955   fCent = new TH1F("fCent","Centrality",200,0,100) ;
956   fOutputList->Add(fCent);
957  
958   // Make common binning
959   const Double_t kMinP = 0.;
960   const Double_t kMaxP = 50.;
961   //const Double_t kTPCSigMim = 40.;
962   //const Double_t kTPCSigMax = 140.;
963
964   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
965   Int_t nBins[16] =  {  250,    10,   60,    20,   100,  300,  50,   40,   200, 200,  250, 200,    3,    5, 100,    8};
966   Double_t min[16] = {kMinP,  -0.5, 1.0,  -1.0,  -6.0,    0,   0,    0,   0.0, 0.0,  0.0,   0, -1.5, -0.5,  80, -1.5};
967   Double_t max[16] = {kMaxP,   9.5, 4.0,   1.0,   4.0,  3.0, 0.1,   40,   200, 2.0, 50.0, 100,  1.5,  4.5, 180,  6.5};
968   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);
969   if(fqahist==1)fOutputList->Add(fEleInfo);
970
971   //<---  Trigger info
972   /*
973   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
974   fOutputList->Add(fClsEBftTrigCut);
975
976   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
977   fOutputList->Add(fClsEAftTrigCut);
978
979   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
980   fOutputList->Add(fClsEAftTrigCut1);
981
982   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
983   fOutputList->Add(fClsEAftTrigCut2);
984
985   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
986   fOutputList->Add(fClsEAftTrigCut3);
987
988   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
989   fOutputList->Add(fClsEAftTrigCut4);
990
991   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
992   fOutputList->Add(fClsETime);
993
994   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
995   fOutputList->Add(fClsETime1);
996
997   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
998   fOutputList->Add(fTrigTimes);
999
1000   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1001   fOutputList->Add(fCellCheck);
1002   */
1003   //<---------- MC
1004
1005   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1006   fOutputList->Add(fInputHFEMC);
1007
1008   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
1009   fOutputList->Add(fInputAlle);
1010
1011   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1012   fOutputList->Add(fIncpTMChfe);
1013
1014   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
1015   fOutputList->Add(fIncpTMChfeAll);
1016
1017   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1018   fOutputList->Add(fIncpTMCM20hfe);
1019
1020   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
1021   fOutputList->Add(fIncpTMCM20hfeAll);
1022
1023   fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1024   fOutputList->Add(fIncpTMCM20hfeCheck);
1025
1026   Int_t nBinspho2[5] =  { 200, 100,    7,    3, 700};
1027   Double_t minpho2[5] = {  0.,  0., -2.5, -0.5, 0.};   
1028   Double_t maxpho2[5] = {100., 50.,  4.5,  2.5, 70.};   
1029
1030   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
1031   fOutputList->Add(fIncpTMCpho);
1032
1033   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1034   fOutputList->Add(fIncpTMCM20pho);
1035
1036   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1037   fOutputList->Add(fPhoElecPtMC);
1038   
1039   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1040   fOutputList->Add(fPhoElecPtMCM20);
1041
1042   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1043   fOutputList->Add(fSameElecPtMC);
1044
1045   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1046   fOutputList->Add(fSameElecPtMCM20);
1047
1048   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1049   fIncpTMCM20pho_pi0e->Sumw2();
1050   fOutputList->Add(fIncpTMCM20pho_pi0e);
1051
1052   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1053   fPhoElecPtMCM20_pi0e->Sumw2();
1054   fOutputList->Add(fPhoElecPtMCM20_pi0e);
1055
1056   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1057   fSameElecPtMCM20_pi0e->Sumw2();
1058   fOutputList->Add(fSameElecPtMCM20_pi0e);
1059  //
1060   fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1061   fIncpTMCM20pho_eta->Sumw2();
1062   fOutputList->Add(fIncpTMCM20pho_eta);
1063
1064   fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1065   fPhoElecPtMCM20_eta->Sumw2();
1066   fOutputList->Add(fPhoElecPtMCM20_eta);
1067
1068   fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1069   fSameElecPtMCM20_eta->Sumw2();
1070   fOutputList->Add(fSameElecPtMCM20_eta);
1071
1072
1073   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1074   fOutputList->Add(CheckNclust);
1075
1076   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1077   fOutputList->Add(CheckNits);
1078
1079   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1080   fOutputList->Add(Hpi0pTcheck);
1081
1082   HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1083   fOutputList->Add(HETApTcheck);
1084
1085   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1086   fOutputList->Add(HphopTcheck);
1087
1088   fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1089   fOutputList->Add(fpTCheck);
1090
1091   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1092   fOutputList->Add(fMomDtoE);
1093
1094   PostData(1,fOutputList);
1095 }
1096
1097 //________________________________________________________________________
1098 void AliAnalysisTaskHFECal::Terminate(Option_t *)
1099 {
1100   // Info("Terminate");
1101         AliAnalysisTaskSE::Terminate();
1102 }
1103
1104 //________________________________________________________________________
1105 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1106 {
1107   // Check single track cuts for a given cut step
1108   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1109   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1110   return kTRUE;
1111 }
1112 //_________________________________________
1113 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
1114 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
1115 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)
1116 {
1117   //Identify non-heavy flavour electrons using Invariant mass method
1118   
1119   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1120   fTrackCuts->SetRequireTPCRefit(kTRUE);
1121   fTrackCuts->SetRequireITSRefit(kTRUE);
1122   fTrackCuts->SetEtaRange(-0.9,0.9);
1123   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1124   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1125   fTrackCuts->SetMinNClustersTPC(90);
1126   
1127   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1128   
1129   double ptEle = track->Pt();  //add
1130   if(ibgevent==0 && w > 0.0)
1131      {
1132       fpTCheck->Fill(ptEle,w);
1133      }
1134
1135   Bool_t flagPhotonicElec = kFALSE;
1136   Bool_t flagConvinatElec = kFALSE;
1137   
1138   int p1 = 0;
1139   if(mce==3)
1140      {
1141        Int_t label = TMath::Abs(track->GetLabel());
1142        TParticle* particle = stack->Particle(label);
1143        p1 = particle->GetFirstMother();
1144      }
1145
1146   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1147   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1148     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1149     if (!trackAsso) {
1150       printf("ERROR: Could not receive track %d\n", jTracks);
1151       continue;
1152     }
1153     if(itrack==jTracks)continue;    
1154     int jbgevent = 0;    
1155
1156     int p2 = 0;
1157     if(mce==3)
1158     {
1159       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1160       TParticle* particle2 = stack->Particle(label2);
1161       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1162       if(MChijing_ass)jbgevent =1;
1163       if(particle2->GetFirstMother()>-1)
1164          p2 = particle2->GetFirstMother();
1165     }
1166
1167     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1168     Double_t mass=999., width = -999;
1169     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1170     
1171     //ptPrim = track->Pt();
1172     ptPrim = ptEle;
1173
1174     dEdxAsso = trackAsso->GetTPCsignal();
1175     ptAsso = trackAsso->Pt();
1176     Int_t chargeAsso = trackAsso->Charge();
1177     Int_t charge = track->Charge();
1178     
1179
1180     if(ptAsso <0.5) continue;
1181     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1182     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1183     
1184     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1185     if(charge>0) fPDGe1 = -11;
1186     if(chargeAsso>0) fPDGe2 = -11;
1187  
1188     //printf("chargeAsso = %d\n",chargeAsso);
1189     //printf("charge = %d\n",charge);
1190     if(charge == chargeAsso) fFlagLS = kTRUE;
1191     if(charge != chargeAsso) fFlagULS = kTRUE;
1192     
1193     //printf("fFlagLS = %d\n",fFlagLS);
1194     //printf("fFlagULS = %d\n",fFlagULS);
1195     //printf("\n");
1196
1197     AliKFParticle ge1(*track, fPDGe1);
1198     AliKFParticle ge2(*trackAsso, fPDGe2);
1199     AliKFParticle recg(ge1, ge2);
1200     
1201     if(recg.GetNDF()<1) continue;
1202     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1203     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1204     
1205     //  for v5-03-70-AN comment
1206     AliKFVertex primV(*pVtx);
1207     primV += recg;
1208     recg.SetProductionVertex(primV);
1209     
1210     recg.SetMassConstraint(0,0.0001);
1211     
1212     openingAngle = ge1.GetAngle(ge2);
1213     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1214     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1215     
1216     
1217     recg.GetMass(mass,width);
1218     
1219     double ishower = 0;
1220     if(shower>0.0 && shower<0.3)ishower = 1;
1221
1222     double phoinfo[9];
1223     phoinfo[0] = cent;
1224     phoinfo[1] = ptPrim;
1225     phoinfo[2] = mass;
1226     phoinfo[3] = nSig;
1227     phoinfo[4] = openingAngle;
1228     phoinfo[5] = ishower;
1229     phoinfo[6] = ep;
1230     phoinfo[7] = mce;
1231     phoinfo[8] = ptAsso;
1232
1233     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1234     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1235     if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1236     if(fFlagULS && ibgevent==0 && jbgevent==0)
1237        {
1238          fInvmassULSmc->Fill(phoinfo,w);
1239        }
1240     //printf("fInvmassCut %f\n",fInvmassCut);
1241     //printf("openingAngle %f\n",fOpeningAngleCut);
1242
1243     if(openingAngle > fOpeningAngleCut) continue;
1244     
1245     // for real data  
1246     //printf("mce =%f\n",mce);
1247     if(mce<-0.5) // mce==-1. is real
1248        {
1249          //printf("Real data\n");
1250          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1251          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1252                flagPhotonicElec = kTRUE;
1253               }
1254          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1255                flagConvinatElec = kTRUE;
1256               }
1257         }
1258     // for MC data  
1259     else
1260        {
1261          //printf("MC data\n");
1262
1263          if(w>0.0)
1264            {
1265            //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1266            if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass,w);
1267            if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass,w);
1268            if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass,w);
1269            if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass,w);
1270            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass,w);
1271            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass,w);
1272            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass,w);
1273            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass,w);
1274           }
1275          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1276          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1277                flagPhotonicElec = kTRUE;
1278               }
1279          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1280                flagConvinatElec = kTRUE;
1281               }
1282         }
1283
1284   }
1285   fFlagPhotonicElec = flagPhotonicElec;
1286   fFlagConvinatElec = flagConvinatElec;
1287   
1288 }
1289 //-------------------------------------------
1290
1291 void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1292 {
1293  //int label = 99999;
1294  //int pid = 99999;
1295
1296  if(part->GetFirstMother()>-1)
1297    {
1298     label = part->GetFirstMother();
1299     pid = stack->Particle(label)->GetPdgCode();
1300    }
1301    cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
1302 }
1303
1304 double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1305 {
1306         double weight = 1.0;
1307
1308         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1309         {
1310                 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1311         }
1312         else
1313         {
1314                 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1315         }
1316   return weight;
1317 }
1318
1319 double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1320 {
1321   double weight = 1.0;
1322
1323   weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1324   return weight;
1325 }
1326
1327
1328 //_________________________________________
1329 void AliAnalysisTaskHFECal::FindTriggerClusters()
1330 {
1331   // constants
1332   const int nModuleCols = 2;
1333   const int nModuleRows = 5;
1334   const int nColsFeeModule = 48;
1335   const int nRowsFeeModule = 24;
1336   const int nColsFaltroModule = 24;
1337   const int nRowsFaltroModule = 12;
1338   //const int faltroWidthMax = 20;
1339
1340   // part 1, trigger extraction -------------------------------------
1341   Int_t globCol, globRow;
1342   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1343   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1344
1345   //Int_t trigtimes[faltroWidthMax];
1346   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1347   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1348   //Double_t fTrigCutLow = 6;
1349   //Double_t fTrigCutHigh = 10;
1350   Double_t fTimeCutLow =  469e-09;
1351   Double_t fTimeCutHigh = 715e-09;
1352
1353   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1354
1355   // erase trigger maps
1356   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1357   {
1358     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1359     {
1360       ftriggersCut[i][j] = 0;
1361       ftriggers[i][j] = 0;
1362       ftriggersTime[i][j] = 0;
1363     }
1364   }
1365
1366   Int_t iglobCol=0, iglobRow=0;
1367   // go through triggers
1368   if( fCaloTrigger->GetEntries() > 0 )
1369   {
1370     // needs reset
1371     fCaloTrigger->Reset();
1372     while( fCaloTrigger->Next() )
1373     {
1374       fCaloTrigger->GetPosition( globCol, globRow );
1375       fCaloTrigger->GetNL0Times( ntimes );
1376       /*
1377       // no L0s
1378       if( ntimes < 1 )   continue;
1379       // get precise timings
1380       fCaloTrigger->GetL0Times( trigtimes );
1381       trigInCut = 0;
1382       for(Int_t i = 0; i < ntimes; i++ )
1383       {
1384         // save the first trigger time in channel
1385         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1386           triggersTime[globCol][globRow] = trigtimes[i];
1387         //printf("trigger times: %d\n",trigtimes[i]);
1388         // check if in cut
1389         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1390           trigInCut = 1;
1391
1392         fTrigTimes->Fill(trigtimes[i]);
1393       }
1394       */
1395    
1396       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1397       Int_t bit = 0;
1398       fCaloTrigger->GetTriggerBits(bit);
1399
1400       Int_t ts = 0;
1401       fCaloTrigger->GetL1TimeSum(ts);
1402       if (ts > 0)ftriggers[globCol][globRow] = 1;
1403       // number of triggered channels in event
1404       nTrigChannel++;
1405       // ... inside cut
1406       if(ts>0 && (bit >> 6 & 0x1))
1407       {
1408         iglobCol = globCol;
1409         iglobRow = globRow;
1410         nTrigChannelCut++;
1411         ftriggersCut[globCol][globRow] = 1;
1412       }
1413
1414     } // calo trigger entries
1415   } // has calo trigger entries
1416
1417   // part 2 go through the clusters here -----------------------------------
1418   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1419   Short_t cellAddr, nSACell, mclabel;
1420   //Int_t nSACell, iSACell, mclabel;
1421   Int_t iSACell;
1422   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
1423   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1424
1425   TRefArray *fCaloClusters = new TRefArray();
1426   fESD->GetEMCALClusters( fCaloClusters ); 
1427   nCluster = fCaloClusters->GetEntries();
1428
1429
1430   // save all cells times if there are clusters  
1431   if( nCluster > 0 ){
1432     // erase time map
1433     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){ 
1434       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1435         cellTime[i][j] = 0.;
1436         cellEnergy[i][j] = 0.;
1437       }
1438     }
1439
1440     // get cells
1441     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1442     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1443     nSACell = fCaloCells->GetNumberOfCells();
1444     for(iSACell = 0; iSACell < nSACell; iSACell++ ){ 
1445       // get the cell info *fCal
1446       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1447
1448       // get cell position 
1449       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta ); 
1450       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1451
1452       // convert co global phi eta  
1453       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1454       geta = ieta + nColsFeeModule*(nSupMod%2);
1455
1456       // save cell time and energy
1457       cellTime[geta][gphi] = cellTimeT;
1458       cellEnergy[geta][gphi] = cellAmp;
1459
1460     }
1461   }
1462
1463   Int_t nClusterTrig, nClusterTrigCut;
1464   UShort_t *cellAddrs;
1465   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1466   Float_t clsPos[3] = {0.,0.,0.};
1467
1468   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1469   {
1470     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1471     if(!cluster || !cluster->IsEMCAL()) continue;
1472
1473     // get cluster cells
1474     nCell = cluster->GetNCells();
1475
1476     // get cluster energy
1477     clsE = cluster->E();
1478
1479     // get cluster position
1480     cluster->GetPosition(clsPos);
1481     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1482     clsEta = clsPosVec.Eta();
1483     clsPhi = clsPosVec.Phi();
1484
1485     // get the cell addresses
1486     cellAddrs = cluster->GetCellsAbsId();
1487
1488     // check if the cluster contains cell, that was marked as triggered
1489     nClusterTrig = 0;
1490     nClusterTrigCut = 0;
1491
1492     // loop the cells to check, if cluser in acceptance
1493     // any cluster with a cell outside acceptance is not considered
1494     for( iCell = 0; iCell < nCell; iCell++ )
1495     {
1496      // check hot cell
1497      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]); 
1498
1499       // get cell position
1500       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1501       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1502
1503       // convert co global phi eta
1504       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1505       geta = ieta + nColsFeeModule*(nSupMod%2);
1506
1507       if( cellTime[geta][gphi] > 0. ){ 
1508         clusterTime += cellTime[geta][gphi];
1509         gCell++;
1510       }
1511
1512       // get corresponding FALTRO
1513       fphi = gphi / 2;
1514       feta = geta / 2;
1515
1516       // try to match with a triggered
1517       if( ftriggers[feta][fphi]==1)
1518       {  nClusterTrig++;
1519       }
1520       if( ftriggersCut[feta][fphi]==1)
1521       { nClusterTrigCut++;
1522       }
1523
1524     } // cells
1525
1526
1527     if( gCell > 0 ) 
1528       clusterTime = clusterTime / (Double_t)gCell;
1529     // fix the reconstriction code time 100ns jumps
1530     if( fESD->GetBunchCrossNumber() % 4 < 2 )
1531       clusterTime -= 0.0000001;
1532
1533     //fClsETime->Fill(clsE,clusterTime);
1534     //fClsEBftTrigCut->Fill(clsE);
1535
1536     if(nClusterTrig>0){
1537       //fClsETime1->Fill(clsE,clusterTime);
1538     }
1539
1540     if(nClusterTrig>0){
1541       cluster->SetChi2(1);
1542       //fClsEAftTrigCut1->Fill(clsE);                                               
1543     }
1544
1545     if(nClusterTrigCut>0){
1546       cluster->SetChi2(2);
1547       //fClsEAftTrigCut2->Fill(clsE);
1548     }
1549
1550     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1551     {
1552       cluster->SetChi2(3);
1553       //fClsEAftTrigCut3->Fill(clsE);
1554     }
1555
1556     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1557     {
1558       // cluster->SetChi2(4);
1559       //fClsEAftTrigCut4->Fill(clsE);
1560     }
1561     if(nClusterTrigCut<1)
1562     {
1563       cluster->SetChi2(0);
1564
1565       //fClsEAftTrigCut->Fill(clsE);
1566     }
1567
1568   } // clusters
1569 }
1570
1571
1572
1573
1574
1575
1576