7a1c203a69eee64891272eb24cd33361618869c5
[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 pi0->g->e
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             {
514              //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
515              //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
516              FindMother(particle_parent, grand_parentlabel, grand_parentPID);
517              double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
518              if(grand_parentPID==111 && mcGrandpT>0.0)
519                 {
520                   // check eta->pi0 decay !
521                    int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
522                    TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
523                    FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
524                    if(grand2_parentPID==221)
525                      {
526                       double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
527                       mcOrgEta = kTRUE;
528                       mcWeight = GetMCweight(mcGrandpT2);  
529                      }
530                    else
531                      {
532                       mcOrgPi0 = kTRUE;
533                       mcWeight = GetMCweight(mcGrandpT);  
534                      }
535                 }
536              if(grand_parentPID==221 && mcGrandpT>0.0)
537                 {
538                    mcOrgEta = kTRUE;
539                    mcWeight = GetMCweightEta(mcGrandpT); 
540                 }
541             }
542          } 
543        if(fabs(mcpid)==11 && mcDtoE)mcele= 1.; 
544        if(fabs(mcpid)==11 && mcBtoE)mcele= 2.; 
545        if(fabs(mcpid)==11 && mcPho)mcele= 3.; 
546       }
547  
548     //printf("weight = %f\n",mcWeight);
549
550     if(TMath::Abs(track->Eta())>0.6) continue;
551     if(TMath::Abs(track->Pt()<2.0)) continue;
552     
553     fTrackPtBefTrkCuts->Fill(track->Pt());              
554     // RecKine: ITSTPC cuts  
555     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
556     
557     //RecKink
558     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
559       if(track->GetKinkIndex(0) != 0) continue;
560     } 
561     
562     // RecPrim
563     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
564     
565     // HFEcuts: ITS layers cuts
566     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
567     
568     // HFE cuts: TPC PID cleanup
569     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
570
571     int nTPCcl = track->GetTPCNcls();
572     int nTPCclF = track->GetTPCNclsF();
573     int nITS = track->GetNcls(0);
574     
575     fTrackPtAftTrkCuts->Fill(track->Pt());              
576     
577     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
578     pt = track->Pt();
579     if(pt<2.0)continue;
580     
581     // Track extrapolation
582     
583     Int_t charge = track->Charge();
584     fTrkpt->Fill(pt);
585     mom = track->P();
586     phi = track->Phi();
587     eta = track->Eta();
588     dEdx = track->GetTPCsignal();
589     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
590     
591         double ncells = -1.0;
592         double m20 = -1.0;
593         double m02 = -1.0;
594         double disp = -1.0;
595         double rmatch = -1.0;  
596         double nmatch = -1.0;
597         double oppstatus = 0.0;
598
599     Bool_t fFlagPhotonicElec = kFALSE;
600     Bool_t fFlagConvinatElec = kFALSE;
601
602     Int_t clsId = track->GetEMCALcluster();
603     if (clsId>0){
604       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
605       if(clust && clust->IsEMCAL()){
606
607                 double clustE = clust->E();
608                 eop = clustE/fabs(mom);
609                 //double clustT = clust->GetTOF();
610                 ncells = clust->GetNCells();
611                 m02 = clust->GetM02();
612                 m20 = clust->GetM20();
613                 disp = clust->GetDispersion();
614                 double delphi = clust->GetTrackDx(); 
615                 double deleta = clust->GetTrackDz(); 
616                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
617                 nmatch = clust->GetNTracksMatched();
618
619                 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
620                 {
621                   SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing);
622                 }
623                 if(fFlagPhotonicElec)oppstatus = 1.0;
624                 if(fFlagConvinatElec)oppstatus = 2.0;
625                 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
626
627                   double valdedx[16];
628                   valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
629                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
630                   valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
631                   valdedx[15] = mcele;
632                   if(fqahist==1)fEleInfo->Fill(valdedx);
633                  
634
635       }
636     }
637      
638     if(nITS<2.5)continue;
639     if(nTPCcl<100)continue;
640    
641     CheckNclust->Fill(nTPCcl); 
642     CheckNits->Fill(nITS); 
643
644     fdEdxBef->Fill(mom,fTPCnSigma);
645     fTPCnsigma->Fill(mom,fTPCnSigma);
646     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
647
648     Int_t pidpassed = 1;
649     
650
651     //--- track accepted
652     AliHFEpidObject hfetrack;
653     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
654     hfetrack.SetRecTrack(track);
655     double binct = 10.5;
656     if((0.0< cent) && (cent<5.0)) binct = 0.5;
657     if((5.0< cent) && (cent<10.0)) binct = 1.5;
658     if((10.0< cent) && (cent<20.0)) binct = 2.5;
659     if((20.0< cent) && (cent<30.0)) binct = 3.5;
660     if((30.0< cent) && (cent<40.0)) binct = 4.5;
661     if((40.0< cent) && (cent<50.0)) binct = 5.5;
662     if((50.0< cent) && (cent<60.0)) binct = 6.5;
663     if((60.0< cent) && (cent<70.0)) binct = 7.5;
664     if((70.0< cent) && (cent<80.0)) binct = 8.5;
665     if((80.0< cent) && (cent<90.0)) binct = 9.5;
666     if((90.0< cent) && (cent<100.0)) binct = 10.5;
667
668     hfetrack.SetCentrality((int)binct); //added
669     hfetrack.SetPbPb();
670     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
671
672     if(pidpassed==0) continue;
673     
674     fTrkEovPAft->Fill(pt,eop);
675     fdEdxAft->Fill(mom,fTPCnSigma);
676
677     // Fill real data
678     fIncpT->Fill(cent,pt);    
679     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
680     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
681
682     if(m20>0.0 && m20<0.3)
683       { 
684        fIncpTM20->Fill(cent,pt);    
685        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
686        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
687      }
688     
689     // MC
690     if(mcele>0) // select MC electrons
691       {
692
693           fIncpTMChfeAll->Fill(cent,pt);    
694           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);    
695
696        if(mcBtoE || mcDtoE) // select B->e & D->e
697          {
698           fIncpTMChfe->Fill(cent,pt);    
699           if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);    
700           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT);    
701          }
702       
703        if(mcPho) // select photonic electrons
704         {
705          double phoval[5];
706          phoval[0] = cent;
707          phoval[1] = pt;
708          phoval[2] = fTPCnSigma;
709          phoval[3] = iHijing;
710          phoval[4] = mcMompT;
711
712          fIncpTMCpho->Fill(phoval);    
713          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
714          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
715
716          if(m20>0.0 && m20<0.3) 
717            {
718             fIncpTMCM20pho->Fill(phoval);    
719             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
720             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
721             // pi0->g->e
722             if(mcWeight>-1)
723               {
724                if(iHijing==1)mcWeight = 1.0; 
725                if(mcOrgPi0)
726                  {
727                   fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);    
728                   if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
729                   if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
730                  }
731                if(mcOrgEta)
732                  {
733                   fIncpTMCM20pho_eta->Fill(phoval,mcWeight);    
734                   if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
735                   if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
736                  }
737                // --- eta
738               }
739            }
740         } 
741       } 
742  }
743  PostData(1, fOutputList);
744 }
745 //_________________________________________
746 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
747 {
748   //--- Check MC
749  
750   //Bool_t mcData = kFALSE;
751   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
752     {
753      fmcData = kTRUE;
754      printf("+++++ MC Data available");
755     }
756   if(fmcData)
757     {
758      printf("++++++++= MC analysis \n");
759     }
760   else
761    {
762      printf("++++++++ real data analysis \n");
763    }
764
765    printf("+++++++ QA hist %d",fqahist);
766
767   //---- Geometry
768   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
769
770   //--------Initialize PID
771   //fPID->SetHasMCData(kFALSE);
772   fPID->SetHasMCData(fmcData);
773   if(!fPID->GetNumberOfPIDdetectors()) 
774     {
775       fPID->AddDetector("TPC", 0);
776       fPID->AddDetector("EMCAL", 1);
777     }
778   
779   Double_t params[4];
780   const char *cutmodel;
781   cutmodel = "pol0";
782   params[0] = -1.0; //sigma min
783   double maxnSig = 3.0;
784   if(fmcData)
785     {
786      params[0] = -5.0; //sigma min
787      maxnSig = 5.0; 
788     } 
789
790   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
791
792   fPID->SortDetectors(); 
793   fPIDqa = new AliHFEpidQAmanager();
794   fPIDqa->Initialize(fPID);
795  
796   //------- fcut --------------  
797   fCuts = new AliHFEcuts();
798   fCuts->CreateStandardCuts();
799   //fCuts->SetMinNClustersTPC(100);
800   fCuts->SetMinNClustersTPC(90);
801   fCuts->SetMinRatioTPCclusters(0.6);
802   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
803   //fCuts->SetMinNClustersITS(3);
804   fCuts->SetMinNClustersITS(2);
805   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
806   fCuts->SetCheckITSLayerStatus(kFALSE);
807   fCuts->SetVertexRange(10.);
808   fCuts->SetTOFPIDStep(kFALSE);
809   fCuts->SetPtRange(2, 50);
810   fCuts->SetMaxImpactParam(3.,3.);
811
812   //--------Initialize correction Framework and Cuts
813   fCFM = new AliCFManager;
814   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
815   fCFM->SetNStepParticle(kNcutSteps);
816   for(Int_t istep = 0; istep < kNcutSteps; istep++)
817     fCFM->SetParticleCutsList(istep, NULL);
818   
819   if(!fCuts){
820     AliWarning("Cuts not available. Default cuts will be used");
821     fCuts = new AliHFEcuts;
822     fCuts->CreateStandardCuts();
823   }
824   fCuts->Initialize(fCFM);
825   
826   //---------Output Tlist
827   fOutputList = new TList();
828   fOutputList->SetOwner();
829   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
830   
831   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
832   fOutputList->Add(fNoEvents);
833   
834   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
835   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5}; 
836   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5}; 
837   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
838   if(fqahist==1)fOutputList->Add(fEMCAccE);
839
840   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
841   fOutputList->Add(fTrkpt);
842   
843   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
844   fOutputList->Add(fTrackPtBefTrkCuts);
845   
846   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
847   fOutputList->Add(fTrackPtAftTrkCuts);
848   
849   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
850   fOutputList->Add(fTPCnsigma);
851   
852   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
853   fOutputList->Add(fTrkEovPBef);
854   
855   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
856   fOutputList->Add(fTrkEovPAft);
857   
858   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
859   fOutputList->Add(fdEdxBef);
860   
861   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
862   fOutputList->Add(fdEdxAft);
863   
864   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
865   fOutputList->Add(fIncpT);
866
867   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
868   fOutputList->Add(fIncpTM20);
869   
870   Int_t nBinspho[9] =  { 200, 100, 500, 12,   50,    4, 200,   8, 100};
871   Double_t minpho[9] = {  0.,  0.,  0., -2.5,  0, -0.5,   0,-1.5,   0};   
872   Double_t maxpho[9] = {100., 50., 0.5, 3.5,   1,  3.5,   2, 6.5,  50};   
873
874   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);
875   if(fqahist==1)fOutputList->Add(fInvmassLS);
876   
877   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);
878   if(fqahist==1)fOutputList->Add(fInvmassULS);
879   
880   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);
881   if(fqahist==1)fOutputList->Add(fInvmassLSmc);
882   
883   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);
884   if(fqahist==1)fOutputList->Add(fInvmassULSmc);
885
886   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 );
887   fInvmassLSmc0->Sumw2();
888   fOutputList->Add(fInvmassLSmc0);
889   
890   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 );
891   fInvmassLSmc1->Sumw2();
892   fOutputList->Add(fInvmassLSmc1);
893
894   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 );
895   fInvmassLSmc2->Sumw2();
896   fOutputList->Add(fInvmassLSmc2);
897
898   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 );
899   fInvmassLSmc3->Sumw2();
900   fOutputList->Add(fInvmassLSmc3);
901
902   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 );
903   fInvmassULSmc0->Sumw2();
904   fOutputList->Add(fInvmassULSmc0);
905
906   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 );
907   fInvmassULSmc1->Sumw2();
908   fOutputList->Add(fInvmassULSmc1);
909
910   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 );
911   fInvmassULSmc2->Sumw2();
912   fOutputList->Add(fInvmassULSmc2);
913
914   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 );
915   fInvmassULSmc3->Sumw2();
916   fOutputList->Add(fInvmassULSmc3);
917
918   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
919   fOutputList->Add(fOpeningAngleLS);
920   
921   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
922   fOutputList->Add(fOpeningAngleULS);
923   
924   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
925   fOutputList->Add(fPhotoElecPt);
926   
927   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
928   fOutputList->Add(fPhoElecPt);
929   
930   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
931   fOutputList->Add(fPhoElecPtM20);
932
933   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
934   fOutputList->Add(fSameElecPt);
935
936   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
937   fOutputList->Add(fSameElecPtM20);
938
939   fCent = new TH1F("fCent","Centrality",200,0,100) ;
940   fOutputList->Add(fCent);
941  
942   // Make common binning
943   const Double_t kMinP = 0.;
944   const Double_t kMaxP = 50.;
945   //const Double_t kTPCSigMim = 40.;
946   //const Double_t kTPCSigMax = 140.;
947
948   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
949   Int_t nBins[16] =  {  250,    10,   60,    20,   100,  300,  50,   40,   200, 200,  250, 200,    3,    5, 100,    8};
950   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};
951   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};
952   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);
953   if(fqahist==1)fOutputList->Add(fEleInfo);
954
955   //<---  Trigger info
956   /*
957   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
958   fOutputList->Add(fClsEBftTrigCut);
959
960   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
961   fOutputList->Add(fClsEAftTrigCut);
962
963   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
964   fOutputList->Add(fClsEAftTrigCut1);
965
966   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
967   fOutputList->Add(fClsEAftTrigCut2);
968
969   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
970   fOutputList->Add(fClsEAftTrigCut3);
971
972   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
973   fOutputList->Add(fClsEAftTrigCut4);
974
975   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
976   fOutputList->Add(fClsETime);
977
978   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
979   fOutputList->Add(fClsETime1);
980
981   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
982   fOutputList->Add(fTrigTimes);
983
984   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
985   fOutputList->Add(fCellCheck);
986   */
987   //<---------- MC
988
989   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
990   fOutputList->Add(fInputHFEMC);
991
992   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
993   fOutputList->Add(fInputAlle);
994
995   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
996   fOutputList->Add(fIncpTMChfe);
997
998   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
999   fOutputList->Add(fIncpTMChfeAll);
1000
1001   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1002   fOutputList->Add(fIncpTMCM20hfe);
1003
1004   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
1005   fOutputList->Add(fIncpTMCM20hfeAll);
1006
1007   fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1008   fOutputList->Add(fIncpTMCM20hfeCheck);
1009
1010   Int_t nBinspho2[5] =  { 200, 100,    7,    3, 100};
1011   Double_t minpho2[5] = {  0.,  0., -2.5, -0.5, 0.};   
1012   Double_t maxpho2[5] = {100., 50.,  4.5,  2.5, 50.};   
1013
1014   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
1015   fOutputList->Add(fIncpTMCpho);
1016
1017   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1018   fOutputList->Add(fIncpTMCM20pho);
1019
1020   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1021   fOutputList->Add(fPhoElecPtMC);
1022   
1023   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1024   fOutputList->Add(fPhoElecPtMCM20);
1025
1026   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1027   fOutputList->Add(fSameElecPtMC);
1028
1029   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1030   fOutputList->Add(fSameElecPtMCM20);
1031
1032   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1033   fIncpTMCM20pho_pi0e->Sumw2();
1034   fOutputList->Add(fIncpTMCM20pho_pi0e);
1035
1036   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1037   fPhoElecPtMCM20_pi0e->Sumw2();
1038   fOutputList->Add(fPhoElecPtMCM20_pi0e);
1039
1040   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1041   fSameElecPtMCM20_pi0e->Sumw2();
1042   fOutputList->Add(fSameElecPtMCM20_pi0e);
1043  //
1044   fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1045   fIncpTMCM20pho_eta->Sumw2();
1046   fOutputList->Add(fIncpTMCM20pho_eta);
1047
1048   fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1049   fPhoElecPtMCM20_eta->Sumw2();
1050   fOutputList->Add(fPhoElecPtMCM20_eta);
1051
1052   fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1053   fSameElecPtMCM20_eta->Sumw2();
1054   fOutputList->Add(fSameElecPtMCM20_eta);
1055
1056
1057   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1058   fOutputList->Add(CheckNclust);
1059
1060   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1061   fOutputList->Add(CheckNits);
1062
1063   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1064   fOutputList->Add(Hpi0pTcheck);
1065
1066   HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1067   fOutputList->Add(HETApTcheck);
1068
1069   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1070   fOutputList->Add(HphopTcheck);
1071
1072   fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1073   fOutputList->Add(fpTCheck);
1074
1075   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1076   fOutputList->Add(fMomDtoE);
1077
1078   PostData(1,fOutputList);
1079 }
1080
1081 //________________________________________________________________________
1082 void AliAnalysisTaskHFECal::Terminate(Option_t *)
1083 {
1084   // Info("Terminate");
1085         AliAnalysisTaskSE::Terminate();
1086 }
1087
1088 //________________________________________________________________________
1089 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1090 {
1091   // Check single track cuts for a given cut step
1092   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1093   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1094   return kTRUE;
1095 }
1096 //_________________________________________
1097 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
1098 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
1099 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)
1100 {
1101   //Identify non-heavy flavour electrons using Invariant mass method
1102   
1103   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1104   fTrackCuts->SetRequireTPCRefit(kTRUE);
1105   fTrackCuts->SetRequireITSRefit(kTRUE);
1106   fTrackCuts->SetEtaRange(-0.9,0.9);
1107   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1108   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1109   fTrackCuts->SetMinNClustersTPC(90);
1110   
1111   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1112   
1113   double ptEle = track->Pt();  //add
1114   if(ibgevent==0 && w > 0.0)
1115      {
1116       fpTCheck->Fill(ptEle,w);
1117      }
1118
1119   Bool_t flagPhotonicElec = kFALSE;
1120   Bool_t flagConvinatElec = kFALSE;
1121   
1122   int p1 = 0;
1123   if(mce==3)
1124      {
1125        Int_t label = TMath::Abs(track->GetLabel());
1126        TParticle* particle = stack->Particle(label);
1127        p1 = particle->GetFirstMother();
1128      }
1129
1130   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1131   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1132     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1133     if (!trackAsso) {
1134       printf("ERROR: Could not receive track %d\n", jTracks);
1135       continue;
1136     }
1137     if(itrack==jTracks)continue;    
1138     int jbgevent = 0;    
1139
1140     int p2 = 0;
1141     if(mce==3)
1142     {
1143       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1144       TParticle* particle2 = stack->Particle(label2);
1145       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1146       if(MChijing_ass)jbgevent =1;
1147       if(particle2->GetFirstMother()>-1)
1148          p2 = particle2->GetFirstMother();
1149     }
1150
1151     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1152     Double_t mass=999., width = -999;
1153     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1154     
1155     //ptPrim = track->Pt();
1156     ptPrim = ptEle;
1157
1158     dEdxAsso = trackAsso->GetTPCsignal();
1159     ptAsso = trackAsso->Pt();
1160     Int_t chargeAsso = trackAsso->Charge();
1161     Int_t charge = track->Charge();
1162     
1163
1164     if(ptAsso <0.5) continue;
1165     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1166     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1167     
1168     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1169     if(charge>0) fPDGe1 = -11;
1170     if(chargeAsso>0) fPDGe2 = -11;
1171  
1172     //printf("chargeAsso = %d\n",chargeAsso);
1173     //printf("charge = %d\n",charge);
1174     if(charge == chargeAsso) fFlagLS = kTRUE;
1175     if(charge != chargeAsso) fFlagULS = kTRUE;
1176     
1177     //printf("fFlagLS = %d\n",fFlagLS);
1178     //printf("fFlagULS = %d\n",fFlagULS);
1179     //printf("\n");
1180
1181     AliKFParticle ge1(*track, fPDGe1);
1182     AliKFParticle ge2(*trackAsso, fPDGe2);
1183     AliKFParticle recg(ge1, ge2);
1184     
1185     if(recg.GetNDF()<1) continue;
1186     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1187     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1188     
1189     //  for v5-03-70-AN comment
1190     AliKFVertex primV(*pVtx);
1191     primV += recg;
1192     recg.SetProductionVertex(primV);
1193     
1194     recg.SetMassConstraint(0,0.0001);
1195     
1196     openingAngle = ge1.GetAngle(ge2);
1197     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1198     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1199     
1200     
1201     recg.GetMass(mass,width);
1202     
1203     double ishower = 0;
1204     if(shower>0.0 && shower<0.3)ishower = 1;
1205
1206     double phoinfo[9];
1207     phoinfo[0] = cent;
1208     phoinfo[1] = ptPrim;
1209     phoinfo[2] = mass;
1210     phoinfo[3] = nSig;
1211     phoinfo[4] = openingAngle;
1212     phoinfo[5] = ishower;
1213     phoinfo[6] = ep;
1214     phoinfo[7] = mce;
1215     phoinfo[8] = ptAsso;
1216
1217     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1218     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1219     if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1220     if(fFlagULS && ibgevent==0 && jbgevent==0)
1221        {
1222          fInvmassULSmc->Fill(phoinfo,w);
1223        }
1224     //printf("fInvmassCut %f\n",fInvmassCut);
1225     //printf("openingAngle %f\n",fOpeningAngleCut);
1226
1227     if(openingAngle > fOpeningAngleCut) continue;
1228     
1229     // for real data  
1230     //printf("mce =%f\n",mce);
1231     if(mce<-0.5) // mce==-1. is real
1232        {
1233          //printf("Real data\n");
1234          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1235          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1236                flagPhotonicElec = kTRUE;
1237               }
1238          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1239                flagConvinatElec = kTRUE;
1240               }
1241         }
1242     // for MC data  
1243     else
1244        {
1245          //printf("MC data\n");
1246
1247          if(w>0.0)
1248            {
1249            if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc0->Fill(ptPrim,mass,w);
1250            if(fFlagULS && ibgevent==0 && jbgevent==0) fInvmassULSmc0->Fill(ptPrim,mass,w);
1251            if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc1->Fill(ptPrim,mass);
1252            if(fFlagULS && ibgevent==0 && jbgevent==0) fInvmassULSmc1->Fill(ptPrim,mass);
1253            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassLSmc2->Fill(ptPrim,mass,w);
1254            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassULSmc2->Fill(ptPrim,mass,w);
1255            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassLSmc3->Fill(ptPrim,mass);
1256            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassULSmc3->Fill(ptPrim,mass);
1257           }
1258          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1259          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1260                flagPhotonicElec = kTRUE;
1261               }
1262          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1263                flagConvinatElec = kTRUE;
1264               }
1265         }
1266
1267   }
1268   fFlagPhotonicElec = flagPhotonicElec;
1269   fFlagConvinatElec = flagConvinatElec;
1270   
1271 }
1272 //-------------------------------------------
1273
1274 void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1275 {
1276  //int label = 99999;
1277  //int pid = 99999;
1278
1279  if(part->GetFirstMother()>-1)
1280    {
1281     label = part->GetFirstMother();
1282     pid = stack->Particle(label)->GetPdgCode();
1283    }
1284    //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
1285 }
1286
1287 double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1288 {
1289         double weight = 1.0;
1290
1291         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1292         {
1293                 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1294         }
1295         else
1296         {
1297                 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1298         }
1299   return weight;
1300 }
1301
1302 double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1303 {
1304   double weight = 1.0;
1305
1306   weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1307   return weight;
1308 }
1309
1310
1311 //_________________________________________
1312 void AliAnalysisTaskHFECal::FindTriggerClusters()
1313 {
1314   // constants
1315   const int nModuleCols = 2;
1316   const int nModuleRows = 5;
1317   const int nColsFeeModule = 48;
1318   const int nRowsFeeModule = 24;
1319   const int nColsFaltroModule = 24;
1320   const int nRowsFaltroModule = 12;
1321   //const int faltroWidthMax = 20;
1322
1323   // part 1, trigger extraction -------------------------------------
1324   Int_t globCol, globRow;
1325   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1326   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1327
1328   //Int_t trigtimes[faltroWidthMax];
1329   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1330   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1331   //Double_t fTrigCutLow = 6;
1332   //Double_t fTrigCutHigh = 10;
1333   Double_t fTimeCutLow =  469e-09;
1334   Double_t fTimeCutHigh = 715e-09;
1335
1336   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1337
1338   // erase trigger maps
1339   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1340   {
1341     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1342     {
1343       ftriggersCut[i][j] = 0;
1344       ftriggers[i][j] = 0;
1345       ftriggersTime[i][j] = 0;
1346     }
1347   }
1348
1349   Int_t iglobCol=0, iglobRow=0;
1350   // go through triggers
1351   if( fCaloTrigger->GetEntries() > 0 )
1352   {
1353     // needs reset
1354     fCaloTrigger->Reset();
1355     while( fCaloTrigger->Next() )
1356     {
1357       fCaloTrigger->GetPosition( globCol, globRow );
1358       fCaloTrigger->GetNL0Times( ntimes );
1359       /*
1360       // no L0s
1361       if( ntimes < 1 )   continue;
1362       // get precise timings
1363       fCaloTrigger->GetL0Times( trigtimes );
1364       trigInCut = 0;
1365       for(Int_t i = 0; i < ntimes; i++ )
1366       {
1367         // save the first trigger time in channel
1368         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1369           triggersTime[globCol][globRow] = trigtimes[i];
1370         //printf("trigger times: %d\n",trigtimes[i]);
1371         // check if in cut
1372         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1373           trigInCut = 1;
1374
1375         fTrigTimes->Fill(trigtimes[i]);
1376       }
1377       */
1378    
1379       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1380       Int_t bit = 0;
1381       fCaloTrigger->GetTriggerBits(bit);
1382
1383       Int_t ts = 0;
1384       fCaloTrigger->GetL1TimeSum(ts);
1385       if (ts > 0)ftriggers[globCol][globRow] = 1;
1386       // number of triggered channels in event
1387       nTrigChannel++;
1388       // ... inside cut
1389       if(ts>0 && (bit >> 6 & 0x1))
1390       {
1391         iglobCol = globCol;
1392         iglobRow = globRow;
1393         nTrigChannelCut++;
1394         ftriggersCut[globCol][globRow] = 1;
1395       }
1396
1397     } // calo trigger entries
1398   } // has calo trigger entries
1399
1400   // part 2 go through the clusters here -----------------------------------
1401   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1402   Short_t cellAddr, nSACell, mclabel;
1403   //Int_t nSACell, iSACell, mclabel;
1404   Int_t iSACell;
1405   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
1406   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1407
1408   TRefArray *fCaloClusters = new TRefArray();
1409   fESD->GetEMCALClusters( fCaloClusters ); 
1410   nCluster = fCaloClusters->GetEntries();
1411
1412
1413   // save all cells times if there are clusters  
1414   if( nCluster > 0 ){
1415     // erase time map
1416     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){ 
1417       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1418         cellTime[i][j] = 0.;
1419         cellEnergy[i][j] = 0.;
1420       }
1421     }
1422
1423     // get cells
1424     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1425     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1426     nSACell = fCaloCells->GetNumberOfCells();
1427     for(iSACell = 0; iSACell < nSACell; iSACell++ ){ 
1428       // get the cell info *fCal
1429       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1430
1431       // get cell position 
1432       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta ); 
1433       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1434
1435       // convert co global phi eta  
1436       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1437       geta = ieta + nColsFeeModule*(nSupMod%2);
1438
1439       // save cell time and energy
1440       cellTime[geta][gphi] = cellTimeT;
1441       cellEnergy[geta][gphi] = cellAmp;
1442
1443     }
1444   }
1445
1446   Int_t nClusterTrig, nClusterTrigCut;
1447   UShort_t *cellAddrs;
1448   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1449   Float_t clsPos[3] = {0.,0.,0.};
1450
1451   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1452   {
1453     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1454     if(!cluster || !cluster->IsEMCAL()) continue;
1455
1456     // get cluster cells
1457     nCell = cluster->GetNCells();
1458
1459     // get cluster energy
1460     clsE = cluster->E();
1461
1462     // get cluster position
1463     cluster->GetPosition(clsPos);
1464     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1465     clsEta = clsPosVec.Eta();
1466     clsPhi = clsPosVec.Phi();
1467
1468     // get the cell addresses
1469     cellAddrs = cluster->GetCellsAbsId();
1470
1471     // check if the cluster contains cell, that was marked as triggered
1472     nClusterTrig = 0;
1473     nClusterTrigCut = 0;
1474
1475     // loop the cells to check, if cluser in acceptance
1476     // any cluster with a cell outside acceptance is not considered
1477     for( iCell = 0; iCell < nCell; iCell++ )
1478     {
1479      // check hot cell
1480      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]); 
1481
1482       // get cell position
1483       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1484       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1485
1486       // convert co global phi eta
1487       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1488       geta = ieta + nColsFeeModule*(nSupMod%2);
1489
1490       if( cellTime[geta][gphi] > 0. ){ 
1491         clusterTime += cellTime[geta][gphi];
1492         gCell++;
1493       }
1494
1495       // get corresponding FALTRO
1496       fphi = gphi / 2;
1497       feta = geta / 2;
1498
1499       // try to match with a triggered
1500       if( ftriggers[feta][fphi]==1)
1501       {  nClusterTrig++;
1502       }
1503       if( ftriggersCut[feta][fphi]==1)
1504       { nClusterTrigCut++;
1505       }
1506
1507     } // cells
1508
1509
1510     if( gCell > 0 ) 
1511       clusterTime = clusterTime / (Double_t)gCell;
1512     // fix the reconstriction code time 100ns jumps
1513     if( fESD->GetBunchCrossNumber() % 4 < 2 )
1514       clusterTime -= 0.0000001;
1515
1516     //fClsETime->Fill(clsE,clusterTime);
1517     //fClsEBftTrigCut->Fill(clsE);
1518
1519     if(nClusterTrig>0){
1520       //fClsETime1->Fill(clsE,clusterTime);
1521     }
1522
1523     if(nClusterTrig>0){
1524       cluster->SetChi2(1);
1525       //fClsEAftTrigCut1->Fill(clsE);                                               
1526     }
1527
1528     if(nClusterTrigCut>0){
1529       cluster->SetChi2(2);
1530       //fClsEAftTrigCut2->Fill(clsE);
1531     }
1532
1533     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1534     {
1535       cluster->SetChi2(3);
1536       //fClsEAftTrigCut3->Fill(clsE);
1537     }
1538
1539     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1540     {
1541       // cluster->SetChi2(4);
1542       //fClsEAftTrigCut4->Fill(clsE);
1543     }
1544     if(nClusterTrigCut<1)
1545     {
1546       cluster->SetChi2(0);
1547
1548       //fClsEAftTrigCut->Fill(clsE);
1549     }
1550
1551   } // clusters
1552 }
1553
1554
1555
1556
1557
1558
1559