updated
[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   ,fGeom(0)
90   ,fOutputList(0)
91   ,fqahist(1) 
92   ,fTrackCuts(0)
93   ,fCuts(0)
94   ,fIdentifiedAsOutInz(kFALSE)
95   ,fPassTheEventCut(kFALSE)
96   ,fRejectKinkMother(kFALSE)
97   ,fmcData(kFALSE)
98   ,fVz(0.0)
99   ,fCFM(0)      
100   ,fPID(0)
101   ,fPIDqa(0)           
102   ,fOpeningAngleCut(0.1)
103   ,fInvmassCut(0.01)    
104   ,fNoEvents(0)
105   ,fEMCAccE(0)
106   ,fTrkpt(0)
107   ,fTrkEovPBef(0)        
108   ,fTrkEovPAft(0)       
109   ,fdEdxBef(0)   
110   ,fdEdxAft(0)   
111   ,fIncpT(0)    
112   ,fIncpTM20(0) 
113   ,fInvmassLS(0)                
114   ,fInvmassULS(0)               
115   ,fOpeningAngleLS(0)   
116   ,fOpeningAngleULS(0)  
117   ,fPhotoElecPt(0)
118   ,fPhoElecPt(0)
119   ,fPhoElecPtM20(0)
120   ,fSameElecPt(0)
121   ,fSameElecPtM20(0)
122   ,fTrackPtBefTrkCuts(0)         
123   ,fTrackPtAftTrkCuts(0)
124   ,fTPCnsigma(0)
125   ,fCent(0)
126   ,fEleInfo(0)
127   /*
128   ,fClsEBftTrigCut(0)    
129   ,fClsEAftTrigCut(0)    
130   ,fClsEAftTrigCut1(0)   
131   ,fClsEAftTrigCut2(0)   
132   ,fClsEAftTrigCut3(0)   
133   ,fClsEAftTrigCut4(0)   
134   ,fClsETime(0)       
135   ,fClsETime1(0)       
136   ,fTrigTimes(0)
137   ,fCellCheck(0)
138   */
139   ,fInputHFEMC(0)
140   ,fInputAlle(0)
141   ,fIncpTMChfe(0)       
142   ,fIncpTMChfeAll(0)    
143   ,fIncpTMCM20hfe(0)    
144   ,fIncpTMCM20hfeAll(0) 
145   ,fIncpTMCpho(0)       
146   ,fIncpTMCM20pho(0)    
147   ,fPhoElecPtMC(0)
148   ,fPhoElecPtMCM20(0)
149   ,fSameElecPtMC(0)
150   ,fSameElecPtMCM20(0)
151 {
152   //Named constructor
153   
154   fPID = new AliHFEpid("hfePid");
155   fTrackCuts = new AliESDtrackCuts();
156   
157   // Define input and output slots here
158   // Input slot #0 works with a TChain
159   DefineInput(0, TChain::Class());
160   // Output slot #0 id reserved by the base class for AOD
161   // Output slot #1 writes into a TH1 container
162   // DefineOutput(1, TH1I::Class());
163   DefineOutput(1, TList::Class());
164   //  DefineOutput(3, TTree::Class());
165 }
166
167 //________________________________________________________________________
168 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() 
169   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
170   ,fESD(0)
171   ,fMC(0)
172   ,fGeom(0)
173   ,fOutputList(0)
174   ,fqahist(1)
175   ,fTrackCuts(0)
176   ,fCuts(0)
177   ,fIdentifiedAsOutInz(kFALSE)
178   ,fPassTheEventCut(kFALSE)
179   ,fRejectKinkMother(kFALSE)
180   ,fmcData(kFALSE)
181   ,fVz(0.0)
182   ,fCFM(0)      
183   ,fPID(0)       
184   ,fPIDqa(0)           
185   ,fOpeningAngleCut(0.1)
186   ,fInvmassCut(0.01)    
187   ,fNoEvents(0)
188   ,fEMCAccE(0)
189   ,fTrkpt(0)
190   ,fTrkEovPBef(0)        
191   ,fTrkEovPAft(0)        
192   ,fdEdxBef(0)   
193   ,fdEdxAft(0)   
194   ,fIncpT(0)     
195   ,fIncpTM20(0)  
196   ,fInvmassLS(0)                
197   ,fInvmassULS(0)               
198   ,fOpeningAngleLS(0)   
199   ,fOpeningAngleULS(0)  
200   ,fPhotoElecPt(0)
201   ,fPhoElecPt(0)
202   ,fPhoElecPtM20(0)
203   ,fSameElecPt(0)
204   ,fSameElecPtM20(0)
205   ,fTrackPtBefTrkCuts(0)         
206   ,fTrackPtAftTrkCuts(0)                  
207   ,fTPCnsigma(0)
208   ,fCent(0)
209   ,fEleInfo(0)
210   /*
211   ,fClsEBftTrigCut(0)    
212   ,fClsEAftTrigCut(0)    
213   ,fClsEAftTrigCut1(0)   
214   ,fClsEAftTrigCut2(0)   
215   ,fClsEAftTrigCut3(0)   
216   ,fClsEAftTrigCut4(0)   
217   ,fClsETime(0)       
218   ,fClsETime1(0)       
219   ,fTrigTimes(0)
220   ,fCellCheck(0)
221   */
222   ,fInputHFEMC(0)
223   ,fInputAlle(0)
224   ,fIncpTMChfe(0)       
225   ,fIncpTMChfeAll(0)    
226   ,fIncpTMCM20hfe(0)    
227   ,fIncpTMCM20hfeAll(0) 
228   ,fIncpTMCpho(0)       
229   ,fIncpTMCM20pho(0)    
230   ,fPhoElecPtMC(0)
231   ,fPhoElecPtMCM20(0)
232   ,fSameElecPtMC(0)
233   ,fSameElecPtMCM20(0)
234 {
235         //Default constructor
236         fPID = new AliHFEpid("hfePid");
237
238         fTrackCuts = new AliESDtrackCuts();
239         
240         // Constructor
241         // Define input and output slots here
242         // Input slot #0 works with a TChain
243         DefineInput(0, TChain::Class());
244         // Output slot #0 id reserved by the base class for AOD
245         // Output slot #1 writes into a TH1 container
246         // DefineOutput(1, TH1I::Class());
247         DefineOutput(1, TList::Class());
248         //DefineOutput(3, TTree::Class());
249 }
250 //_________________________________________
251
252 AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
253 {
254   //Destructor 
255   
256   delete fOutputList;
257   delete fGeom;
258   delete fPID;
259   delete fCFM;
260   delete fPIDqa;
261   delete fTrackCuts;
262 }
263 //_________________________________________
264
265 void AliAnalysisTaskHFECal::UserExec(Option_t*)
266 {
267   //Main loop
268   //Called for each event
269   
270   // create pointer to event
271   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
272   if (!fESD) {
273     printf("ERROR: fESD not available\n");
274     return;
275   }
276   
277   if(!fCuts){
278     AliError("HFE cuts not available");
279     return;
280   }
281   
282   if(!fPID->IsInitialized()){ 
283     // Initialize PID with the given run number
284     AliWarning("PID not initialised, get from Run no");
285     fPID->InitializePID(fESD->GetRunNumber());
286   }
287  
288   if(fmcData)fMC = MCEvent();
289   AliStack* stack = NULL;
290   if(fmcData && fMC)stack = fMC->Stack();
291
292   Float_t cent = -1.;
293   AliCentrality *centrality = fESD->GetCentrality(); 
294   cent = centrality->GetCentralityPercentile("V0M");
295
296   //---- fill MC track info
297   if(fmcData && fMC)
298     {
299     Int_t nParticles = stack->GetNtrack();
300     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
301       TParticle* particle = stack->Particle(iParticle);
302       int fPDG = particle->GetPdgCode(); 
303       double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
304       double pTMC = particle->Pt();
305       double proR = particle->R();
306       double etaMC = particle->Eta();
307       if(fabs(etaMC)>0.6)continue;
308       Bool_t mcInDtoE= kFALSE;
309       Bool_t mcInBtoE= kFALSE;
310
311       if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
312         {
313             int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();  
314             if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
315             if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
316             if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)fInputHFEMC->Fill(cent,pTMC);
317          }
318
319          if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
320
321       }
322     } 
323
324   fNoEvents->Fill(0);
325
326   Int_t fNOtrks =  fESD->GetNumberOfTracks();
327   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
328   
329   Double_t pVtxZ = -999;
330   pVtxZ = pVtx->GetZ();
331   
332   if(TMath::Abs(pVtxZ)>10) return;
333   fNoEvents->Fill(1);
334   
335   if(fNOtrks<2) return;
336   fNoEvents->Fill(2);
337   
338   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
339   if(!pidResponse){
340     AliDebug(1, "Using default PID Response");
341     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); 
342   }
343   
344   fPID->SetPIDResponse(pidResponse);
345   
346   fCFM->SetRecEventInfo(fESD);
347   
348   //Float_t cent = -1.;
349   //AliCentrality *centrality = fESD->GetCentrality(); 
350   //cent = centrality->GetCentralityPercentile("V0M");
351   fCent->Fill(cent);
352   
353   if(cent>90.) return;
354         
355  // Calorimeter info.
356  
357    FindTriggerClusters();
358
359   // make EMCAL array 
360   for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
361      {
362       AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
363       if(clust->IsEMCAL())
364         {
365          double clustE = clust->E();
366          float  emcx[3]; // cluster pos
367          clust->GetPosition(emcx);
368          TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
369          double emcphi = clustpos.Phi(); 
370          double emceta = clustpos.Eta();
371          double calInfo[5];
372          calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2(); 
373          //if(clustE>1.5)fEMCAccE->Fill(calInfo); 
374          //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo); 
375         }
376    }
377
378   // Track loop 
379   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
380     AliESDtrack* track = fESD->GetTrack(iTracks);
381     if (!track) {
382       printf("ERROR: Could not receive track %d\n", iTracks);
383       continue;
384     }
385    
386     Bool_t mcPho = kFALSE;
387     Bool_t mcDtoE= kFALSE;
388     Bool_t mcBtoE= kFALSE;
389     double mcele = -1.;
390     double mcpT = 0.0;
391     if(fmcData && fMC && stack)
392       {
393        Int_t label = TMath::Abs(track->GetLabel());
394        TParticle* particle = stack->Particle(label);
395        int mcpid = particle->GetPdgCode();
396        mcpT = particle->Pt();
397        //printf("MCpid = %d",mcpid);
398        if(particle->GetFirstMother()>-1)
399          {
400           int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
401           if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
402           if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
403           if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
404          } 
405        if(fabs(mcpid)==11 && mcDtoE)mcele= 1.; 
406        if(fabs(mcpid)==11 && mcBtoE)mcele= 2.; 
407        if(fabs(mcpid)==11 && mcPho)mcele= 3.; 
408       }
409  
410     if(TMath::Abs(track->Eta())>0.6) continue;
411     if(TMath::Abs(track->Pt()<2.0)) continue;
412     
413     fTrackPtBefTrkCuts->Fill(track->Pt());              
414     // RecKine: ITSTPC cuts  
415     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
416     
417     //RecKink
418     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
419       if(track->GetKinkIndex(0) != 0) continue;
420     } 
421     
422     // RecPrim
423     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
424     
425     // HFEcuts: ITS layers cuts
426     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
427     
428     // HFE cuts: TPC PID cleanup
429     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
430
431     
432     fTrackPtAftTrkCuts->Fill(track->Pt());              
433     
434     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
435     pt = track->Pt();
436     if(pt<2.0)continue;
437     
438     // Track extrapolation
439     
440     Int_t charge = track->Charge();
441     fTrkpt->Fill(pt);
442     mom = track->P();
443     phi = track->Phi();
444     eta = track->Eta();
445     dEdx = track->GetTPCsignal();
446     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
447     
448         double ncells = -1.0;
449         double m20 = -1.0;
450         double m02 = -1.0;
451         double disp = -1.0;
452         double rmatch = -1.0;  
453         double nmatch = -1.0;
454         double oppstatus = 0.0;
455
456     Bool_t fFlagPhotonicElec = kFALSE;
457     Bool_t fFlagConvinatElec = kFALSE;
458
459     Int_t clsId = track->GetEMCALcluster();
460     if (clsId>0){
461       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
462       if(clust && clust->IsEMCAL()){
463
464                 double clustE = clust->E();
465                 eop = clustE/fabs(mom);
466                 //double clustT = clust->GetTOF();
467                 ncells = clust->GetNCells();
468                 m02 = clust->GetM02();
469                 m20 = clust->GetM20();
470                 disp = clust->GetDispersion();
471                 double delphi = clust->GetTrackDx(); 
472                 double deleta = clust->GetTrackDz(); 
473                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
474                 nmatch = clust->GetNTracksMatched();
475
476                 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
477                 {
478                   SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele);
479                 }
480                 if(fFlagPhotonicElec)oppstatus = 1.0;
481                 if(fFlagConvinatElec)oppstatus = 2.0;
482                 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
483
484                   double valdedx[16];
485                   valdedx[0] = pt; valdedx[1] = dEdx; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
486                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = m02; valdedx[9] = m20; valdedx[10] = mcpT;
487                   valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = clust->Chi2();
488                   valdedx[15] = mcele;
489                   if(fqahist==1)fEleInfo->Fill(valdedx);
490                  
491
492       }
493     }
494         
495     fdEdxBef->Fill(mom,fTPCnSigma);
496     fTPCnsigma->Fill(mom,fTPCnSigma);
497     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
498
499     Int_t pidpassed = 1;
500     
501
502     //--- track accepted
503     AliHFEpidObject hfetrack;
504     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
505     hfetrack.SetRecTrack(track);
506     double binct = 10.5;
507     if((0.0< cent) && (cent<5.0)) binct = 0.5;
508     if((5.0< cent) && (cent<10.0)) binct = 1.5;
509     if((10.0< cent) && (cent<20.0)) binct = 2.5;
510     if((20.0< cent) && (cent<30.0)) binct = 3.5;
511     if((30.0< cent) && (cent<40.0)) binct = 4.5;
512     if((40.0< cent) && (cent<50.0)) binct = 5.5;
513     if((50.0< cent) && (cent<60.0)) binct = 6.5;
514     if((60.0< cent) && (cent<70.0)) binct = 7.5;
515     if((70.0< cent) && (cent<80.0)) binct = 8.5;
516     if((80.0< cent) && (cent<90.0)) binct = 9.5;
517     if((90.0< cent) && (cent<100.0)) binct = 10.5;
518
519     hfetrack.SetCentrality((int)binct); //added
520     hfetrack.SetPbPb();
521     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
522
523     if(pidpassed==0) continue;
524     
525     fTrkEovPAft->Fill(pt,eop);
526     fdEdxAft->Fill(mom,fTPCnSigma);
527
528     // Fill real data
529     fIncpT->Fill(cent,pt);    
530     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
531     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
532
533     if(m20>0.0 && m20<0.3)
534       { 
535        fIncpTM20->Fill(cent,pt);    
536        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
537        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
538      }
539     
540     // MC
541     if(mcele>0)
542       {
543
544           fIncpTMChfeAll->Fill(cent,pt);    
545           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);    
546
547        if(mcBtoE || mcDtoE)
548          {
549           fIncpTMChfe->Fill(cent,pt);    
550           if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);    
551          }
552        if(mcPho)
553         {
554          double phoval[3];
555          phoval[0] = cent;
556          phoval[1] = pt;
557          phoval[2] = fTPCnSigma;
558
559          fIncpTMCpho->Fill(phoval);    
560          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
561          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
562
563          if(m20>0.0 && m20<0.3) 
564            {
565             fIncpTMCM20pho->Fill(phoval);    
566             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
567             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
568            }
569         } 
570       } 
571  }
572  PostData(1, fOutputList);
573 }
574 //_________________________________________
575 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
576 {
577   //--- Check MC
578  
579   //Bool_t mcData = kFALSE;
580   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
581     {
582      fmcData = kTRUE;
583      printf("+++++ MC Data available");
584     }
585   if(fmcData)
586     {
587      printf("++++++++= MC analysis \n");
588     }
589   else
590    {
591      printf("++++++++ real data analysis \n");
592    }
593
594    printf("+++++++ QA hist %d",fqahist);
595
596   //---- Geometry
597   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
598
599   //--------Initialize PID
600   //fPID->SetHasMCData(kFALSE);
601   fPID->SetHasMCData(fmcData);
602   if(!fPID->GetNumberOfPIDdetectors()) 
603     {
604       fPID->AddDetector("TPC", 0);
605       fPID->AddDetector("EMCAL", 1);
606     }
607   
608   Double_t params[4];
609   const char *cutmodel;
610   cutmodel = "pol0";
611   params[0] = -1.0; //sigma min
612   double maxnSig = 3.0;
613   if(fmcData)
614     {
615      params[0] = -5.0; //sigma min
616      maxnSig = 5.0; 
617     } 
618
619   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
620
621   fPID->SortDetectors(); 
622   fPIDqa = new AliHFEpidQAmanager();
623   fPIDqa->Initialize(fPID);
624   
625   //--------Initialize correction Framework and Cuts
626   fCFM = new AliCFManager;
627   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
628   fCFM->SetNStepParticle(kNcutSteps);
629   for(Int_t istep = 0; istep < kNcutSteps; istep++)
630     fCFM->SetParticleCutsList(istep, NULL);
631   
632   if(!fCuts){
633     AliWarning("Cuts not available. Default cuts will be used");
634     fCuts = new AliHFEcuts;
635     fCuts->CreateStandardCuts();
636   }
637   fCuts->Initialize(fCFM);
638   
639   //---------Output Tlist
640   fOutputList = new TList();
641   fOutputList->SetOwner();
642   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
643   
644   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
645   fOutputList->Add(fNoEvents);
646   
647   Int_t binsE[5] =    {250, 100,  1000, 100,   10};
648   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5}; 
649   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5}; 
650   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
651   if(fqahist==1)fOutputList->Add(fEMCAccE);
652
653   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
654   fOutputList->Add(fTrkpt);
655   
656   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
657   fOutputList->Add(fTrackPtBefTrkCuts);
658   
659   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
660   fOutputList->Add(fTrackPtAftTrkCuts);
661   
662   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
663   fOutputList->Add(fTPCnsigma);
664   
665   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
666   fOutputList->Add(fTrkEovPBef);
667   
668   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
669   fOutputList->Add(fTrkEovPAft);
670   
671   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
672   fOutputList->Add(fdEdxBef);
673   
674   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
675   fOutputList->Add(fdEdxAft);
676   
677   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50);
678   fOutputList->Add(fIncpT);
679
680   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",100,0,100,100,0,50);
681   fOutputList->Add(fIncpTM20);
682   
683   Int_t nBinspho[8] =  { 100, 100, 500, 12,   50,    4, 200,   8};
684   Double_t minpho[8] = {  0.,  0.,  0., -2.5,  0, -0.5,   0,-1.5};   
685   Double_t maxpho[8] = {100., 50., 0.5, 3.5,   1,  3.5,   2, 6.5};   
686
687   fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 8, nBinspho,minpho, maxpho);
688   fOutputList->Add(fInvmassLS);
689   
690   fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 8, nBinspho,minpho, maxpho);
691   fOutputList->Add(fInvmassULS);
692   
693   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
694   fOutputList->Add(fOpeningAngleLS);
695   
696   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
697   fOutputList->Add(fOpeningAngleULS);
698   
699   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
700   fOutputList->Add(fPhotoElecPt);
701   
702   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",100,0,100,100,0,50);
703   fOutputList->Add(fPhoElecPt);
704   
705   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",100,0,100,100,0,50);
706   fOutputList->Add(fPhoElecPtM20);
707
708   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",100,0,100,100,0,50);
709   fOutputList->Add(fSameElecPt);
710
711   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",100,0,100,100,0,50);
712   fOutputList->Add(fSameElecPtM20);
713
714   fCent = new TH1F("fCent","Centrality",100,0,100) ;
715   fOutputList->Add(fCent);
716  
717   // Make common binning
718   const Double_t kMinP = 0.;
719   const Double_t kMaxP = 50.;
720   const Double_t kTPCSigMim = 40.;
721   const Double_t kTPCSigMax = 140.;
722
723   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
724   Int_t nBins[16] =  {  250,        200,   60,    20,   100,  300,  50,   40,   200, 200,  250, 100,    3,    5,   10,    8};
725   Double_t min[16] = {kMinP,  kTPCSigMim, 1.0,  -1.0,  -6.0,    0,   0,    0,   0.0, 0.0,  0.0,   0, -1.5, -0.5, -0.5, -1.5};
726   Double_t max[16] = {kMaxP,  kTPCSigMax, 4.0,   1.0,   4.0,  3.0, 0.1,   40,   2.0, 2.0, 50.0, 100,  1.5,  4.5,  9.5,  6.5};
727   fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;M02;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
728   if(fqahist==1)fOutputList->Add(fEleInfo);
729
730   //<---  Trigger info
731   /*
732   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
733   fOutputList->Add(fClsEBftTrigCut);
734
735   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
736   fOutputList->Add(fClsEAftTrigCut);
737
738   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
739   fOutputList->Add(fClsEAftTrigCut1);
740
741   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
742   fOutputList->Add(fClsEAftTrigCut2);
743
744   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
745   fOutputList->Add(fClsEAftTrigCut3);
746
747   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
748   fOutputList->Add(fClsEAftTrigCut4);
749
750   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
751   fOutputList->Add(fClsETime);
752
753   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
754   fOutputList->Add(fClsETime1);
755
756   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
757   fOutputList->Add(fTrigTimes);
758
759   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
760   fOutputList->Add(fCellCheck);
761   */
762   //<---------- MC
763
764   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",100,0,100,100,0,50);
765   fOutputList->Add(fInputHFEMC);
766
767   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",100,0,100,100,0,50);
768   fOutputList->Add(fInputAlle);
769
770   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",100,0,100,100,0,50);
771   fOutputList->Add(fIncpTMChfe);
772
773   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",100,0,100,100,0,50);
774   fOutputList->Add(fIncpTMChfeAll);
775
776   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",100,0,100,100,0,50);
777   fOutputList->Add(fIncpTMCM20hfe);
778
779   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",100,0,100,100,0,50);
780   fOutputList->Add(fIncpTMCM20hfeAll);
781
782
783   Int_t nBinspho2[3] =  { 100, 100,    7};
784   Double_t minpho2[3] = {  0.,  0., -2.5};   
785   Double_t maxpho2[3] = {100., 50.,  4.5};   
786
787   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",3,nBinspho2,minpho2,maxpho2);
788   fOutputList->Add(fIncpTMCpho);
789
790   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",3,nBinspho2,minpho2,maxpho2);
791   fOutputList->Add(fIncpTMCM20pho);
792
793   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",3,nBinspho2,minpho2,maxpho2);
794   fOutputList->Add(fPhoElecPtMC);
795   
796   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",3,nBinspho2,minpho2,maxpho2);
797   fOutputList->Add(fPhoElecPtMCM20);
798
799   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",3,nBinspho2,minpho2,maxpho2);
800   fOutputList->Add(fSameElecPtMC);
801
802   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",3,nBinspho2,minpho2,maxpho2);
803   fOutputList->Add(fSameElecPtMCM20);
804
805
806   PostData(1,fOutputList);
807 }
808
809 //________________________________________________________________________
810 void AliAnalysisTaskHFECal::Terminate(Option_t *)
811 {
812   // Info("Terminate");
813         AliAnalysisTaskSE::Terminate();
814 }
815
816 //________________________________________________________________________
817 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
818 {
819   // Check single track cuts for a given cut step
820   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
821   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
822   return kTRUE;
823 }
824 //_________________________________________
825 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
826 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
827 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)
828 {
829   //Identify non-heavy flavour electrons using Invariant mass method
830   
831   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
832   fTrackCuts->SetRequireTPCRefit(kTRUE);
833   fTrackCuts->SetRequireITSRefit(kTRUE);
834   fTrackCuts->SetEtaRange(-0.9,0.9);
835   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
836   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
837   fTrackCuts->SetMinNClustersTPC(90);
838   
839   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
840   
841   Bool_t flagPhotonicElec = kFALSE;
842   Bool_t flagConvinatElec = kFALSE;
843   
844   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
845   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
846     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
847     if (!trackAsso) {
848       printf("ERROR: Could not receive track %d\n", jTracks);
849       continue;
850     }
851     if(itrack==jTracks)continue;    
852
853     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
854     Double_t mass=999., width = -999;
855     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
856     
857     ptPrim = track->Pt();
858
859     dEdxAsso = trackAsso->GetTPCsignal();
860     ptAsso = trackAsso->Pt();
861     Int_t chargeAsso = trackAsso->Charge();
862     Int_t charge = track->Charge();
863     
864     //if(ptAsso <0.3) continue;
865     if(ptAsso <0.5) continue;
866     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
867     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
868     
869     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
870     if(charge>0) fPDGe1 = -11;
871     if(chargeAsso>0) fPDGe2 = -11;
872     
873     if(charge == chargeAsso) fFlagLS = kTRUE;
874     if(charge != chargeAsso) fFlagULS = kTRUE;
875     
876     AliKFParticle ge1(*track, fPDGe1);
877     AliKFParticle ge2(*trackAsso, fPDGe2);
878     AliKFParticle recg(ge1, ge2);
879     
880     if(recg.GetNDF()<1) continue;
881     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
882     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
883     
884     AliKFVertex primV(*pVtx);
885     primV += recg;
886     recg.SetProductionVertex(primV);
887     
888     recg.SetMassConstraint(0,0.0001);
889     
890     openingAngle = ge1.GetAngle(ge2);
891     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
892     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
893     
894     
895     recg.GetMass(mass,width);
896     
897     double ishower = 0;
898     if(shower>0.0 && shower<0.3)ishower = 1;
899
900     double phoinfo[8];
901     phoinfo[0] = cent;
902     phoinfo[1] = ptPrim;
903     phoinfo[2] = mass;
904     phoinfo[3] = nSig;
905     phoinfo[4] = openingAngle;
906     phoinfo[5] = ishower;
907     phoinfo[6] = ep;
908     phoinfo[7] = mce;
909
910     if(fFlagLS) fInvmassLS->Fill(phoinfo);
911     if(fFlagULS) fInvmassULS->Fill(phoinfo);
912     
913     if(openingAngle > fOpeningAngleCut) continue;
914       
915     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
916       flagPhotonicElec = kTRUE;
917     }
918     if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
919       flagConvinatElec = kTRUE;
920     }
921     
922   }
923   fFlagPhotonicElec = flagPhotonicElec;
924   fFlagConvinatElec = flagConvinatElec;
925   
926 }
927
928
929 //_________________________________________
930 void AliAnalysisTaskHFECal::FindTriggerClusters()
931 {
932   // constants
933   const int nModuleCols = 2;
934   const int nModuleRows = 5;
935   const int nColsFeeModule = 48;
936   const int nRowsFeeModule = 24;
937   const int nColsFaltroModule = 24;
938   const int nRowsFaltroModule = 12;
939   //const int faltroWidthMax = 20;
940
941   // part 1, trigger extraction -------------------------------------
942   Int_t globCol, globRow;
943   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
944   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
945
946   //Int_t trigtimes[faltroWidthMax];
947   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
948   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
949   //Double_t fTrigCutLow = 6;
950   //Double_t fTrigCutHigh = 10;
951   Double_t fTimeCutLow =  469e-09;
952   Double_t fTimeCutHigh = 715e-09;
953
954   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
955
956   // erase trigger maps
957   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
958   {
959     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
960     {
961       ftriggersCut[i][j] = 0;
962       ftriggers[i][j] = 0;
963       ftriggersTime[i][j] = 0;
964     }
965   }
966
967   Int_t iglobCol=0, iglobRow=0;
968   // go through triggers
969   if( fCaloTrigger->GetEntries() > 0 )
970   {
971     // needs reset
972     fCaloTrigger->Reset();
973     while( fCaloTrigger->Next() )
974     {
975       fCaloTrigger->GetPosition( globCol, globRow );
976       fCaloTrigger->GetNL0Times( ntimes );
977       /*
978       // no L0s
979       if( ntimes < 1 )   continue;
980       // get precise timings
981       fCaloTrigger->GetL0Times( trigtimes );
982       trigInCut = 0;
983       for(Int_t i = 0; i < ntimes; i++ )
984       {
985         // save the first trigger time in channel
986         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
987           triggersTime[globCol][globRow] = trigtimes[i];
988         //printf("trigger times: %d\n",trigtimes[i]);
989         // check if in cut
990         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
991           trigInCut = 1;
992
993         fTrigTimes->Fill(trigtimes[i]);
994       }
995       */
996    
997       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
998       Int_t bit = 0;
999       fCaloTrigger->GetTriggerBits(bit);
1000
1001       Int_t ts = 0;
1002       fCaloTrigger->GetL1TimeSum(ts);
1003       if (ts > 0)ftriggers[globCol][globRow] = 1;
1004       // number of triggered channels in event
1005       nTrigChannel++;
1006       // ... inside cut
1007       if(ts>0 && (bit >> 6 & 0x1))
1008       {
1009         iglobCol = globCol;
1010         iglobRow = globRow;
1011         nTrigChannelCut++;
1012         ftriggersCut[globCol][globRow] = 1;
1013       }
1014
1015     } // calo trigger entries
1016   } // has calo trigger entries
1017
1018   // part 2 go through the clusters here -----------------------------------
1019   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1020   Short_t cellAddr, nSACell, mclabel;
1021   //Int_t nSACell, iSACell, mclabel;
1022   Int_t iSACell;
1023   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
1024   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1025
1026   TRefArray *fCaloClusters = new TRefArray();
1027   fESD->GetEMCALClusters( fCaloClusters ); 
1028   nCluster = fCaloClusters->GetEntries();
1029
1030
1031   // save all cells times if there are clusters  
1032   if( nCluster > 0 ){
1033     // erase time map
1034     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){ 
1035       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1036         cellTime[i][j] = 0.;
1037         cellEnergy[i][j] = 0.;
1038       }
1039     }
1040
1041     // get cells
1042     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1043     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1044     nSACell = fCaloCells->GetNumberOfCells();
1045     for(iSACell = 0; iSACell < nSACell; iSACell++ ){ 
1046       // get the cell info *fCal
1047       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1048
1049       // get cell position 
1050       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta ); 
1051       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1052
1053       // convert co global phi eta  
1054       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1055       geta = ieta + nColsFeeModule*(nSupMod%2);
1056
1057       // save cell time and energy
1058       cellTime[geta][gphi] = cellTimeT;
1059       cellEnergy[geta][gphi] = cellAmp;
1060
1061     }
1062   }
1063
1064   Int_t nClusterTrig, nClusterTrigCut;
1065   UShort_t *cellAddrs;
1066   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1067   Float_t clsPos[3] = {0.,0.,0.};
1068
1069   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1070   {
1071     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1072     if(!cluster || !cluster->IsEMCAL()) continue;
1073
1074     // get cluster cells
1075     nCell = cluster->GetNCells();
1076
1077     // get cluster energy
1078     clsE = cluster->E();
1079
1080     // get cluster position
1081     cluster->GetPosition(clsPos);
1082     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1083     clsEta = clsPosVec.Eta();
1084     clsPhi = clsPosVec.Phi();
1085
1086     // get the cell addresses
1087     cellAddrs = cluster->GetCellsAbsId();
1088
1089     // check if the cluster contains cell, that was marked as triggered
1090     nClusterTrig = 0;
1091     nClusterTrigCut = 0;
1092
1093     // loop the cells to check, if cluser in acceptance
1094     // any cluster with a cell outside acceptance is not considered
1095     for( iCell = 0; iCell < nCell; iCell++ )
1096     {
1097      // check hot cell
1098      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]); 
1099
1100       // get cell position
1101       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1102       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1103
1104       // convert co global phi eta
1105       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1106       geta = ieta + nColsFeeModule*(nSupMod%2);
1107
1108       if( cellTime[geta][gphi] > 0. ){ 
1109         clusterTime += cellTime[geta][gphi];
1110         gCell++;
1111       }
1112
1113       // get corresponding FALTRO
1114       fphi = gphi / 2;
1115       feta = geta / 2;
1116
1117       // try to match with a triggered
1118       if( ftriggers[feta][fphi]==1)
1119       {  nClusterTrig++;
1120       }
1121       if( ftriggersCut[feta][fphi]==1)
1122       { nClusterTrigCut++;
1123       }
1124
1125     } // cells
1126
1127
1128     if( gCell > 0 ) 
1129       clusterTime = clusterTime / (Double_t)gCell;
1130     // fix the reconstriction code time 100ns jumps
1131     if( fESD->GetBunchCrossNumber() % 4 < 2 )
1132       clusterTime -= 0.0000001;
1133
1134     //fClsETime->Fill(clsE,clusterTime);
1135     //fClsEBftTrigCut->Fill(clsE);
1136
1137     if(nClusterTrig>0){
1138       //fClsETime1->Fill(clsE,clusterTime);
1139     }
1140
1141     if(nClusterTrig>0){
1142       cluster->SetChi2(1);
1143       //fClsEAftTrigCut1->Fill(clsE);                                               
1144     }
1145
1146     if(nClusterTrigCut>0){
1147       cluster->SetChi2(2);
1148       //fClsEAftTrigCut2->Fill(clsE);
1149     }
1150
1151     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1152     {
1153       cluster->SetChi2(3);
1154       //fClsEAftTrigCut3->Fill(clsE);
1155     }
1156
1157     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1158     {
1159       // cluster->SetChi2(4);
1160       //fClsEAftTrigCut4->Fill(clsE);
1161     }
1162     if(nClusterTrigCut<1)
1163     {
1164       cluster->SetChi2(0);
1165
1166       //fClsEAftTrigCut->Fill(clsE);
1167     }
1168
1169   } // clusters
1170 }
1171
1172
1173
1174
1175
1176
1177