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