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 "AliAnalysisTaskFlowTPCTOFQCSP.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 "AliAODPid.h"
97 #include "AliFlowTrack.h"
98 #include "AliAnalysisTaskVnV0.h"
99 #include "AliSelectNonHFE.h"
102 class AliFlowTrackCuts;
106 ClassImp(AliAnalysisTaskFlowTPCTOFQCSP)
107 //________________________________________________________________________
108 AliAnalysisTaskFlowTPCTOFQCSP::AliAnalysisTaskFlowTPCTOFQCSP(const char *name)
109 : AliAnalysisTaskSE(name)
114 ,fIdentifiedAsOutInz(kFALSE)
115 ,fPassTheEventCut(kFALSE)
119 ,fCutsRP(0) // track cuts for reference particles
120 ,fNullCuts(0) // dummy cuts for flow event tracks
121 ,fFlowEvent(0) //! flow events (one for each inv mass band)
122 ,fkCentralityMethod(0)
142 ,fCentralityNoPass(0)
149 ,fMultCorAfterCuts(0)
157 ,fMultCorBeforeCuts(0)
162 ,fQAPIDSparse(kFALSE)
169 ,fOpeningAngleCut(0.1)
173 ,fNonHFE(new AliSelectNonHFE)
178 // fPID = new AliHFEpid("hfePid");
179 // Define input and output slots here
180 // Input slot #0 works with a TChain
181 DefineInput(0, TChain::Class());
182 // Output slot #0 id reserved by the base class for AOD
183 // Output slot #1 writes into a TH1 container
184 // DefineOutput(1, TH1I::Class());
185 DefineOutput(1, TList::Class());
186 DefineOutput(2, AliFlowEventSimple::Class());
189 //________________________________________________________________________
190 AliAnalysisTaskFlowTPCTOFQCSP::AliAnalysisTaskFlowTPCTOFQCSP()
191 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElectFlow")
196 ,fIdentifiedAsOutInz(kFALSE)
197 ,fPassTheEventCut(kFALSE)
201 ,fCutsRP(0) // track cuts for reference particles
202 ,fNullCuts(0) // dummy cuts for flow event tracks
203 ,fFlowEvent(0) //! flow events (one for each inv mass band)
204 ,fkCentralityMethod(0)
224 ,fCentralityNoPass(0)
231 ,fMultCorAfterCuts(0)
239 ,fMultCorBeforeCuts(0)
244 ,fQAPIDSparse(kFALSE)
251 ,fOpeningAngleCut(0.1)
255 ,fNonHFE(new AliSelectNonHFE)
258 //Default constructor
259 // fPID = new AliHFEpid("hfePid");
261 // Define input and output slots here
262 // Input slot #0 works with a TChain
263 DefineInput(0, TChain::Class());
264 // Output slot #0 id reserved by the base class for AOD
265 // Output slot #1 writes into a TH1 container
266 // DefineOutput(1, TH1I::Class());
267 DefineOutput(1, TList::Class());
268 DefineOutput(2, AliFlowEventSimple::Class());
269 //DefineOutput(3, TTree::Class());
271 //_________________________________________
273 AliAnalysisTaskFlowTPCTOFQCSP::~AliAnalysisTaskFlowTPCTOFQCSP()
282 if (fOutputList) delete fOutputList;
283 if (fFlowEvent) delete fFlowEvent;
286 //_________________________________________
288 void AliAnalysisTaskFlowTPCTOFQCSP::UserExec(Option_t*)
291 //Called for each event
293 // create pointer to event
295 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
300 printf("ERROR: fAOD not available\n");
306 AliError("HFE cuts not available");
310 // if(!fPID->IsInitialized())
312 // Initialize PID with the given run number
313 // AliWarning("PID not initialised, get from Run no");
314 // fPID->InitializePID(fAOD->GetRunNumber());
317 // cout << "kTrigger == " << fTrigger <<endl;
320 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ) return;
323 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
326 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA) ) return;
329 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) ) return;
333 //---------------CENTRALITY AND EVENT SELECTION-----------------------
337 Int_t fNOtrks = fAOD->GetNumberOfTracks();
339 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
340 if (!trkVtx || trkVtx->GetNContributors()<=0)return;
341 TString vtxTtl = trkVtx->GetTitle();
342 if (!vtxTtl.Contains("VertexerTracks"))return;
343 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
344 if (!spdVtx || spdVtx->GetNContributors()<=0)return;
345 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
346 vtxz = trkVtx->GetZ();
347 if(TMath::Abs(vtxz)>fVz)return;
350 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
351 if(fNOtrks<2) return;
354 Bool_t pass = kFALSE; //to select centrality
355 CheckCentrality(fAOD,pass);
362 PlotVZeroMultiplcities(fAOD);
365 PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent); //Calculate event plane Qvector and EP resolution for inclusive
368 fCFM->SetRecEventInfo(fAOD);
370 // Look for kink mother
371 Int_t numberofvertices = fAOD->GetNumberOfVertices();
372 Double_t listofmotherkink[numberofvertices];
373 Int_t numberofmotherkink = 0;
374 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
375 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
376 if(!aodvertex) continue;
377 if(aodvertex->GetType()==AliAODVertex::kKink) {
378 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
379 if(!mother) continue;
380 Int_t idmother = mother->GetID();
381 listofmotherkink[numberofmotherkink] = idmother;
382 //printf("ID %d\n",idmother);
383 numberofmotherkink++;
387 //=============================================V0EP from Alex======================================================================
388 Double_t qxEPa = 0, qyEPa = 0;
389 Double_t qxEPc = 0, qyEPc = 0;
391 Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
392 Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
395 Double_t Qx2 = 0, Qy2 = 0;
397 for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
399 AliAODTrack* aodTrack = fAOD->GetTrack(iT);
404 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
407 if (!aodTrack->TestFilterBit(128))
410 Qx2 += TMath::Cos(2*aodTrack->Phi());
411 Qy2 += TMath::Sin(2*aodTrack->Phi());
414 Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
416 EPVzA->Fill(evPlAngV0A);
417 EPVzC->Fill(evPlAngV0C);
418 EPTPC->Fill(evPlAngTPC);
420 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
421 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
422 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
423 //====================================================================================================================
425 AliAODTrack *track = NULL;
428 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
433 track = fAOD->GetTrack(iTracks);
436 printf("ERROR: Could not receive track %d\n", iTracks);
440 if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
442 //--------------------------------------hfe begin-----------------------------------------------------------
443 //==========================================================================================================
444 //======================================track cuts==========================================================
445 if(track->Eta()<-0.8 || track->Eta()>0.8) continue; //eta cuts on candidates
447 // RecKine: ITSTPC cuts
448 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
450 // Reject kink mother
451 Bool_t kinkmotherpass = kTRUE;
452 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
453 if(track->GetID() == listofmotherkink[kinkmother]) {
454 kinkmotherpass = kFALSE;
458 if(!kinkmotherpass) continue;
461 // if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; //deleted for DCA absence
462 // HFEcuts: ITS layers cuts
463 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
464 // HFE cuts: TPC PID cleanup
465 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
466 //==========================================================================================================
467 Double_t eta = track->Eta();
468 Double_t phi = track->Phi();
469 Double_t pt = track->Pt(); //pt track after cuts
470 if(pt<fpTCut) continue;
471 //==========================================================================================================
472 //=========================================PID==============================================================
473 if(track->GetTPCsignalN() < fTPCS) continue;
474 Float_t fTPCnSigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
475 Float_t fTOFnSigma = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
477 // Float_t eDEDX = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE);
479 // fTPCnsigma->Fill(track->P(),fTPCnSigma);
480 // fTOFns->Fill(track->P(),fTOFnSigma);
482 Double_t valPidSparse[4] = {
488 fQAPidSparse->Fill(valPidSparse);
490 /* || track->GetTPCsignal() < 72 || track->GetTPCsignal() > 95*/
492 if(fTOFnSigma < fminTOFnSigma || fTOFnSigma > fmaxTOFnSigma) continue;
493 if(fTPCnSigma < fminTPCnsigma || fTPCnSigma > fmaxTPCnsigma) continue; //cuts on nsigma tpc
494 //==========================================================================================================
495 //=========================================QA PID SPARSE====================================================
496 Float_t timeTOF = track->GetTOFsignal();
497 Double_t intTime[5] = {-99., -99., -99., -99., -99.};
498 track->GetIntegratedTimes(intTime);
499 Float_t timeElec = intTime[0];
500 Float_t intLength = 2.99792458e-2* timeElec;
502 if ((intLength > 0) && (timeTOF > 0))
503 beta = intLength/2.99792458e-2/timeTOF;
506 Double_t valPid[4] = {
508 track->GetTPCsignal(),
512 fQAPid->Fill(valPid);
517 if(phi<1.4 || phi >3.14)continue; //to have same EMCal phi acceptance
520 //==========================================================================================================
521 //============================Event Plane Method with V0====================================================
522 Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
523 Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
524 Double_t v2Phi[3] = {
529 //=========================================================================================================
530 fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
531 fTOFnsAft->Fill(track->P(),fTOFnSigma);
532 fTOFBetaAft->Fill(track->P(),beta);
533 fInclusiveElecPt->Fill(pt);
536 //=========================================================================================================
537 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
538 AliFlowTrack *sTrack = new AliFlowTrack();
540 sTrack->SetID(track->GetID());
541 sTrack->SetForRPSelection(kFALSE);
542 sTrack->SetForPOISelection(kTRUE);
543 sTrack->SetMass(263732);
544 for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
546 // cout << " no of rps " << iRPs << endl;
547 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
549 if (!iRP->InRPSelection()) continue;
550 if( sTrack->GetID() == iRP->GetID())
552 if(fDebug) printf(" was in RP set");
553 // cout << sTrack->GetID() <<" == " << iRP->GetID() << " was in RP set====REMOVED" <<endl;
554 iRP->SetForRPSelection(kFALSE);
555 fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
557 } //end of for loop on RPs
558 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
562 fNonHFE = new AliSelectNonHFE();
563 fNonHFE->SetAODanalysis(kTRUE);
564 fNonHFE->SetInvariantMassCut(fInvmassCut);
565 if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
566 //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
567 //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
568 fNonHFE->SetAlgorithm("DCA"); //KF
569 fNonHFE->SetPIDresponse(fPIDResponse);
570 fNonHFE->SetTrackCuts(-3,3);
572 fNonHFE->SetHistAngleBack(fOpeningAngleLS);
573 fNonHFE->SetHistAngle(fOpeningAngleULS);
574 //fNonHFE->SetHistDCABack(fDCABack);
575 //fNonHFE->SetHistDCA(fDCA);
576 fNonHFE->SetHistMassBack(fInvmassLS1);
577 fNonHFE->SetHistMass(fInvmassULS1);
579 fNonHFE->FindNonHFE(iTracks,track,fAOD);
581 // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
582 // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
583 // Bool_t fUlsIsPartner = kFALSE;
584 // Bool_t fLsIsPartner = kFALSE;
585 if(fNonHFE->IsULS()){
586 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
587 fULSElecPt->Fill(track->Pt());
592 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
593 fLSElecPt->Fill(track->Pt());
598 //=========================================================================================================
599 //----------------------Selection and Flow of Photonic Electrons-----------------------------
600 Bool_t fFlagPhotonicElec = kFALSE;
601 SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
602 if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
603 // Semi inclusive electron
604 if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
606 //-------------------------------------------------------------------------------------------
609 PostData(1, fOutputList);
610 PostData(2, fFlowEvent);
611 //----------hfe end---------
613 //_________________________________________
614 void AliAnalysisTaskFlowTPCTOFQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track, Bool_t &fFlagPhotonicElec)
617 //Identify non-heavy flavour electrons using Invariant mass method
618 Bool_t flagPhotonicElec = kFALSE;
620 for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
621 AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
623 printf("ERROR: Could not receive track %d\n", jTracks);
626 // if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; // TESTBIT FOR AOD double Counting
627 if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
628 if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)|| (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
631 if(jTracks == itrack) continue;
632 Double_t ptAsso=-999., nsigma=-999.0;
633 Double_t mass=-999., width = -999;
634 Double_t openingAngle = -999.;
635 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
638 ptAsso = trackAsso->Pt();
639 Short_t chargeAsso = trackAsso->Charge();
640 Short_t charge = track->Charge();
641 nsigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
643 if(trackAsso->GetTPCNcls() < 80) continue;
644 if(nsigma < -3 || nsigma > 3) continue;
645 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
646 if(ptAsso <0.3) continue;
648 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
649 if(charge>0) fPDGe1 = -11;
650 if(chargeAsso>0) fPDGe2 = -11;
652 if(charge == chargeAsso) fFlagLS = kTRUE;
653 if(charge != chargeAsso) fFlagULS = kTRUE;
655 AliKFParticle::SetField(fAOD->GetMagneticField());
656 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
657 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
658 AliKFParticle recg(ge1, ge2);
660 if(recg.GetNDF()<1) continue;
661 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
662 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
663 recg.GetMass(mass,width);
665 openingAngle = ge1.GetAngle(ge2);
666 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
667 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
668 if(fOP_angle){if(openingAngle > fOpeningAngleCut) continue;}
671 if(fFlagLS) fInvmassLS1->Fill(mass);
672 if(fFlagULS) fInvmassULS1->Fill(mass);
674 if(mass<fInvmassCut){
675 if(fFlagULS){fULSElecPt->Fill(track->Pt());}
676 if(fFlagLS){fLSElecPt->Fill(track->Pt());}
679 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
680 flagPhotonicElec = kTRUE;
684 fFlagPhotonicElec = flagPhotonicElec;
686 //___________________________________________
687 void AliAnalysisTaskFlowTPCTOFQCSP::UserCreateOutputObjects()
691 //----------hfe initialising begin---------
692 fNullCuts = new AliFlowTrackCuts("null_cuts");
694 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
695 cc->SetNbinsMult(10000);
697 cc->SetMultMax(10000);
703 cc->SetNbinsPhi(180);
705 cc->SetPhiMax(TMath::TwoPi());
716 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
717 AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
719 AliFatal("Input handler needed");
722 fPIDResponse=inputHandler->GetPIDResponse();
724 //pid response object
725 if (!fPIDResponse)AliError("PIDResponse object was not created");
728 //--------Initialize PID
729 // fPID->SetHasMCData(kFALSE);
730 // if(!fPID->GetNumberOfPIDdetectors())
732 // fPID->AddDetector("TPC", 0);
733 // fPID->AddDetector("EMCAL", 1);
736 // fPID->SortDetectors();
737 // fPIDqa = new AliHFEpidQAmanager();
738 // fPIDqa->Initialize(fPID);
742 //--------Initialize correction Framework and Cuts
743 fCFM = new AliCFManager;
744 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
745 fCFM->SetNStepParticle(kNcutSteps);
746 for(Int_t istep = 0; istep < kNcutSteps; istep++)
747 fCFM->SetParticleCutsList(istep, NULL);
750 AliWarning("Cuts not available. Default cuts will be used");
751 fCuts = new AliHFEcuts;
752 fCuts->CreateStandardCuts();
756 fCuts->Initialize(fCFM);
757 //----------hfe initialising end--------
758 //---------Output Tlist
759 fOutputList = new TList();
760 fOutputList->SetOwner();
761 // fOutputList->Add(fPIDqa->MakeList("PIDQA"));
763 fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
764 fOutputList->Add(fNoEvents);
766 // fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
767 // fOutputList->Add(fTPCnsigma);
769 fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
770 fOutputList->Add(fTPCnsigmaAft);
772 // fTOFns = new TH2F("fTOFns","track TOFnSigma",1000,0,50,200,-10,10);
773 // fOutputList->Add(fTOFns);
775 fTOFnsAft = new TH2F("fTOFnsAft","track TOFnSigma",1000,0,50,200,-10,10);
776 fOutputList->Add(fTOFnsAft);
778 fTOFBetaAft = new TH2F("fTOFBetaAft","track TOFBeta",1000,0,50,120,0,1.2);
779 fOutputList->Add(fTOFBetaAft);
781 fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
782 fOutputList->Add(fInclusiveElecPt);
784 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
785 fOutputList->Add(fPhotoElecPt);
787 fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
788 fOutputList->Add(fSemiInclElecPt);
790 fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
791 fOutputList->Add(fULSElecPt);
793 fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
794 fOutputList->Add(fLSElecPt);
796 fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
797 fOutputList->Add(fInvmassLS1);
799 fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
800 fOutputList->Add(fInvmassULS1);
802 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
803 fOutputList->Add(fCentralityPass);
805 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
806 fOutputList->Add(fCentralityNoPass);
808 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
809 fOutputList->Add(fPhi);
811 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
812 fOutputList->Add(fEta);
814 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
815 fOutputList->Add(fVZEROA);
817 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
818 fOutputList->Add(fVZEROC);
820 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
821 fOutputList->Add(fTPCM);
823 fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
824 fOutputList->Add(fvertex);
826 fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
827 fOutputList->Add(fMultCorBeforeCuts);
829 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
830 fOutputList->Add(fMultCorAfterCuts);
832 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
833 fOutputList->Add(fMultvsCentr);
835 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
836 fOutputList->Add(fOpeningAngleLS);
838 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
839 fOutputList->Add(fOpeningAngleULS);
843 //----------------------------------------------------------------------------
844 EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
845 fOutputList->Add(EPVzA);
846 EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
847 fOutputList->Add(EPVzC);
848 EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
849 fOutputList->Add(EPTPC);
850 //----------------------------------------------------------------------------
851 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
852 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
853 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
854 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
855 fOutputList->Add(fSubEventDPhiv2);
856 //================================Event Plane with VZERO=====================
857 const Int_t nPtBins = 10;
858 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};
860 Int_t bins[3] = { 50, 50, nPtBins};
861 Double_t xmin[3] = { -1., -1., 0};
862 Double_t xmax[3] = { 1., 1., 8};
863 fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
864 // Set bin limits for axes which are not standard binned
865 fV2Phi->SetBinEdges(2, binsPt);
867 fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
868 fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
869 fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
870 fOutputList->Add(fV2Phi);
871 //----------------------------------------------------------------------------
872 //----------------------------------------------------------------------------
874 Int_t binsQA[4] = { 150, 100, 120, 3};
875 Double_t xminQA[4] = { 0., 50, 0, -1.5};
876 Double_t xmaxQA[4] = { 15., 150, 1.2, 1.5};
877 fQAPid = new THnSparseF("fQAPid", "p:dEdx:beta:ch", 4, binsQA, xminQA, xmaxQA);
878 fQAPid->GetAxis(0)->SetTitle("p (Gev/c");
879 fQAPid->GetAxis(1)->SetTitle("dE/dx");
880 fQAPid->GetAxis(2)->SetTitle("#beta (TOF)");
881 fQAPid->GetAxis(3)->SetTitle("charge");
882 fOutputList->Add(fQAPid);
884 //===========================================================================
885 Int_t binsQA2[5] = { 200, 200, 200, 100};
886 Double_t xminQA2[5] = { 0., -10, -4, 0};
887 Double_t xmaxQA2[5] = { 20., 10, 4, 6.28};
888 fQAPidSparse = new THnSparseF("fQAPidSparse", "pt:p:tpcnsigma:tofnsigma:phi", 4, binsQA2, xminQA2, xmaxQA2);
889 fQAPidSparse->GetAxis(0)->SetTitle("pt (Gev/c");
890 fQAPidSparse->GetAxis(1)->SetTitle("tpc nsigma");
891 fQAPidSparse->GetAxis(2)->SetTitle("tofnsigma");
892 fQAPidSparse->GetAxis(3)->SetTitle("phi");
893 fOutputList->Add(fQAPidSparse);
894 //===========================================================================
895 PostData(1,fOutputList);
896 // create and post flowevent
897 fFlowEvent = new AliFlowEvent(10000);
898 PostData(2, fFlowEvent);
900 //________________________________________________________________________
901 void AliAnalysisTaskFlowTPCTOFQCSP::Terminate(Option_t *)
903 // Info("Terminate");
904 AliAnalysisTaskSE::Terminate();
906 //_____________________________________________________________________________
907 template <typename T> void AliAnalysisTaskFlowTPCTOFQCSP::PlotVZeroMultiplcities(const T* event) const
909 // QA multiplicity plots
910 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
911 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
913 //_____________________________________________________________________________
914 template <typename T> void AliAnalysisTaskFlowTPCTOFQCSP::SetNullCuts(T* event)
917 if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
918 fCutsRP->SetEvent(event, MCEvent());
919 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
920 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
921 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
922 fNullCuts->SetEvent(event, MCEvent());
924 //_____________________________________________________________________________
925 void AliAnalysisTaskFlowTPCTOFQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
927 //Prepare flow events
929 FlowEv->Fill(fCutsRP, fNullCuts);
930 FlowEv->SetReferenceMultiplicity(iMulti);
931 FlowEv->DefineDeadZone(0, 0, 0, 0);
932 // FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
934 //_____________________________________________________________________________
935 Bool_t AliAnalysisTaskFlowTPCTOFQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
937 // Check single track cuts for a given cut step
938 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
939 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
942 //_________________________________________
943 void AliAnalysisTaskFlowTPCTOFQCSP::CheckCentrality(AliAODEvent* event, Bool_t ¢ralitypass)
945 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
946 if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
947 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
948 // cout << "--------------Centrality evaluated-------------------------"<<endl;
949 if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
951 fCentralityNoPass->Fill(fCentrality);
952 //cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
953 centralitypass = kFALSE;
956 // cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
957 centralitypass = kTRUE;
959 //to remove the bias introduced by multeplicity outliers---------------------
960 Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
961 Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
963 if (TMath::Abs(centv0 - centTrk) > 5.0){
964 centralitypass = kFALSE;
965 fCentralityNoPass->Fill(fCentrality);
967 const Int_t nGoodTracks = event->GetNumberOfTracks();
969 Float_t multTPC(0.); // tpc mult estimate
970 Float_t multGlob(0.); // global multiplicity
971 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
972 AliAODTrack* trackAOD = event->GetTrack(iTracks);
973 if (!trackAOD) continue;
974 if (!(trackAOD->TestFilterBit(1))) continue;
975 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;
978 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
979 AliAODTrack* trackAOD = event->GetTrack(iTracks);
980 if (!trackAOD) continue;
981 if (!(trackAOD->TestFilterBit(16))) continue;
982 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;
983 Double_t b[2] = {-99., -99.};
984 Double_t bCov[3] = {-99., -99., -99.};
985 if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
986 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
989 // printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
990 // if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))){ 2010
991 if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
992 centralitypass = kFALSE;
993 fCentralityNoPass->Fill(fCentrality);
995 fMultCorBeforeCuts->Fill(multGlob, multTPC);
998 fCentralityPass->Fill(fCentrality);
999 fMultCorAfterCuts->Fill(multGlob, multTPC);
1000 fMultvsCentr->Fill(fCentrality, multTPC);
1003 //_____________________________________________________________________________
1004 void AliAnalysisTaskFlowTPCTOFQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
1006 // Set a centrality range ]min, max] and define the method to use for centrality selection
1007 fCentralityMin = CentralityMin;
1008 fCentralityMax = CentralityMax;
1009 fkCentralityMethod = CentralityMethod;
1011 //_____________________________________________________________________________
1012 void AliAnalysisTaskFlowTPCTOFQCSP::SetIDCuts(Double_t minTPCnsigma, Double_t maxTPCnsigma, Double_t minTOFnSigma, Double_t maxTOFnSigma)
1015 fminTPCnsigma = minTPCnsigma;
1016 fmaxTPCnsigma = maxTPCnsigma;
1017 fminTOFnSigma = minTOFnSigma;
1018 fmaxTOFnSigma = maxTOFnSigma;
1020 //_____________________________________________________________________________
1022 //_____________________________________________________________________________