updated
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskElecV2.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 v2 with EMCal triggered events
17 // Author: Denise Godoy
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 "AliAnalysisTaskElecV2.h"
39 #include "TGeoGlobalMagField.h"
40 #include "AliLog.h"
41 #include "AliAnalysisTaskSE.h"
42 #include "TRefArray.h"
43 #include "TVector.h"
44 #include "AliESDInputHandler.h"
45 #include "AliESDpid.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliPhysicsSelection.h"
48 #include "AliESDCaloCluster.h"
49 #include "AliAODCaloCluster.h"
50 #include "AliEMCALRecoUtils.h"
51 #include "AliEMCALGeometry.h"
52 #include "AliGeomManager.h"
53 #include "stdio.h"
54 #include "TGeoManager.h"
55 #include "iostream"
56 #include "fstream"
57
58 #include "AliEMCALTrack.h"
59 #include "AliMagF.h"
60
61 #include "AliKFParticle.h"
62 #include "AliKFVertex.h"
63
64 #include "AliMCEventHandler.h"
65 #include "AliMCEvent.h"
66 #include "AliMCParticle.h"
67 #include "AliStack.h"
68
69 #include "AliPID.h"
70 #include "AliPIDResponse.h"
71 #include "AliHFEcontainer.h"
72 #include "AliHFEcuts.h"
73 #include "AliHFEpid.h"
74 #include "AliHFEpidBase.h"
75 #include "AliHFEpidQAmanager.h"
76 #include "AliHFEtools.h"
77 #include "AliCFContainer.h"
78 #include "AliCFManager.h"
79
80 #include "AliEventplane.h"
81 #include "AliCentrality.h"
82
83 ClassImp(AliAnalysisTaskElecV2)
84 //________________________________________________________________________
85 AliAnalysisTaskElecV2::AliAnalysisTaskElecV2(const char *name) 
86   : AliAnalysisTaskSE(name)
87   ,fESD(0)
88   ,fMC(0)
89   ,fOutputList(0)
90   ,fTrackCuts(0)
91   ,fCuts(0)
92   ,fIdentifiedAsOutInz(kFALSE)
93   ,fPassTheEventCut(kFALSE)
94   ,fRejectKinkMother(kFALSE)
95   ,fIsMC(kFALSE)
96   ,fVz(0.0)
97   ,fCFM(0)      
98   ,fPID(0)
99   ,fPIDqa(0)           
100   ,fOpeningAngleCut(0.1)
101   ,fInvmassCut(0.01)    
102   ,fNoEvents(0)
103   ,fTrkpt(0)
104   ,fTrkEovPBef(0)        
105   ,fTrkEovPAft(0)       
106   ,fdEdxBef(0)   
107   ,fdEdxAft(0)   
108   ,fInvmassLS(0)                
109   ,fInvmassULS(0)               
110   ,fOpeningAngleLS(0)   
111   ,fOpeningAngleULS(0)  
112   ,fPhotoElecPt(0)
113   ,fSemiInclElecPt(0)
114   ,fMCphotoElecPt(0)
115   ,fTrackPtBefTrkCuts(0)         
116   ,fTrackPtAftTrkCuts(0)
117   ,fTPCnsigma(0)
118   ,fCent(0)
119   ,fevPlaneV0A(0)
120   ,fevPlaneV0C(0)
121   ,fevPlaneTPC(0)
122   ,fTPCsubEPres(0)
123   ,fEPres(0)
124   ,fCorr(0)
125   ,feTPCV2(0)
126   ,feV2(0)
127   ,fphoteV2(0)
128   ,fChargPartV2(0)
129 {
130   //Named constructor
131   
132   fPID = new AliHFEpid("hfePid");
133   fTrackCuts = new AliESDtrackCuts();
134   
135   // Define input and output slots here
136   // Input slot #0 works with a TChain
137   DefineInput(0, TChain::Class());
138   // Output slot #0 id reserved by the base class for AOD
139   // Output slot #1 writes into a TH1 container
140   // DefineOutput(1, TH1I::Class());
141   DefineOutput(1, TList::Class());
142   //  DefineOutput(3, TTree::Class());
143 }
144
145 //________________________________________________________________________
146 AliAnalysisTaskElecV2::AliAnalysisTaskElecV2() 
147   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
148   ,fESD(0)
149   ,fMC(0)
150   ,fOutputList(0)
151   ,fTrackCuts(0)
152   ,fCuts(0)
153   ,fIdentifiedAsOutInz(kFALSE)
154   ,fPassTheEventCut(kFALSE)
155   ,fRejectKinkMother(kFALSE)
156   ,fIsMC(kFALSE)
157   ,fVz(0.0)
158   ,fCFM(0)      
159   ,fPID(0)       
160   ,fPIDqa(0)           
161   ,fOpeningAngleCut(0.1)
162   ,fInvmassCut(0.01)    
163   ,fNoEvents(0)
164   ,fTrkpt(0)
165   ,fTrkEovPBef(0)        
166   ,fTrkEovPAft(0)        
167   ,fdEdxBef(0)   
168   ,fdEdxAft(0)   
169   ,fInvmassLS(0)                
170   ,fInvmassULS(0)               
171   ,fOpeningAngleLS(0)   
172   ,fOpeningAngleULS(0)  
173   ,fPhotoElecPt(0)
174   ,fSemiInclElecPt(0)
175   ,fMCphotoElecPt(0)
176   ,fTrackPtBefTrkCuts(0)         
177   ,fTrackPtAftTrkCuts(0)                  
178   ,fTPCnsigma(0)
179   ,fCent(0)
180   ,fevPlaneV0A(0)
181   ,fevPlaneV0C(0)
182   ,fevPlaneTPC(0)
183   ,fTPCsubEPres(0)
184   ,fEPres(0)
185   ,fCorr(0)
186   ,feTPCV2(0)
187   ,feV2(0)
188   ,fphoteV2(0)
189   ,fChargPartV2(0)
190 {
191         //Default constructor
192         fPID = new AliHFEpid("hfePid");
193
194         fTrackCuts = new AliESDtrackCuts();
195         
196         // Constructor
197         // Define input and output slots here
198         // Input slot #0 works with a TChain
199         DefineInput(0, TChain::Class());
200         // Output slot #0 id reserved by the base class for AOD
201         // Output slot #1 writes into a TH1 container
202         // DefineOutput(1, TH1I::Class());
203         DefineOutput(1, TList::Class());
204         //DefineOutput(3, TTree::Class());
205 }
206 //_________________________________________
207
208 AliAnalysisTaskElecV2::~AliAnalysisTaskElecV2()
209 {
210   //Destructor 
211   
212   delete fOutputList;
213   delete fPID;
214   delete fCFM;
215   delete fPIDqa;
216   delete fTrackCuts;
217 }
218 //_________________________________________
219
220 void AliAnalysisTaskElecV2::UserExec(Option_t*)
221 {
222   //Main loop
223   //Called for each event
224   
225   // create pointer to event
226   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
227   if (!fESD) {
228     printf("ERROR: fESD not available\n");
229     return;
230   }
231   
232   if(!fCuts){
233     AliError("HFE cuts not available");
234     return;
235   }
236   
237   if(!fPID->IsInitialized()){ 
238     // Initialize PID with the given run number
239     AliWarning("PID not initialised, get from Run no");
240     fPID->InitializePID(fESD->GetRunNumber());
241   }
242  
243   if(fIsMC)fMC = MCEvent();
244   AliStack* stack = NULL;
245   if(fIsMC && fMC) stack = fMC->Stack();
246  
247   Int_t fNOtrks =  fESD->GetNumberOfTracks();
248   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
249   
250   Double_t pVtxZ = -999;
251   pVtxZ = pVtx->GetZ();
252   
253   if(TMath::Abs(pVtxZ)>10) return;
254   fNoEvents->Fill(0);
255   
256   if(fNOtrks<2) return;
257   
258   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
259   if(!pidResponse){
260     AliDebug(1, "Using default PID Response");
261     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); 
262   }
263   
264   fPID->SetPIDResponse(pidResponse);
265   
266   fCFM->SetRecEventInfo(fESD);
267   
268   Float_t cent = -1.;
269   AliCentrality *centrality = fESD->GetCentrality(); 
270   cent = centrality->GetCentralityPercentile("V0M");
271   fCent->Fill(cent);
272   
273   if(cent>90.) return;
274         
275   //Event planes
276   
277   Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
278   if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
279   fevPlaneV0A->Fill(evPlaneV0A,cent);
280   
281   Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
282   if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
283   fevPlaneV0C->Fill(evPlaneV0C,cent);
284   
285   AliEventplane* esdTPCep = fESD->GetEventplane();
286   TVector2 *standardQ = 0x0;
287   Double_t qx = -999., qy = -999.;
288   standardQ = esdTPCep->GetQVector(); 
289   if(!standardQ)return;
290  
291   qx = standardQ->X();
292   qy = standardQ->Y();
293   
294   TVector2 qVectorfortrack;
295   qVectorfortrack.Set(qx,qy);
296   Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
297   fevPlaneTPC->Fill(evPlaneTPC,cent);
298   
299   TVector2 *qsub1a = esdTPCep->GetQsub1();
300   TVector2 *qsub2a = esdTPCep->GetQsub2();
301   Double_t evPlaneResTPC = -999.;
302   if(qsub1a && qsub2a)
303   {
304           evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
305   }
306     
307   fTPCsubEPres->Fill(evPlaneResTPC,cent);
308   
309   Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
310   fEPres->Fill(evPlaneRes);
311   
312   // Track loop 
313   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
314     AliESDtrack* track = fESD->GetTrack(iTracks);
315     if (!track) {
316       printf("ERROR: Could not receive track %d\n", iTracks);
317       continue;
318     }
319     
320     if(TMath::Abs(track->Eta())>0.7) continue;
321     
322     fTrackPtBefTrkCuts->Fill(track->Pt());              
323     
324     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
325     
326     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
327       if(track->GetKinkIndex(0) != 0) continue;
328     } 
329     
330     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
331     
332     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
333     
334     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
335     
336     fTrackPtAftTrkCuts->Fill(track->Pt());              
337     
338     Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999., clsPhi=-999., clsEta=-999., wclsE = -999., wEovP = -999.;
339     
340     
341     pt = track->Pt();
342     if(pt<2) continue;
343     fTrkpt->Fill(pt);
344         
345     Int_t clsId = track->GetEMCALcluster();
346     if (clsId>0){
347       AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
348       if(cluster && cluster->IsEMCAL()){
349         clsE = cluster->E();
350         Float_t clusterPosition[3]={0,0,0};
351         cluster->GetPosition(clusterPosition);
352         TVector3 clsPosVec(clusterPosition[0],clusterPosition[1],clusterPosition[2]);
353         clsPhi = clsPosVec.Phi();
354         clsEta = clsPosVec.Eta();
355         
356         wclsE = GetclusterE(iTracks, clsPhi, clsEta);
357         
358       }
359     }
360
361
362     
363     p = track->P();
364     phi = track->Phi();
365     dEdx = track->GetTPCsignal();
366     EovP = clsE/p;
367     wEovP = wclsE/p;
368     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
369     fdEdxBef->Fill(p,dEdx);
370     fTPCnsigma->Fill(p,fTPCnSigma);
371        
372     Bool_t fFlagPhotonicElec = kFALSE, fFlagLS=kFALSE, fFlagULS=kFALSE;
373     Double_t openingAngle = -999., mass=999., width = -999;
374             
375     SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec,fFlagLS,fFlagULS,openingAngle,mass,width);
376     
377     Double_t corr[14]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagLS,fFlagULS,openingAngle,mass,width};
378     fCorr->Fill(corr);
379        
380     if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
381     Int_t pidpassed = 1;
382     
383     //--- track accepted
384     AliHFEpidObject hfetrack;
385     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
386     hfetrack.SetRecTrack(track);
387     hfetrack.SetPbPb();
388     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
389     
390     Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
391     fChargPartV2->Fill(corrV2); 
392
393     if(fTPCnSigma >= -0.5){
394         Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track); 
395         Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track); 
396         TVector2 newQVectorfortrack;
397         newQVectorfortrack.Set(qX,qY);
398         Double_t corrV2TPC = -999.;
399         corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; 
400
401         Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
402
403         feTPCV2->Fill(correctedV2);
404     }
405     
406     if(pidpassed==0) continue;
407     
408     Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track); 
409     Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track); 
410     TVector2 newQVectorfortrack;
411     newQVectorfortrack.Set(qX,qY);
412     Double_t corrV2TPC = -999.;
413     corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; 
414         
415     Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
416     
417     feV2->Fill(correctedV2);
418     
419     fTrkEovPAft->Fill(pt,EovP);
420     fdEdxAft->Fill(p,dEdx);
421     
422     openingAngle = -999.;
423     mass=999.;
424     width = -999;
425     fFlagPhotonicElec = kFALSE;
426     fFlagLS=kFALSE;
427     fFlagULS=kFALSE;
428     SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec,fFlagLS,fFlagULS,openingAngle,mass,width);
429     
430     if(fIsMC && fMC && stack){
431       Int_t label = track->GetLabel();
432       if(label>0){
433         TParticle *particle = stack->Particle(label);
434         if(particle){
435             Int_t partPDG = particle->GetPdgCode();
436             Double_t partPt = particle->Pt();
437             Int_t IsElec = 0;
438             if (TMath::Abs(partPDG)==11) IsElec = 1;
439           
440             Int_t idMother = particle->GetFirstMother();
441             if (idMother>0){
442             TParticle *mother = stack->Particle(idMother);
443             Int_t motherPDG = mother->GetPdgCode();
444             
445             Double_t mc[9]={partPt,fFlagPhotonicElec,fFlagLS,fFlagULS,IsElec,openingAngle,mass,cent,pt};
446             
447             if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// gamma, pi0, eta
448           }
449         }
450       }
451     }
452     
453     
454   if(fFlagPhotonicElec){
455     fphoteV2->Fill(correctedV2);
456     fPhotoElecPt->Fill(pt);
457   }
458     
459   if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
460
461  }
462  PostData(1, fOutputList);
463 }
464 //_________________________________________
465 void AliAnalysisTaskElecV2::UserCreateOutputObjects()
466 {
467   //--- Check MC
468   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
469     fIsMC = kTRUE;
470     printf("+++++ MC Data available");
471   }
472   //--------Initialize PID
473   fPID->SetHasMCData(fIsMC);
474   
475   if(!fPID->GetNumberOfPIDdetectors()) 
476     {
477       fPID->AddDetector("TPC", 0);
478       fPID->AddDetector("EMCAL", 1);
479     }
480   
481   fPID->SortDetectors(); 
482   fPIDqa = new AliHFEpidQAmanager();
483   fPIDqa->Initialize(fPID);
484   
485   //--------Initialize correction Framework and Cuts
486   fCFM = new AliCFManager;
487   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
488   fCFM->SetNStepParticle(kNcutSteps);
489   for(Int_t istep = 0; istep < kNcutSteps; istep++)
490     fCFM->SetParticleCutsList(istep, NULL);
491   
492   if(!fCuts){
493     AliWarning("Cuts not available. Default cuts will be used");
494     fCuts = new AliHFEcuts;
495     fCuts->CreateStandardCuts();
496   }
497   fCuts->Initialize(fCFM);
498   
499   //---------Output Tlist
500   fOutputList = new TList();
501   fOutputList->SetOwner();
502   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
503   
504   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
505   fOutputList->Add(fNoEvents);
506   
507   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
508   fOutputList->Add(fTrkpt);
509   
510   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
511   fOutputList->Add(fTrackPtBefTrkCuts);
512   
513   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
514   fOutputList->Add(fTrackPtAftTrkCuts);
515   
516   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
517   fOutputList->Add(fTPCnsigma);
518   
519   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
520   fOutputList->Add(fTrkEovPBef);
521   
522   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
523   fOutputList->Add(fTrkEovPAft);
524   
525   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
526   fOutputList->Add(fdEdxBef);
527   
528   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
529   fOutputList->Add(fdEdxAft);
530   
531   fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
532   fOutputList->Add(fInvmassLS);
533   
534   fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
535   fOutputList->Add(fInvmassULS);
536   
537   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
538   fOutputList->Add(fOpeningAngleLS);
539   
540   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
541   fOutputList->Add(fOpeningAngleULS);
542   
543   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
544   fOutputList->Add(fPhotoElecPt);
545   
546   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
547   fOutputList->Add(fSemiInclElecPt);
548   
549   fCent = new TH1F("fCent","Centrality",100,0,100) ;
550   fOutputList->Add(fCent);
551   
552   fevPlaneV0A = new TH2F("fevPlaneV0A","V0A EP",100,0,TMath::Pi(),90,0,90);
553   fOutputList->Add(fevPlaneV0A);
554   
555   fevPlaneV0C = new TH2F("fevPlaneV0C","V0C EP",100,0,TMath::Pi(),90,0,90);
556   fOutputList->Add(fevPlaneV0C);
557   
558   fevPlaneTPC = new TH2F("fevPlaneTPC","TPC EP",100,0,TMath::Pi(),90,0,90);
559   fOutputList->Add(fevPlaneTPC);
560     
561   fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
562   fOutputList->Add(fTPCsubEPres);
563   
564   Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
565   Double_t xminv1[4]={-1,-1,-1,0};
566   Double_t xmaxv1[4]={1,1,1,90}; 
567   fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
568   fOutputList->Add(fEPres);
569         
570   //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagLS,fFlagULS,openingAngle,mass,width
571   Int_t binsv2[14]={100,100,90,100,100,100,100,100,3,3,3,100,100,100}; 
572   Double_t xminv2[14]={0,-3.5,0,0,0,0,0,0,-1,-1,-1,0,0,-5};
573   Double_t xmaxv2[14]={2*TMath::Pi(),3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi(),2,2,2,1,0.5,5}; 
574   fCorr = new THnSparseD ("fCorr","Correlations",14,binsv2,xminv2,xmaxv2);
575   fOutputList->Add(fCorr);
576     
577   Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
578   Double_t xminv3[5]={0,0,-1,-1,-1};
579   Double_t xmaxv3[5]={90,50,1,1,1}; 
580   feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
581   fOutputList->Add(feV2);
582   
583   Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
584   Double_t xminv4[5]={0,0,-1,-1,-1};
585   Double_t xmaxv4[5]={90,50,1,1,1}; 
586   fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
587   fOutputList->Add(fphoteV2);
588   
589   Int_t binsv5[7]={100,90,100,100,100,100,100}; // phi, cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
590   Double_t xminv5[7]={0,0,0,0,-1,-1,-1};
591   Double_t xmaxv5[7]={2*TMath::Pi(),90,50,3,1,1,1}; 
592   fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",7,binsv5,xminv5,xmaxv5);
593   fOutputList->Add(fChargPartV2);
594   
595   Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
596   Double_t xminv6[5]={0,0,-1,-1,-1};
597   Double_t xmaxv6[5]={90,50,1,1,1}; 
598   feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
599   fOutputList->Add(feTPCV2);
600   
601   Int_t binsv7[9]={100,3,3,3,3,100,100,90,100}; // partPt,  fFlagPhotonicElec,  fFlagLS,  fFlagULS,  IsElec,  openingAngle,  mass,  cent,  pt
602   Double_t xminv7[9]={0,-1,-1,-1,-1,0,0,0,0};
603   Double_t xmaxv7[9]={50,2,2,2,2,1,0.5,90,50}; 
604   fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",9,binsv7,xminv7,xmaxv7);
605   fOutputList->Add(fMCphotoElecPt);
606    
607   PostData(1,fOutputList);
608 }
609
610 //________________________________________________________________________
611 void AliAnalysisTaskElecV2::Terminate(Option_t *)
612 {
613   // Info("Terminate");
614         AliAnalysisTaskSE::Terminate();
615 }
616
617 //________________________________________________________________________
618 Bool_t AliAnalysisTaskElecV2::ProcessCutStep(Int_t cutStep, AliVParticle *track)
619 {
620   // Check single track cuts for a given cut step
621   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
622   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
623   return kTRUE;
624 }
625 //_________________________________________
626 void AliAnalysisTaskElecV2::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec,Bool_t &fFlagLS,Bool_t &fFlagULS, Double_t &openingAngle,Double_t &mass,Double_t &width)
627 {
628   //Identify non-heavy flavour electrons using Invariant mass method
629   
630   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
631   fTrackCuts->SetRequireTPCRefit(kTRUE);
632   fTrackCuts->SetRequireITSRefit(kTRUE);
633   fTrackCuts->SetEtaRange(-0.7,0.7);
634   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
635   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
636   fTrackCuts->SetMinNClustersTPC(100);
637   
638   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
639   
640   Bool_t flagPhotonicElec = kFALSE;
641   
642   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
643     
644     if(jTracks==iTracks) continue;
645     
646     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
647     if (!trackAsso) {
648       printf("ERROR: Could not receive track %d\n", jTracks);
649       continue;
650     }
651     
652     Double_t dEdxAsso = -999., ptAsso=-999.;
653     
654     dEdxAsso = trackAsso->GetTPCsignal();
655     ptAsso = trackAsso->Pt();
656     Int_t chargeAsso = trackAsso->Charge();
657     Int_t charge = track->Charge();
658     
659     if(ptAsso <0.3) continue;
660     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
661     if(dEdxAsso <70 || dEdxAsso>100) continue;
662     
663     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
664     if(charge>0) fPDGe1 = -11;
665     if(chargeAsso>0) fPDGe2 = -11;
666     
667     if(charge == chargeAsso) fFlagLS = kTRUE;
668     if(charge != chargeAsso) fFlagULS = kTRUE;
669     
670     AliKFParticle ge1(*track, fPDGe1);
671     AliKFParticle ge2(*trackAsso, fPDGe2);
672     AliKFParticle recg(ge1, ge2);
673     
674     if(recg.GetNDF()<1) continue;
675     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
676     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
677     
678     AliKFVertex primV(*pVtx);
679     primV += recg;
680     recg.SetProductionVertex(primV);
681     
682     recg.SetMassConstraint(0,0.0001);
683     
684     openingAngle = ge1.GetAngle(ge2);
685     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
686     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
687     
688     recg.GetMass(mass,width);
689     
690     if(openingAngle > fOpeningAngleCut) continue;
691         
692     if(fFlagLS) fInvmassLS->Fill(mass);
693     if(fFlagULS) fInvmassULS->Fill(mass);
694           
695     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
696     
697   }
698   fFlagPhotonicElec = flagPhotonicElec;
699   
700 }
701 //_________________________________________
702 Double_t AliAnalysisTaskElecV2::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
703 {
704   //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
705   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB); 
706   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
707   Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
708   
709   return cos2DeltaPhi;
710 }
711
712 //_________________________________________
713 Double_t AliAnalysisTaskElecV2::GetDeltaPhi(Double_t phiA,Double_t phiB) const
714 {
715   //Get phi-psi_EP
716   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB); 
717   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
718   
719   return dPhi;
720 }
721 //_________________________________________
722 Double_t AliAnalysisTaskElecV2::GetclusterE(Int_t iTrack, Double_t clsPhi, Double_t clsEta) const
723 {
724   //Return E
725   for (Int_t jTracks = 0; jTracks < fESD->GetNumberOfTracks(); jTracks++){
726
727     if(jTracks==iTrack) continue;
728
729     AliESDtrack* wtrack = fESD->GetTrack(jTracks);
730     if (!wtrack) continue;
731
732     Double_t wclsPhi=-999., wclsEta=-999., dPhi=-999., dEta=-999., dR=-999., wclsE=-999.;
733
734     Int_t wclsId = wtrack->GetEMCALcluster();
735     if (wclsId>0){
736       AliESDCaloCluster *wcluster = fESD->GetCaloCluster(wclsId);
737       if(wcluster && wcluster->IsEMCAL()){
738         Float_t wclusterPosition[3]={0,0,0};
739         wcluster->GetPosition(wclusterPosition);
740         TVector3 clsPosVec(wclusterPosition[0],wclusterPosition[1],wclusterPosition[2]);
741         wclsPhi = clsPosVec.Phi();
742         wclsEta = clsPosVec.Eta();
743
744         dPhi = TMath::Abs(wclsPhi - clsPhi);
745         dEta = TMath::Abs(wclsEta - clsEta);
746         dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
747         if(dR>0.15){
748           wclsE = wcluster->E();
749           return wclsE;
750         }
751       }
752     }
753   }
754   return -999.;
755 }