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