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 TOF //
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 "AliAnalysisTaskFlowITSTPCTOFQCSP.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 "AliGeomManager.h"
62 #include "TGeoManager.h"
65 #include "AliEMCALTrack.h"
66 //#include "AliEMCALTracker.h"
68 #include "AliKFParticle.h"
69 #include "AliKFVertex.h"
70 #include "AliHFEcontainer.h"
71 #include "AliHFEcuts.h"
72 #include "AliHFEpid.h"
73 #include "AliHFEpidBase.h"
74 #include "AliHFEpidQAmanager.h"
75 #include "AliHFEtools.h"
76 #include "AliCFContainer.h"
77 #include "AliCFManager.h"
78 #include "AliKFParticle.h"
79 #include "AliKFVertex.h"
80 #include "AliCentrality.h"
81 #include "AliVEvent.h"
83 #include "AliMCEvent.h"
85 #include "AliFlowCandidateTrack.h"
86 #include "AliFlowTrackCuts.h"
87 #include "AliFlowEventSimple.h"
88 #include "AliFlowCommonConstants.h"
89 #include "AliFlowEvent.h"
92 #include "AliESDVZERO.h"
93 #include "AliAODVZERO.h"
95 #include "AliPIDResponse.h"
96 #include "AliFlowTrack.h"
97 #include "AliAnalysisTaskVnV0.h"
98 #include "AliSelectNonHFE.h"
101 class AliFlowTrackCuts;
105 ClassImp(AliAnalysisTaskFlowITSTPCTOFQCSP)
106 //________________________________________________________________________
107 AliAnalysisTaskFlowITSTPCTOFQCSP::AliAnalysisTaskFlowITSTPCTOFQCSP(const char *name)
108 : AliAnalysisTaskSE(name)
113 ,fIdentifiedAsOutInz(kFALSE)
114 ,fPassTheEventCut(kFALSE)
118 ,fCutsRP(0) // track cuts for reference particles
119 ,fNullCuts(0) // dummy cuts for flow event tracks
120 ,fFlowEvent(0) //! flow events (one for each inv mass band)
121 ,fkCentralityMethod(0)
139 ,fTPCnsigmaVSptAft(0)
144 ,fCentralityNoPass(0)
151 ,fMultCorAfterCuts(0)
159 ,fMultCorBeforeCuts(0)
161 ,fminITSnsigmaLowpT(-1)
162 ,fmaxITSnsigmaLowpT(1)
163 ,fminITSnsigmaHighpT(-2)
164 ,fmaxITSnsigmaHighpT(2)
165 ,fminTPCnsigmaLowpT(-1)
166 ,fmaxTPCnsigmaLowpT(3)
167 ,fminTPCnsigmaHighpT(0)
168 ,fmaxTPCnsigmaHighpT(3)
169 //,fQAPIDSparse(kFALSE)
175 ,fOpeningAngleCut(0.1)
179 ,fNonHFE(new AliSelectNonHFE)
182 ,fTPCnsigmaAftITSTOF(0)
187 ,fTPCvsITSafterTOF(0)
190 ,fMultCorAfterCentrBeforeCuts(0)
191 ,fMultCorAfterVZTRKComp(0)
192 ,fCentralityBeforePileup(0)
193 ,fCentralityAfterVZTRK(0)
194 ,fCentralityAfterCorrCut(0)
195 ,fMultCorAfterCorrCut(0)
199 ,fSubEventDPhiv2new(0)
201 ,fHistCentrDistr(0x0)
202 ,fCentralityNoPassForFlattening(0)
214 fPID = new AliHFEpid("hfePid");
215 // Define input and output slots here
216 // Input slot #0 works with a TChain
217 DefineInput(0, TChain::Class());
218 // Output slot #0 id reserved by the base class for AOD
219 // Output slot #1 writes into a TH1 container
220 // DefineOutput(1, TH1I::Class());
221 DefineOutput(1, TList::Class());
222 DefineOutput(2, AliFlowEventSimple::Class());
225 //________________________________________________________________________
226 AliAnalysisTaskFlowITSTPCTOFQCSP::AliAnalysisTaskFlowITSTPCTOFQCSP()
227 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElectFlow")
232 ,fIdentifiedAsOutInz(kFALSE)
233 ,fPassTheEventCut(kFALSE)
237 ,fCutsRP(0) // track cuts for reference particles
238 ,fNullCuts(0) // dummy cuts for flow event tracks
239 ,fFlowEvent(0) //! flow events (one for each inv mass band)
240 ,fkCentralityMethod(0)
258 ,fTPCnsigmaVSptAft(0)
263 ,fCentralityNoPass(0)
270 ,fMultCorAfterCuts(0)
278 ,fMultCorBeforeCuts(0)
280 ,fminITSnsigmaLowpT(-1)
281 ,fmaxITSnsigmaLowpT(1)
282 ,fminITSnsigmaHighpT(-2)
283 ,fmaxITSnsigmaHighpT(2)
284 ,fminTPCnsigmaLowpT(-1)
285 ,fmaxTPCnsigmaLowpT(3)
286 ,fminTPCnsigmaHighpT(0)
287 ,fmaxTPCnsigmaHighpT(3)
288 //,fQAPIDSparse(kFALSE)
294 ,fOpeningAngleCut(0.1)
298 ,fNonHFE(new AliSelectNonHFE)
301 ,fTPCnsigmaAftITSTOF(0)
306 ,fTPCvsITSafterTOF(0)
309 ,fMultCorAfterCentrBeforeCuts(0)
310 ,fMultCorAfterVZTRKComp(0)
311 ,fCentralityBeforePileup(0)
312 ,fCentralityAfterVZTRK(0)
313 ,fCentralityAfterCorrCut(0)
314 ,fMultCorAfterCorrCut(0)
318 ,fSubEventDPhiv2new(0)
320 ,fHistCentrDistr(0x0)
321 ,fCentralityNoPassForFlattening(0)
331 //Default constructor
332 fPID = new AliHFEpid("hfePid");
334 // Define input and output slots here
335 // Input slot #0 works with a TChain
336 DefineInput(0, TChain::Class());
337 // Output slot #0 id reserved by the base class for AOD
338 // Output slot #1 writes into a TH1 container
339 // DefineOutput(1, TH1I::Class());
340 DefineOutput(1, TList::Class());
341 DefineOutput(2, AliFlowEventSimple::Class());
342 //DefineOutput(3, TTree::Class());
344 //_________________________________________
346 AliAnalysisTaskFlowITSTPCTOFQCSP::~AliAnalysisTaskFlowITSTPCTOFQCSP()
352 // delete fPIDResponse;
355 if (fOutputList) delete fOutputList;
356 if (fFlowEvent) delete fFlowEvent;
359 //_________________________________________
361 void AliAnalysisTaskFlowITSTPCTOFQCSP::UserExec(Option_t*)
364 //Called for each event
366 // create pointer to event
368 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
373 printf("ERROR: fAOD not available\n");
379 AliError("HFE cuts not available");
383 if(!fPID->IsInitialized())
385 // Initialize PID with the given run number
386 AliWarning("PID not initialised, get from Run no");
387 fPID->InitializePID(fAOD->GetRunNumber());
390 // cout << "kTrigger == " << fTrigger <<endl;
393 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral)) return;
397 if ( !(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny) ) return;
399 TString firedTriggerClasses = static_cast<const AliAODEvent*>(InputEvent())->GetFiredTriggerClasses();
401 if ( ! ( firedTriggerClasses.Contains("CVLN_B2-B-NOPF-ALLNOTRD") || firedTriggerClasses.Contains("CVLN_R1-B-NOPF-ALLNOTRD") || firedTriggerClasses.Contains("CSEMI_R1-B-NOPF-ALLNOTRD") ) ) return;
404 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA)) return;
407 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return;
410 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kCentral | AliVEvent::kSemiCentral))) return;
413 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
419 //---------------CENTRALITY AND EVENT SELECTION-----------------------
423 Int_t fNOtrks = fAOD->GetNumberOfTracks();
425 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
426 if (!trkVtx || trkVtx->GetNContributors()<=0)return;
427 TString vtxTtl = trkVtx->GetTitle();
428 if (!vtxTtl.Contains("VertexerTracks"))return;
429 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
430 if (!spdVtx || spdVtx->GetNContributors()<=0)return;
431 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
432 vtxz = trkVtx->GetZ();
433 if(TMath::Abs(vtxz)>fVz)return;
436 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
437 if(fNOtrks<2) return;
440 Bool_t pass = kFALSE; //to select centrality
441 CheckCentrality(fAOD,pass);
448 PlotVZeroMultiplcities(fAOD);
451 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent); //Calculate event plane Qvector and EP resolution for inclusive
453 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
456 AliDebug(1, "Using default PID Response");
457 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
460 fPID->SetPIDResponse(pidResponse);
462 fCFM->SetRecEventInfo(fAOD);
464 // Look for kink mother
465 Int_t numberofvertices = fAOD->GetNumberOfVertices();
466 Double_t listofmotherkink[numberofvertices];
467 Int_t numberofmotherkink = 0;
468 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
469 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
470 if(!aodvertex) continue;
471 if(aodvertex->GetType()==AliAODVertex::kKink) {
472 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
473 if(!mother) continue;
474 Int_t idmother = mother->GetID();
475 listofmotherkink[numberofmotherkink] = idmother;
476 //printf("ID %d\n",idmother);
477 numberofmotherkink++;
481 //=============================================V0EP from Alex======================================================================
482 Double_t qxEPa = 0, qyEPa = 0;
483 Double_t qxEPc = 0, qyEPc = 0;
484 Double_t qxEP = 0, qyEP = 0;
486 Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
487 Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
488 Double_t evPlAngV0 = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 10, 2, qxEP, qyEP);
491 Double_t Qx2 = 0, Qy2 = 0;
492 Double_t Qx2p = 0, Qy2p = 0;
493 Double_t Qx2n = 0, Qy2n = 0;
495 for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
497 AliAODTrack* aodTrack = fAOD->GetTrack(iT);
502 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
505 if (!aodTrack->TestFilterBit(128))
509 if(aodTrack->Eta()>0 && aodTrack->Eta()<0.8){
511 Qx2p += TMath::Cos(2*aodTrack->Phi());
512 Qy2p += TMath::Sin(2*aodTrack->Phi());
514 if(aodTrack->Eta()<0 && aodTrack->Eta()> -0.8){
516 Qx2n += TMath::Cos(2*aodTrack->Phi());
517 Qy2n += TMath::Sin(2*aodTrack->Phi());
521 Qx2 += TMath::Cos(2*aodTrack->Phi());
522 Qy2 += TMath::Sin(2*aodTrack->Phi());
529 Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
530 Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
531 Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
533 EPVzA->Fill(evPlAngV0A);
534 EPVzC->Fill(evPlAngV0C);
535 EPTPC->Fill(evPlAngTPC);
537 EPTPCn->Fill(evPlAngTPCn);
538 EPTPCp->Fill(evPlAngTPCp);
539 EPVz->Fill(evPlAngV0);
541 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
542 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
543 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
546 fSubEventDPhiv2new->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))); // vzero - tpcp
547 fSubEventDPhiv2new->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))); // vzero - tpcn
548 fSubEventDPhiv2new->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
551 //====================================================================================================================
553 AliAODTrack *track = NULL;
556 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
558 track = fAOD->GetTrack(iTracks);
561 printf("ERROR: Could not receive track %d\n", iTracks);
564 if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
566 //--------------------------------------hfe begin-----------------------------------------------------------
567 //==========================================================================================================
568 //======================================track cuts==========================================================
569 if(track->Eta()<-0.8 || track->Eta()>0.8) continue; //eta cuts on candidates
571 // RecKine: ITSTPC cuts
572 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
574 // Reject kink mother
575 Bool_t kinkmotherpass = kTRUE;
576 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
577 if(track->GetID() == listofmotherkink[kinkmother]) {
578 kinkmotherpass = kFALSE;
582 if(!kinkmotherpass) continue;
585 // if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; //deleted for DCA absence
586 // HFEcuts: ITS layers cuts
587 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
588 // HFE cuts: TPC PID cleanup
589 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
590 //==========================================================================================================
591 Double_t eta = track->Eta();
592 Double_t phi = track->Phi();
595 if(phi<1.4 || phi >3.14)continue; //to have same EMCal phi acceptance
600 Double_t pt = track->Pt(); //pt track after cuts
601 if(pt<fpTCutmin || pt>fpTCutmax) continue;
602 //==========================================================================================================
603 //=========================================PID==============================================================
604 if(track->GetTPCsignalN() < fTPCS) continue;
605 Float_t fITSnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasITS(track, AliPID::kElectron) : 1000;
606 Float_t fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
607 Float_t fTOFnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTOF(track, AliPID::kElectron) : 1000;
608 // Float_t eDEDX = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE);
609 fITSnsigma->Fill(track->P(),fITSnSigma);
610 fTPCnsigma->Fill(track->P(),fTPCnSigma);
611 fTOFns->Fill(track->P(),fTOFnSigma);
612 fITSvsTOF->Fill(fTOFnSigma,fITSnSigma);
613 fTPCvsITS->Fill(fTPCnSigma,fITSnSigma);
614 fTPCvsTOF->Fill(fTPCnSigma,fTOFnSigma);
617 if(fTOFnSigma < fminTOFnSigma || fTOFnSigma > fmaxTOFnSigma) continue;
618 }//cuts on nsigma tof full pt range
620 fITSnsigmaAftTOF->Fill(track->P(),fITSnSigma);
621 fTPCnsigmaAftTOF->Fill(track->P(),fTPCnSigma);
622 fTPCvsITSafterTOF->Fill(fTPCnSigma,fITSnSigma);
624 Double_t valPidSparse[3] = {
629 fQAPidSparse->Fill(valPidSparse);
633 if(fITSnSigma < fminITSnsigmaLowpT || fITSnSigma > fmaxITSnsigmaLowpT)continue;
634 }//cuts on nsigma its low pt
636 if(fITSnSigma < fminITSnsigmaHighpT || fITSnSigma > fmaxITSnsigmaHighpT)continue;
637 }//cuts on nsigma its high pt
638 fTPCnsigmaAftITSTOF->Fill(track->P(),fTPCnSigma);
639 if(pt >= 0.25 && pt < 1.5){
640 if(fTPCnSigma < fminTPCnsigmaLowpT || fTPCnSigma > fmaxTPCnsigmaLowpT) continue;
641 }//cuts on nsigma tpc lowpt
643 if(fTPCnSigma < fminTPCnsigmaHighpT || fTPCnSigma > fmaxTPCnsigmaHighpT) continue;
644 }//cuts on nsigma tpc high pt
645 fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
646 fTPCnsigmaVSptAft->Fill(pt,fTPCnSigma);
648 //==========================================================================================================
649 //=========================================QA PID SPARSE====================================================
650 Float_t timeTOF = track->GetTOFsignal();
651 Double_t intTime[5] = {-99., -99., -99., -99., -99.};
652 track->GetIntegratedTimes(intTime);
653 Float_t timeElec = intTime[0];
654 Float_t intLength = 2.99792458e-2* timeElec;
656 if ((intLength > 0) && (timeTOF > 0))
657 beta = intLength/2.99792458e-2/timeTOF;
660 // Double_t valPid[4] = {
662 // track->GetTPCsignal(),
666 // fQAPid->Fill(valPid);
670 fITSnsigmaAft->Fill(track->P(),fITSnSigma);
671 fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
672 fTOFnsAft->Fill(track->P(),fTOFnSigma);
673 fTOFBetaAft->Fill(track->P(),beta);
674 fInclusiveElecPt->Fill(pt);
677 //=========================================================================================================
678 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
679 AliFlowTrack *sTrack = new AliFlowTrack();
681 sTrack->SetID(track->GetID());
682 sTrack->SetForRPSelection(kTRUE);
683 sTrack->SetForPOISelection(kTRUE);
684 sTrack->SetMass(263732);
685 for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
687 // cout << " no of rps " << iRPs << endl;
688 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
690 if (!iRP->InRPSelection()) continue;
691 if( sTrack->GetID() == iRP->GetID())
693 if(fDebug) printf(" was in RP set");
694 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set====REMOVED" <<endl;
695 iRP->SetForRPSelection(kFALSE);
696 // fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
698 } //end of for loop on RPs
699 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
700 fFlowEvent->SetNumberOfPOIs(fFlowEvent->GetNumberOfPOIs()+1);
701 //============================Event Plane Method with V0====================================================
702 Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
703 Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
704 Double_t v2Phi[3] = {
710 Double_t v2PhiVz = TMath::Cos(2*(phi - evPlAngV0));
711 Double_t v2PhiV0tot[2] = {
714 fV2Phivzerotot->Fill(v2PhiV0tot);
716 //==========================================================================================================
717 //=========================================================================================================
721 fNonHFE = new AliSelectNonHFE();
722 fNonHFE->SetAODanalysis(kTRUE);
723 fNonHFE->SetInvariantMassCut(fInvmassCut);
724 if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
725 //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
726 //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
727 fNonHFE->SetAlgorithm("DCA"); //KF
728 fNonHFE->SetPIDresponse(pidResponse);
729 fNonHFE->SetTrackCuts(-3,3);
731 fNonHFE->SetHistAngleBack(fOpeningAngleLS);
732 fNonHFE->SetHistAngle(fOpeningAngleULS);
733 //fNonHFE->SetHistDCABack(fDCABack);
734 //fNonHFE->SetHistDCA(fDCA);
735 fNonHFE->SetHistMassBack(fInvmassLS1);
736 fNonHFE->SetHistMass(fInvmassULS1);
738 fNonHFE->FindNonHFE(iTracks,track,fAOD);
740 // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
741 // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
742 // Bool_t fUlsIsPartner = kFALSE;
743 // Bool_t fLsIsPartner = kFALSE;
744 if(fNonHFE->IsULS()){
745 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
746 fULSElecPt->Fill(track->Pt());
751 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
752 fLSElecPt->Fill(track->Pt());
757 //=========================================================================================================
758 //----------------------Selection and Flow of Photonic Electrons-----------------------------
759 Bool_t fFlagPhotonicElec = kFALSE;
760 SelectPhotonicElectron(iTracks,track,fTPCnSigma,evPlAngV0,fFlagPhotonicElec);
761 if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
762 // Semi inclusive electron
763 if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
765 //-------------------------------------------------------------------------------------------
768 PostData(1, fOutputList);
769 PostData(2, fFlowEvent);
771 //----------hfe end---------
773 //_________________________________________
774 void AliAnalysisTaskFlowITSTPCTOFQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track,Float_t fTPCnSigma,Double_t evPlAngV0, Bool_t &fFlagPhotonicElec)
777 //Identify non-heavy flavour electrons using Invariant mass method
778 Bool_t flagPhotonicElec = kFALSE;
780 for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
781 AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
783 printf("ERROR: Could not receive track %d\n", jTracks);
786 // if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
787 if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
790 if(!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
793 if(!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
795 // if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)|| (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
798 if(jTracks == itrack) continue;
799 Double_t ptAsso=-999., nsigma=-999.0;
800 Double_t mass=-999., width = -999;
801 Double_t openingAngle = -999.;
802 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
805 ptAsso = trackAsso->Pt();
806 Short_t chargeAsso = trackAsso->Charge();
807 Short_t charge = track->Charge();
808 // nsigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
809 nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
812 //if(trackAsso->GetTPCNcls() < 80) continue;
813 if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
814 if(nsigma < -3 || nsigma > 3) continue;
815 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
816 if(ptAsso < fptminAsso) continue;
818 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
819 if(charge>0) fPDGe1 = -11;
820 if(chargeAsso>0) fPDGe2 = -11;
822 if(charge == chargeAsso) fFlagLS = kTRUE;
823 if(charge != chargeAsso) fFlagULS = kTRUE;
825 AliKFParticle::SetField(fAOD->GetMagneticField());
826 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
827 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
828 AliKFParticle recg(ge1, ge2);
830 if(recg.GetNDF()<1) continue;
831 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
832 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
833 recg.GetMass(mass,width);
835 openingAngle = ge1.GetAngle(ge2);
836 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
837 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
838 if(fOP_angle){if(openingAngle > fOpeningAngleCut) continue;}
841 if(fFlagLS) fInvmassLS1->Fill(mass);
842 if(fFlagULS) fInvmassULS1->Fill(mass);
845 Double_t MassSparseULS[3] = {
849 fSparseMassULS->Fill(MassSparseULS);
852 Double_t MassSparseLS[3] = {
856 fSparseMassLS->Fill(MassSparseLS);
860 if(mass<fInvmassCut){
861 if(fFlagULS){fULSElecPt->Fill(track->Pt());}
862 if(fFlagLS){fLSElecPt->Fill(track->Pt());}
866 Double_t phi = track->Phi();
867 Float_t DeltaPhi_eEP = TVector2::Phi_0_2pi(phi - evPlAngV0);
868 if(DeltaPhi_eEP > TMath::Pi()) {DeltaPhi_eEP = DeltaPhi_eEP - TMath::Pi();}
871 if(mass<fInvmassCut){
873 Double_t ulsSparse[3] = {
878 fSparsephipsiULS->Fill(ulsSparse);
881 Double_t lsSparse[3] = {
886 fSparsephipsiLS->Fill(lsSparse);
892 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
893 flagPhotonicElec = kTRUE;
897 fFlagPhotonicElec = flagPhotonicElec;
899 //___________________________________________
900 void AliAnalysisTaskFlowITSTPCTOFQCSP::UserCreateOutputObjects()
904 //----------hfe initialising begin---------
905 fNullCuts = new AliFlowTrackCuts("null_cuts");
907 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
908 cc->SetNbinsMult(10000);
910 cc->SetMultMax(10000);
916 cc->SetNbinsPhi(180);
918 cc->SetPhiMax(TMath::TwoPi());
929 // AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
930 // AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
931 // if (!inputHandler){
932 // AliFatal("Input handler needed");
935 // fPIDResponse=inputHandler->GetPIDResponse();
937 //pid response object
938 // if (!fPIDResponse)AliError("PIDResponse object was not created");
941 //--------Initialize PID
942 fPID->SetHasMCData(kFALSE);
943 if(!fPID->GetNumberOfPIDdetectors())
945 fPID->AddDetector("ITS", 0);
946 fPID->AddDetector("TOF", 1);
947 fPID->AddDetector("TPC", 2);
951 fPID->SortDetectors();
952 fPIDqa = new AliHFEpidQAmanager();
953 fPIDqa->Initialize(fPID);
957 //--------Initialize correction Framework and Cuts
958 fCFM = new AliCFManager;
959 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
960 fCFM->SetNStepParticle(kNcutSteps);
961 for(Int_t istep = 0; istep < kNcutSteps; istep++)
962 fCFM->SetParticleCutsList(istep, NULL);
965 AliWarning("Cuts not available. Default cuts will be used");
966 fCuts = new AliHFEcuts;
967 fCuts->CreateStandardCuts();
971 fCuts->Initialize(fCFM);
972 //----------hfe initialising end--------
973 //---------Output Tlist
974 fOutputList = new TList();
975 fOutputList->SetOwner();
976 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
978 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
979 fOutputList->Add(fNoEvents);
981 fITSnsigma = new TH2F("fITSnsigma", "ITS - n sigma before HFE pid",600,0,6,400,-20,20);
982 fOutputList->Add(fITSnsigma);
984 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",600,0,6,400,-20,20);
985 fOutputList->Add(fTPCnsigma);
987 fITSnsigmaAft = new TH2F("fITSnsigmaAft", "ITS - n sigma after HFE pid",1000,0,10,300,-10,20);
988 fOutputList->Add(fITSnsigmaAft);
989 fITSvsTOF = new TH2F("fITSvsTOF", "ITS tof",400,-20,20,400,-20,20);
990 fOutputList->Add(fITSvsTOF);
991 fTPCvsITS = new TH2F("TPCvsITS", "TPC ITS",400,-20,20,400,-20,20);
992 fOutputList->Add(fTPCvsITS);
993 fTPCvsTOF = new TH2F("TPCvsTOF", "TPC TOF",400,-20,20,400,-20,20);
994 fOutputList->Add(fTPCvsTOF);
995 fTPCvsITSafterTOF = new TH2F("TPCvsITSafterTOF", "TPC ITS",400,-20,20,400,-20,20);
996 fOutputList->Add(fTPCvsITSafterTOF);
999 fITSnsigmaAftTOF = new TH2F("fITSnsigmaAftTOF", "ITS - n sigma after HFE pid",600,0,6,400,-20,20);
1000 fOutputList->Add(fITSnsigmaAftTOF);
1002 fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
1003 fOutputList->Add(fTPCnsigmaAft);
1005 fTPCnsigmaVSptAft = new TH2F("fTPCnsigmaVSptAft", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
1006 fOutputList->Add(fTPCnsigmaVSptAft);
1008 fTPCnsigmaAftITSTOF = new TH2F("fTPCnsigmaAftITSTOF", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
1009 fOutputList->Add(fTPCnsigmaAftITSTOF);
1011 fTPCnsigmaAftTOF = new TH2F("fTPCnsigmaAftTOF", "TPC - n sigma after HFE pid",600,0,6,400,-20,20);
1012 fOutputList->Add(fTPCnsigmaAftTOF);
1014 fTOFns = new TH2F("fTOFns","track TOFnSigma",600,0,6,400,-20,20);
1015 fOutputList->Add(fTOFns);
1017 fTOFnsAft = new TH2F("fTOFnsAft","track TOFnSigma",600,0,6,400,-20,20);
1018 fOutputList->Add(fTOFnsAft);
1020 fTOFBetaAft = new TH2F("fTOFBetaAft","track TOFBeta",600,0,6,120,0,1.2);
1021 fOutputList->Add(fTOFBetaAft);
1023 fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",100,0,5);
1024 fOutputList->Add(fInclusiveElecPt);
1026 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,5);
1027 fOutputList->Add(fPhotoElecPt);
1029 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,5);
1030 fOutputList->Add(fSemiInclElecPt);
1032 fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",100,0,5);
1033 fOutputList->Add(fULSElecPt);
1035 fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",100,0,5);
1036 fOutputList->Add(fLSElecPt);
1038 fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
1039 fOutputList->Add(fInvmassLS1);
1041 fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
1042 fOutputList->Add(fInvmassULS1);
1044 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
1045 fOutputList->Add(fCentralityPass);
1047 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
1048 fOutputList->Add(fCentralityNoPass);
1050 fCentralityNoPassForFlattening = new TH1F("fCentralityNoPassForFlattening", "Centrality No Pass for flattening", 101, -1, 100);
1051 fOutputList->Add(fCentralityNoPassForFlattening);
1053 fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
1054 fOutputList->Add(fCentralityBeforePileup);
1056 fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
1057 fOutputList->Add(fCentralityAfterVZTRK);
1059 fCentralityAfterCorrCut = new TH1F("fCentralityAfterCorrCut", "fCentralityAfterCorrCut Pass", 101, -1, 100);
1060 fOutputList->Add(fCentralityAfterCorrCut);
1062 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
1063 fOutputList->Add(fPhi);
1065 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
1066 fOutputList->Add(fEta);
1068 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
1069 fOutputList->Add(fVZEROA);
1071 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
1072 fOutputList->Add(fVZEROC);
1074 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
1075 fOutputList->Add(fTPCM);
1077 fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
1078 fOutputList->Add(fvertex);
1080 fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1081 fOutputList->Add(fMultCorBeforeCuts);
1083 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1084 fOutputList->Add(fMultCorAfterCuts);
1086 fMultCorAfterCentrBeforeCuts = new TH2F("fMultCorAfterCentrBeforeCuts", "TPC vs Global multiplicity (After CC before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1087 fOutputList->Add(fMultCorAfterCentrBeforeCuts);
1089 fMultCorAfterVZTRKComp = new TH2F("fMultCorAfterVZTRKComp", "TPC vs Global multiplicity (After V0-TRK); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1090 fOutputList->Add(fMultCorAfterVZTRKComp);
1092 fMultCorAfterCorrCut = new TH2F("fMultCorAfterCorrCut", "TPC vs Global multiplicity (After CorrCut); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1093 fOutputList->Add(fMultCorAfterCorrCut);
1095 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
1096 fOutputList->Add(fMultvsCentr);
1098 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1099 fOutputList->Add(fOpeningAngleLS);
1101 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1102 fOutputList->Add(fOpeningAngleULS);
1106 //----------------------------------------------------------------------------
1107 EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
1108 fOutputList->Add(EPVzA);
1109 EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
1110 fOutputList->Add(EPVzC);
1111 EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
1112 fOutputList->Add(EPTPC);
1113 EPVz = new TH1D("EPVz", "EPVz", 80, -2, 2);
1114 fOutputList->Add(EPVz);
1115 EPTPCp = new TH1D("EPTPCp", "EPTPCp", 80, -2, 2);
1116 fOutputList->Add(EPTPCp);
1117 EPTPCn = new TH1D("EPTPCn", "EPTPCn", 80, -2, 2);
1118 fOutputList->Add(EPTPCn);
1121 //----------------------------------------------------------------------------
1122 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
1123 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1124 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1125 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1126 fOutputList->Add(fSubEventDPhiv2);
1128 fSubEventDPhiv2new = new TProfile("fSubEventDPhiv2new", "fSubEventDPhiv2new", 3, 0, 3);
1129 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1130 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1131 fSubEventDPhiv2new->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1132 fOutputList->Add(fSubEventDPhiv2new);
1134 //================================Event Plane with VZERO=====================
1135 const Int_t nPtBins = 12;
1136 Double_t binsPt[nPtBins+1] = {0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5};
1138 Int_t bins[3] = { 50, 50, nPtBins};
1139 Double_t xmin[3] = { -1., -1., 0};
1140 Double_t xmax[3] = { 1., 1., 5};
1141 fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
1142 // Set bin limits for axes which are not standard binned
1143 fV2Phi->SetBinEdges(2, binsPt);
1145 fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
1146 fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
1147 fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1148 fOutputList->Add(fV2Phi);
1149 //================================Event Plane with VZERO=====================
1150 Int_t binsV[2] = { 50, 100};
1151 Double_t xminV[2] = { -1., 0};
1152 Double_t xmaxV[2] = { 1., 5};
1153 fV2Phivzerotot = new THnSparseF("fV2Phivzerotot", "v2:pt", 2, binsV, xminV, xmaxV);
1154 // Set bin limits for axes which are not standard binned
1155 //fV2Phivzerotot->SetBinEdges(1, binsPt);
1157 fV2Phivzerotot->GetAxis(0)->SetTitle("v_{2} (V0)");
1158 fV2Phivzerotot->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
1159 fOutputList->Add(fV2Phivzerotot);
1163 //----------------------------------------------------------------------------
1164 //----------------------------------------------------------------------------
1165 // if(fQAPIDSparse){
1166 // Int_t binsQA[4] = { 150, 100, 120, 3};
1167 // Double_t xminQA[4] = { 0., 50, 0, -1.5};
1168 // Double_t xmaxQA[4] = { 15., 150, 1.2, 1.5};
1169 // fQAPid = new THnSparseF("fQAPid", "p:dEdx:beta:ch", 4, binsQA, xminQA, xmaxQA);
1170 // fQAPid->GetAxis(0)->SetTitle("p (Gev/c");
1171 // fQAPid->GetAxis(1)->SetTitle("dE/dx");
1172 // fQAPid->GetAxis(2)->SetTitle("#beta (TOF)");
1173 // fQAPid->GetAxis(3)->SetTitle("charge");
1174 // fOutputList->Add(fQAPid);
1176 //===========================================================================
1177 Int_t binsQA2[3] = { 100, 40, 150/*, 60*/};
1178 Double_t xminQA2[3] = { 0., -2, -15/*, -3*/};
1179 Double_t xmaxQA2[3] = { 5., 2, 15/*, 3*/};
1180 fQAPidSparse = new THnSparseF("fQAPidSparse", "pt:itsnsigma:tpcnsigma", 3, binsQA2, xminQA2, xmaxQA2);
1181 fQAPidSparse->GetAxis(0)->SetTitle("pt (Gev/c)");
1182 fQAPidSparse->GetAxis(1)->SetTitle("itsnsigma");
1183 fQAPidSparse->GetAxis(2)->SetTitle("tpcnsigma");
1184 fOutputList->Add(fQAPidSparse);
1185 //===========================================================================
1189 Int_t binsphipsi[3] = { 100, 150, 6};
1190 Double_t xminphipsi[3] = { 0., -15, 0};
1191 Double_t xmaxphipsi[3] = { 5., 15, TMath::Pi()};
1192 fSparsephipsiULS = new THnSparseF("fSparsephipsiULS", "pt:tpcnsigma:DeltaPhiULS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1193 fSparsephipsiULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1194 fSparsephipsiULS->GetAxis(1)->SetTitle("tpcnsigma");
1195 fSparsephipsiULS->GetAxis(2)->SetTitle("DeltaPhiULS");
1196 fOutputList->Add(fSparsephipsiULS);
1198 fSparsephipsiLS = new THnSparseF("fSparsephipsiLS", "pt:tpcnsigma:DeltaPhiLS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1199 fSparsephipsiLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1200 fSparsephipsiLS->GetAxis(1)->SetTitle("tpcnsigma");
1201 fSparsephipsiLS->GetAxis(2)->SetTitle("DeltaPhiLS");
1202 fOutputList->Add(fSparsephipsiLS);
1205 Int_t binsmass[2] = { 100, 200};
1206 Double_t xminmass[2] = { 0., 0};
1207 Double_t xmaxmass[2] = { 5., 1.};
1208 fSparseMassULS = new THnSparseF("fSparseMassULS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1209 fSparseMassULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1210 fSparseMassULS->GetAxis(1)->SetTitle("mass");
1211 fOutputList->Add(fSparseMassULS);
1213 fSparseMassLS = new THnSparseF("fSparseMassLS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1214 fSparseMassLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1215 fSparseMassLS->GetAxis(1)->SetTitle("mass");
1216 fOutputList->Add(fSparseMassLS);
1219 PostData(1,fOutputList);
1220 // create and post flowevent
1221 fFlowEvent = new AliFlowEvent(10000);
1222 PostData(2, fFlowEvent);
1225 //________________________________________________________________________
1226 void AliAnalysisTaskFlowITSTPCTOFQCSP::Terminate(Option_t *)
1228 // Info("Terminate");
1229 AliAnalysisTaskSE::Terminate();
1231 //_____________________________________________________________________________
1232 template <typename T> void AliAnalysisTaskFlowITSTPCTOFQCSP::PlotVZeroMultiplcities(const T* event) const
1234 // QA multiplicity plots
1235 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
1236 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
1238 //_____________________________________________________________________________
1239 template <typename T> void AliAnalysisTaskFlowITSTPCTOFQCSP::SetNullCuts(T* event)
1242 if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
1243 fCutsRP->SetEvent(event, MCEvent());
1244 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
1245 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
1246 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
1247 fNullCuts->SetEvent(event, MCEvent());
1249 //_____________________________________________________________________________
1250 void AliAnalysisTaskFlowITSTPCTOFQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
1252 //Prepare flow events
1253 FlowEv->ClearFast();
1254 FlowEv->Fill(fCutsRP, fNullCuts);
1255 FlowEv->SetReferenceMultiplicity(iMulti);
1256 FlowEv->DefineDeadZone(0, 0, 0, 0);
1257 // FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
1259 //_____________________________________________________________________________
1260 Bool_t AliAnalysisTaskFlowITSTPCTOFQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1262 // Check single track cuts for a given cut step
1263 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1264 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1267 //_________________________________________
1268 void AliAnalysisTaskFlowITSTPCTOFQCSP::CheckCentrality(AliAODEvent* event, Bool_t ¢ralitypass)
1270 //============================Multiplicity TPV vs Global===============================================================================
1271 const Int_t nGoodTracks = event->GetNumberOfTracks();
1272 Float_t multTPC(0.); // tpc mult estimate
1273 Float_t multGlob(0.); // global multiplicity
1274 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1275 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1276 if (!trackAOD) continue;
1277 if (!(trackAOD->TestFilterBit(1))) continue;
1278 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;
1281 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1282 AliAODTrack* trackAOD = event->GetTrack(iTracks);
1283 if (!trackAOD) continue;
1284 if (!(trackAOD->TestFilterBit(16))) continue;
1285 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;
1286 Double_t b[2] = {-99., -99.};
1287 Double_t bCov[3] = {-99., -99., -99.};
1288 if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
1289 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1292 fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1293 //============================================================================================================================
1294 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
1295 if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
1296 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
1297 // cout << "--------------Centrality evaluated-------------------------"<<endl;
1298 if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
1300 fCentralityNoPass->Fill(fCentrality);
1301 // cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
1302 centralitypass = kFALSE;
1305 // cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
1306 centralitypass = kTRUE;
1308 if (centralitypass){
1309 fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1310 fCentralityBeforePileup->Fill(fCentrality);
1311 }//...after centrality selectrion
1312 //============================================================================================================================
1313 //to remove the bias introduced by multeplicity outliers---------------------
1314 Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
1315 Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
1316 if (TMath::Abs(centv0 - centTrk) > 5.0){
1317 centralitypass = kFALSE;
1318 fCentralityNoPass->Fill(fCentrality);
1320 if (centralitypass){
1321 fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1322 fCentralityAfterVZTRK->Fill(fCentrality);
1323 }//...after centrality selectrion
1324 //============================================================================================================================
1326 if(fTrigger==1 || fTrigger==4 || fTrigger==5){
1327 if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
1328 // cout <<" Trigger ==" <<fTrigger<< endl;
1329 centralitypass = kFALSE;
1330 fCentralityNoPass->Fill(fCentrality);
1334 if(! (multTPC > (77.9 + 1.395*multGlob) && multTPC < (187.3 + 1.665*multGlob))){
1335 // cout <<" Trigger ==" <<fTrigger<< endl;
1336 centralitypass = kFALSE;
1337 fCentralityNoPass->Fill(fCentrality);
1341 if (centralitypass){
1342 fMultCorAfterCorrCut->Fill(multGlob, multTPC);
1343 fCentralityAfterCorrCut->Fill(fCentrality);
1344 }//...after CORR CUT
1345 //=================================All cuts are passed==================++++==================================================
1346 //=================================Now Centrality flattening for central trigger==================++++==================================================
1347 if(fTrigger==0 || fTrigger==4){
1348 if(!IsEventSelectedForCentrFlattening(fCentrality)){
1349 centralitypass = kFALSE;
1350 fCentralityNoPassForFlattening->Fill(fCentrality);
1353 //==============================fill histo after all cuts==============================++++==================================================
1355 fCentralityPass->Fill(fCentrality);
1356 fMultCorAfterCuts->Fill(multGlob, multTPC);
1357 fMultvsCentr->Fill(fCentrality, multTPC);
1360 //_____________________________________________________________________________
1361 void AliAnalysisTaskFlowITSTPCTOFQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
1363 // Set a centrality range ]min, max] and define the method to use for centrality selection
1364 fCentralityMin = CentralityMin;
1365 fCentralityMax = CentralityMax;
1366 fkCentralityMethod = CentralityMethod;
1368 //_____________________________________________________________________________
1369 void AliAnalysisTaskFlowITSTPCTOFQCSP::SetIDCuts(Double_t minTOFnSigma, Double_t maxTOFnSigma, Double_t minITSnsigmalowpt, Double_t maxITSnsigmalowpt, Double_t minITSnsigmahighpt, Double_t maxITSnsigmahighpt, Double_t minTPCnsigmalowpt, Double_t maxTPCnsigmalowpt, Double_t minTPCnsigmahighpt, Double_t maxTPCnsigmahighpt)
1372 fminTOFnSigma = minTOFnSigma;
1373 fmaxTOFnSigma = maxTOFnSigma;
1374 fminITSnsigmaLowpT = minITSnsigmalowpt;
1375 fmaxITSnsigmaLowpT = maxITSnsigmalowpt;
1376 fminITSnsigmaHighpT = minITSnsigmahighpt;
1377 fmaxITSnsigmaHighpT = maxITSnsigmahighpt;
1378 fminTPCnsigmaLowpT = minTPCnsigmalowpt;
1379 fmaxTPCnsigmaLowpT = maxTPCnsigmalowpt;
1380 fminTPCnsigmaHighpT = minTPCnsigmahighpt;
1381 fmaxTPCnsigmaHighpT = maxTPCnsigmahighpt;
1384 //_____________________________________________________________________________
1385 //_____________________________________________________________________________
1386 void AliAnalysisTaskFlowITSTPCTOFQCSP::SetpTCuttrack(Double_t ptmin, Double_t ptmax)
1392 //_____________________________________________________________________________
1393 //_____________________________________________________________________________
1394 void AliAnalysisTaskFlowITSTPCTOFQCSP::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
1395 // set the histo for centrality flattening
1396 // the centrality is flatten in the range minCentr,maxCentr
1397 // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
1398 // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
1399 // 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)
1400 // 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
1402 if(maxCentr<minCentr){
1403 AliWarning("AliAnalysisCheckCorrdist::Wrong centralities values while setting the histogram for centrality flattening");
1406 if(fHistCentrDistr)delete fHistCentrDistr;
1407 fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
1408 fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
1409 Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
1410 Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
1411 fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
1412 Double_t ref=0.,bincont=0.,binrefwidth=1.;
1414 if(TMath::Abs(centrRef)<0.0001){
1415 binref=fHistCentrDistr->GetMinimumBin();
1416 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1417 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1419 else if(centrRef>0.){
1420 binref=h->FindBin(centrRef);
1421 if(binref<1||binref>h->GetNbinsX()){
1422 AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
1424 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1425 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1428 if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
1429 binref=fHistCentrDistr->GetMaximumBin();
1430 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1431 ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
1434 for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
1435 if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
1436 bincont=h->GetBinContent(j);
1437 fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
1438 fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
1441 h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
1445 fHistCentrDistr->SetBinContent(0,switchTRand);
1450 //-------------------------------------------------
1451 Bool_t AliAnalysisTaskFlowITSTPCTOFQCSP::IsEventSelectedForCentrFlattening(Float_t centvalue){
1453 // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
1454 // Can be faster if it was required that fHistCentrDistr covers
1455 // exactly the desired centrality range (e.g. part of the lines below should be done during the
1456 // setting of the histo) and TH1::SetMinimum called
1459 if(!fHistCentrDistr) return kTRUE;
1460 // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
1461 // if(maxbin>fHistCentrDistr->GetNbinsX()){
1462 // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
1465 Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
1466 Double_t bincont=fHistCentrDistr->GetBinContent(bin);
1467 Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
1469 if(fHistCentrDistr->GetBinContent(0)<-0.9999){
1470 if(gRandom->Uniform(1.)<bincont)return kTRUE;
1474 if(centDigits*100.<bincont)return kTRUE;
1478 //---------------------------------------------------------------------------
1481 //_____________________________________________________________________________