387617c9fd1b4359b137200eac7d28e4a1e3c5f8
[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[4];
648          phoval[0] = cent;
649          phoval[1] = pt;
650          phoval[2] = fTPCnSigma;
651          phoval[3] = iHijing;
652
653          fIncpTMCpho->Fill(phoval);    
654          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
655          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
656
657          if(m20>0.0 && m20<0.3) 
658            {
659             fIncpTMCM20pho->Fill(phoval);    
660             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
661             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
662             // pi0->g->e
663             if(mcWeight>-1)
664               {
665                fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);    
666                if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
667                if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
668               }
669            }
670         } 
671       } 
672  }
673  PostData(1, fOutputList);
674 }
675 //_________________________________________
676 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
677 {
678   //--- Check MC
679  
680   //Bool_t mcData = kFALSE;
681   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
682     {
683      fmcData = kTRUE;
684      printf("+++++ MC Data available");
685     }
686   if(fmcData)
687     {
688      printf("++++++++= MC analysis \n");
689     }
690   else
691    {
692      printf("++++++++ real data analysis \n");
693    }
694
695    printf("+++++++ QA hist %d",fqahist);
696
697   //---- Geometry
698   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
699
700   //--------Initialize PID
701   //fPID->SetHasMCData(kFALSE);
702   fPID->SetHasMCData(fmcData);
703   if(!fPID->GetNumberOfPIDdetectors()) 
704     {
705       fPID->AddDetector("TPC", 0);
706       fPID->AddDetector("EMCAL", 1);
707     }
708   
709   Double_t params[4];
710   const char *cutmodel;
711   cutmodel = "pol0";
712   params[0] = -1.0; //sigma min
713   double maxnSig = 3.0;
714   if(fmcData)
715     {
716      params[0] = -5.0; //sigma min
717      maxnSig = 5.0; 
718     } 
719
720   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
721
722   fPID->SortDetectors(); 
723   fPIDqa = new AliHFEpidQAmanager();
724   fPIDqa->Initialize(fPID);
725  
726   //------- fcut --------------  
727   fCuts = new AliHFEcuts();
728   fCuts->CreateStandardCuts();
729   //fCuts->SetMinNClustersTPC(100);
730   fCuts->SetMinNClustersTPC(90);
731   fCuts->SetMinRatioTPCclusters(0.6);
732   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
733   //fCuts->SetMinNClustersITS(3);
734   fCuts->SetMinNClustersITS(2);
735   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
736   fCuts->SetCheckITSLayerStatus(kFALSE);
737   fCuts->SetVertexRange(10.);
738   fCuts->SetTOFPIDStep(kFALSE);
739   fCuts->SetPtRange(2, 50);
740   fCuts->SetMaxImpactParam(3.,3.);
741
742   //--------Initialize correction Framework and Cuts
743   fCFM = new AliCFManager;
744   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
745   fCFM->SetNStepParticle(kNcutSteps);
746   for(Int_t istep = 0; istep < kNcutSteps; istep++)
747     fCFM->SetParticleCutsList(istep, NULL);
748   
749   if(!fCuts){
750     AliWarning("Cuts not available. Default cuts will be used");
751     fCuts = new AliHFEcuts;
752     fCuts->CreateStandardCuts();
753   }
754   fCuts->Initialize(fCFM);
755   
756   //---------Output Tlist
757   fOutputList = new TList();
758   fOutputList->SetOwner();
759   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
760   
761   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
762   fOutputList->Add(fNoEvents);
763   
764   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
765   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5}; 
766   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5}; 
767   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
768   if(fqahist==1)fOutputList->Add(fEMCAccE);
769
770   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
771   fOutputList->Add(fTrkpt);
772   
773   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
774   fOutputList->Add(fTrackPtBefTrkCuts);
775   
776   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
777   fOutputList->Add(fTrackPtAftTrkCuts);
778   
779   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
780   fOutputList->Add(fTPCnsigma);
781   
782   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
783   fOutputList->Add(fTrkEovPBef);
784   
785   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
786   fOutputList->Add(fTrkEovPAft);
787   
788   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
789   fOutputList->Add(fdEdxBef);
790   
791   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
792   fOutputList->Add(fdEdxAft);
793   
794   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
795   fOutputList->Add(fIncpT);
796
797   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
798   fOutputList->Add(fIncpTM20);
799   
800   Int_t nBinspho[9] =  { 200, 100, 500, 12,   50,    4, 200,   8, 100};
801   Double_t minpho[9] = {  0.,  0.,  0., -2.5,  0, -0.5,   0,-1.5,   0};   
802   Double_t maxpho[9] = {100., 50., 0.5, 3.5,   1,  3.5,   2, 6.5,  50};   
803
804   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);
805   fOutputList->Add(fInvmassLS);
806   
807   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);
808   fOutputList->Add(fInvmassULS);
809   
810   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
811   fOutputList->Add(fOpeningAngleLS);
812   
813   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
814   fOutputList->Add(fOpeningAngleULS);
815   
816   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
817   fOutputList->Add(fPhotoElecPt);
818   
819   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
820   fOutputList->Add(fPhoElecPt);
821   
822   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
823   fOutputList->Add(fPhoElecPtM20);
824
825   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
826   fOutputList->Add(fSameElecPt);
827
828   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
829   fOutputList->Add(fSameElecPtM20);
830
831   fCent = new TH1F("fCent","Centrality",200,0,100) ;
832   fOutputList->Add(fCent);
833  
834   // Make common binning
835   const Double_t kMinP = 0.;
836   const Double_t kMaxP = 50.;
837   //const Double_t kTPCSigMim = 40.;
838   //const Double_t kTPCSigMax = 140.;
839
840   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
841   Int_t nBins[16] =  {  250,    10,   60,    20,   100,  300,  50,   40,   200, 200,  250, 200,    3,    5, 100,    8};
842   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};
843   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};
844   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);
845   if(fqahist==1)fOutputList->Add(fEleInfo);
846
847   //<---  Trigger info
848   /*
849   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
850   fOutputList->Add(fClsEBftTrigCut);
851
852   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
853   fOutputList->Add(fClsEAftTrigCut);
854
855   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
856   fOutputList->Add(fClsEAftTrigCut1);
857
858   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
859   fOutputList->Add(fClsEAftTrigCut2);
860
861   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
862   fOutputList->Add(fClsEAftTrigCut3);
863
864   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
865   fOutputList->Add(fClsEAftTrigCut4);
866
867   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
868   fOutputList->Add(fClsETime);
869
870   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
871   fOutputList->Add(fClsETime1);
872
873   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
874   fOutputList->Add(fTrigTimes);
875
876   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
877   fOutputList->Add(fCellCheck);
878   */
879   //<---------- MC
880
881   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
882   fOutputList->Add(fInputHFEMC);
883
884   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
885   fOutputList->Add(fInputAlle);
886
887   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
888   fOutputList->Add(fIncpTMChfe);
889
890   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
891   fOutputList->Add(fIncpTMChfeAll);
892
893   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
894   fOutputList->Add(fIncpTMCM20hfe);
895
896   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
897   fOutputList->Add(fIncpTMCM20hfeAll);
898
899
900   Int_t nBinspho2[4] =  { 200, 100,    7, 3};
901   Double_t minpho2[4] = {  0.,  0., -2.5, -0.5};   
902   Double_t maxpho2[4] = {100., 50.,  4.5, 2.5};   
903
904   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",4,nBinspho2,minpho2,maxpho2);
905   fOutputList->Add(fIncpTMCpho);
906
907   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",4,nBinspho2,minpho2,maxpho2);
908   fOutputList->Add(fIncpTMCM20pho);
909
910   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",4,nBinspho2,minpho2,maxpho2);
911   fOutputList->Add(fPhoElecPtMC);
912   
913   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",4,nBinspho2,minpho2,maxpho2);
914   fOutputList->Add(fPhoElecPtMCM20);
915
916   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",4,nBinspho2,minpho2,maxpho2);
917   fOutputList->Add(fSameElecPtMC);
918
919   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",4,nBinspho2,minpho2,maxpho2);
920   fOutputList->Add(fSameElecPtMCM20);
921
922   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",4,nBinspho2,minpho2,maxpho2);
923   fOutputList->Add(fIncpTMCM20pho_pi0e);
924
925   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",4,nBinspho2,minpho2,maxpho2);
926   fOutputList->Add(fPhoElecPtMCM20_pi0e);
927
928   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",4,nBinspho2,minpho2,maxpho2);
929   fOutputList->Add(fSameElecPtMCM20_pi0e);
930
931   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
932   fOutputList->Add(CheckNclust);
933
934   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
935   fOutputList->Add(CheckNits);
936
937   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
938   fOutputList->Add(Hpi0pTcheck);
939
940   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
941   fOutputList->Add(HphopTcheck);
942
943   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
944   fOutputList->Add(fMomDtoE);
945
946   PostData(1,fOutputList);
947 }
948
949 //________________________________________________________________________
950 void AliAnalysisTaskHFECal::Terminate(Option_t *)
951 {
952   // Info("Terminate");
953         AliAnalysisTaskSE::Terminate();
954 }
955
956 //________________________________________________________________________
957 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
958 {
959   // Check single track cuts for a given cut step
960   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
961   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
962   return kTRUE;
963 }
964 //_________________________________________
965 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
966 //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
967 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)
968 {
969   //Identify non-heavy flavour electrons using Invariant mass method
970   
971   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
972   fTrackCuts->SetRequireTPCRefit(kTRUE);
973   fTrackCuts->SetRequireITSRefit(kTRUE);
974   fTrackCuts->SetEtaRange(-0.9,0.9);
975   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
976   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
977   fTrackCuts->SetMinNClustersTPC(90);
978   
979   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
980   
981   Bool_t flagPhotonicElec = kFALSE;
982   Bool_t flagConvinatElec = kFALSE;
983   
984   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
985   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
986     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
987     if (!trackAsso) {
988       printf("ERROR: Could not receive track %d\n", jTracks);
989       continue;
990     }
991     if(itrack==jTracks)continue;    
992
993     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
994     Double_t mass=999., width = -999;
995     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
996     
997     ptPrim = track->Pt();
998
999     dEdxAsso = trackAsso->GetTPCsignal();
1000     ptAsso = trackAsso->Pt();
1001     Int_t chargeAsso = trackAsso->Charge();
1002     Int_t charge = track->Charge();
1003     
1004
1005     if(ptAsso <0.5) continue;
1006     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1007     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1008     
1009     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1010     if(charge>0) fPDGe1 = -11;
1011     if(chargeAsso>0) fPDGe2 = -11;
1012  
1013     //printf("chargeAsso = %d\n",chargeAsso);
1014     //printf("charge = %d\n",charge);
1015     if(charge == chargeAsso) fFlagLS = kTRUE;
1016     if(charge != chargeAsso) fFlagULS = kTRUE;
1017     
1018     //printf("fFlagLS = %d\n",fFlagLS);
1019     //printf("fFlagULS = %d\n",fFlagULS);
1020     //printf("\n");
1021
1022     AliKFParticle ge1(*track, fPDGe1);
1023     AliKFParticle ge2(*trackAsso, fPDGe2);
1024     AliKFParticle recg(ge1, ge2);
1025     
1026     if(recg.GetNDF()<1) continue;
1027     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1028     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1029     
1030     AliKFVertex primV(*pVtx);
1031     primV += recg;
1032     recg.SetProductionVertex(primV);
1033     
1034     recg.SetMassConstraint(0,0.0001);
1035     
1036     openingAngle = ge1.GetAngle(ge2);
1037     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1038     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1039     
1040     
1041     recg.GetMass(mass,width);
1042     
1043     double ishower = 0;
1044     if(shower>0.0 && shower<0.3)ishower = 1;
1045
1046     double phoinfo[9];
1047     phoinfo[0] = cent;
1048     phoinfo[1] = ptPrim;
1049     phoinfo[2] = mass;
1050     phoinfo[3] = nSig;
1051     phoinfo[4] = openingAngle;
1052     phoinfo[5] = ishower;
1053     phoinfo[6] = ep;
1054     phoinfo[7] = mce;
1055     phoinfo[8] = ptAsso;
1056
1057     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1058     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1059     
1060     //printf("fInvmassCut %f\n",fInvmassCut);
1061     //printf("openingAngle %f\n",fOpeningAngleCut);
1062
1063     if(openingAngle > fOpeningAngleCut) continue;
1064       
1065     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1066       flagPhotonicElec = kTRUE;
1067     }
1068     if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1069       flagConvinatElec = kTRUE;
1070     }
1071     
1072   }
1073   fFlagPhotonicElec = flagPhotonicElec;
1074   fFlagConvinatElec = flagConvinatElec;
1075   
1076 }
1077
1078
1079 //_________________________________________
1080 void AliAnalysisTaskHFECal::FindTriggerClusters()
1081 {
1082   // constants
1083   const int nModuleCols = 2;
1084   const int nModuleRows = 5;
1085   const int nColsFeeModule = 48;
1086   const int nRowsFeeModule = 24;
1087   const int nColsFaltroModule = 24;
1088   const int nRowsFaltroModule = 12;
1089   //const int faltroWidthMax = 20;
1090
1091   // part 1, trigger extraction -------------------------------------
1092   Int_t globCol, globRow;
1093   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1094   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1095
1096   //Int_t trigtimes[faltroWidthMax];
1097   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1098   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1099   //Double_t fTrigCutLow = 6;
1100   //Double_t fTrigCutHigh = 10;
1101   Double_t fTimeCutLow =  469e-09;
1102   Double_t fTimeCutHigh = 715e-09;
1103
1104   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1105
1106   // erase trigger maps
1107   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1108   {
1109     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1110     {
1111       ftriggersCut[i][j] = 0;
1112       ftriggers[i][j] = 0;
1113       ftriggersTime[i][j] = 0;
1114     }
1115   }
1116
1117   Int_t iglobCol=0, iglobRow=0;
1118   // go through triggers
1119   if( fCaloTrigger->GetEntries() > 0 )
1120   {
1121     // needs reset
1122     fCaloTrigger->Reset();
1123     while( fCaloTrigger->Next() )
1124     {
1125       fCaloTrigger->GetPosition( globCol, globRow );
1126       fCaloTrigger->GetNL0Times( ntimes );
1127       /*
1128       // no L0s
1129       if( ntimes < 1 )   continue;
1130       // get precise timings
1131       fCaloTrigger->GetL0Times( trigtimes );
1132       trigInCut = 0;
1133       for(Int_t i = 0; i < ntimes; i++ )
1134       {
1135         // save the first trigger time in channel
1136         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1137           triggersTime[globCol][globRow] = trigtimes[i];
1138         //printf("trigger times: %d\n",trigtimes[i]);
1139         // check if in cut
1140         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1141           trigInCut = 1;
1142
1143         fTrigTimes->Fill(trigtimes[i]);
1144       }
1145       */
1146    
1147       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1148       Int_t bit = 0;
1149       fCaloTrigger->GetTriggerBits(bit);
1150
1151       Int_t ts = 0;
1152       fCaloTrigger->GetL1TimeSum(ts);
1153       if (ts > 0)ftriggers[globCol][globRow] = 1;
1154       // number of triggered channels in event
1155       nTrigChannel++;
1156       // ... inside cut
1157       if(ts>0 && (bit >> 6 & 0x1))
1158       {
1159         iglobCol = globCol;
1160         iglobRow = globRow;
1161         nTrigChannelCut++;
1162         ftriggersCut[globCol][globRow] = 1;
1163       }
1164
1165     } // calo trigger entries
1166   } // has calo trigger entries
1167
1168   // part 2 go through the clusters here -----------------------------------
1169   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1170   Short_t cellAddr, nSACell, mclabel;
1171   //Int_t nSACell, iSACell, mclabel;
1172   Int_t iSACell;
1173   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
1174   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1175
1176   TRefArray *fCaloClusters = new TRefArray();
1177   fESD->GetEMCALClusters( fCaloClusters ); 
1178   nCluster = fCaloClusters->GetEntries();
1179
1180
1181   // save all cells times if there are clusters  
1182   if( nCluster > 0 ){
1183     // erase time map
1184     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){ 
1185       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1186         cellTime[i][j] = 0.;
1187         cellEnergy[i][j] = 0.;
1188       }
1189     }
1190
1191     // get cells
1192     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1193     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1194     nSACell = fCaloCells->GetNumberOfCells();
1195     for(iSACell = 0; iSACell < nSACell; iSACell++ ){ 
1196       // get the cell info *fCal
1197       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1198
1199       // get cell position 
1200       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta ); 
1201       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1202
1203       // convert co global phi eta  
1204       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1205       geta = ieta + nColsFeeModule*(nSupMod%2);
1206
1207       // save cell time and energy
1208       cellTime[geta][gphi] = cellTimeT;
1209       cellEnergy[geta][gphi] = cellAmp;
1210
1211     }
1212   }
1213
1214   Int_t nClusterTrig, nClusterTrigCut;
1215   UShort_t *cellAddrs;
1216   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1217   Float_t clsPos[3] = {0.,0.,0.};
1218
1219   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1220   {
1221     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1222     if(!cluster || !cluster->IsEMCAL()) continue;
1223
1224     // get cluster cells
1225     nCell = cluster->GetNCells();
1226
1227     // get cluster energy
1228     clsE = cluster->E();
1229
1230     // get cluster position
1231     cluster->GetPosition(clsPos);
1232     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1233     clsEta = clsPosVec.Eta();
1234     clsPhi = clsPosVec.Phi();
1235
1236     // get the cell addresses
1237     cellAddrs = cluster->GetCellsAbsId();
1238
1239     // check if the cluster contains cell, that was marked as triggered
1240     nClusterTrig = 0;
1241     nClusterTrigCut = 0;
1242
1243     // loop the cells to check, if cluser in acceptance
1244     // any cluster with a cell outside acceptance is not considered
1245     for( iCell = 0; iCell < nCell; iCell++ )
1246     {
1247      // check hot cell
1248      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]); 
1249
1250       // get cell position
1251       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1252       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1253
1254       // convert co global phi eta
1255       gphi = iphi + nRowsFeeModule*(nSupMod/2);
1256       geta = ieta + nColsFeeModule*(nSupMod%2);
1257
1258       if( cellTime[geta][gphi] > 0. ){ 
1259         clusterTime += cellTime[geta][gphi];
1260         gCell++;
1261       }
1262
1263       // get corresponding FALTRO
1264       fphi = gphi / 2;
1265       feta = geta / 2;
1266
1267       // try to match with a triggered
1268       if( ftriggers[feta][fphi]==1)
1269       {  nClusterTrig++;
1270       }
1271       if( ftriggersCut[feta][fphi]==1)
1272       { nClusterTrigCut++;
1273       }
1274
1275     } // cells
1276
1277
1278     if( gCell > 0 ) 
1279       clusterTime = clusterTime / (Double_t)gCell;
1280     // fix the reconstriction code time 100ns jumps
1281     if( fESD->GetBunchCrossNumber() % 4 < 2 )
1282       clusterTime -= 0.0000001;
1283
1284     //fClsETime->Fill(clsE,clusterTime);
1285     //fClsEBftTrigCut->Fill(clsE);
1286
1287     if(nClusterTrig>0){
1288       //fClsETime1->Fill(clsE,clusterTime);
1289     }
1290
1291     if(nClusterTrig>0){
1292       cluster->SetChi2(1);
1293       //fClsEAftTrigCut1->Fill(clsE);                                               
1294     }
1295
1296     if(nClusterTrigCut>0){
1297       cluster->SetChi2(2);
1298       //fClsEAftTrigCut2->Fill(clsE);
1299     }
1300
1301     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1302     {
1303       cluster->SetChi2(3);
1304       //fClsEAftTrigCut3->Fill(clsE);
1305     }
1306
1307     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1308     {
1309       // cluster->SetChi2(4);
1310       //fClsEAftTrigCut4->Fill(clsE);
1311     }
1312     if(nClusterTrigCut<1)
1313     {
1314       cluster->SetChi2(0);
1315
1316       //fClsEAftTrigCut->Fill(clsE);
1317     }
1318
1319   } // clusters
1320 }
1321
1322
1323
1324
1325
1326
1327