1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Class for heavy-flavour electron v2 with EMCal triggered events
17 // Author: Denise Godoy
25 #include "THnSparse.h"
26 #include "TLorentzVector.h"
30 #include "AliAnalysisTask.h"
31 #include "AliAnalysisManager.h"
33 #include "AliESDEvent.h"
34 #include "AliESDHandler.h"
35 #include "AliAODEvent.h"
36 #include "AliAODHandler.h"
38 #include "AliAnalysisTaskFlowTPCEMCalEP.h"
39 #include "TGeoGlobalMagField.h"
41 #include "AliAnalysisTaskSE.h"
42 #include "TRefArray.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"
54 #include "TGeoManager.h"
58 #include "AliEMCALTrack.h"
61 #include "AliKFParticle.h"
62 #include "AliKFVertex.h"
64 #include "AliMCEventHandler.h"
65 #include "AliMCEvent.h"
66 #include "AliMCParticle.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"
80 #include "AliEventplane.h"
81 #include "AliCentrality.h"
83 ClassImp(AliAnalysisTaskFlowTPCEMCalEP)
84 //________________________________________________________________________
85 AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP(const char *name)
86 : AliAnalysisTaskSE(name)
92 ,fIdentifiedAsOutInz(kFALSE)
93 ,fPassTheEventCut(kFALSE)
94 ,fRejectKinkMother(kFALSE)
100 ,fOpeningAngleCut(0.1)
115 ,fTrackPtBefTrkCuts(0)
116 ,fTrackPtAftTrkCuts(0)
146 for(Int_t k = 0; k < 6; k++) {
153 fPID = new AliHFEpid("hfePid");
154 fTrackCuts = new AliESDtrackCuts();
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());
166 //________________________________________________________________________
167 AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
168 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
174 ,fIdentifiedAsOutInz(kFALSE)
175 ,fPassTheEventCut(kFALSE)
176 ,fRejectKinkMother(kFALSE)
182 ,fOpeningAngleCut(0.1)
197 ,fTrackPtBefTrkCuts(0)
198 ,fTrackPtAftTrkCuts(0)
226 //Default constructor
228 for(Int_t k = 0; k < 6; k++) {
235 fPID = new AliHFEpid("hfePid");
237 fTrackCuts = new AliESDtrackCuts();
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());
249 //_________________________________________
251 AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
261 //_________________________________________
263 void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
266 //Called for each event
268 // create pointer to event
269 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
271 printf("ERROR: fESD not available\n");
276 AliError("HFE cuts not available");
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());
286 if(fIsMC)fMC = MCEvent();
287 AliStack* stack = NULL;
288 if(fIsMC && fMC) stack = fMC->Stack();
290 Int_t fNOtrks = fESD->GetNumberOfTracks();
291 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
293 Double_t pVtxZ = -999;
294 pVtxZ = pVtx->GetZ();
296 if(TMath::Abs(pVtxZ)>10) return;
299 if(fNOtrks<2) return;
301 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
303 AliDebug(1, "Using default PID Response");
304 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
307 fPID->SetPIDResponse(pidResponse);
309 fCFM->SetRecEventInfo(fESD);
312 AliCentrality *centrality = fESD->GetCentrality();
313 cent = centrality->GetCentralityPercentile("V0M");
317 if(cent>=20 && cent<=40) kCent = 0;
318 if(cent>=30 && cent<=50) kCent = 1;
322 Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
323 if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
324 fevPlaneV0A->Fill(evPlaneV0A,cent);
326 Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
327 if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
328 fevPlaneV0C->Fill(evPlaneV0C,cent);
330 AliEventplane* esdTPCep = fESD->GetEventplane();
331 TVector2 *standardQ = 0x0;
332 Double_t qx = -999., qy = -999.;
333 standardQ = esdTPCep->GetQVector();
334 if(!standardQ)return;
339 TVector2 qVectorfortrack;
340 qVectorfortrack.Set(qx,qy);
341 Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
342 fevPlaneTPC->Fill(evPlaneTPC,cent);
344 TVector2 *qsub1a = esdTPCep->GetQsub1();
345 TVector2 *qsub2a = esdTPCep->GetQsub2();
346 Double_t evPlaneResTPC = -999.;
349 evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
352 fTPCsubEPres->Fill(evPlaneResTPC,cent);
354 Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
355 fEPres->Fill(evPlaneRes);
357 if(fIsMC && fMC && stack && cent>20 && cent<40){
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;
366 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
368 if(!MChijing)iHijing = 0;
370 if(fPDG==111)fPi0Weight->Fill(pTMC,iHijing);//pi0
371 if(fPDG==221)fEtaWeight->Fill(pTMC,iHijing);//eta
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-
377 Int_t idMother = particle->GetFirstMother();
379 TParticle *mother = stack->Particle(idMother);
380 int motherPDG = mother->GetPdgCode();
381 if(fPDG==22 && motherPDG!=111 && motherPDG!=221)fGammaWeight->Fill(pTMC,iHijing);//gamma
389 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
390 AliESDtrack* track = fESD->GetTrack(iTracks);
392 printf("ERROR: Could not receive track %d\n", iTracks);
396 if(TMath::Abs(track->Eta())>0.7) continue;
398 fTrackPtBefTrkCuts->Fill(track->Pt());
400 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
402 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
403 if(track->GetKinkIndex(0) != 0) continue;
406 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
408 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
410 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
412 fTrackPtAftTrkCuts->Fill(track->Pt());
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.;
420 Int_t clsId = track->GetEMCALcluster();
422 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
423 if(cluster && cluster->IsEMCAL()){
425 m20 = cluster->GetM20();
426 m02 = cluster->GetM02();
432 dEdx = track->GetTPCsignal();
435 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
436 fdEdxBef->Fill(p,dEdx);
437 fTPCnsigma->Fill(p,fTPCnSigma);
439 //Remove electron candidate from the event plane
440 Float_t evPlaneCorrTPC = evPlaneTPC;
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;
449 Bool_t fFlagPhotonicElec = kFALSE;
450 Bool_t fFlagPhotonicElecBCG = kFALSE;
452 SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG);
454 Double_t corr[12]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneCorrTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C), static_cast<Double_t>(fFlagPhotonicElec), static_cast<Double_t>(fFlagPhotonicElecBCG),m02,m20};
458 Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
459 Int_t whichPart = -99;
460 Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
461 Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
462 Double_t weight = 1.;
463 Double_t Dweight = 1.;
466 Bool_t pi0Decay= kFALSE;
467 Bool_t etaDecay= kFALSE;
470 Double_t phiD = -999.,phie = -999.,phieRec = -999.,ptD = -999.,pte = -999.,pteRec = -999.;
472 if(fIsMC && fMC && stack){
473 Int_t label = track->GetLabel();
475 TParticle *particle = stack->Particle(label);
477 partPDG = particle->GetPdgCode();
478 partPt = particle->Pt();
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
485 MChijing = fMC->IsFromBGEvent(label);
488 if(!MChijing) iHijing = 0; // 0 if enhanced sample
490 Int_t idMother = particle->GetFirstMother();
492 TParticle *mother = stack->Particle(idMother);
493 motherPt = mother->Pt();
494 motherPDG = mother->GetPdgCode();
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();
505 Double_t cocktail[9]={phiD,phie,phieRec,ptD,pte,pteRec,evPlaneV0A, static_cast<Double_t>(iHijing), static_cast<Double_t>(kCent)};
506 fCocktail->Fill(cocktail);
508 //*************************************************************
509 if (motherPDG==22) whichFirstMother = 3; //gamma
510 if (motherPDG==111) whichFirstMother = 2; //pi0
511 if (motherPDG==221) whichFirstMother = 1; //eta
513 Int_t idSecondMother = particle->GetSecondMother();
514 if (idSecondMother>0){
515 TParticle *secondMother = stack->Particle(idSecondMother);
516 secondMotherPt = secondMother->Pt();
517 secondMotherPDG = secondMother->GetPdgCode();
519 if (secondMotherPDG==111) whichSecondMother = 2; //pi0
520 if (secondMotherPDG==221) whichSecondMother = 1; //eta
522 Int_t idThirdMother = secondMother->GetFirstMother();
523 if (idThirdMother>0){
524 TParticle *thirdMother = stack->Particle(idThirdMother);
525 thirdMotherPt = thirdMother->Pt();
526 thirdMotherPDG = thirdMother->GetPdgCode();
528 if (thirdMotherPDG==221) whichThirdMother = 1; //eta
533 if (cent>30 && cent<50 && TMath::Abs(partPDG)==11){
535 if (motherPDG==421 || TMath::Abs(motherPDG)==411){ // D
536 Double_t phi_D = -99., Deltaphi_De=-99.;
537 phi_D = mother->Phi();
538 Deltaphi_De = phi_D - phi;
540 fD0_e->Fill(pt,Deltaphi_De);
542 Dweight = GetDweight(0,motherPt);
543 if(iHijing==1) Dweight = 1.0;
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);
552 if (motherPDG==421){ // D0
554 Dweight = GetDweight(1,motherPt);
555 if(iHijing==1) Dweight = 1.0;
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);
564 if (motherPDG==411){ // D+
565 Dweight = GetDweight(2,motherPt);
566 if(iHijing==1) Dweight = 1.0;
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);
575 if (motherPDG==-411){ //D-
576 Dweight = GetDweight(3,motherPt);
577 if(iHijing==1) Dweight = 1.0;
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);
589 if (TMath::Abs(partPDG)==11 && motherPDG==111 && secondMotherPDG!=221){ //not eta -> pi0 -> e
590 weight = GetPi0weight(motherPt);
593 if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG!=221){ //not eta -> pi0 -> gamma -> e
594 weight = GetPi0weight(secondMotherPt);
599 if (TMath::Abs(partPDG)==11 && motherPDG==221){ //eta -> e
600 weight = GetEtaweight(motherPt);
603 if (TMath::Abs(partPDG)==11 && motherPDG==111 && secondMotherPDG==221){ //eta -> pi0 -> e
604 weight = GetEtaweight(secondMotherPt);
607 if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==221){ //eta -> gamma -> e
608 weight = GetEtaweight(secondMotherPt);
611 if (TMath::Abs(partPDG)==11 && motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG==221){ //eta -> pi0 -> gamma -> e
612 weight = GetEtaweight(thirdMotherPt);
616 if (cent>30 && cent<50 && fTPCnSigma>-1 && fTPCnSigma<3 && EovP>1 && EovP<1.3 && (motherPDG==22 || motherPDG==111 || motherPDG==221)){
617 if(iHijing==1) weight = 1.0;
620 fTot_pi0e->Fill(partPt,weight);
621 if(fFlagPhotonicElec) fPhot_pi0e->Fill(partPt,weight);
622 if(fFlagPhotonicElecBCG) fPhotBCG_pi0e->Fill(partPt,weight);
626 fTot_etae->Fill(partPt,weight);
627 if(fFlagPhotonicElec) fPhot_etae->Fill(partPt,weight);
628 if(fFlagPhotonicElecBCG) fPhotBCG_etae->Fill(partPt,weight);
634 Double_t mc[15]={EovP,fTPCnSigma,partPt, static_cast<Double_t>(fFlagPhotonicElec), static_cast<Double_t>(fFlagPhotonicElecBCG), static_cast<Double_t>(whichPart),cent,pt, static_cast<Double_t>(whichFirstMother), static_cast<Double_t>(whichSecondMother), static_cast<Double_t>(whichThirdMother), static_cast<Double_t>(iHijing),motherPt,secondMotherPt,thirdMotherPt};
636 // fMCphotoElecPt->Fill(mc);
637 if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// mother = gamma, pi0, eta
644 if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
648 AliHFEpidObject hfetrack;
649 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
650 hfetrack.SetRecTrack(track);
652 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
654 Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
655 fChargPartV2->Fill(corrV2);
657 if(fTPCnSigma >= -0.5){
658 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
659 feTPCV2->Fill(correctedV2);
662 if(pidpassed==0) continue;
664 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
666 feV2->Fill(correctedV2);
668 fTrkEovPAft->Fill(pt,EovP);
669 fdEdxAft->Fill(p,dEdx);
671 if(fFlagPhotonicElec){
672 fphoteV2->Fill(correctedV2);
673 fPhotoElecPt->Fill(pt);
676 if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
679 PostData(1, fOutputList);
681 //_________________________________________
682 void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
685 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
687 printf("+++++ MC Data available");
689 //--------Initialize PID
690 fPID->SetHasMCData(fIsMC);
692 if(!fPID->GetNumberOfPIDdetectors())
694 fPID->AddDetector("TPC", 0);
695 fPID->AddDetector("EMCAL", 1);
698 fPID->SortDetectors();
699 fPIDqa = new AliHFEpidQAmanager();
700 fPIDqa->Initialize(fPID);
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);
710 AliWarning("Cuts not available. Default cuts will be used");
711 fCuts = new AliHFEcuts;
712 fCuts->CreateStandardCuts();
714 fCuts->Initialize(fCFM);
716 //---------Output Tlist
717 fOutputList = new TList();
718 fOutputList->SetOwner();
719 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
721 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
722 fOutputList->Add(fNoEvents);
724 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
725 fOutputList->Add(fTrkpt);
727 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
728 fOutputList->Add(fTrackPtBefTrkCuts);
730 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
731 fOutputList->Add(fTrackPtAftTrkCuts);
733 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
734 fOutputList->Add(fTPCnsigma);
736 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
737 fOutputList->Add(fTrkEovPBef);
739 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
740 fOutputList->Add(fTrkEovPAft);
742 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
743 fOutputList->Add(fdEdxBef);
745 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
746 fOutputList->Add(fdEdxAft);
748 fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
749 fOutputList->Add(fInvmassLS);
751 fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
752 fOutputList->Add(fInvmassULS);
754 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
755 fOutputList->Add(fOpeningAngleLS);
757 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
758 fOutputList->Add(fOpeningAngleULS);
760 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
761 fOutputList->Add(fPhotoElecPt);
763 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
764 fOutputList->Add(fSemiInclElecPt);
766 fCent = new TH1F("fCent","Centrality",100,0,100) ;
767 fOutputList->Add(fCent);
769 fevPlaneV0A = new TH2F("fevPlaneV0A","V0A EP",100,0,TMath::Pi(),90,0,90);
770 fOutputList->Add(fevPlaneV0A);
772 fevPlaneV0C = new TH2F("fevPlaneV0C","V0C EP",100,0,TMath::Pi(),90,0,90);
773 fOutputList->Add(fevPlaneV0C);
775 fevPlaneTPC = new TH2F("fevPlaneTPC","TPC EP",100,0,TMath::Pi(),90,0,90);
776 fOutputList->Add(fevPlaneTPC);
778 fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
779 fOutputList->Add(fTPCsubEPres);
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);
787 //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG,m20,m02
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);
792 fOutputList->Add(fCorr);
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);
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);
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);
810 fOutputList->Add(fChargPartV2);
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);
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};
820 Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
821 Double_t xmaxv7[15]={3,3.5,50,2,2,4,90,50,4,4,4,2,50,50,50};
822 fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
823 fOutputList->Add(fMCphotoElecPt);
825 fGammaWeight = new TH2F("fGammaWeight", "Gamma weight",100,0,50,3,-1,2);
826 fOutputList->Add(fGammaWeight);
828 fPi0Weight = new TH2F("fPi0Weight", "Pi0 weight",100,0,50,3,-1,2);
829 fOutputList->Add(fPi0Weight);
831 fEtaWeight = new TH2F("fEtaWeight", "Eta weight",100,0,50,3,-1,2);
832 fOutputList->Add(fEtaWeight);
834 fD0Weight = new TH2F("fD0Weight", "D0 weight",100,0,50,3,-1,2);
835 fOutputList->Add(fD0Weight);
837 fDplusWeight = new TH2F("fDplusWeight", "D+ weight",100,0,50,3,-1,2);
838 fOutputList->Add(fDplusWeight);
840 fDminusWeight = new TH2F("fDminusWeight", "D- weight",100,0,50,3,-1,2);
841 fOutputList->Add(fDminusWeight);
843 fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
844 fOutputList->Add(fD0_e);
846 for(Int_t k = 0; k < 6; k++) {
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);
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);
858 fOutputList->Add(fDe[k]);
859 fOutputList->Add(fD0e[k]);
860 fOutputList->Add(fDpluse[k]);
861 fOutputList->Add(fDminuse[k]);
866 double bin_v2[9] = {2,2.5,3,4,6,8,10,13,18};
868 fTot_pi0e = new TH1F("fTot_pi0e","fTot_pi0e",nbin_v2,bin_v2);
869 fOutputList->Add(fTot_pi0e);
871 fPhot_pi0e = new TH1F("fPhot_pi0e","fPhot_pi0e",nbin_v2,bin_v2);
872 fOutputList->Add(fPhot_pi0e);
874 fPhotBCG_pi0e = new TH1F("fPhotBCG_pi0e","fPhotBCG_pi0e",nbin_v2,bin_v2);
875 fOutputList->Add(fPhotBCG_pi0e);
877 fTot_etae = new TH1F("fTot_etae","fTot_etae",nbin_v2,bin_v2);
878 fOutputList->Add(fTot_etae);
880 fPhot_etae = new TH1F("fPhot_etae","fPhot_etae",nbin_v2,bin_v2);
881 fOutputList->Add(fPhot_etae);
883 fPhotBCG_etae = new TH1F("fPhotBCG_etae","fPhotBCG_etae",nbin_v2,bin_v2);
884 fOutputList->Add(fPhotBCG_etae);
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);
894 PostData(1,fOutputList);
897 //________________________________________________________________________
898 void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
900 // Info("Terminate");
901 AliAnalysisTaskSE::Terminate();
904 //________________________________________________________________________
905 Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
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;
912 //_________________________________________
913 void AliAnalysisTaskFlowTPCEMCalEP::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec, Bool_t &fFlagPhotonicElecBCG)
915 //Identify non-heavy flavour electrons using Invariant mass method
917 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
918 fTrackCuts->SetRequireTPCRefit(kTRUE);
919 fTrackCuts->SetRequireITSRefit(kTRUE);
920 fTrackCuts->SetEtaRange(-0.7,0.7);
921 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
922 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
923 fTrackCuts->SetMinNClustersTPC(100);
925 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
927 Bool_t flagPhotonicElec = kFALSE;
928 Bool_t flagPhotonicElecBCG = kFALSE;
930 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
932 if(jTracks==iTracks) continue;
934 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
936 printf("ERROR: Could not receive track %d\n", jTracks);
940 Double_t dEdxAsso = -999., ptAsso=-999.;
941 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
942 Double_t openingAngle = -999., mass=999., width = -999;
944 dEdxAsso = trackAsso->GetTPCsignal();
945 ptAsso = trackAsso->Pt();
946 Int_t chargeAsso = trackAsso->Charge();
947 Int_t charge = track->Charge();
949 if(ptAsso <0.5) continue;
950 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
951 if(dEdxAsso <65 || dEdxAsso>100) continue;
953 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
954 if(charge>0) fPDGe1 = -11;
955 if(chargeAsso>0) fPDGe2 = -11;
957 if(charge == chargeAsso) fFlagLS = kTRUE;
958 if(charge != chargeAsso) fFlagULS = kTRUE;
960 AliKFParticle ge1(*track, fPDGe1);
961 AliKFParticle ge2(*trackAsso, fPDGe2);
962 AliKFParticle recg(ge1, ge2);
964 if(recg.GetNDF()<1) continue;
965 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
966 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
968 AliKFVertex primV(*pVtx);
970 recg.SetProductionVertex(primV);
972 recg.SetMassConstraint(0,0.0001);
974 openingAngle = ge1.GetAngle(ge2);
975 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
976 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
978 if(openingAngle > fOpeningAngleCut) continue;
980 recg.GetMass(mass,width);
982 if(fFlagLS) fInvmassLS->Fill(mass);
983 if(fFlagULS) fInvmassULS->Fill(mass);
985 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
986 if(mass<fInvmassCut && fFlagLS && !flagPhotonicElecBCG) flagPhotonicElecBCG = kTRUE;
989 fFlagPhotonicElec = flagPhotonicElec;
990 fFlagPhotonicElecBCG = flagPhotonicElecBCG;
993 //_________________________________________
994 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
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);
1001 return cos2DeltaPhi;
1003 //_________________________________________
1004 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
1007 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
1008 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
1012 //_________________________________________
1013 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT) const
1016 double weight = 1.0;
1018 if(mcPi0pT>0.0 && mcPi0pT<5.0)
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));
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);
1031 //_________________________________________
1032 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT) const
1035 double weight = 1.0;
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));
1042 //_________________________________________
1043 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDweight(Int_t whichD, Double_t mcDpT) const
1046 double weight = 1.0;
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-