]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalEP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix of sigmaZ for crossing tracklets from Alex
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskFlowTPCEMCalEP.cxx
CommitLineData
1d08b6e8 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
914042c2 38#include "AliAnalysisTaskFlowTPCEMCalEP.h"
1d08b6e8 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
3a94d6fe 64#include "AliMCEventHandler.h"
65#include "AliMCEvent.h"
66#include "AliMCParticle.h"
67#include "AliStack.h"
68
1d08b6e8 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
914042c2 83ClassImp(AliAnalysisTaskFlowTPCEMCalEP)
1d08b6e8 84//________________________________________________________________________
914042c2 85AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP(const char *name)
1d08b6e8 86 : AliAnalysisTaskSE(name)
87 ,fESD(0)
3a94d6fe 88 ,fMC(0)
1d08b6e8 89 ,fOutputList(0)
90 ,fTrackCuts(0)
91 ,fCuts(0)
92 ,fIdentifiedAsOutInz(kFALSE)
93 ,fPassTheEventCut(kFALSE)
94 ,fRejectKinkMother(kFALSE)
3a94d6fe 95 ,fIsMC(kFALSE)
1d08b6e8 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)
3a94d6fe 114 ,fMCphotoElecPt(0)
1d08b6e8 115 ,fTrackPtBefTrkCuts(0)
116 ,fTrackPtAftTrkCuts(0)
117 ,fTPCnsigma(0)
118 ,fCent(0)
44fdf7ac 119 ,fevPlaneV0A(0)
120 ,fevPlaneV0C(0)
121 ,fevPlaneTPC(0)
1d08b6e8 122 ,fTPCsubEPres(0)
123 ,fEPres(0)
124 ,fCorr(0)
125 ,feTPCV2(0)
126 ,feV2(0)
127 ,fphoteV2(0)
128 ,fChargPartV2(0)
6522b3b5 129 ,fGammaWeight(0)
130 ,fPi0Weight(0)
131 ,fEtaWeight(0)
def484dc 132 ,fD0Weight(0)
133 ,fDplusWeight(0)
134 ,fDminusWeight(0)
937b2ff2 135 ,fD0_e(0)
1c424962 136 ,fTot_pi0e(0)
137 ,fPhot_pi0e(0)
138 ,fPhotBCG_pi0e(0)
139 ,fTot_etae(0)
140 ,fPhot_etae(0)
141 ,fPhotBCG_etae(0)
535b11a0 142 ,fCocktail(0)
1d08b6e8 143{
144 //Named constructor
937b2ff2 145
146 for(Int_t k = 0; k < 6; k++) {
147 fDe[k]= NULL;
148 fD0e[k]= NULL;
149 fDpluse[k]= NULL;
150 fDminuse[k]= NULL;
151 }
152
1d08b6e8 153 fPID = new AliHFEpid("hfePid");
154 fTrackCuts = new AliESDtrackCuts();
155
156 // Define input and output slots here
157 // Input slot #0 works with a TChain
158 DefineInput(0, TChain::Class());
159 // Output slot #0 id reserved by the base class for AOD
160 // Output slot #1 writes into a TH1 container
161 // DefineOutput(1, TH1I::Class());
162 DefineOutput(1, TList::Class());
163 // DefineOutput(3, TTree::Class());
164}
165
166//________________________________________________________________________
914042c2 167AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
1d08b6e8 168 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
169 ,fESD(0)
3a94d6fe 170 ,fMC(0)
1d08b6e8 171 ,fOutputList(0)
172 ,fTrackCuts(0)
173 ,fCuts(0)
174 ,fIdentifiedAsOutInz(kFALSE)
175 ,fPassTheEventCut(kFALSE)
176 ,fRejectKinkMother(kFALSE)
3a94d6fe 177 ,fIsMC(kFALSE)
1d08b6e8 178 ,fVz(0.0)
179 ,fCFM(0)
180 ,fPID(0)
181 ,fPIDqa(0)
182 ,fOpeningAngleCut(0.1)
183 ,fInvmassCut(0.01)
184 ,fNoEvents(0)
185 ,fTrkpt(0)
186 ,fTrkEovPBef(0)
187 ,fTrkEovPAft(0)
188 ,fdEdxBef(0)
189 ,fdEdxAft(0)
190 ,fInvmassLS(0)
191 ,fInvmassULS(0)
192 ,fOpeningAngleLS(0)
193 ,fOpeningAngleULS(0)
194 ,fPhotoElecPt(0)
195 ,fSemiInclElecPt(0)
3a94d6fe 196 ,fMCphotoElecPt(0)
1d08b6e8 197 ,fTrackPtBefTrkCuts(0)
198 ,fTrackPtAftTrkCuts(0)
199 ,fTPCnsigma(0)
200 ,fCent(0)
44fdf7ac 201 ,fevPlaneV0A(0)
202 ,fevPlaneV0C(0)
203 ,fevPlaneTPC(0)
1d08b6e8 204 ,fTPCsubEPres(0)
205 ,fEPres(0)
206 ,fCorr(0)
207 ,feTPCV2(0)
208 ,feV2(0)
209 ,fphoteV2(0)
210 ,fChargPartV2(0)
6522b3b5 211 ,fGammaWeight(0)
212 ,fPi0Weight(0)
213 ,fEtaWeight(0)
def484dc 214 ,fD0Weight(0)
215 ,fDplusWeight(0)
216 ,fDminusWeight(0)
937b2ff2 217 ,fD0_e(0)
1c424962 218 ,fTot_pi0e(0)
219 ,fPhot_pi0e(0)
220 ,fPhotBCG_pi0e(0)
221 ,fTot_etae(0)
222 ,fPhot_etae(0)
223 ,fPhotBCG_etae(0)
535b11a0 224 ,fCocktail(0)
1d08b6e8 225{
226 //Default constructor
937b2ff2 227
228 for(Int_t k = 0; k < 6; k++) {
229 fDe[k]= NULL;
230 fD0e[k]= NULL;
231 fDpluse[k]= NULL;
232 fDminuse[k]= NULL;
233 }
234
1d08b6e8 235 fPID = new AliHFEpid("hfePid");
236
237 fTrackCuts = new AliESDtrackCuts();
238
239 // Constructor
240 // Define input and output slots here
241 // Input slot #0 works with a TChain
242 DefineInput(0, TChain::Class());
243 // Output slot #0 id reserved by the base class for AOD
244 // Output slot #1 writes into a TH1 container
245 // DefineOutput(1, TH1I::Class());
246 DefineOutput(1, TList::Class());
247 //DefineOutput(3, TTree::Class());
248}
249//_________________________________________
250
914042c2 251AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
1d08b6e8 252{
253 //Destructor
254
255 delete fOutputList;
256 delete fPID;
257 delete fCFM;
258 delete fPIDqa;
259 delete fTrackCuts;
260}
261//_________________________________________
262
914042c2 263void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
1d08b6e8 264{
265 //Main loop
266 //Called for each event
267
268 // create pointer to event
269 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
270 if (!fESD) {
271 printf("ERROR: fESD not available\n");
272 return;
273 }
274
275 if(!fCuts){
276 AliError("HFE cuts not available");
277 return;
278 }
279
280 if(!fPID->IsInitialized()){
281 // Initialize PID with the given run number
282 AliWarning("PID not initialised, get from Run no");
283 fPID->InitializePID(fESD->GetRunNumber());
284 }
285
3a94d6fe 286 if(fIsMC)fMC = MCEvent();
287 AliStack* stack = NULL;
288 if(fIsMC && fMC) stack = fMC->Stack();
1d08b6e8 289
290 Int_t fNOtrks = fESD->GetNumberOfTracks();
291 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
292
293 Double_t pVtxZ = -999;
294 pVtxZ = pVtx->GetZ();
295
296 if(TMath::Abs(pVtxZ)>10) return;
297 fNoEvents->Fill(0);
298
299 if(fNOtrks<2) return;
300
301 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
302 if(!pidResponse){
303 AliDebug(1, "Using default PID Response");
304 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
305 }
306
307 fPID->SetPIDResponse(pidResponse);
308
309 fCFM->SetRecEventInfo(fESD);
310
311 Float_t cent = -1.;
312 AliCentrality *centrality = fESD->GetCentrality();
313 cent = centrality->GetCentralityPercentile("V0M");
314 fCent->Fill(cent);
535b11a0 315
316 int kCent = 3;
317 if(cent>=20 && cent<=40) kCent = 0;
318 if(cent>=30 && cent<=50) kCent = 1;
6522b3b5 319
1d08b6e8 320 //Event planes
321
322 Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
323 if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
5945f829 324 fevPlaneV0A->Fill(evPlaneV0A,cent);
1d08b6e8 325
326 Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
327 if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
5945f829 328 fevPlaneV0C->Fill(evPlaneV0C,cent);
1d08b6e8 329
330 AliEventplane* esdTPCep = fESD->GetEventplane();
44fdf7ac 331 TVector2 *standardQ = 0x0;
1d08b6e8 332 Double_t qx = -999., qy = -999.;
44fdf7ac 333 standardQ = esdTPCep->GetQVector();
ee89e744 334 if(!standardQ)return;
335
336 qx = standardQ->X();
337 qy = standardQ->Y();
338
1d08b6e8 339 TVector2 qVectorfortrack;
340 qVectorfortrack.Set(qx,qy);
341 Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
5945f829 342 fevPlaneTPC->Fill(evPlaneTPC,cent);
44fdf7ac 343
1d08b6e8 344 TVector2 *qsub1a = esdTPCep->GetQsub1();
345 TVector2 *qsub2a = esdTPCep->GetQsub2();
346 Double_t evPlaneResTPC = -999.;
347 if(qsub1a && qsub2a)
348 {
349 evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
350 }
351
352 fTPCsubEPres->Fill(evPlaneResTPC,cent);
353
354 Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
355 fEPres->Fill(evPlaneRes);
356
535b11a0 357 if(fIsMC && fMC && stack && cent>20 && cent<40){
6522b3b5 358 Int_t nParticles = stack->GetNtrack();
359 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
360 TParticle* particle = stack->Particle(iParticle);
361 int fPDG = particle->GetPdgCode();
362 double pTMC = particle->Pt();
363 double etaMC = particle->Eta();
364 if(fabs(etaMC)>0.7)continue;
365
366 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
367 int iHijing = 1;
368 if(!MChijing)iHijing = 0;
369
370 if(fPDG==111)fPi0Weight->Fill(pTMC,iHijing);//pi0
371 if(fPDG==221)fEtaWeight->Fill(pTMC,iHijing);//eta
def484dc 372 if(fPDG==421)fD0Weight->Fill(pTMC,iHijing);//D0
373 if(fPDG==411)fDplusWeight->Fill(pTMC,iHijing);//D+
374 if(fPDG==-411)fDminusWeight->Fill(pTMC,iHijing);//D-
375
376
6522b3b5 377 Int_t idMother = particle->GetFirstMother();
378 if (idMother>0){
379 TParticle *mother = stack->Particle(idMother);
380 int motherPDG = mother->GetPdgCode();
381 if(fPDG==22 && motherPDG!=111 && motherPDG!=221)fGammaWeight->Fill(pTMC,iHijing);//gamma
382 }
383
384 }
385 }
386
387
1d08b6e8 388 // Track loop
389 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
390 AliESDtrack* track = fESD->GetTrack(iTracks);
391 if (!track) {
392 printf("ERROR: Could not receive track %d\n", iTracks);
393 continue;
394 }
395
396 if(TMath::Abs(track->Eta())>0.7) continue;
397
398 fTrackPtBefTrkCuts->Fill(track->Pt());
b1a1875f 399
1d08b6e8 400 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
401
1d08b6e8 402 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
403 if(track->GetKinkIndex(0) != 0) continue;
404 }
405
1d08b6e8 406 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
407
1d08b6e8 408 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
409
1d08b6e8 410 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
411
412 fTrackPtAftTrkCuts->Fill(track->Pt());
413
1c424962 414 Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999., wclsE = -999., wEovP = -999., m02= -999., m20= -999.;
32d22164 415
1d08b6e8 416 pt = track->Pt();
44fdf7ac 417 if(pt<2) continue;
1d08b6e8 418 fTrkpt->Fill(pt);
3a94d6fe 419
1d08b6e8 420 Int_t clsId = track->GetEMCALcluster();
421 if (clsId>0){
422 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
423 if(cluster && cluster->IsEMCAL()){
424 clsE = cluster->E();
1c424962 425 m20 = cluster->GetM20();
426 m02 = cluster->GetM02();
1d08b6e8 427 }
428 }
b1a1875f 429
1d08b6e8 430 p = track->P();
431 phi = track->Phi();
432 dEdx = track->GetTPCsignal();
433 EovP = clsE/p;
b1a1875f 434 wEovP = wclsE/p;
1d08b6e8 435 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
436 fdEdxBef->Fill(p,dEdx);
437 fTPCnsigma->Fill(p,fTPCnSigma);
32d22164 438
439 //Remove electron candidate from the event plane
6522b3b5 440 Float_t evPlaneCorrTPC = evPlaneTPC;
32d22164 441 if(dEdx>70 && dEdx<90){
442 Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
443 Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
444 TVector2 newQVectorfortrack;
445 newQVectorfortrack.Set(qX,qY);
446 evPlaneCorrTPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
447 }
448
06d74a7b 449 Bool_t fFlagPhotonicElec = kFALSE;
450 Bool_t fFlagPhotonicElecBCG = kFALSE;
5945f829 451
06d74a7b 452 SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG);
453
1c424962 454 Double_t corr[12]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneCorrTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG,m02,m20};
1d08b6e8 455 fCorr->Fill(corr);
06d74a7b 456
32d22164 457
458 Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
6522b3b5 459 Int_t whichPart = -99;
32d22164 460 Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
1c424962 461 Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
462 Double_t weight = 1.;
def484dc 463 Double_t Dweight = 1.;
32d22164 464 Bool_t MChijing;
465
1c424962 466 Bool_t pi0Decay= kFALSE;
467 Bool_t etaDecay= kFALSE;
535b11a0 468
469
470 Double_t phiD = -999.,phie = -999.,phieRec = -999.,ptD = -999.,pte = -999.,pteRec = -999.;
471
32d22164 472 if(fIsMC && fMC && stack){
06d74a7b 473 Int_t label = track->GetLabel();
474 if(label>0){
475 TParticle *particle = stack->Particle(label);
476 if(particle){
32d22164 477 partPDG = particle->GetPdgCode();
478 partPt = particle->Pt();
479
6522b3b5 480 if (TMath::Abs(partPDG)==11) whichPart = 0; //electron
481 if (partPDG==22) whichPart = 3; //gamma
482 if (partPDG==111) whichPart = 2; //pi0
483 if (partPDG==221) whichPart = 1; //eta
484
32d22164 485 MChijing = fMC->IsFromBGEvent(label);
06d74a7b 486
32d22164 487 int iHijing = 1;
488 if(!MChijing) iHijing = 0; // 0 if enhanced sample
489
06d74a7b 490 Int_t idMother = particle->GetFirstMother();
491 if (idMother>0){
32d22164 492 TParticle *mother = stack->Particle(idMother);
493 motherPt = mother->Pt();
494 motherPDG = mother->GetPdgCode();
495
535b11a0 496 //*************** for cocktail ********************************
497 if (TMath::Abs(partPDG)==11 && (motherPDG==421 || TMath::Abs(motherPDG)==411)){
498 pteRec = track->Pt();
499 phieRec = track->Phi();
500 phie = particle->Phi();
501 pte = particle->Pt();
502 phiD = mother->Phi();
503 ptD = mother->Pt();
504
505 Double_t cocktail[9]={phiD,phie,phieRec,ptD,pte,pteRec,evPlaneV0A,iHijing,kCent};
506 fCocktail->Fill(cocktail);
507 }
508 //*************************************************************
32d22164 509 if (motherPDG==22) whichFirstMother = 3; //gamma
510 if (motherPDG==111) whichFirstMother = 2; //pi0
511 if (motherPDG==221) whichFirstMother = 1; //eta
512
513 Int_t idSecondMother = particle->GetSecondMother();
514 if (idSecondMother>0){
515 TParticle *secondMother = stack->Particle(idSecondMother);
516 secondMotherPt = secondMother->Pt();
517 secondMotherPDG = secondMother->GetPdgCode();
518
519 if (secondMotherPDG==111) whichSecondMother = 2; //pi0
520 if (secondMotherPDG==221) whichSecondMother = 1; //eta
521
522 Int_t idThirdMother = secondMother->GetFirstMother();
523 if (idThirdMother>0){
524 TParticle *thirdMother = stack->Particle(idThirdMother);
525 thirdMotherPt = thirdMother->Pt();
526 thirdMotherPDG = thirdMother->GetPdgCode();
527
528 if (thirdMotherPDG==221) whichThirdMother = 1; //eta
529 }
530 }
def484dc 531
532
533 if (cent>30 && cent<50 && TMath::Abs(partPDG)==11){
937b2ff2 534
535 if (motherPDG==421 || TMath::Abs(motherPDG)==411){ // D
908103bf 536 Double_t phi_D = -99., Deltaphi_De=-99.;
537 phi_D = mother->Phi();
538 Deltaphi_De = phi_D - phi;
539
540 fD0_e->Fill(pt,Deltaphi_De);
541
937b2ff2 542 Dweight = GetDweight(0,motherPt);
543 if(iHijing==1) Dweight = 1.0;
544
545 if (motherPt>=2 && motherPt<3) fDe[0]->Fill(pt,Dweight);
546 if (motherPt>=3 && motherPt<4) fDe[1]->Fill(pt,Dweight);
547 if (motherPt>=4 && motherPt<6) fDe[2]->Fill(pt,Dweight);
548 if (motherPt>=6 && motherPt<8) fDe[3]->Fill(pt,Dweight);
549 if (motherPt>=8 && motherPt<12) fDe[4]->Fill(pt,Dweight);
550 if (motherPt>=12 && motherPt<16) fDe[5]->Fill(pt,Dweight);
551 }
552 if (motherPDG==421){ // D0
908103bf 553
937b2ff2 554 Dweight = GetDweight(1,motherPt);
555 if(iHijing==1) Dweight = 1.0;
556
557 if (motherPt>=2 && motherPt<3) fD0e[0]->Fill(pt,Dweight);
558 if (motherPt>=3 && motherPt<4) fD0e[1]->Fill(pt,Dweight);
559 if (motherPt>=4 && motherPt<6) fD0e[2]->Fill(pt,Dweight);
560 if (motherPt>=6 && motherPt<8) fD0e[3]->Fill(pt,Dweight);
561 if (motherPt>=8 && motherPt<12) fD0e[4]->Fill(pt,Dweight);
562 if (motherPt>=12 && motherPt<16) fD0e[5]->Fill(pt,Dweight);
563 }
564 if (motherPDG==411){ // D+
565 Dweight = GetDweight(2,motherPt);
566 if(iHijing==1) Dweight = 1.0;
567
568 if (motherPt>=2 && motherPt<3) fDpluse[0]->Fill(pt,Dweight);
569 if (motherPt>=3 && motherPt<4) fDpluse[1]->Fill(pt,Dweight);
570 if (motherPt>=4 && motherPt<6) fDpluse[2]->Fill(pt,Dweight);
571 if (motherPt>=6 && motherPt<8) fDpluse[3]->Fill(pt,Dweight);
572 if (motherPt>=8 && motherPt<12) fDpluse[4]->Fill(pt,Dweight);
573 if (motherPt>=12 && motherPt<16) fDpluse[5]->Fill(pt,Dweight);
574 }
575 if (motherPDG==-411){ //D-
576 Dweight = GetDweight(3,motherPt);
577 if(iHijing==1) Dweight = 1.0;
578
579 if (motherPt>=2 && motherPt<3) fDminuse[0]->Fill(pt,Dweight);
580 if (motherPt>=3 && motherPt<4) fDminuse[1]->Fill(pt,Dweight);
581 if (motherPt>=4 && motherPt<6) fDminuse[2]->Fill(pt,Dweight);
582 if (motherPt>=6 && motherPt<8) fDminuse[3]->Fill(pt,Dweight);
583 if (motherPt>=8 && motherPt<12) fDminuse[4]->Fill(pt,Dweight);
584 if (motherPt>=12 && motherPt<16) fDminuse[5]->Fill(pt,Dweight);
585 }
def484dc 586 }
587
1c424962 588 //pi0 decay
589 if (TMath::Abs(partPDG)==11 && motherPDG==111 && secondMotherPDG!=221){ //not eta -> pi0 -> e
590 weight = GetPi0weight(motherPt);
591 pi0Decay = kTRUE;
592 }
593 if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG!=221){ //not eta -> pi0 -> gamma -> e
594 weight = GetPi0weight(secondMotherPt);
595 pi0Decay = kTRUE;
596 }
597
598 //eta decay
599 if (TMath::Abs(partPDG)==11 && motherPDG==221){ //eta -> e
600 weight = GetEtaweight(motherPt);
601 etaDecay = kTRUE;
602 }
603 if (TMath::Abs(partPDG)==11 && motherPDG==111 && secondMotherPDG==221){ //eta -> pi0 -> e
604 weight = GetEtaweight(secondMotherPt);
605 etaDecay = kTRUE;
606 }
607 if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==221){ //eta -> gamma -> e
608 weight = GetEtaweight(secondMotherPt);
609 etaDecay = kTRUE;
610 }
611 if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG==221){ //eta -> pi0 -> gamma -> e
612 weight = GetEtaweight(thirdMotherPt);
613 etaDecay = kTRUE;
614 }
615
def484dc 616 if (cent>30 && cent<50 && fTPCnSigma>-1 && fTPCnSigma<3 && EovP>1 && EovP<1.3 && (motherPDG==22 || motherPDG==111 || motherPDG==221)){
1c424962 617 if(iHijing==1) weight = 1.0;
618
619 if (pi0Decay){
620 fTot_pi0e->Fill(partPt,weight);
621 if(fFlagPhotonicElec) fPhot_pi0e->Fill(partPt,weight);
622 if(fFlagPhotonicElecBCG) fPhotBCG_pi0e->Fill(partPt,weight);
623 }
624
625 if (etaDecay){
626 fTot_etae->Fill(partPt,weight);
627 if(fFlagPhotonicElec) fPhot_etae->Fill(partPt,weight);
628 if(fFlagPhotonicElecBCG) fPhotBCG_etae->Fill(partPt,weight);
629 }
630
631 }
632
633
6522b3b5 634 Double_t mc[15]={EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,whichFirstMother,whichSecondMother,whichThirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt};
1c424962 635
636// fMCphotoElecPt->Fill(mc);
637 if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// mother = gamma, pi0, eta
06d74a7b 638 }
639 }
640 }
641 }
06d74a7b 642
1d08b6e8 643
644 if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
645 Int_t pidpassed = 1;
646
647 //--- track accepted
648 AliHFEpidObject hfetrack;
649 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
650 hfetrack.SetRecTrack(track);
651 hfetrack.SetPbPb();
652 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
653
44fdf7ac 654 Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
1d08b6e8 655 fChargPartV2->Fill(corrV2);
b3bd01b3 656
657 if(fTPCnSigma >= -0.5){
32d22164 658 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
b3bd01b3 659 feTPCV2->Fill(correctedV2);
660 }
661
662 if(pidpassed==0) continue;
663
32d22164 664 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
1d08b6e8 665
1d08b6e8 666 feV2->Fill(correctedV2);
667
668 fTrkEovPAft->Fill(pt,EovP);
669 fdEdxAft->Fill(p,dEdx);
670
3a94d6fe 671 if(fFlagPhotonicElec){
672 fphoteV2->Fill(correctedV2);
673 fPhotoElecPt->Fill(pt);
674 }
675
676 if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
677
1d08b6e8 678 }
679 PostData(1, fOutputList);
680}
681//_________________________________________
914042c2 682void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
1d08b6e8 683{
3a94d6fe 684 //--- Check MC
685 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
686 fIsMC = kTRUE;
687 printf("+++++ MC Data available");
688 }
1d08b6e8 689 //--------Initialize PID
3a94d6fe 690 fPID->SetHasMCData(fIsMC);
691
1d08b6e8 692 if(!fPID->GetNumberOfPIDdetectors())
693 {
694 fPID->AddDetector("TPC", 0);
695 fPID->AddDetector("EMCAL", 1);
696 }
697
698 fPID->SortDetectors();
699 fPIDqa = new AliHFEpidQAmanager();
700 fPIDqa->Initialize(fPID);
701
702 //--------Initialize correction Framework and Cuts
703 fCFM = new AliCFManager;
704 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
705 fCFM->SetNStepParticle(kNcutSteps);
706 for(Int_t istep = 0; istep < kNcutSteps; istep++)
707 fCFM->SetParticleCutsList(istep, NULL);
708
709 if(!fCuts){
710 AliWarning("Cuts not available. Default cuts will be used");
711 fCuts = new AliHFEcuts;
712 fCuts->CreateStandardCuts();
713 }
714 fCuts->Initialize(fCFM);
715
716 //---------Output Tlist
717 fOutputList = new TList();
718 fOutputList->SetOwner();
719 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
720
721 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
722 fOutputList->Add(fNoEvents);
723
724 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
725 fOutputList->Add(fTrkpt);
726
727 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
728 fOutputList->Add(fTrackPtBefTrkCuts);
729
730 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
731 fOutputList->Add(fTrackPtAftTrkCuts);
732
733 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
734 fOutputList->Add(fTPCnsigma);
735
736 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
737 fOutputList->Add(fTrkEovPBef);
738
739 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
740 fOutputList->Add(fTrkEovPAft);
741
742 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
743 fOutputList->Add(fdEdxBef);
744
745 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
746 fOutputList->Add(fdEdxAft);
747
748 fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
749 fOutputList->Add(fInvmassLS);
750
751 fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
752 fOutputList->Add(fInvmassULS);
753
754 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
755 fOutputList->Add(fOpeningAngleLS);
756
757 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
758 fOutputList->Add(fOpeningAngleULS);
759
760 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
761 fOutputList->Add(fPhotoElecPt);
762
763 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
764 fOutputList->Add(fSemiInclElecPt);
765
766 fCent = new TH1F("fCent","Centrality",100,0,100) ;
767 fOutputList->Add(fCent);
768
5945f829 769 fevPlaneV0A = new TH2F("fevPlaneV0A","V0A EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 770 fOutputList->Add(fevPlaneV0A);
771
5945f829 772 fevPlaneV0C = new TH2F("fevPlaneV0C","V0C EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 773 fOutputList->Add(fevPlaneV0C);
774
5945f829 775 fevPlaneTPC = new TH2F("fevPlaneTPC","TPC EP",100,0,TMath::Pi(),90,0,90);
44fdf7ac 776 fOutputList->Add(fevPlaneTPC);
777
1d08b6e8 778 fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
779 fOutputList->Add(fTPCsubEPres);
780
781 Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
782 Double_t xminv1[4]={-1,-1,-1,0};
783 Double_t xmaxv1[4]={1,1,1,90};
784 fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
785 fOutputList->Add(fEPres);
786
06d74a7b 787 //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG,m20,m02
1c424962 788 Int_t binsv2[12]={100,100,90,100,100,100,100,100,3,3,100,100};
789 Double_t xminv2[12]={0,-3.5,0,0,0,0,0,0,-1,-1,0,0};
790 Double_t xmaxv2[12]={2*TMath::Pi(),3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi(),2,2,1,1};
791 fCorr = new THnSparseD ("fCorr","Correlations",12,binsv2,xminv2,xmaxv2);
1d08b6e8 792 fOutputList->Add(fCorr);
793
794 Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
795 Double_t xminv3[5]={0,0,-1,-1,-1};
796 Double_t xmaxv3[5]={90,50,1,1,1};
797 feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
798 fOutputList->Add(feV2);
799
800 Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
801 Double_t xminv4[5]={0,0,-1,-1,-1};
802 Double_t xmaxv4[5]={90,50,1,1,1};
803 fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
804 fOutputList->Add(fphoteV2);
805
44fdf7ac 806 Int_t binsv5[7]={100,90,100,100,100,100,100}; // phi, cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
807 Double_t xminv5[7]={0,0,0,0,-1,-1,-1};
808 Double_t xmaxv5[7]={2*TMath::Pi(),90,50,3,1,1,1};
809 fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",7,binsv5,xminv5,xmaxv5);
1d08b6e8 810 fOutputList->Add(fChargPartV2);
811
812 Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
813 Double_t xminv6[5]={0,0,-1,-1,-1};
814 Double_t xmaxv6[5]={90,50,1,1,1};
815 feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
816 fOutputList->Add(feTPCV2);
817
6522b3b5 818 //EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,firstMother,secondMother,thirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt
819 Int_t binsv7[15]={100,100,100,3,3,5,90,100,5,5,5,3,100,100,100};
32d22164 820 Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
6522b3b5 821 Double_t xmaxv7[15]={3,3.5,50,2,2,4,90,50,4,4,4,2,50,50,50};
32d22164 822 fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
3a94d6fe 823 fOutputList->Add(fMCphotoElecPt);
6522b3b5 824
825 fGammaWeight = new TH2F("fGammaWeight", "Gamma weight",100,0,50,3,-1,2);
826 fOutputList->Add(fGammaWeight);
827
828 fPi0Weight = new TH2F("fPi0Weight", "Pi0 weight",100,0,50,3,-1,2);
829 fOutputList->Add(fPi0Weight);
830
831 fEtaWeight = new TH2F("fEtaWeight", "Eta weight",100,0,50,3,-1,2);
832 fOutputList->Add(fEtaWeight);
1c424962 833
def484dc 834 fD0Weight = new TH2F("fD0Weight", "D0 weight",100,0,50,3,-1,2);
835 fOutputList->Add(fD0Weight);
836
837 fDplusWeight = new TH2F("fDplusWeight", "D+ weight",100,0,50,3,-1,2);
838 fOutputList->Add(fDplusWeight);
839
840 fDminusWeight = new TH2F("fDminusWeight", "D- weight",100,0,50,3,-1,2);
841 fOutputList->Add(fDminusWeight);
842
908103bf 843 fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
937b2ff2 844 fOutputList->Add(fD0_e);
def484dc 845
937b2ff2 846 for(Int_t k = 0; k < 6; k++) {
847
848 TString De_name = Form("fDe%d",k);
849 TString D0e_name = Form("fD0e%d",k);
850 TString Dpluse_name = Form("fDpluse%d",k);
851 TString Dminuse_name = Form("fDminuse%d",k);
852
853 fDe[k] = new TH1F((const char*)De_name,"",100,0,50);
854 fD0e[k] = new TH1F((const char*)D0e_name,"",100,0,50);
855 fDpluse[k] = new TH1F((const char*)Dpluse_name,"",100,0,50);
856 fDminuse[k] = new TH1F((const char*)Dminuse_name,"",100,0,50);
857
858 fOutputList->Add(fDe[k]);
859 fOutputList->Add(fD0e[k]);
860 fOutputList->Add(fDpluse[k]);
861 fOutputList->Add(fDminuse[k]);
862
863 }
535b11a0 864
1c424962 865 int nbin_v2 = 8;
866 double bin_v2[9] = {2,2.5,3,4,6,8,10,13,18};
867
868 fTot_pi0e = new TH1F("fTot_pi0e","fTot_pi0e",nbin_v2,bin_v2);
869 fOutputList->Add(fTot_pi0e);
6522b3b5 870
1c424962 871 fPhot_pi0e = new TH1F("fPhot_pi0e","fPhot_pi0e",nbin_v2,bin_v2);
872 fOutputList->Add(fPhot_pi0e);
873
874 fPhotBCG_pi0e = new TH1F("fPhotBCG_pi0e","fPhotBCG_pi0e",nbin_v2,bin_v2);
875 fOutputList->Add(fPhotBCG_pi0e);
876
877 fTot_etae = new TH1F("fTot_etae","fTot_etae",nbin_v2,bin_v2);
878 fOutputList->Add(fTot_etae);
879
880 fPhot_etae = new TH1F("fPhot_etae","fPhot_etae",nbin_v2,bin_v2);
881 fOutputList->Add(fPhot_etae);
882
883 fPhotBCG_etae = new TH1F("fPhotBCG_etae","fPhotBCG_etae",nbin_v2,bin_v2);
884 fOutputList->Add(fPhotBCG_etae);
535b11a0 885
886 //phiD,phie,phieRec,ptD,pte,pteRec,evPlaneV0,iHijing,kCent
887 Int_t binsv8[9]={100,100,100,100,100,100,100,4,4};
888 Double_t xminv8[9]={0,0,0,0,0,0,0,-2,-2};
889 Double_t xmaxv8[9]={2*TMath::Pi(),2*TMath::Pi(),2*TMath::Pi(),50,50,50,TMath::Pi(),2,2};
890 fCocktail = new THnSparseD ("fCocktail","Correlations",9,binsv8,xminv8,xmaxv8);
891 fOutputList->Add(fCocktail);
892
1c424962 893
1d08b6e8 894 PostData(1,fOutputList);
895}
896
897//________________________________________________________________________
914042c2 898void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
1d08b6e8 899{
900 // Info("Terminate");
901 AliAnalysisTaskSE::Terminate();
902}
903
904//________________________________________________________________________
914042c2 905Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1d08b6e8 906{
907 // Check single track cuts for a given cut step
908 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
909 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
910 return kTRUE;
911}
912//_________________________________________
914042c2 913void AliAnalysisTaskFlowTPCEMCalEP::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec, Bool_t &fFlagPhotonicElecBCG)
1d08b6e8 914{
915 //Identify non-heavy flavour electrons using Invariant mass method
916
917 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
918 fTrackCuts->SetRequireTPCRefit(kTRUE);
5945f829 919 fTrackCuts->SetRequireITSRefit(kTRUE);
1d08b6e8 920 fTrackCuts->SetEtaRange(-0.7,0.7);
921 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
922 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
923 fTrackCuts->SetMinNClustersTPC(100);
924
925 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
926
927 Bool_t flagPhotonicElec = kFALSE;
06d74a7b 928 Bool_t flagPhotonicElecBCG = kFALSE;
1d08b6e8 929
b1a1875f 930 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
931
5945f829 932 if(jTracks==iTracks) continue;
b1a1875f 933
1d08b6e8 934 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
935 if (!trackAsso) {
936 printf("ERROR: Could not receive track %d\n", jTracks);
937 continue;
938 }
939
5945f829 940 Double_t dEdxAsso = -999., ptAsso=-999.;
06d74a7b 941 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
942 Double_t openingAngle = -999., mass=999., width = -999;
1d08b6e8 943
944 dEdxAsso = trackAsso->GetTPCsignal();
945 ptAsso = trackAsso->Pt();
946 Int_t chargeAsso = trackAsso->Charge();
947 Int_t charge = track->Charge();
948
06d74a7b 949 if(ptAsso <0.5) continue;
1d08b6e8 950 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
06d74a7b 951 if(dEdxAsso <65 || dEdxAsso>100) continue;
1d08b6e8 952
953 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
954 if(charge>0) fPDGe1 = -11;
955 if(chargeAsso>0) fPDGe2 = -11;
956
957 if(charge == chargeAsso) fFlagLS = kTRUE;
958 if(charge != chargeAsso) fFlagULS = kTRUE;
959
960 AliKFParticle ge1(*track, fPDGe1);
961 AliKFParticle ge2(*trackAsso, fPDGe2);
962 AliKFParticle recg(ge1, ge2);
963
964 if(recg.GetNDF()<1) continue;
965 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
966 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
967
968 AliKFVertex primV(*pVtx);
969 primV += recg;
970 recg.SetProductionVertex(primV);
971
972 recg.SetMassConstraint(0,0.0001);
973
974 openingAngle = ge1.GetAngle(ge2);
975 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
976 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
977
5945f829 978 if(openingAngle > fOpeningAngleCut) continue;
06d74a7b 979
980 recg.GetMass(mass,width);
5945f829 981
1d08b6e8 982 if(fFlagLS) fInvmassLS->Fill(mass);
983 if(fFlagULS) fInvmassULS->Fill(mass);
984
5945f829 985 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
06d74a7b 986 if(mass<fInvmassCut && fFlagLS && !flagPhotonicElecBCG) flagPhotonicElecBCG = kTRUE;
1d08b6e8 987
988 }
989 fFlagPhotonicElec = flagPhotonicElec;
06d74a7b 990 fFlagPhotonicElecBCG = flagPhotonicElecBCG;
1d08b6e8 991
992}
993//_________________________________________
914042c2 994Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
1d08b6e8 995{
996 //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
997 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
998 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
999 Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
1000
1001 return cos2DeltaPhi;
1002}
1d08b6e8 1003//_________________________________________
914042c2 1004Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
1d08b6e8 1005{
1006 //Get phi-psi_EP
1007 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
1008 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
1009
1010 return dPhi;
1011}
1c424962 1012//_________________________________________
1013Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT) const
1014{
937b2ff2 1015 //Get Pi0 weight
1c424962 1016 double weight = 1.0;
1017
1018 if(mcPi0pT>0.0 && mcPi0pT<5.0)
1019 {
1020
def484dc 1021// weight = (2.877091*mcPi0pT)/(TMath::Power(0.706963+mcPi0pT/3.179309,17.336628)*exp(-mcPi0pT));
1022 weight = (2.392024*mcPi0pT)/(TMath::Power(0.688810+mcPi0pT/3.005145,16.811845)*exp(-mcPi0pT));
1c424962 1023 }
1024 else
1025 {
def484dc 1026// weight = (0.0004*mcPi0pT)/TMath::Power(-0.176181+mcPi0pT/3.989747,5.629235);
1027 weight = (0.000186*mcPi0pT)/TMath::Power(-0.606279+mcPi0pT/3.158642,4.365540);
1c424962 1028 }
1029 return weight;
1030}
1031//_________________________________________
1032Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT) const
1033{
937b2ff2 1034 //Get eta weight
1c424962 1035 double weight = 1.0;
1036
def484dc 1037// weight = (0.818052*mcEtapT)/(TMath::Power(0.358651+mcEtapT/2.878631,9.494043));
1038 weight = (0.622703*mcEtapT)/(TMath::Power(0.323045+mcEtapT/2.736407,9.180356));
1c424962 1039
1040 return weight;
1041}
def484dc 1042//_________________________________________
937b2ff2 1043Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDweight(Int_t whichD, Double_t mcDpT) const
def484dc 1044{
937b2ff2 1045 //get D weights
def484dc 1046 double weight = 1.0;
937b2ff2 1047
1048 if (whichD == 0) weight = 0.271583*TMath::Landau(mcDpT,3.807103,1.536753,0); // D
1049 if (whichD == 1) weight = 0.300771*TMath::Landau(mcDpT,3.725771,1.496980,0); // D0
1050 if (whichD == 2) weight = 0.247280*TMath::Landau(mcDpT,3.746811,1.607551,0); // D+
1051 if (whichD == 3) weight = 0.249410*TMath::Landau(mcDpT,3.611508,1.632196,0); //D-
1052
def484dc 1053 return weight;
1054}