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"
37 #include "AliVEvent.h"
39 #include "AliAnalysisTaskFlowTPCEMCalEP.h"
40 #include "TGeoGlobalMagField.h"
42 #include "AliAnalysisTaskSE.h"
43 #include "TRefArray.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"
55 #include "TGeoManager.h"
59 #include "AliEMCALTrack.h"
62 #include "AliKFParticle.h"
63 #include "AliKFVertex.h"
65 #include "AliMCEventHandler.h"
66 #include "AliMCEvent.h"
67 #include "AliMCParticle.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"
81 #include "AliEventplane.h"
82 #include "AliCentrality.h"
84 #include "AliSelectNonHFE.h"
87 ClassImp(AliAnalysisTaskFlowTPCEMCalEP)
88 //________________________________________________________________________
89 AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP(const char *name)
90 : AliAnalysisTaskSE(name)
100 ,fIdentifiedAsOutInz(kFALSE)
101 ,fPassTheEventCut(kFALSE)
102 ,fRejectKinkMother(kFALSE)
109 ,fOpeningAngleCut(0.1)
115 ,fnonHFEalgorithm("KF")
125 ,fTrackPtBefTrkCuts(0)
126 ,fTrackPtAftTrkCuts(0)
162 for(Int_t k = 0; k < 6; k++) {
169 fPID = new AliHFEpid("hfePid");
170 fTrackCuts = new AliESDtrackCuts();
171 fNonHFE = new AliSelectNonHFE();
175 // Define input and output slots here
176 // Input slot #0 works with a TChain
177 DefineInput(0, TChain::Class());
178 // Output slot #0 id reserved by the base class for AOD
179 // Output slot #1 writes into a TH1 container
180 // DefineOutput(1, TH1I::Class());
181 DefineOutput(1, TList::Class());
182 // DefineOutput(3, TTree::Class());
185 //________________________________________________________________________
186 AliAnalysisTaskFlowTPCEMCalEP::AliAnalysisTaskFlowTPCEMCalEP()
187 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
197 ,fIdentifiedAsOutInz(kFALSE)
198 ,fPassTheEventCut(kFALSE)
199 ,fRejectKinkMother(kFALSE)
206 ,fOpeningAngleCut(0.1)
212 ,fnonHFEalgorithm("KF")
222 ,fTrackPtBefTrkCuts(0)
223 ,fTrackPtAftTrkCuts(0)
258 //Default constructor
260 for(Int_t k = 0; k < 6; k++) {
267 fPID = new AliHFEpid("hfePid");
268 fTrackCuts = new AliESDtrackCuts();
269 fNonHFE = new AliSelectNonHFE();
275 // Define input and output slots here
276 // Input slot #0 works with a TChain
277 DefineInput(0, TChain::Class());
278 // Output slot #0 id reserved by the base class for AOD
279 // Output slot #1 writes into a TH1 container
280 // DefineOutput(1, TH1I::Class());
281 DefineOutput(1, TList::Class());
282 //DefineOutput(3, TTree::Class());
284 //_________________________________________
286 AliAnalysisTaskFlowTPCEMCalEP::~AliAnalysisTaskFlowTPCEMCalEP()
296 //_________________________________________
298 void AliAnalysisTaskFlowTPCEMCalEP::UserExec(Option_t*)
301 //Called for each event
303 // create pointer to event
304 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
306 printf("ERROR: fESD not available\n");
310 fVevent = dynamic_cast<AliVEvent*>(InputEvent());
312 printf("ERROR: fVEvent not available\n");
317 AliError("HFE cuts not available");
321 if(!fPID->IsInitialized()){
322 // Initialize PID with the given run number
323 AliWarning("PID not initialised, get from Run no");
324 if(fIsAOD)fPID->InitializePID(fAOD->GetRunNumber());
325 if(!fIsAOD) fPID->InitializePID(fESD->GetRunNumber());
328 if(fIsMC)fMC = MCEvent();
329 AliStack* stack = NULL;
330 if(fIsMC && fMC) stack = fMC->Stack();
332 Int_t fNOtrks =fESD->GetNumberOfTracks();
333 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
335 Double_t pVtxZ = -999;
336 pVtxZ = pVtx->GetZ();
338 if(TMath::Abs(pVtxZ)>10) return;
341 if(fNOtrks<2) return;
343 fpidResponse = fInputHandler->GetPIDResponse();
345 AliDebug(1, "Using default PID Response");
346 fpidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
349 fPID->SetPIDResponse(fpidResponse);
351 fCFM->SetRecEventInfo(fVevent);
354 AliCentrality *centrality = fESD->GetCentrality();
355 cent = centrality->GetCentralityPercentile("V0M");
358 // cout<<"TEST: "<<fInvmassCut<< " "<< fOpeningAngleCut<<" "<<fnonHFEalgorithm<<" "<<fminCent<<" "<<fmaxCent<<" here!!!!!!!!!!!" <<endl;
362 Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
363 if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
364 fevPlaneV0A->Fill(evPlaneV0A);
366 Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
367 if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
368 fevPlaneV0C->Fill(evPlaneV0C);
370 Double_t evPlaneV0 = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0",fESD,2));
371 if(evPlaneV0 > TMath::Pi()) evPlaneV0 = evPlaneV0 - TMath::Pi();
372 fevPlaneV0->Fill(evPlaneV0);
374 AliEventplane* esdTPCep = fESD->GetEventplane();
375 TVector2 *standardQ = 0x0;
376 Double_t qx = -999., qy = -999.;
377 standardQ = esdTPCep->GetQVector();
378 if(!standardQ)return;
383 TVector2 qVectorfortrack;
384 qVectorfortrack.Set(qx,qy);
385 Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
386 fevPlaneTPC->Fill(evPlaneTPC);
388 //Event plane resolutions
390 // --> 2 subevent method (only for TPC EP)
392 TVector2 *qsub1a = esdTPCep->GetQsub1();
393 TVector2 *qsub2a = esdTPCep->GetQsub2();
394 Double_t evPlaneResTPC = -999.;
395 if(qsub1a && qsub2a){
396 evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
399 fTPCsubEPres->Fill(evPlaneResTPC);
401 // --> 3 event method (V0, V0A, and V0C EP)
403 Double_t Qx2pos = -999., Qy2pos = -999., Qx2neg = -999., Qy2neg = -999., Qweight = 1;
405 for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
407 AliVParticle* vparticle = fVevent->GetTrack(iTracks);
409 printf("ERROR: Could not receive track %d\n", iTracks);
413 AliESDtrack *trackEP = dynamic_cast<AliESDtrack*>(vparticle);
416 printf("ERROR: Could not receive track %d\n", iTracks);
420 if(TMath::Abs(trackEP->Eta())>0.8 || trackEP->Pt() < 0.15 || trackEP->Pt() > 4) continue;
422 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, trackEP)) continue;
424 if (trackEP->Pt() < 2) Qweight = trackEP->Pt()/2;
426 if(trackEP->Eta()>0 && trackEP->Eta()<0.8){
427 Qx2pos += Qweight*TMath::Cos(2*trackEP->Phi());
428 Qy2pos += Qweight*TMath::Sin(2*trackEP->Phi());
430 if(trackEP->Eta()<0 && trackEP->Eta()>-0.8){
431 Qx2neg += Qweight*TMath::Cos(2*trackEP->Phi());
432 Qy2neg += Qweight*TMath::Sin(2*trackEP->Phi());
434 }//track loop only for EP
436 Double_t evPlaneTPCneg = TMath::ATan2(Qy2neg, Qx2neg)/2;
437 Double_t evPlaneTPCpos = TMath::ATan2(Qy2pos, Qx2pos)/2;
439 Double_t evPlaneRes[7]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),
440 GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0,evPlaneTPCpos),
441 GetCos2DeltaPhi(evPlaneV0,evPlaneTPCneg),GetCos2DeltaPhi(evPlaneTPCpos,evPlaneTPCneg),cent};
442 fEPres->Fill(evPlaneRes);
445 if(fIsMC && fMC && stack){
446 Int_t nParticles = stack->GetNtrack();
447 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
448 TParticle* particle = stack->Particle(iParticle);
449 int fPDG = particle->GetPdgCode();
450 double pTMC = particle->Pt();
451 double etaMC = particle->Eta();
452 if(fabs(etaMC)>0.7)continue;
454 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
456 if(!MChijing)iHijing = 0;
458 if(fPDG==111)fPi0Weight->Fill(pTMC,iHijing);//pi0
459 if(fPDG==221)fEtaWeight->Fill(pTMC,iHijing);//eta
460 if(fPDG==421)fD0Weight->Fill(pTMC,iHijing);//D0
461 if(fPDG==411)fDplusWeight->Fill(pTMC,iHijing);//D+
462 if(fPDG==-411)fDminusWeight->Fill(pTMC,iHijing);//D-
464 Int_t idMother = particle->GetFirstMother();
466 TParticle *mother = stack->Particle(idMother);
467 int motherPDG = mother->GetPdgCode();
468 if(fPDG==22 && motherPDG!=111 && motherPDG!=221)fGammaWeight->Fill(pTMC,iHijing);//gamma
473 Double_t ptRange[8] = {2, 2.5, 3, 4, 6, 8, 10, 13};
474 Double_t ptDmeson[7] = {2,3,4,6,8,12,16};
475 Double_t deltaPhiRange[7];
476 for(Int_t j=0;j<7;j++){
477 deltaPhiRange[j] = j*(TMath::Pi()/6);
481 for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++) {
483 AliVParticle* vparticle = fVevent->GetTrack(iTracks);
485 printf("ERROR: Could not receive track %d\n", iTracks);
489 AliVTrack *vtrack = dynamic_cast<AliVTrack*>(vparticle);
490 AliESDtrack *track = dynamic_cast<AliESDtrack*>(vparticle);
492 if (TMath::Abs(track->Eta())>0.7) continue;
494 fTrackPtBefTrkCuts->Fill(track->Pt());
496 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
498 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
499 if(track->GetKinkIndex(0) != 0) continue;
502 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
504 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
506 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
508 fTrackPtAftTrkCuts->Fill(track->Pt());
510 Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999.,
511 wclsE = -999., wEovP = -999., m02= -999., m20= -999.;
517 Int_t clsId = track->GetEMCALcluster();
519 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
520 if(cluster && cluster->IsEMCAL()){
522 m20 = cluster->GetM20();
523 m02 = cluster->GetM02();
529 dEdx = track->GetTPCsignal();
532 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
533 fdEdxBef->Fill(p,dEdx);
534 fTPCnsigma->Fill(p,fTPCnSigma);
536 //Remove electron candidate from the event plane
537 Float_t evPlaneCorrTPC = evPlaneTPC;
538 if(dEdx>70 && dEdx<90){
539 Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
540 Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
541 TVector2 newQVectorfortrack;
542 newQVectorfortrack.Set(qX,qY);
543 evPlaneCorrTPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
546 Bool_t fFlagPhotonicElec = kFALSE;
547 Bool_t fFlagPhotonicElecBCG = kFALSE;
549 //Non-HFE reconstruction
550 fNonHFE->SetPIDresponse(fpidResponse);
551 fNonHFE->FindNonHFE(iTracks,vparticle,fVevent);
552 Int_t *fUlsPartner = fNonHFE->GetPartnersULS(); // Pointer to the ULS partners index
553 Int_t *fLsPartner = fNonHFE->GetPartnersLS(); // Pointer to the LS partners index
555 if (fNonHFE->IsULS()) fFlagPhotonicElec=kTRUE;
556 if (fNonHFE->IsLS()) fFlagPhotonicElecBCG=kTRUE;
558 fNonHFE->SetHistAngleBack(fOpAngleBack);
559 fNonHFE->SetHistAngle(fOpAngle);
560 fNonHFE->SetHistMassBack(fInvMassBack);
561 fNonHFE->SetHistMass(fInvMass);
562 if (fnonHFEalgorithm == "DCA")fNonHFE->SetHistDCABack(fDCABack);
563 if (fnonHFEalgorithm == "DCA")fNonHFE->SetHistDCA(fDCA);
566 Double_t corr[11]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneV0),GetCos2DeltaPhi(phi,evPlaneV0), static_cast<Double_t>(fFlagPhotonicElec),
567 static_cast<Double_t>(fFlagPhotonicElecBCG),m02,m20};
570 Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
571 Int_t whichPart = -99;
572 Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
573 Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
574 Double_t weight = 1.;
575 Double_t Dweight = 1.;
578 Bool_t pi0Decay= kFALSE;
579 Bool_t etaDecay= kFALSE;
582 Double_t phiD = -999.,phie = -999.,phieRec = -999.,ptD = -999.,pte = -999.,pteRec = -999.;
584 if(fIsMC && fMC && stack){
585 Int_t label = track->GetLabel();
587 TParticle *particle = stack->Particle(label);
589 partPDG = particle->GetPdgCode();
590 partPt = particle->Pt();
592 if (TMath::Abs(partPDG)==11) whichPart = 0; //electron
593 if (partPDG==22) whichPart = 3; //gamma
594 if (partPDG==111) whichPart = 2; //pi0
595 if (partPDG==221) whichPart = 1; //eta
597 MChijing = fMC->IsFromBGEvent(label);
600 if(!MChijing) iHijing = 0; // 0 if enhanced sample
602 Int_t idMother = particle->GetFirstMother();
604 TParticle *mother = stack->Particle(idMother);
605 motherPt = mother->Pt();
606 motherPDG = mother->GetPdgCode();
608 if (motherPDG==22) whichFirstMother = 3; //gamma
609 if (motherPDG==111) whichFirstMother = 2; //pi0
610 if (motherPDG==221) whichFirstMother = 1; //eta
612 Int_t idSecondMother = particle->GetSecondMother();
613 if (idSecondMother>0){
614 TParticle *secondMother = stack->Particle(idSecondMother);
615 secondMotherPt = secondMother->Pt();
616 secondMotherPDG = secondMother->GetPdgCode();
618 if (secondMotherPDG==111) whichSecondMother = 2; //pi0
619 if (secondMotherPDG==221) whichSecondMother = 1; //eta
621 Int_t idThirdMother = secondMother->GetFirstMother();
622 if (idThirdMother>0){
623 TParticle *thirdMother = stack->Particle(idThirdMother);
624 thirdMotherPt = thirdMother->Pt();
625 thirdMotherPDG = thirdMother->GetPdgCode();
627 if (thirdMotherPDG==221) whichThirdMother = 1; //eta
632 if (TMath::Abs(partPDG)==11){
635 if (motherPDG==421 || TMath::Abs(motherPDG)==411){ // D
636 Double_t phi_D = -99., Deltaphi_De=-99.;
637 phi_D = mother->Phi();
638 Deltaphi_De = phi_D - phi;
640 fD0_e->Fill(pt,Deltaphi_De);
642 Dweight= GetDweight(0,motherPt,cent);
643 if(iHijing==1) Dweight = 1.0;
644 for(Int_t i=0;i<6;i++){
645 if (motherPt>=ptDmeson[i] && motherPt<ptDmeson[i+1]) fDe[i]->Fill(pt,Dweight);
648 if (motherPDG==421){ // D0
649 Dweight= GetDweight(1,motherPt,cent);
650 if(iHijing==1) Dweight = 1.0;
651 for(Int_t i=0;i<6;i++){
652 if (motherPt>=ptDmeson[i] && motherPt<ptDmeson[i+1]) fD0e[i]->Fill(pt,Dweight);
655 if (motherPDG==411){ // D+
656 Dweight= GetDweight(2,motherPt,cent);
657 if(iHijing==1) Dweight = 1.0;
658 for(Int_t i=0;i<6;i++){
659 if (motherPt>=ptDmeson[i] && motherPt<ptDmeson[i+1]) fDpluse[i]->Fill(pt,Dweight);
662 if (motherPDG==-411){ //D-
663 Dweight= GetDweight(3,motherPt,cent);
664 if(iHijing==1) Dweight = 1.0;
665 for(Int_t i=0;i<6;i++){
666 if (motherPt>=ptDmeson[i] && motherPt<ptDmeson[i+1]) fDminuse[i]->Fill(pt,Dweight);
671 if (motherPDG==111 && secondMotherPDG!=221){ //not eta -> pi0 -> e
672 weight = GetPi0weight(motherPt,cent);
675 if (motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG!=221){ //not eta -> pi0 -> gamma -> e
676 weight = GetPi0weight(secondMotherPt,cent);
681 if (motherPDG==221){ //eta -> e
682 weight = GetEtaweight(motherPt,cent);
685 if (motherPDG==111 && secondMotherPDG==221){ //eta -> pi0 -> e
686 weight = GetEtaweight(secondMotherPt,cent);
689 if (motherPDG==22 && secondMotherPDG==221){ //eta -> gamma -> e
690 weight = GetEtaweight(secondMotherPt,cent);
693 if (motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG==221){ //eta -> pi0 -> gamma -> e
694 weight = GetEtaweight(thirdMotherPt,cent);
699 if (fTPCnSigma>-1 && fTPCnSigma<3 && EovP>1 && EovP<1.3 && (motherPDG==22 || motherPDG==111 || motherPDG==221)){
700 if(iHijing==1) weight = 1.0;
703 fTot_pi0e->Fill(partPt,weight);
704 if(fFlagPhotonicElec) fPhot_pi0e->Fill(partPt,weight);
705 if(fFlagPhotonicElecBCG) fPhotBCG_pi0e->Fill(partPt,weight);
708 fTot_etae->Fill(partPt,weight);
709 if(fFlagPhotonicElec) fPhot_etae->Fill(partPt,weight);
710 if(fFlagPhotonicElecBCG) fPhotBCG_etae->Fill(partPt,weight);
714 Double_t mc[15]={EovP,fTPCnSigma,partPt, static_cast<Double_t>(fFlagPhotonicElec),
715 static_cast<Double_t>(fFlagPhotonicElecBCG), static_cast<Double_t>(whichPart),cent,pt,
716 static_cast<Double_t>(whichFirstMother), static_cast<Double_t>(whichSecondMother),
717 static_cast<Double_t>(whichThirdMother), static_cast<Double_t>(iHijing),motherPt,secondMotherPt,thirdMotherPt};
719 //fMCphotoElecPt->Fill(mc);
720 if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// mother = gamma, pi0, eta
726 if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
730 AliHFEpidObject hfetrack;
731 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
732 hfetrack.SetRecTrack(track);
734 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
736 Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
737 GetCos2DeltaPhi(phi,evPlaneV0C)};
738 fChargPartV2->Fill(corrV2);
740 if(fTPCnSigma >= -0.5){
741 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
742 GetCos2DeltaPhi(phi,evPlaneV0C)};
743 feTPCV2->Fill(correctedV2);
746 if(pidpassed==0) continue;
748 Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),
749 GetCos2DeltaPhi(phi,evPlaneV0C)};
751 feV2->Fill(correctedV2);
752 fTrkEovPAft->Fill(pt,EovP);
753 fdEdxAft->Fill(p,dEdx);
755 if(fFlagPhotonicElec){
756 fphoteV2->Fill(correctedV2);
757 fPhotoElecPt->Fill(pt);
760 if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
763 PostData(1, fOutputList);
765 //_________________________________________
766 void AliAnalysisTaskFlowTPCEMCalEP::UserCreateOutputObjects()
769 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
771 printf("+++++ MC Data available");
773 //--------Initialize PID
774 fPID->SetHasMCData(fIsMC);
776 if(!fPID->GetNumberOfPIDdetectors())
778 fPID->AddDetector("TPC", 0);
779 fPID->AddDetector("EMCAL", 1);
782 fPID->SortDetectors();
783 fPIDqa = new AliHFEpidQAmanager();
784 fPIDqa->Initialize(fPID);
786 //--------Initialize correction Framework and Cuts
787 fCFM = new AliCFManager;
788 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
789 fCFM->SetNStepParticle(kNcutSteps);
790 for(Int_t istep = 0; istep < kNcutSteps; istep++)
791 fCFM->SetParticleCutsList(istep, NULL);
794 AliWarning("Cuts not available. Default cuts will be used");
795 fCuts = new AliHFEcuts;
796 fCuts->CreateStandardCuts();
798 fCuts->Initialize(fCFM);
800 //---------Output Tlist
801 fOutputList = new TList();
802 fOutputList->SetOwner();
803 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
805 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
806 fOutputList->Add(fNoEvents);
808 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
809 fOutputList->Add(fTrkpt);
811 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
812 fOutputList->Add(fTrackPtBefTrkCuts);
814 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
815 fOutputList->Add(fTrackPtAftTrkCuts);
817 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
818 fOutputList->Add(fTPCnsigma);
820 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
821 fOutputList->Add(fTrkEovPBef);
823 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
824 fOutputList->Add(fTrkEovPAft);
826 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
827 fOutputList->Add(fdEdxBef);
829 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
830 fOutputList->Add(fdEdxAft);
832 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
833 fOutputList->Add(fPhotoElecPt);
835 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
836 fOutputList->Add(fSemiInclElecPt);
838 fCent = new TH1F("fCent","Centrality",100,0,100) ;
839 fOutputList->Add(fCent);
841 fevPlaneV0A = new TH1F("fevPlaneV0A","V0A EP",100,0,TMath::Pi());
842 fOutputList->Add(fevPlaneV0A);
844 fevPlaneV0C = new TH1F("fevPlaneV0C","V0C EP",100,0,TMath::Pi());
845 fOutputList->Add(fevPlaneV0C);
847 fevPlaneV0 = new TH1F("fevPlaneV0","V0 EP",100,0,TMath::Pi());
848 fOutputList->Add(fevPlaneV0);
850 fevPlaneTPC = new TH1F("fevPlaneTPC","TPC EP",100,0,TMath::Pi());
851 fOutputList->Add(fevPlaneTPC);
853 fTPCsubEPres = new TH1F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1);
854 fOutputList->Add(fTPCsubEPres);
856 Int_t binsv1[7]={100,100,100,100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, V0-TPCpos, V0-TPCneg, TPCpos-TPCneg, cent
857 Double_t xminv1[7]={-1,-1,-1,-1,-1,-1,0};
858 Double_t xmaxv1[7]={1,1,1,1,1,1,90};
859 fEPres = new THnSparseD ("fEPres","EP resolution",7,binsv1,xminv1,xmaxv1);
860 fOutputList->Add(fEPres);
862 //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneV0),GetCos2DeltaPhi(phi,evPlaneV0), fFlagPhotonicElec, fFlagPhotonicElecBCG,m02,m20
863 Int_t binsv2[11]={100,200,90,100,100,100,100,3,3,100,100};
864 Double_t xminv2[11]={0,-5,0,0,0,0,-1,-1,-1,0,0};
865 Double_t xmaxv2[11]={2*TMath::Pi(),5,90,50,3,TMath::Pi(),1,2,2,1,1};
866 fCorr = new THnSparseD ("fCorr","Correlations",11,binsv2,xminv2,xmaxv2);
867 fOutputList->Add(fCorr);
869 Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
870 Double_t xminv3[5]={0,0,-1,-1,-1};
871 Double_t xmaxv3[5]={90,50,1,1,1};
872 feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
873 fOutputList->Add(feV2);
875 Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
876 Double_t xminv4[5]={0,0,-1,-1,-1};
877 Double_t xmaxv4[5]={90,50,1,1,1};
878 fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
879 fOutputList->Add(fphoteV2);
881 Int_t binsv5[7]={100,90,100,100,100,100,100}; // phi, cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
882 Double_t xminv5[7]={0,0,0,0,-1,-1,-1};
883 Double_t xmaxv5[7]={2*TMath::Pi(),90,50,3,1,1,1};
884 fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",7,binsv5,xminv5,xmaxv5);
885 fOutputList->Add(fChargPartV2);
887 Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
888 Double_t xminv6[5]={0,0,-1,-1,-1};
889 Double_t xmaxv6[5]={90,50,1,1,1};
890 feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
891 fOutputList->Add(feTPCV2);
893 //EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,firstMother,secondMother,thirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt
894 Int_t binsv7[15]={100,100,100,3,3,5,90,100,5,5,5,3,100,100,100};
895 Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
896 Double_t xmaxv7[15]={3,3.5,50,2,2,4,90,50,4,4,4,2,50,50,50};
897 fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
898 fOutputList->Add(fMCphotoElecPt);
900 fGammaWeight = new TH2F("fGammaWeight", "Gamma weight",100,0,50,3,-1,2);
901 fOutputList->Add(fGammaWeight);
903 fPi0Weight = new TH2F("fPi0Weight", "Pi0 weight",100,0,50,3,-1,2);
904 fOutputList->Add(fPi0Weight);
906 fEtaWeight = new TH2F("fEtaWeight", "Eta weight",100,0,50,3,-1,2);
907 fOutputList->Add(fEtaWeight);
909 fD0Weight = new TH2F("fD0Weight", "D0 weight",100,0,50,3,-1,2);
910 fOutputList->Add(fD0Weight);
912 fDplusWeight = new TH2F("fDplusWeight", "D+ weight",100,0,50,3,-1,2);
913 fOutputList->Add(fDplusWeight);
915 fDminusWeight = new TH2F("fDminusWeight", "D- weight",100,0,50,3,-1,2);
916 fOutputList->Add(fDminusWeight);
918 fD0_e = new TH2F("fD0_e", "D0 vs e",100,0,50,200,-6.3,6.3);
919 fOutputList->Add(fD0_e);
921 for(Int_t k = 0; k < 6; k++) {
923 TString De_name = Form("fDe%d",k);
924 TString D0e_name = Form("fD0e%d",k);
925 TString Dpluse_name = Form("fDpluse%d",k);
926 TString Dminuse_name = Form("fDminuse%d",k);
928 fDe[k] = new TH1F((const char*)De_name,"",100,0,50);
929 fD0e[k] = new TH1F((const char*)D0e_name,"",100,0,50);
930 fDpluse[k] = new TH1F((const char*)Dpluse_name,"",100,0,50);
931 fDminuse[k] = new TH1F((const char*)Dminuse_name,"",100,0,50);
933 fOutputList->Add(fDe[k]);
934 fOutputList->Add(fD0e[k]);
935 fOutputList->Add(fDpluse[k]);
936 fOutputList->Add(fDminuse[k]);
941 double bin_v2[9] = {2,2.5,3,4,6,8,10,13,18};
943 fTot_pi0e = new TH1F("fTot_pi0e","fTot_pi0e",nbin_v2,bin_v2);
944 fOutputList->Add(fTot_pi0e);
946 fPhot_pi0e = new TH1F("fPhot_pi0e","fPhot_pi0e",nbin_v2,bin_v2);
947 fOutputList->Add(fPhot_pi0e);
949 fPhotBCG_pi0e = new TH1F("fPhotBCG_pi0e","fPhotBCG_pi0e",nbin_v2,bin_v2);
950 fOutputList->Add(fPhotBCG_pi0e);
952 fTot_etae = new TH1F("fTot_etae","fTot_etae",nbin_v2,bin_v2);
953 fOutputList->Add(fTot_etae);
955 fPhot_etae = new TH1F("fPhot_etae","fPhot_etae",nbin_v2,bin_v2);
956 fOutputList->Add(fPhot_etae);
958 fPhotBCG_etae = new TH1F("fPhotBCG_etae","fPhotBCG_etae",nbin_v2,bin_v2);
959 fOutputList->Add(fPhotBCG_etae);
961 fInvMass = new TH1F("fInvMass","",200,0,0.3);
962 fOutputList->Add(fInvMass);
964 fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
965 fOutputList->Add(fInvMassBack);
967 fDCA = new TH1F("fDCA","",200,0,1);
968 fOutputList->Add(fDCA);
970 fDCABack = new TH1F("fDCABack","",200,0,1);
971 fOutputList->Add(fDCABack);
973 fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
974 fOutputList->Add(fOpAngle);
976 fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
977 fOutputList->Add(fOpAngleBack);
979 PostData(1,fOutputList);
982 //________________________________________________________________________
983 void AliAnalysisTaskFlowTPCEMCalEP::Terminate(Option_t *)
985 // Info("Terminate");
986 AliAnalysisTaskSE::Terminate();
989 //________________________________________________________________________
990 Bool_t AliAnalysisTaskFlowTPCEMCalEP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
992 // Check single track cuts for a given cut step
993 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
994 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
997 //_________________________________________
998 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
1000 //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
1001 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
1002 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
1003 Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
1005 return cos2DeltaPhi;
1007 //_________________________________________
1008 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetDeltaPhi(Double_t phiA,Double_t phiB) const
1011 Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
1012 if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
1016 //_________________________________________
1017 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetPi0weight(Double_t mcPi0pT, Float_t cent) const
1020 double weight = 1.0;
1022 if(mcPi0pT>0.0 && mcPi0pT<5.0){
1023 if (cent>20.0 && cent<40.0) weight = (2.877091*mcPi0pT)/(TMath::Power(0.706963+mcPi0pT/3.179309,17.336628)*exp(-mcPi0pT));
1024 if (cent>30.0 && cent<50.0) weight = (2.392024*mcPi0pT)/(TMath::Power(0.688810+mcPi0pT/3.005145,16.811845)*exp(-mcPi0pT));
1027 if (cent>20.0 && cent<40.0) weight = (0.0004*mcPi0pT)/TMath::Power(-0.176181+mcPi0pT/3.989747,5.629235);
1028 if (cent>30.0 && cent<50.0) weight = (0.000186*mcPi0pT)/TMath::Power(-0.606279+mcPi0pT/3.158642,4.365540);
1032 //_________________________________________
1033 Double_t AliAnalysisTaskFlowTPCEMCalEP::GetEtaweight(Double_t mcEtapT, Float_t cent) const
1036 double weight = 1.0;
1038 if (cent>20.0 && cent<40.0) weight = (0.818052*mcEtapT)/(TMath::Power(0.358651+mcEtapT/2.878631,9.494043));
1039 if (cent>30.0 && cent<50.0) 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, Float_t cent) const
1046 double weight = 1.0;
1048 if (cent>30.0 && cent<50.0){
1049 if (whichD == 0) weight = 0.271583*TMath::Landau(mcDpT,3.807103,1.536753,0); // D
1050 if (whichD == 1) weight = 0.300771*TMath::Landau(mcDpT,3.725771,1.496980,0); // D0
1051 if (whichD == 2) weight = 0.247280*TMath::Landau(mcDpT,3.746811,1.607551,0); // D+
1052 if (whichD == 3) weight = 0.249410*TMath::Landau(mcDpT,3.611508,1.632196,0); //D-
1056 //_________________________________________
1057 void AliAnalysisTaskFlowTPCEMCalEP::InitParameters()
1061 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1062 fTrackCuts->SetRequireTPCRefit(kTRUE);
1063 fTrackCuts->SetRequireITSRefit(kTRUE);
1064 fTrackCuts->SetEtaRange(-0.7,0.7);
1065 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1066 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1067 fTrackCuts->SetMinNClustersTPC(100);
1068 fTrackCuts->SetPtRange(0.5,100);
1070 fNonHFE->SetAODanalysis(kFALSE);
1071 fNonHFE->SetInvariantMassCut(fInvmassCut);
1072 fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
1073 fNonHFE->SetChi2OverNDFCut(fChi2Cut);
1074 fNonHFE->SetAlgorithm(fnonHFEalgorithm); //KF or DCA
1075 if (fnonHFEalgorithm == "DCA") fNonHFE->SetDCACut(fDCAcut);
1076 fNonHFE->SetTrackCuts(-3.5,3.5,fTrackCuts);