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 **************************************************************************/
17 ////////////////////////////////////////////////////////////////////////
19 // Task for Heavy Flavour Electron Flow with TPC plus EMCal //
20 // Non-Photonic Electron identified with Invariant mass //
21 // analysis methos in function SelectPhotonicElectron //
24 // Author: Andrea Dubla (Utrecht University) //
27 ////////////////////////////////////////////////////////////////////////
34 #include "THnSparse.h"
35 #include "TLorentzVector.h"
38 #include "AliAnalysisTask.h"
39 #include "AliAnalysisManager.h"
40 #include "AliESDEvent.h"
41 #include "AliESDHandler.h"
42 #include "AliAODEvent.h"
43 #include "AliAODHandler.h"
44 #include "AliAnalysisTaskFlowTPCEMCalQCSP.h"
45 #include "TGeoGlobalMagField.h"
47 #include "AliAnalysisTaskSE.h"
48 #include "TRefArray.h"
50 #include "AliESDInputHandler.h"
51 #include "AliESDpid.h"
52 #include "AliAODInputHandler.h"
53 #include "AliAODPid.h"
54 #include "AliESDtrackCuts.h"
55 #include "AliPhysicsSelection.h"
56 #include "AliCentralitySelectionTask.h"
57 #include "AliESDCaloCluster.h"
58 #include "AliAODCaloCluster.h"
59 #include "AliESDCaloTrigger.h"
60 #include "AliEMCALRecoUtils.h"
61 #include "AliEMCALGeometry.h"
62 #include "AliGeomManager.h"
64 #include "TGeoManager.h"
67 #include "AliEMCALTrack.h"
68 //#include "AliEMCALTracker.h"
70 #include "AliKFParticle.h"
71 #include "AliKFVertex.h"
73 #include "AliPIDResponse.h"
74 #include "AliHFEcontainer.h"
75 #include "AliHFEcuts.h"
76 #include "AliHFEpid.h"
77 #include "AliHFEpidBase.h"
78 #include "AliHFEpidQAmanager.h"
79 #include "AliHFEtools.h"
80 #include "AliCFContainer.h"
81 #include "AliCFManager.h"
82 #include "AliKFParticle.h"
83 #include "AliKFVertex.h"
84 #include "AliCentrality.h"
85 #include "AliVEvent.h"
87 #include "AliMCEvent.h"
89 #include "AliFlowCandidateTrack.h"
90 #include "AliFlowTrackCuts.h"
91 #include "AliFlowEventSimple.h"
92 #include "AliFlowCommonConstants.h"
93 #include "AliFlowEvent.h"
96 #include "AliESDVZERO.h"
97 #include "AliAODVZERO.h"
99 #include "AliPIDResponse.h"
100 #include "AliFlowTrack.h"
101 #include "AliAnalysisTaskVnV0.h"
102 #include "AliSelectNonHFE.h"
105 class AliFlowTrackCuts;
109 ClassImp(AliAnalysisTaskFlowTPCEMCalQCSP)
110 //________________________________________________________________________
111 AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP(const char *name)
112 : AliAnalysisTaskSE(name)
118 ,fIdentifiedAsOutInz(kFALSE)
119 ,fPassTheEventCut(kFALSE)
123 ,fCutsRP(0) // track cuts for reference particles
124 ,fNullCuts(0) // dummy cuts for flow event tracks
125 ,fFlowEvent(0) //! flow events (one for each inv mass band)
126 ,fkCentralityMethod(0)
145 ,fCentralityNoPass(0)
161 ,fMultCorAfterCuts(0)
168 ,fSparseElectronHadron(0)
170 ,fMultCorBeforeCuts(0)
171 ,fSideBandsFlow(kFALSE)
172 ,fPhiminusPsi(kFALSE)
173 ,fFlowEventCont(0) //! flow events (one for each inv mass band)
175 ,fSparseElectronpurity(0)
178 ,fNonHFE(new AliSelectNonHFE)
185 ,fMultCorAfterCentrBeforeCuts(0)
186 ,fMultCorAfterVZTRKComp(0)
187 ,fCentralityBeforePileup(0)
188 ,fCentralityAfterVZTRK(0)
189 ,fCentralityAfterCorrCut(0)
190 ,fMultCorAfterCorrCut(0)
194 ,fSubEventDPhiv2new(0)
196 ,fHistCentrDistr(0x0)
197 ,fCentralityNoPassForFlattening(0)
198 ,fInvmassLS1highpt(0)
199 ,fInvmassULS1highpt(0)
207 fPID = new AliHFEpid("hfePid");
208 // Define input and output slots here
209 // Input slot #0 works with a TChain
210 DefineInput(0, TChain::Class());
211 // Output slot #0 id reserved by the base class for AOD
212 // Output slot #1 writes into a TH1 container
213 // DefineOutput(1, TH1I::Class());
214 DefineOutput(1, TList::Class());
215 DefineOutput(2, AliFlowEventSimple::Class());
217 DefineOutput(3, AliFlowEventSimple::Class());
219 // DefineOutput(3, TTree::Class());
222 //________________________________________________________________________
223 AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP()
224 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskFlowTPCEMCalQCSP")
230 ,fIdentifiedAsOutInz(kFALSE)
231 ,fPassTheEventCut(kFALSE)
235 ,fCutsRP(0) // track cuts for reference particles
236 ,fNullCuts(0) // dummy cuts for flow event tracks
237 ,fFlowEvent(0) //! flow events (one for each inv mass band)
238 ,fkCentralityMethod(0)
257 ,fCentralityNoPass(0)
273 ,fMultCorAfterCuts(0)
280 ,fSparseElectronHadron(0)
282 ,fMultCorBeforeCuts(0)
283 ,fSideBandsFlow(kFALSE)
284 ,fPhiminusPsi(kFALSE)
285 ,fFlowEventCont(0) //! flow events (one for each inv mass band)
287 ,fSparseElectronpurity(0)
290 ,fNonHFE(new AliSelectNonHFE)
297 ,fMultCorAfterCentrBeforeCuts(0)
298 ,fMultCorAfterVZTRKComp(0)
299 ,fCentralityBeforePileup(0)
300 ,fCentralityAfterVZTRK(0)
301 ,fCentralityAfterCorrCut(0)
302 ,fMultCorAfterCorrCut(0)
306 ,fSubEventDPhiv2new(0)
308 ,fHistCentrDistr(0x0)
309 ,fCentralityNoPassForFlattening(0)
310 ,fInvmassLS1highpt(0)
311 ,fInvmassULS1highpt(0)
317 //Default constructor
318 fPID = new AliHFEpid("hfePid");
320 // Define input and output slots here
321 // Input slot #0 works with a TChain
322 DefineInput(0, TChain::Class());
323 // Output slot #0 id reserved by the base class for AOD
324 // Output slot #1 writes into a TH1 container
325 // DefineOutput(1, TH1I::Class());
326 DefineOutput(1, TList::Class());
327 DefineOutput(2, AliFlowEventSimple::Class());
328 // DefineOutput(3, TTree::Class());
330 DefineOutput(3, AliFlowEventSimple::Class());
332 //DefineOutput(3, TTree::Class());
334 //_________________________________________
336 AliAnalysisTaskFlowTPCEMCalQCSP::~AliAnalysisTaskFlowTPCEMCalQCSP()
345 if (fDCA) delete fNonHFE;
346 if (fOutputList) delete fOutputList;
347 if (fFlowEvent) delete fFlowEvent;
348 if (fFlowEventCont) delete fFlowEventCont;
351 //_________________________________________
353 void AliAnalysisTaskFlowTPCEMCalQCSP::UserExec(Option_t*)
356 //Called for each event
358 // create pointer to event
360 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
363 printf("ERROR: fAOD not available\n");
369 AliError("HFE cuts not available");
373 if(!fPID->IsInitialized())
375 // Initialize PID with the given run number
376 AliWarning("PID not initialised, get from Run no");
377 fPID->InitializePID(fAOD->GetRunNumber());
380 // cout << "kTrigger == " << fTrigger <<endl;
382 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral)) return;
386 if ( !(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny) ) return;
388 TString firedTriggerClasses = static_cast<const AliAODEvent*>(InputEvent())->GetFiredTriggerClasses();
390 if ( ! ( firedTriggerClasses.Contains("CVLN_B2-B-NOPF-ALLNOTRD") || firedTriggerClasses.Contains("CVLN_R1-B-NOPF-ALLNOTRD") || firedTriggerClasses.Contains("CSEMI_R1-B-NOPF-ALLNOTRD") ) ) return;
394 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA)) return;
397 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return;
400 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kCentral | AliVEvent::kSemiCentral))) return;
403 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
407 //---------------CENTRALITY AND EVENT SELECTION-----------------------
408 Int_t fNOtrks = fAOD->GetNumberOfTracks();
410 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
411 if (!trkVtx || trkVtx->GetNContributors()<=0)return;
412 TString vtxTtl = trkVtx->GetTitle();
413 if (!vtxTtl.Contains("VertexerTracks"))return;
414 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
415 if (!spdVtx || spdVtx->GetNContributors()<=0)return;
416 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
417 vtxz = trkVtx->GetZ();
418 if(TMath::Abs(vtxz)>10)return;
421 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
422 if(fNOtrks<2) return;
424 Bool_t pass = kFALSE; //to select centrality and pile up protection
425 CheckCentrality(fAOD,pass);
430 PlotVZeroMultiplcities(fAOD);
433 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent); //Calculate event plane Qvector and EP resolution for inclusive
436 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEventCont); //Calculate event plane Qvector and EP resolution for inclusive
439 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
442 AliDebug(1, "Using default PID Response");
443 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
446 fPID->SetPIDResponse(pidResponse);
447 fCFM->SetRecEventInfo(fAOD);
449 // Look for kink mother
450 Int_t numberofvertices = fAOD->GetNumberOfVertices();
451 Double_t listofmotherkink[numberofvertices];
452 Int_t numberofmotherkink = 0;
453 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
454 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
455 if(!aodvertex) continue;
456 if(aodvertex->GetType()==AliAODVertex::kKink) {
457 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
458 if(!mother) continue;
459 Int_t idmother = mother->GetID();
460 listofmotherkink[numberofmotherkink] = idmother;
461 //printf("ID %d\n",idmother);
462 numberofmotherkink++;
467 //=============================================V0EP from Alex======================================================================
468 Double_t qxEPa = 0, qyEPa = 0;
469 Double_t qxEPc = 0, qyEPc = 0;
470 Double_t qxEP = 0, qyEP = 0;
472 Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
473 Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
474 Double_t evPlAngV0 = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 10, 2, qxEP, qyEP);
477 Double_t Qx2 = 0, Qy2 = 0;
478 Double_t Qx2p = 0, Qy2p = 0;
479 Double_t Qx2n = 0, Qy2n = 0;
481 for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
483 AliAODTrack* aodTrack = fAOD->GetTrack(iT);
488 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
491 if (!aodTrack->TestFilterBit(128))
495 if(aodTrack->Eta()>0 && aodTrack->Eta()<0.8){
497 Qx2p += TMath::Cos(2*aodTrack->Phi());
498 Qy2p += TMath::Sin(2*aodTrack->Phi());
500 if(aodTrack->Eta()<0 && aodTrack->Eta()> -0.8){
502 Qx2n += TMath::Cos(2*aodTrack->Phi());
503 Qy2n += TMath::Sin(2*aodTrack->Phi());
507 Qx2 += TMath::Cos(2*aodTrack->Phi());
508 Qy2 += TMath::Sin(2*aodTrack->Phi());
515 Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
516 Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
517 Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
519 EPVzA->Fill(evPlAngV0A);
520 EPVzC->Fill(evPlAngV0C);
521 EPTPC->Fill(evPlAngTPC);
523 EPTPCn->Fill(evPlAngTPCn);
524 EPTPCp->Fill(evPlAngTPCp);
525 EPVz->Fill(evPlAngV0);
530 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
531 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
532 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
535 fSubEventDPhiv2new->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))); // vzero - tpcp
536 fSubEventDPhiv2new->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))); // vzero - tpcn
537 fSubEventDPhiv2new->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
542 //====================================================================================================================
545 AliAODTrack *track = NULL;
548 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
550 track = fAOD->GetTrack(iTracks);
553 printf("ERROR: Could not receive track %d\n", iTracks);
557 if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
558 //----------hfe begin---------
559 if(track->Eta()<-0.7 || track->Eta()>0.7) continue; //eta cuts on candidates
561 // RecKine: ITSTPC cuts
562 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
564 // Reject kink mother
565 Bool_t kinkmotherpass = kTRUE;
566 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
567 if(track->GetID() == listofmotherkink[kinkmother]) {
568 kinkmotherpass = kFALSE;
572 if(!kinkmotherpass) continue;
575 // if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; //deleted for DCA absence
576 // HFEcuts: ITS layers cuts
577 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
578 // HFE cuts: TPC PID cleanup
579 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
581 Double_t fClsE = -999, p = -999, fEovP=-999, pt = -999, fTPCnSigma=0;
582 // Track extrapolation
583 Int_t fClsId = track->GetEMCALcluster();
584 if(fClsId < 0) continue;
585 AliAODCaloCluster *cluster = fAOD->GetCaloCluster(fClsId);
586 if(TMath::Abs(cluster->GetTrackDx()) > 0.05 || TMath::Abs(cluster->GetTrackDz()) > 0.05) continue;
588 pt = track->Pt(); //pt track after cuts
589 if(pt<fpTCut) continue;
590 fClsE = cluster->E();
592 // dEdx = track->GetTPCsignal();
594 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
595 Double_t m20 =cluster->GetM20();
596 Double_t m02 =cluster->GetM02();
597 Double_t disp=cluster->GetDispersion();
598 if(fTPCnSigma >= -1 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,fEovP);
599 fTPCnsigma->Fill(p,fTPCnSigma);
600 // fdEdxBef->Fill(p,dEdx);
601 Double_t eta = track->Eta();
602 Double_t phi = track->Phi();
603 //-----------------------Phiminupsi method to remove the contamination-----------------------------------------------
604 //-----------------------fTPCnSigma < -3.5 hadrons will be selected from this region--------------------------
605 Float_t dPhi_aeh = TVector2::Phi_0_2pi(phi - evPlAngV0A);
606 if(dPhi_aeh > TMath::Pi()) dPhi_aeh = dPhi_aeh - TMath::Pi();
607 Float_t dPhi_ceh = TVector2::Phi_0_2pi(phi - evPlAngV0C);
608 if(dPhi_ceh > TMath::Pi()) dPhi_ceh = dPhi_ceh - TMath::Pi();
611 Double_t valueElh[8] = {
620 fSparseElectronHadron->Fill(valueElh);
622 //----------------------------------------------------------------------------------------------------------
623 //---------------------------From here usual electron selection---------------------------------------------
624 //----------------------------------------------------------------------------------------------------------
625 if(m20 < fminM20 || m20 > fmaxM20) continue;
626 if(m02 < fminM02 || m02 > fmaxM02) continue;
627 if(disp > fDispersion ) continue;
628 //---------------------------------for purity---------------------------------------------------------------
630 Double_t valuepurity[3] = {
634 fSparseElectronpurity->Fill(valuepurity);
636 //----------------------------------------------------------------------------------------------------------
637 //----------------------------------------------------------------------------------------------------------
638 if(fTPCnSigma < fminTPC || fTPCnSigma > fmaxTPC) continue; //cuts on nsigma tpc and EoP
639 //===============================Flow Event for Contamination=============================================
641 if(fEovP>0 && fEovP<0.6){
642 AliFlowTrack *sTrackCont = new AliFlowTrack();
643 sTrackCont->Set(track);
644 sTrackCont->SetID(track->GetID());
645 sTrackCont->SetForRPSelection(kTRUE);
646 sTrackCont->SetForPOISelection(kTRUE);
647 sTrackCont->SetMass(2637);
648 for(int iRPs=0; iRPs!=fFlowEventCont->NumberOfTracks(); ++iRPs)
650 // cout << " no of rps " << iRPs << endl;
651 AliFlowTrack *iRPCont = dynamic_cast<AliFlowTrack*>(fFlowEventCont->GetTrack( iRPs ));
652 if (!iRPCont) continue;
653 if (!iRPCont->InRPSelection()) continue;
654 if( sTrackCont->GetID() == iRPCont->GetID())
656 if(fDebug) printf(" was in RP set");
657 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set" <<endl;
658 iRPCont->SetForRPSelection(kFALSE);
659 // fFlowEventCont->SetNumberOfRPs(fFlowEventCont->GetNumberOfRPs() - 1);
661 } //end of for loop on RPs
662 fFlowEventCont->InsertTrack(((AliFlowTrack*) sTrackCont));
663 fFlowEventCont->SetNumberOfPOIs(fFlowEventCont->GetNumberOfPOIs()+1);
667 //==========================================================================================================
668 //===============================From here eovP cut is used fro QC, SP and EPV0=============================
669 if(fEovP < fminEovP || fEovP >fmaxEovP) continue;
670 //==========================================================================================================
671 //============================Event Plane Method with V0====================================================
672 Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
673 Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
674 Double_t v2Phi[3] = {
680 Double_t v2PhiVz = TMath::Cos(2*(phi - evPlAngV0));
681 Double_t v2PhiV0tot[2] = {
684 fV2Phivzerotot->Fill(v2PhiV0tot);
685 //=========================================================================================================
686 fTPCnsigmaAft->Fill(p,fTPCnSigma);
687 fInclusiveElecPt->Fill(pt);
690 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
691 AliFlowTrack *sTrack = new AliFlowTrack();
693 sTrack->SetID(track->GetID());
694 sTrack->SetForRPSelection(kTRUE);
695 sTrack->SetForPOISelection(kTRUE);
696 sTrack->SetMass(263732);
697 for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
699 // cout << " no of rps " << iRPs << endl;
700 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
702 if (!iRP->InRPSelection()) continue;
703 if( sTrack->GetID() == iRP->GetID())
705 if(fDebug) printf(" was in RP set");
706 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set" <<endl;
707 iRP->SetForRPSelection(kFALSE);
708 // fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
710 } //end of for loop on RPs
711 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
712 fFlowEvent->SetNumberOfPOIs(fFlowEvent->GetNumberOfPOIs()+1);
716 //----------------------Selection of Photonic Electrons DCA-----------------------------
717 fNonHFE = new AliSelectNonHFE();
718 fNonHFE->SetAODanalysis(kTRUE);
719 fNonHFE->SetInvariantMassCut(fInvmassCut);
720 if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
721 //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
722 //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
723 fNonHFE->SetAlgorithm("DCA"); //KF
724 fNonHFE->SetPIDresponse(pidResponse);
725 fNonHFE->SetTrackCuts(-3,3);
727 fNonHFE->SetHistAngleBack(fOpeningAngleLS);
728 fNonHFE->SetHistAngle(fOpeningAngleULS);
729 //fNonHFE->SetHistDCABack(fDCABack);
730 //fNonHFE->SetHistDCA(fDCA);
731 fNonHFE->SetHistMassBack(fInvmassLS1);
732 fNonHFE->SetHistMass(fInvmassULS1);
734 fNonHFE->FindNonHFE(iTracks,track,fAOD);
736 // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
737 // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
738 // Bool_t fUlsIsPartner = kFALSE;
739 // Bool_t fLsIsPartner = kFALSE;
740 if(fNonHFE->IsULS()){
741 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
742 fULSElecPt->Fill(track->Pt());
747 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
748 fLSElecPt->Fill(track->Pt());
754 //----------------------Selection of Photonic Electrons KFParticle-----------------------------
755 Bool_t fFlagPhotonicElec = kFALSE;
756 SelectPhotonicElectron(iTracks,track,fEovP, evPlAngV0, fFlagPhotonicElec);
757 if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
758 // Semi inclusive electron
759 if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
766 PostData(1, fOutputList);
767 PostData(2, fFlowEvent);
769 PostData(3, fFlowEventCont);
772 //----------hfe end---------
774 //_________________________________________
775 void AliAnalysisTaskFlowTPCEMCalQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track,Double_t fEovP,Double_t evPlAngV0, Bool_t &fFlagPhotonicElec)
777 //Identify non-heavy flavour electrons using Invariant mass method KF
779 Bool_t flagPhotonicElec = kFALSE;
781 for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
782 AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
784 printf("ERROR: Could not receive track %d\n", jTracks);
787 // if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
788 if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
789 // if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit) || (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
792 if(!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
795 if(!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
797 if(jTracks == itrack) continue;
798 Double_t ptAsso=-999., nsigma=-999.0;
799 Double_t mass=-999., width = -999;
800 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
801 Double_t openingAngle = -999.;
802 Double_t ptcutonmasshighpt = track->Pt();
804 ptAsso = trackAsso->Pt();
805 Short_t chargeAsso = trackAsso->Charge();
806 Short_t charge = track->Charge();
808 nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
811 if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
812 if(nsigma < -3 || nsigma > 3) continue;
813 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
814 if(ptAsso <0.3) continue;
816 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
817 if(charge>0) fPDGe1 = -11;
818 if(chargeAsso>0) fPDGe2 = -11;
820 if(charge == chargeAsso) fFlagLS = kTRUE;
821 if(charge != chargeAsso) fFlagULS = kTRUE;
823 AliKFParticle::SetField(fAOD->GetMagneticField());
824 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
825 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
826 AliKFParticle recg(ge1, ge2);
828 if(recg.GetNDF()<1) continue;
829 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
830 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
831 recg.GetMass(mass,width);
833 openingAngle = ge1.GetAngle(ge2);
834 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
835 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
836 if(fOP_angle)if(openingAngle > fOpeningAngleCut) continue;
839 if(fFlagLS) fInvmassLS1->Fill(mass);
840 if(fFlagULS) fInvmassULS1->Fill(mass);
843 Double_t MassSparseULS[3] = {
847 fSparseMassULS->Fill(MassSparseULS);
850 Double_t MassSparseLS[3] = {
854 fSparseMassLS->Fill(MassSparseLS);
858 if(ptcutonmasshighpt >= 8.){
859 if(fFlagLS) fInvmassLS1highpt->Fill(mass);
860 if(fFlagULS) fInvmassULS1highpt->Fill(mass);
864 if(mass<fInvmassCut){
865 if(fFlagULS){fULSElecPt->Fill(track->Pt());}
866 if(fFlagLS){fLSElecPt->Fill(track->Pt());}
871 Double_t phi = track->Phi();
872 Float_t DeltaPhi_eEP = TVector2::Phi_0_2pi(phi - evPlAngV0);
873 if(DeltaPhi_eEP > TMath::Pi()) {DeltaPhi_eEP = DeltaPhi_eEP - TMath::Pi();}
876 if(mass<fInvmassCut){
878 Double_t ulsSparse[3] = {
883 fSparsephipsiULS->Fill(ulsSparse);
886 Double_t lsSparse[3] = {
891 fSparsephipsiLS->Fill(lsSparse);
897 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
898 flagPhotonicElec = kTRUE;
901 fFlagPhotonicElec = flagPhotonicElec;
903 //__________________________________________________________________________________
904 void AliAnalysisTaskFlowTPCEMCalQCSP::UserCreateOutputObjects()
908 //----------hfe initialising begin---------
909 fNullCuts = new AliFlowTrackCuts("null_cuts");
911 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
912 cc->SetNbinsMult(10000);
914 cc->SetMultMax(10000);
920 cc->SetNbinsPhi(180);
922 cc->SetPhiMax(TMath::TwoPi());
932 //--------Initialize PID
933 fPID->SetHasMCData(kFALSE);
934 if(!fPID->GetNumberOfPIDdetectors())
936 fPID->AddDetector("TPC", 0);
937 fPID->AddDetector("EMCAL", 1);
940 fPID->SortDetectors();
941 fPIDqa = new AliHFEpidQAmanager();
942 fPIDqa->Initialize(fPID);
944 //--------Initialize correction Framework and Cuts
945 fCFM = new AliCFManager;
946 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
947 fCFM->SetNStepParticle(kNcutSteps);
948 for(Int_t istep = 0; istep < kNcutSteps; istep++)
949 fCFM->SetParticleCutsList(istep, NULL);
952 AliWarning("Cuts not available. Default cuts will be used");
953 fCuts = new AliHFEcuts;
954 fCuts->CreateStandardCuts();
958 fCuts->Initialize(fCFM);
959 //----------hfe initialising end--------
960 //---------Output Tlist
961 fOutputList = new TList();
962 fOutputList->SetOwner();
963 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
965 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
966 fOutputList->Add(fNoEvents);
968 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
969 fOutputList->Add(fTPCnsigma);
971 fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
972 fOutputList->Add(fTPCnsigmaAft);
974 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",1000,0,50,100,0,2);
975 fOutputList->Add(fTrkEovPBef);
977 // fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",1000,0,50,150,0,150);
978 // fOutputList->Add(fdEdxBef);
980 fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
981 fOutputList->Add(fInclusiveElecPt);
983 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
984 fOutputList->Add(fPhotoElecPt);
986 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
987 fOutputList->Add(fSemiInclElecPt);
989 fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
990 fOutputList->Add(fULSElecPt);
992 fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
993 fOutputList->Add(fLSElecPt);
995 fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
996 fOutputList->Add(fInvmassLS1);
998 fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
999 fOutputList->Add(fInvmassULS1);
1001 fInvmassLS1highpt = new TH1F("fInvmassLS1highpt", "Inv mass of LS (e,e); mass(GeV/c^2) highpt; counts;", 1000,0,1.0);
1002 fOutputList->Add(fInvmassLS1highpt);
1004 fInvmassULS1highpt = new TH1F("fInvmassULS1highpt", "Inv mass of ULS (e,e); mass(GeV/c^2) highpt; counts;", 1000,0,1.0);
1005 fOutputList->Add(fInvmassULS1highpt);
1007 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
1008 fOutputList->Add(fCentralityPass);
1010 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
1011 fOutputList->Add(fCentralityNoPass);
1013 fCentralityNoPassForFlattening = new TH1F("fCentralityNoPassForFlattening", "Centrality No Pass for flattening", 101, -1, 100);
1014 fOutputList->Add(fCentralityNoPassForFlattening);
1016 fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
1017 fOutputList->Add(fCentralityBeforePileup);
1019 fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
1020 fOutputList->Add(fCentralityAfterVZTRK);
1022 fCentralityAfterCorrCut = new TH1F("fCentralityAfterCorrCut", "fCentralityAfterCorrCut Pass", 101, -1, 100);
1023 fOutputList->Add(fCentralityAfterCorrCut);
1025 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
1026 fOutputList->Add(fPhi);
1028 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
1029 fOutputList->Add(fEta);
1031 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
1032 fOutputList->Add(fVZEROA);
1034 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
1035 fOutputList->Add(fVZEROC);
1037 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
1038 fOutputList->Add(fTPCM);
1040 fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
1041 fOutputList->Add(fvertex);
1043 fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1044 fOutputList->Add(fMultCorBeforeCuts);
1046 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1047 fOutputList->Add(fMultCorAfterCuts);
1049 fMultCorAfterCentrBeforeCuts = new TH2F("fMultCorAfterCentrBeforeCuts", "TPC vs Global multiplicity (After CC before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1050 fOutputList->Add(fMultCorAfterCentrBeforeCuts);
1052 fMultCorAfterVZTRKComp = new TH2F("fMultCorAfterVZTRKComp", "TPC vs Global multiplicity (After V0-TRK); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1053 fOutputList->Add(fMultCorAfterVZTRKComp);
1055 fMultCorAfterCorrCut = new TH2F("fMultCorAfterCorrCut", "TPC vs Global multiplicity (After CorrCut); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1056 fOutputList->Add(fMultCorAfterCorrCut);
1058 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
1059 fOutputList->Add(fMultvsCentr);
1061 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1062 fOutputList->Add(fOpeningAngleLS);
1064 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1065 fOutputList->Add(fOpeningAngleULS);
1068 //----------------------------------------------------------------------------
1069 EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
1070 fOutputList->Add(EPVzA);
1071 EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
1072 fOutputList->Add(EPVzC);
1073 EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
1074 fOutputList->Add(EPTPC);
1077 EPVz = new TH1D("EPVz", "EPVz", 80, -2, 2);
1078 fOutputList->Add(EPVz);
1079 EPTPCp = new TH1D("EPTPCp", "EPTPCp", 80, -2, 2);
1080 fOutputList->Add(EPTPCp);
1081 EPTPCn = new TH1D("EPTPCn", "EPTPCn", 80, -2, 2);
1082 fOutputList->Add(EPTPCn);
1084 //----------------------------------------------------------------------------
1085 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
1086 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1087 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1088 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1089 fOutputList->Add(fSubEventDPhiv2);
1093 fSubEventDPhiv2new = new TProfile("fSubEventDPhiv2new", "fSubEventDPhiv2new", 3, 0, 3);
1094 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1095 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1096 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1097 fOutputList->Add(fSubEventDPhiv2new);
1098 //================================Event Plane with VZERO A & C=====================
1099 const Int_t nPtBins = 12;
1100 Double_t binsPt[nPtBins+1] = {0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10, 13};
1102 Int_t bins[3] = { 50, 50, nPtBins};
1103 Double_t xmin[3] = { -1., -1., 0};
1104 Double_t xmax[3] = { 1., 1., 13};
1105 fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
1106 // Set bin limits for axes which are not standard binned
1107 fV2Phi->SetBinEdges(2, binsPt);
1109 fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
1110 fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
1111 fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1112 fOutputList->Add(fV2Phi);
1116 //================================Event Plane with VZERO=====================
1117 // const Int_t nPtBins = 10;
1118 // Double_t binsPt[nPtBins+1] = {0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
1120 Int_t binsV[2] = { 50, nPtBins};
1121 Double_t xminV[2] = { -1., 0};
1122 Double_t xmaxV[2] = { 1., 13};
1123 fV2Phivzerotot = new THnSparseF("fV2Phivzerotot", "v2:pt", 2, binsV, xminV, xmaxV);
1124 // Set bin limits for axes which are not standard binned
1125 fV2Phivzerotot->SetBinEdges(1, binsPt);
1127 fV2Phivzerotot->GetAxis(0)->SetTitle("v_{2} (V0)");
1128 fV2Phivzerotot->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
1129 fOutputList->Add(fV2Phivzerotot);
1133 //----------------------------------------------------------------------------
1136 Int_t binsvElectH[8]={ 600, 200, 200 ,100, 100, 100, 10, 10}; //pt, E/p,TPCnSigma,M20,M02,Disp Phi-psiV0A ,Phi-PsiV0C,eta (commented)
1137 Double_t xminvElectH[8]={0, 0, -10 , 0, 0, 0, 0, 0};
1138 Double_t xmaxvElectH[8]={20, 2, 10 , 2, 2, 2, TMath::Pi(), TMath::Pi()};
1139 fSparseElectronHadron = new THnSparseD("ElectronHadron","ElectronHadron",8,binsvElectH,xminvElectH,xmaxvElectH);
1140 fSparseElectronHadron->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1141 fSparseElectronHadron->GetAxis(1)->SetTitle("EovP");
1142 fSparseElectronHadron->GetAxis(2)->SetTitle("TPCnSigma");
1143 fSparseElectronHadron->GetAxis(3)->SetTitle("M20");
1144 fSparseElectronHadron->GetAxis(4)->SetTitle("M02");
1145 fSparseElectronHadron->GetAxis(5)->SetTitle("Disp");
1146 fSparseElectronHadron->GetAxis(6)->SetTitle("phiminuspsi V0A");
1147 fSparseElectronHadron->GetAxis(7)->SetTitle("phiminuspsi V0C");
1148 fOutputList->Add(fSparseElectronHadron);
1150 //----------------------------------------------------------------------------
1151 //----------------------------------------------------------------------------
1153 Int_t binsvpurity[3]={ 600,200, 200}; //pt, E/p,TPCnSigma
1154 Double_t xminvpurity[3]={0, 0, -10};
1155 Double_t xmaxvpurity[3]={30, 2, 10};
1156 fSparseElectronpurity = new THnSparseD("Electronpurity","Electronpurity",3,binsvpurity,xminvpurity,xmaxvpurity);
1157 fSparseElectronpurity->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1158 fSparseElectronpurity->GetAxis(1)->SetTitle("EovP");
1159 fSparseElectronpurity->GetAxis(2)->SetTitle("TPCnSigma");
1160 fOutputList->Add(fSparseElectronpurity);
1162 //----------------------------------------------------------------------------
1165 Int_t binsphipsi[3] = { 100, 200, 6};
1166 Double_t xminphipsi[3] = { 0., 0, 0};
1167 Double_t xmaxphipsi[3] = { 10., 2, TMath::Pi()};
1168 fSparsephipsiULS = new THnSparseF("fSparsephipsiULS", "pt:eop:DeltaPhiULS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1169 fSparsephipsiULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1170 fSparsephipsiULS->GetAxis(1)->SetTitle("eop");
1171 fSparsephipsiULS->GetAxis(2)->SetTitle("DeltaPhiULS");
1172 fOutputList->Add(fSparsephipsiULS);
1174 fSparsephipsiLS = new THnSparseF("fSparsephipsiLS", "pt:eop:DeltaPhiLS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1175 fSparsephipsiLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1176 fSparsephipsiLS->GetAxis(1)->SetTitle("eop");
1177 fSparsephipsiLS->GetAxis(2)->SetTitle("DeltaPhiLS");
1178 fOutputList->Add(fSparsephipsiLS);
1180 Int_t binsmass[2] = { 100, 200};
1181 Double_t xminmass[2] = { 0., 0};
1182 Double_t xmaxmass[2] = { 10., 1.};
1183 fSparseMassULS = new THnSparseF("fSparseMassULS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1184 fSparseMassULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1185 fSparseMassULS->GetAxis(1)->SetTitle("mass");
1186 fOutputList->Add(fSparseMassULS);
1188 fSparseMassLS = new THnSparseF("fSparseMassLS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1189 fSparseMassLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1190 fSparseMassLS->GetAxis(1)->SetTitle("mass");
1191 fOutputList->Add(fSparseMassLS);
1195 PostData(1,fOutputList);
1196 // create and post flowevent
1197 fFlowEvent = new AliFlowEvent(10000);
1198 PostData(2, fFlowEvent);
1201 fFlowEventCont = new AliFlowEvent(10000);
1202 PostData(3, fFlowEventCont);
1206 //________________________________________________________________________
1207 void AliAnalysisTaskFlowTPCEMCalQCSP::Terminate(Option_t *)
1209 // Info("Terminate");
1210 AliAnalysisTaskSE::Terminate();
1212 //_____________________________________________________________________________
1213 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::PlotVZeroMultiplcities(const T* event) const
1215 // QA multiplicity plots
1216 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
1217 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
1219 //_____________________________________________________________________________
1220 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::SetNullCuts(T* event)
1223 if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
1224 fCutsRP->SetEvent(event, MCEvent());
1225 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
1226 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
1227 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
1228 fNullCuts->SetEvent(event, MCEvent());
1230 //_____________________________________________________________________________
1231 void AliAnalysisTaskFlowTPCEMCalQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
1233 //Prepare flow events
1234 FlowEv->ClearFast();
1235 FlowEv->Fill(fCutsRP, fNullCuts);
1236 FlowEv->SetReferenceMultiplicity(iMulti);
1237 FlowEv->DefineDeadZone(0, 0, 0, 0);
1238 // FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
1240 //_____________________________________________________________________________
1241 Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1243 // Check single track cuts for a given cut step
1244 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1245 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1248 //_________________________________________
1249 void AliAnalysisTaskFlowTPCEMCalQCSP::CheckCentrality(AliAODEvent* event, Bool_t ¢ralitypass)
1251 //============================Multiplicity TPV vs Global===============================================================================
1252 const Int_t nGoodTracks = event->GetNumberOfTracks();
1253 Float_t multTPC(0.); // tpc mult estimate
1254 Float_t multGlob(0.); // global multiplicity
1255 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1256 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1257 if (!trackAOD) continue;
1258 if (!(trackAOD->TestFilterBit(1))) continue;
1259 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
1262 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1263 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1264 if (!trackAOD) continue;
1265 if (!(trackAOD->TestFilterBit(16))) continue;
1266 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
1267 Double_t b[2] = {-99., -99.};
1268 Double_t bCov[3] = {-99., -99., -99.};
1269 if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
1270 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1273 fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1274 //============================================================================================================================
1275 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
1276 if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
1277 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
1278 // cout << "--------------Centrality evaluated-------------------------"<<endl;
1279 if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
1281 fCentralityNoPass->Fill(fCentrality);
1282 // cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
1283 centralitypass = kFALSE;
1286 // cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
1287 centralitypass = kTRUE;
1289 if (centralitypass){
1290 fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1291 fCentralityBeforePileup->Fill(fCentrality);
1292 }//...after centrality selectrion
1293 //============================================================================================================================
1294 //to remove the bias introduced by multeplicity outliers---------------------
1295 Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
1296 Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
1297 if (TMath::Abs(centv0 - centTrk) > 5.0){
1298 centralitypass = kFALSE;
1299 fCentralityNoPass->Fill(fCentrality);
1301 if (centralitypass){
1302 fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1303 fCentralityAfterVZTRK->Fill(fCentrality);
1304 }//...after centrality selectrion
1305 //============================================================================================================================
1307 if(fTrigger==1 || fTrigger==4 || fTrigger==5){
1308 if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
1309 // cout <<" Trigger ==" <<fTrigger<< endl;
1310 centralitypass = kFALSE;
1311 fCentralityNoPass->Fill(fCentrality);
1315 if(! (multTPC > (77.9 + 1.395*multGlob) && multTPC < (187.3 + 1.665*multGlob))){
1316 // cout <<" Trigger ==" <<fTrigger<< endl;
1317 centralitypass = kFALSE;
1318 fCentralityNoPass->Fill(fCentrality);
1322 if (centralitypass){
1323 fMultCorAfterCorrCut->Fill(multGlob, multTPC);
1324 fCentralityAfterCorrCut->Fill(fCentrality);
1325 }//...after CORR CUT
1326 //=================================All cuts are passed==================++++==================================================
1327 //=================================Now Centrality flattening for central trigger==================++++==================================================
1328 if(fTrigger==0 || fTrigger==4){
1329 if(!IsEventSelectedForCentrFlattening(fCentrality)){
1330 centralitypass = kFALSE;
1331 fCentralityNoPassForFlattening->Fill(fCentrality);
1334 //==============================fill histo after all cuts==============================++++==================================================
1336 fCentralityPass->Fill(fCentrality);
1337 fMultCorAfterCuts->Fill(multGlob, multTPC);
1338 fMultvsCentr->Fill(fCentrality, multTPC);
1341 //_____________________________________________________________________________
1342 void AliAnalysisTaskFlowTPCEMCalQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
1344 // Set a centrality range ]min, max] and define the method to use for centrality selection
1345 fCentralityMin = CentralityMin;
1346 fCentralityMax = CentralityMax;
1347 fkCentralityMethod = CentralityMethod;
1349 //_____________________________________________________________________________
1350 void AliAnalysisTaskFlowTPCEMCalQCSP::SetIDCuts(Double_t minTPC, Double_t maxTPC, Double_t minEovP, Double_t maxEovP, Double_t minM20, Double_t maxM20, Double_t minM02, Double_t maxM02, Double_t Dispersion)
1361 fDispersion = Dispersion;
1363 //_____________________________________________________________________________
1364 //_____________________________________________________________________________
1365 void AliAnalysisTaskFlowTPCEMCalQCSP::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
1366 // set the histo for centrality flattening
1367 // the centrality is flatten in the range minCentr,maxCentr
1368 // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
1369 // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
1370 // negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat)
1371 // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom
1373 if(maxCentr<minCentr){
1374 AliWarning("AliAnalysisCheckCorrdist::Wrong centralities values while setting the histogram for centrality flattening");
1377 if(fHistCentrDistr)delete fHistCentrDistr;
1378 fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
1379 fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
1380 Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
1381 Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
1382 fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
1383 Double_t ref=0.,bincont=0.,binrefwidth=1.;
1385 if(TMath::Abs(centrRef)<0.0001){
1386 binref=fHistCentrDistr->GetMinimumBin();
1387 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1388 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1390 else if(centrRef>0.){
1391 binref=h->FindBin(centrRef);
1392 if(binref<1||binref>h->GetNbinsX()){
1393 AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
1395 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1396 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1399 if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
1400 binref=fHistCentrDistr->GetMaximumBin();
1401 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1402 ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
1405 for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
1406 if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
1407 bincont=h->GetBinContent(j);
1408 fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
1409 fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
1412 h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
1416 fHistCentrDistr->SetBinContent(0,switchTRand);
1421 //-------------------------------------------------
1422 Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::IsEventSelectedForCentrFlattening(Float_t centvalue){
1424 // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
1425 // Can be faster if it was required that fHistCentrDistr covers
1426 // exactly the desired centrality range (e.g. part of the lines below should be done during the
1427 // setting of the histo) and TH1::SetMinimum called
1430 if(!fHistCentrDistr) return kTRUE;
1431 // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
1432 // if(maxbin>fHistCentrDistr->GetNbinsX()){
1433 // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
1436 Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
1437 Double_t bincont=fHistCentrDistr->GetBinContent(bin);
1438 Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
1440 if(fHistCentrDistr->GetBinContent(0)<-0.9999){
1441 if(gRandom->Uniform(1.)<bincont)return kTRUE;
1445 if(centDigits*100.<bincont)return kTRUE;
1449 //---------------------------------------------------------------------------