1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 #include "AliTwoParticlePIDCorr.h"
17 #include "AliVParticle.h"
30 #include "AliCentrality.h"
31 #include "Riostream.h"
33 #include "AliAnalysisDataSlot.h"
34 #include "AliAnalysisDataContainer.h"
37 #include "AliCFContainer.h"
39 #include "THnSparse.h"
43 #include "AliESDpid.h"
44 #include "AliAODpidUtil.h"
45 #include <AliPIDResponse.h>
46 #include "AliPIDCombined.h"
48 #include <AliInputEventHandler.h>
49 #include "AliAODInputHandler.h"
51 #include "AliAnalysisTaskSE.h"
52 #include "AliAnalysisManager.h"
53 #include "AliCentrality.h"
55 #include "AliVEvent.h"
56 #include "AliAODEvent.h"
57 #include "AliAODTrack.h"
58 #include "AliVTrack.h"
60 #include "AliAODcascade.h"
62 #include "THnSparse.h"
64 #include "AliAODMCHeader.h"
65 #include "AliAODMCParticle.h"
66 #include "AliMCEventHandler.h"
67 #include "AliMCEvent.h"
68 #include "AliMCParticle.h"
69 #include "TParticle.h"
70 #include <TDatabasePDG.h>
71 #include <TParticlePDG.h>
73 #include "AliGenCocktailEventHeader.h"
74 #include "AliGenEventHeader.h"
75 #include "AliCollisionGeometry.h"
76 #include "AliOADBContainer.h"
78 #include "AliEventPoolManager.h"
79 #include "AliAnalysisUtils.h"
80 using namespace AliPIDNameSpace;
83 ClassImp(AliTwoParticlePIDCorr)
84 ClassImp(LRCParticlePID)
86 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
87 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
88 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
89 //Source code::dphicorrelations,VnV0, TaskBFpsi, AliHelperPID,
91 //________________________________________________________________________
92 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
97 fCentralityMethod("V0A"),
98 fPPVsMultUtils(kFALSE),
100 fRequestEventPlane(kFALSE),
101 fRequestEventPlanemixing(kFALSE),
102 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
107 fSharedClusterCut(-1),
108 fSharedTPCmapCut(-1),
109 fSharedfraction_Pair_cut(-1),
111 skipParticlesAbove(0),
116 ffilltrigassoUNID(kFALSE),
117 ffilltrigUNIDassoID(kFALSE),
118 ffilltrigIDassoUNID(kTRUE),
119 ffilltrigIDassoID(kFALSE),
120 ffilltrigIDassoIDMCTRUTH(kFALSE),
121 fMaxNofMixingTracks(50000),
122 fPtOrderMCTruth(kTRUE),
123 fPtOrderDataReco(kTRUE),
124 fWeightPerEvent(kFALSE),
125 fTriggerSpeciesSelection(kFALSE),
126 fAssociatedSpeciesSelection(kFALSE),
127 fRandomizeReactionPlane(kFALSE),
128 fTriggerSpecies(SpPion),
129 fAssociatedSpecies(SpPion),
132 fSelectHighestPtTrig(kFALSE),
133 fcontainPIDtrig(kTRUE),
134 fcontainPIDasso(kFALSE),
136 frejectPileUp(kFALSE),
141 fselectprimaryTruth(kTRUE),
142 fonlyprimarydatareco(kFALSE),
145 ffillhistQAReco(kFALSE),
146 ffillhistQATruth(kFALSE),
147 kTrackVariablesPair(0),
181 fhistJetTrigestimate(0),
182 fTwoTrackDistancePtdip(0x0),
183 fTwoTrackDistancePtdipmix(0x0),
184 fCentralityCorrelation(0x0),
185 fHistVZEROAGainEqualizationMap(0),
186 fHistVZEROCGainEqualizationMap(0),
187 fHistVZEROChannelGainEqualizationMap(0),
188 fCentralityWeights(0),
191 fHistEQVZEROvsTPCmultiplicity(0x0),
192 fHistEQVZEROAvsTPCmultiplicity(0x0),
193 fHistEQVZEROCvsTPCmultiplicity(0x0),
194 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
195 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
196 fHistVZEROCvsVZEROAmultiplicity(0x0),
197 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
198 fHistVZEROSignal(0x0),
199 fHistEventPlaneTruth(0x0),
200 fHistPsiMinusPhi(0x0),
215 gReactionPlane(999.),
243 fControlConvResoncances(0),
258 fCorrelatonTruthPrimary(0),
259 fCorrelatonTruthPrimarymix(0),
265 fTHnCorrIDUNIDmix(0),
267 fTHnTrigcountMCTruthPrim(0),
270 fAnalysisType("AOD"),
272 ftwoTrackEfficiencyCutDataReco(kTRUE),
273 fTwoTrackCutMinRadius(0.8),
274 fTwoTrackCutMaxRadius(2.5),
275 twoTrackEfficiencyCutValue(0.02),
281 fRequestTOFPID(kTRUE),
282 fPIDType(NSigmaTPCTOF),
283 fFIllPIDQAHistos(kTRUE),
286 fdiffPIDcutvalues(kFALSE),
291 fHighPtKaonNSigmaPID(-1),
292 fHighPtKaonSigma(3.5),
293 fUseExclusiveNSigma(kFALSE),
294 fRemoveTracksT0Fill(kFALSE),
296 fTriggerSelectCharge(0),
297 fAssociatedSelectCharge(0),
298 fTriggerRestrictEta(-1),
299 fEtaOrdering(kFALSE),
300 fCutConversions(kFALSE),
301 fCutResonances(kFALSE),
302 fRejectResonanceDaughters(-1),
304 fInjectedSignals(kFALSE),
305 fRemoveWeakDecays(kFALSE),
306 fRemoveDuplicates(kFALSE),
307 fapplyTrigefficiency(kFALSE),
308 fapplyAssoefficiency(kFALSE),
309 ffillefficiency(kFALSE),
310 fmesoneffrequired(kFALSE),
311 fkaonprotoneffrequired(kFALSE),
315 fUsev0DaughterPID(kFALSE),
316 fMinPtDaughter(1.0),// v0 related cut starts here
319 fMaxDCADaughter(1.0),
322 fHistRawPtCentInvK0s(0x0),
323 fHistRawPtCentInvLambda(0x0),
324 fHistRawPtCentInvAntiLambda(0x0),
325 fHistFinalPtCentInvK0s(0x0),
326 fHistFinalPtCentInvLambda(0x0),
327 fHistFinalPtCentInvAntiLambda(0x0),
331 fCutctauAntiLambda(7.8),
338 for ( Int_t i = 0; i < 16; i++) {
342 for ( Int_t i = 0; i < 6; i++ ){
343 fTrackHistEfficiency[i] = NULL;
344 effcorection[i]=NULL;
347 for ( Int_t i = 0; i < 2; i++ ){
348 fTwoTrackDistancePt[i]=NULL;
349 fTwoTrackDistancePtmix[i]=NULL;
352 for(Int_t ipart=0;ipart<NSpecies;ipart++)
353 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
354 fnsigmas[ipart][ipid]=999.;
356 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
358 for(Int_t i = 0; i != 2; ++i)
359 for(Int_t j = 0; j != 2; ++j)
360 for(Int_t iC = 0; iC < 9; iC++){
361 fMeanQ[iC][i][j] = 0.;
362 fWidthQ[iC][i][j] = 1.;
363 fMeanQv3[iC][i][j] = 0.;
364 fWidthQv3[iC][i][j] = 1.;
368 //________________________________________________________________________
369 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
370 :AliAnalysisTaskSE(name),
374 fCentralityMethod("V0A"),
375 fPPVsMultUtils(kFALSE),
377 fRequestEventPlane(kFALSE),
378 fRequestEventPlanemixing(kFALSE),
379 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
384 fSharedClusterCut(-1),
385 fSharedTPCmapCut(-1),
386 fSharedfraction_Pair_cut(-1),
388 skipParticlesAbove(0),
393 ffilltrigassoUNID(kFALSE),
394 ffilltrigUNIDassoID(kFALSE),
395 ffilltrigIDassoUNID(kTRUE),
396 ffilltrigIDassoID(kFALSE),
397 ffilltrigIDassoIDMCTRUTH(kFALSE),
398 fMaxNofMixingTracks(50000),
399 fPtOrderMCTruth(kTRUE),
400 fPtOrderDataReco(kTRUE),
401 fWeightPerEvent(kFALSE),
402 fTriggerSpeciesSelection(kFALSE),
403 fAssociatedSpeciesSelection(kFALSE),
404 fRandomizeReactionPlane(kFALSE),
405 fTriggerSpecies(SpPion),
406 fAssociatedSpecies(SpPion),
409 fSelectHighestPtTrig(kFALSE),
410 fcontainPIDtrig(kTRUE),
411 fcontainPIDasso(kFALSE),
413 frejectPileUp(kFALSE),
418 fselectprimaryTruth(kTRUE),
419 fonlyprimarydatareco(kFALSE),
422 ffillhistQAReco(kFALSE),
423 ffillhistQATruth(kFALSE),
424 kTrackVariablesPair(0),
458 fhistJetTrigestimate(0),
459 fTwoTrackDistancePtdip(0x0),
460 fTwoTrackDistancePtdipmix(0x0),
461 fCentralityCorrelation(0x0),
462 fHistVZEROAGainEqualizationMap(0),
463 fHistVZEROCGainEqualizationMap(0),
464 fHistVZEROChannelGainEqualizationMap(0),
465 fCentralityWeights(0),
468 fHistEQVZEROvsTPCmultiplicity(0x0),
469 fHistEQVZEROAvsTPCmultiplicity(0x0),
470 fHistEQVZEROCvsTPCmultiplicity(0x0),
471 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
472 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
473 fHistVZEROCvsVZEROAmultiplicity(0x0),
474 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
475 fHistVZEROSignal(0x0),
476 fHistEventPlaneTruth(0x0),
477 fHistPsiMinusPhi(0x0),
492 gReactionPlane(999.),
520 fControlConvResoncances(0),
535 fCorrelatonTruthPrimary(0),
536 fCorrelatonTruthPrimarymix(0),
542 fTHnCorrIDUNIDmix(0),
544 fTHnTrigcountMCTruthPrim(0),
547 fAnalysisType("AOD"),
549 ftwoTrackEfficiencyCutDataReco(kTRUE),
550 fTwoTrackCutMinRadius(0.8),
551 fTwoTrackCutMaxRadius(2.5),
552 twoTrackEfficiencyCutValue(0.02),
558 fRequestTOFPID(kTRUE),
559 fPIDType(NSigmaTPCTOF),
560 fFIllPIDQAHistos(kTRUE),
563 fdiffPIDcutvalues(kFALSE),
568 fHighPtKaonNSigmaPID(-1),
569 fHighPtKaonSigma(3.5),
570 fUseExclusiveNSigma(kFALSE),
571 fRemoveTracksT0Fill(kFALSE),
573 fTriggerSelectCharge(0),
574 fAssociatedSelectCharge(0),
575 fTriggerRestrictEta(-1),
576 fEtaOrdering(kFALSE),
577 fCutConversions(kFALSE),
578 fCutResonances(kFALSE),
579 fRejectResonanceDaughters(-1),
581 fInjectedSignals(kFALSE),
582 fRemoveWeakDecays(kFALSE),
583 fRemoveDuplicates(kFALSE),
584 fapplyTrigefficiency(kFALSE),
585 fapplyAssoefficiency(kFALSE),
586 ffillefficiency(kFALSE),
587 fmesoneffrequired(kFALSE),
588 fkaonprotoneffrequired(kFALSE),
592 fUsev0DaughterPID(kFALSE),
593 fMinPtDaughter(1.0),// v0 related cut starts here
596 fMaxDCADaughter(1.0),
599 fHistRawPtCentInvK0s(0x0),
600 fHistRawPtCentInvLambda(0x0),
601 fHistRawPtCentInvAntiLambda(0x0),
602 fHistFinalPtCentInvK0s(0x0),
603 fHistFinalPtCentInvLambda(0x0),
604 fHistFinalPtCentInvAntiLambda(0x0),
608 fCutctauAntiLambda(7.8),
614 // The last in the above list should not have a comma after it
618 for ( Int_t i = 0; i < 16; i++) {
622 for ( Int_t i = 0; i < 6; i++ ){
623 fTrackHistEfficiency[i] = NULL;
624 effcorection[i]=NULL;
628 for ( Int_t i = 0; i < 2; i++ ){
629 fTwoTrackDistancePt[i]=NULL;
630 fTwoTrackDistancePtmix[i]=NULL;
633 for(Int_t ipart=0;ipart<NSpecies;ipart++)
634 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
635 fnsigmas[ipart][ipid]=999.;
637 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
639 for(Int_t i = 0; i != 2; ++i)
640 for(Int_t j = 0; j != 2; ++j)
641 for(Int_t iC = 0; iC < 9; iC++){
642 fMeanQ[iC][i][j] = 0.;
643 fWidthQ[iC][i][j] = 1.;
644 fMeanQv3[iC][i][j] = 0.;
645 fWidthQv3[iC][i][j] = 1.;
649 // Define input and output slots here (never in the dummy constructor)
650 // Input slot #0 works with a TChain - it is connected to the default input container
651 // Output slot #1 writes into a TH1 container
652 DefineInput(0, TChain::Class());
654 DefineOutput(1, TList::Class()); // for output list
655 DefineOutput(2, TList::Class());
656 DefineOutput(3, TList::Class());
660 //________________________________________________________________________
661 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
663 // Destructor. Clean-up the output list, but not the histograms that are put inside
664 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
665 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
670 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
675 if(fRequestEventPlane){
676 if (fList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
681 if (fPID) delete fPID;
682 if (fPIDCombined) delete fPIDCombined;
685 //________________________________________________________________________
687 //////////////////////////////////////////////////////////////////////////////////////////////////
689 TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
690 // returns histo named name
691 return (TH2F*) fOutputList->FindObject(name);
694 //////////////////////////////////////////////////////////////////////////////////////////////////
696 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
700 // Puts the argument in the range [-pi/2,3 pi/2].
703 if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
704 if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
709 //________________________________________________________________________
710 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
713 // Called once (on the worker node)
715 // global switch disabling the reference
716 // (to avoid "Replacing existing TH1" if several wagons are created in train)
717 Bool_t oldStatus = TH1::AddDirectoryStatus();
718 TH1::AddDirectory(kFALSE);
720 const Int_t nPsiTOF = 10;
721 const Int_t nCentrBin = 9;
724 fOutput = new TList();
725 fOutput->SetOwner(); // IMPORTANT!
727 fOutputList = new TList;
728 fOutputList->SetOwner();
729 fOutputList->SetName("PIDQAList");
731 if(fRequestEventPlane){
734 fList->SetName("EPQAList");
736 fEventCounter = new TH1F("fEventCounter","EventCounter", 19, 0.5,19.5);
737 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
738 fEventCounter->GetXaxis()->SetBinLabel(3,"After PileUP Cut");//only for Data
739 fEventCounter->GetXaxis()->SetBinLabel(5,"Have A Vertex");
740 fEventCounter->GetXaxis()->SetBinLabel(7,"After vertex Cut");
741 fEventCounter->GetXaxis()->SetBinLabel(9,"Getting centrality");
742 fEventCounter->GetXaxis()->SetBinLabel(11,"After centrality flattening");
743 fEventCounter->GetXaxis()->SetBinLabel(13,"Within 0-100% centrality");
744 fEventCounter->GetXaxis()->SetBinLabel(15,"Event Analyzed");
745 //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
746 fOutput->Add(fEventCounter);
748 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
749 fOutput->Add(fEtaSpectrasso);
751 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
752 fOutput->Add(fphiSpectraasso);
754 if(fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
755 fOutput->Add(fCentralityCorrelation);
758 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")
760 TString gCentName[8] = {"V0A","V0C","V0M","V0AEq","V0CEq","V0MEq","CL1","ZNA"};
761 fHistCentStats = new TH2F("fHistCentStats",
762 "Centrality statistics;;Cent percentile",
763 8,-0.5,7.5,220,-5,105);
764 for(Int_t i = 1; i <= 8; i++){
765 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
766 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
768 fOutput->Add(fHistCentStats);
771 if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
773 fhistcentrality=new TH1F("fhistcentrality","referencemultiplicity",30001,-0.5,30000.5);
774 fOutput->Add(fhistcentrality);
777 fhistcentrality=new TH1F("fhistcentrality","centrality",220,-5,105);
778 fOutput->Add(fhistcentrality);
780 if(fCentralityMethod=="MC_b"){
781 fhistImpactParm=new TH1F("fhistImpactParm","Impact_Parameter",300,0,30);
782 fOutput->Add(fhistImpactParm);
785 TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
786 fHistRefmult = new TH2F("fHistRefmult",
787 "Reference multiplicity",
788 4,-0.5,3.5,10000,0,20000);
789 for(Int_t i = 1; i <= 4; i++){
790 fHistRefmult->GetXaxis()->SetBinLabel(i,gmultName[i-1].Data());
792 fOutput->Add(fHistRefmult);
794 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
795 //TPC vs EQVZERO multiplicity
796 fHistEQVZEROvsTPCmultiplicity = new TH2F("fHistEQVZEROvsTPCmultiplicity","EqVZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
797 fHistEQVZEROvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO multiplicity (a.u.)");
798 fHistEQVZEROvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
799 fOutput->Add(fHistEQVZEROvsTPCmultiplicity);
802 fHistEQVZEROAvsTPCmultiplicity = new TH2F("fHistEQVZEROAvsTPCmultiplicity","EqVZERO_A vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
803 fHistEQVZEROAvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
804 fHistEQVZEROAvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
805 fOutput->Add(fHistEQVZEROAvsTPCmultiplicity);
808 fHistEQVZEROCvsTPCmultiplicity = new TH2F("fHistEQVZEROCvsTPCmultiplicity","EqVZERO_C vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
809 fHistEQVZEROCvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
810 fHistEQVZEROCvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
811 fOutput->Add(fHistEQVZEROCvsTPCmultiplicity);
813 //EQVZERO vs VZERO multiplicity
814 fHistVZEROCvsEQVZEROCmultiplicity = new TH2F("fHistVZEROCvsEQVZEROCmultiplicity","EqVZERO_C vs VZERO_C multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
815 fHistVZEROCvsEQVZEROCmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
816 fHistVZEROCvsEQVZEROCmultiplicity->GetYaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
817 fOutput->Add(fHistVZEROCvsEQVZEROCmultiplicity);
820 fHistVZEROAvsEQVZEROAmultiplicity = new TH2F("fHistVZEROAvsEQVZEROAmultiplicity","EqVZERO_A vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
821 fHistVZEROAvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
822 fHistVZEROAvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
823 fOutput->Add(fHistVZEROAvsEQVZEROAmultiplicity);
826 //VZEROC vs VZEROA multiplicity
827 fHistVZEROCvsVZEROAmultiplicity = new TH2F("fHistVZEROCvsVZEROAmultiplicity","VZERO_C vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
828 fHistVZEROCvsVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
829 fHistVZEROCvsVZEROAmultiplicity->GetYaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
830 fOutput->Add(fHistVZEROCvsVZEROAmultiplicity);
834 //EQVZEROC vs EQVZEROA multiplicity
835 fHistEQVZEROCvsEQVZEROAmultiplicity = new TH2F("fHistEQVZEROCvsEQVZEROAmultiplicity","EqVZERO_C vs EqVZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
836 fHistEQVZEROCvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
837 fHistEQVZEROCvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
838 fOutput->Add(fHistEQVZEROCvsEQVZEROAmultiplicity);
840 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
841 fOutput->Add(fHistVZEROSignal);
845 if(fRequestEventPlane){
848 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
849 fList->Add(fHistPsiMinusPhi);
851 fEventPlanePID = new TH3F("fEventPlanePID",";centrality;eventplane;PID",20,0.0,100.0,4,-0.5,3.5,4,-0.5,3.5);
852 fList->Add(fEventPlanePID);
856 if(fCutConversions || fCutResonances)
858 fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
859 fOutput->Add(fControlConvResoncances);
862 fHistoTPCdEdx = new TH2F("fHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
863 fOutputList->Add(fHistoTPCdEdx);
864 fHistoTOFbeta = new TH2F(Form("fHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
865 fOutputList->Add(fHistoTOFbeta);
867 fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
868 fOutputList->Add(fTPCTOFPion3d);
870 fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
871 fOutputList->Add(fTPCTOFKaon3d);
873 fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
874 fOutputList->Add(fTPCTOFProton3d);
878 fPionPt = new TH1F("fPionPt","p_{T} distribution",200,0.,10.);
879 fOutputList->Add(fPionPt);
880 fPionEta= new TH1F("fPionEta","#eta distribution",360,-1.8,1.8);
881 fOutputList->Add(fPionEta);
882 fPionPhi = new TH1F("fPionPhi","#phi distribution",340,0,6.8);
883 fOutputList->Add(fPionPhi);
885 fKaonPt = new TH1F("fKaonPt","p_{T} distribution",200,0.,10.);
886 fOutputList->Add(fKaonPt);
887 fKaonEta= new TH1F("fKaonEta","#eta distribution",360,-1.8,1.8);
888 fOutputList->Add(fKaonEta);
889 fKaonPhi = new TH1F("fKaonPhi","#phi distribution",340,0,6.8);
890 fOutputList->Add(fKaonPhi);
892 fProtonPt = new TH1F("fProtonPt","p_{T} distribution",200,0.,10.);
893 fOutputList->Add(fProtonPt);
894 fProtonEta= new TH1F("fProtonEta","#eta distribution",360,-1.8,1.8);
895 fOutputList->Add(fProtonEta);
896 fProtonPhi= new TH1F("fProtonPhi","#phi distribution",340,0,6.8);
897 fOutputList->Add(fProtonPhi);
900 fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
901 fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
902 fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
903 fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);
904 fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
905 fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
906 fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
907 fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
908 fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
909 fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
910 fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
911 fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
912 fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
913 fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
914 fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
915 fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
917 for(Int_t i = 0; i < 16; i++)
919 fOutput->Add(fHistQA[i]);
922 fPriHistShare = new TH1F ("fPriHistShare","Shared clusters, primaries;#shared clusters;counts",160,0,160);
923 fOutput->Add(fPriHistShare);
925 Int_t eventplaneaxis=0;
927 if (fRequestEventPlane) eventplaneaxis=1;
929 kTrackVariablesPair=6+SetChargeAxis+eventplaneaxis;
931 if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
933 if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
935 if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8+SetChargeAxis+eventplaneaxis;
938 // two particle histograms
939 Int_t anaSteps = 1; // analysis steps
940 const char* title = "d^{2}N_{ch}/d#varphid#eta";
942 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
943 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
944 TString* axisTitlePair; // axis titles for track variables
945 axisTitlePair=new TString[kTrackVariablesPair];
947 TString defaultBinningStr;
948 defaultBinningStr = "eta: -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\n"
949 "p_t_assoc: 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 8.0,10.0\n"
950 "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
951 "p_t_eff:0.0,0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0,5.5, 6.0, 7.0, 8.0,9.0,10.0\n"
952 "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
953 "delta_phi: -1.570796, -1.483530, -1.396263, -1.308997, -1.221730, -1.134464, -1.047198, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.087266, 0.0, 0.087266, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.047198, 1.134464, 1.221730, 1.308997, 1.396263, 1.483530, 1.570796, 1.658063, 1.745329, 1.832596, 1.919862, 2.007129, 2.094395, 2.181662, 2.268928, 2.356194, 2.443461, 2.530727, 2.617994, 2.705260, 2.792527, 2.879793, 2.967060, 3.054326, 3.141593, 3.228859, 3.316126, 3.403392, 3.490659, 3.577925, 3.665191, 3.752458, 3.839724, 3.926991, 4.014257, 4.101524, 4.188790, 4.276057, 4.363323, 4.450590, 4.537856, 4.625123, 4.712389\n" // this binning starts at -pi/2 and is modulo 3
954 "delta_eta: -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,2.1, 2.2, 2.3, 2.4\n"
955 "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n"
956 "multiplicity_mixing: 0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1\n";
960 defaultBinningStr += "InvariantMass:0.200,0.300,0.398,0.399,0.4,0.401,0.402,0.403,0.404,0.405,0.406,0.407,0.408,0.409,0.41,0.411,0.412,0.413,0.414,0.415,0.416,0.417,0.418,0.419,0.42,0.421,0.422,0.423,0.424,0.425,0.426,0.427,0.428,0.429,0.43,0.431,0.432,0.433,0.434,0.435,0.436,0.437,0.438,0.439,0.44,0.441,0.442,0.443,0.444,0.445,0.446,0.447,0.448,0.449,0.45,0.451,0.452,0.453,0.454,0.455,0.456,0.457,0.458,0.459,0.46,0.461,0.462,0.463,0.464,0.465,0.466,0.467,0.468,0.469,0.47,0.471,0.472,0.473,0.474,0.475,0.476,0.477,0.478,0.479,0.48,0.481,0.482,0.483,0.484,0.485,0.486,0.487,0.488,0.489,0.49,0.491,0.492,0.493,0.494,0.495,0.496,0.497,0.498,0.499,0.5,0.501,0.502,0.503,0.504,0.505,0.506,0.507,0.508,0.509,0.51,0.511,0.512,0.513,0.514,0.515,0.516,0.517,0.518,0.519,0.52,0.521,0.522,0.523,0.524,0.525,0.526,0.527,0.528,0.529,0.53,0.531,0.532,0.533,0.534,0.535,0.536,0.537,0.538,0.539,0.54,0.541,0.542,0.543,0.544,0.545,0.546,0.547,0.548,0.549,0.55,0.551,0.552,0.553,0.554,0.555,0.556,0.557,0.558,0.559,0.56,0.561,0.562,0.563,0.564,0.565,0.566,0.567,0.568,0.569,0.57,0.571,0.572,0.573,0.574,0.575,0.576,0.577,0.578,0.579,0.58,0.581,0.582,0.583,0.584,0.585,0.586,0.587,0.588,0.589,0.59,0.591,0.592,0.593,0.594,0.595,0.596,0.597,0.598,0.599,0.600,0.700,0.800,0.900,1.000,1.065,1.066,1.067,1.068,1.069,1.07,1.071,1.072,1.073,1.074,1.075,1.076,1.077,1.078,1.079,1.08,1.081,1.082,1.083,1.084,1.085,1.086,1.087,1.088,1.089,1.09,1.091,1.092,1.093,1.094,1.095,1.096,1.097,1.098,1.099,1.1,1.101,1.102,1.103,1.104,1.105,1.106,1.107,1.108,1.109,1.11,1.111,1.112,1.113,1.114,1.115,1.116,1.117,1.118,1.119,1.12,1.121,1.122,1.123,1.124,1.125,1.126,1.127,1.128,1.129,1.13,1.131,1.132,1.133,1.134,1.135,1.136,1.137,1.138,1.139,1.14,1.141,1.142,1.143,1.144,1.145,1.146,1.147,1.148,1.149,1.15,1.151,1.152,1.153,1.154,1.155,1.156,1.157,1.158,1.159,1.16,1.161,1.162,1.163,1.164,1.165\n";
962 if(fRequestEventPlane){
963 defaultBinningStr += "eventPlane: -0.5,0.5,1.5,2.5,3.5\n"; // Event Plane Bins (Psi: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest))
965 if(fRequestEventPlanemixing){
966 defaultBinningStr += "eventPlanemixing: 0.0*TMath::DegToRad(), 30.0*TMath::DegToRad(), 60.0*TMath::DegToRad(), 90.0*TMath::DegToRad(), 120.0*TMath::DegToRad(),150.0*TMath::DegToRad(),180.1*TMath::DegToRad()\n";
969 defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5\n"; // course
972 defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
975 if(SetChargeAxis==2){
976 defaultBinningStr += "TrigCharge: -2.0,0.0,2.0\n"; // course
977 defaultBinningStr += "AssoCharge: -2.0,0.0,2.0\n"; // course
979 // =========================================================
980 // Customization (adopted from AliUEHistograms)
981 // =========================================================
983 TObjArray* lines = defaultBinningStr.Tokenize("\n");
984 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
986 TString line(lines->At(i)->GetName());
987 TString tag = line(0, line.Index(":")+1);
988 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
989 fBinningString += line + "\n";
991 AliInfo(Form("Using custom binning for %s", tag.Data()));
994 fBinningString += fCustomBinning;
996 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
998 // =========================================================
1000 // =========================================================
1002 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
1003 axisTitlePair[0] = " multiplicity";
1005 dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
1006 axisTitlePair[1] = "v_{Z} (cm)";
1008 dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
1009 axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
1011 dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
1012 axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
1014 dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
1015 axisTitlePair[4] = "#Delta#eta";
1017 dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
1018 axisTitlePair[5] = "#Delta#varphi (rad)";
1022 if(fRequestEventPlane){
1023 dBinsPair[dim_val] = GetBinning(fBinningString, "eventPlane", iBinPair[dim_val]);
1024 axisTitlePair[dim_val] = "#varphi - #Psi_{2} (a.u.)";
1028 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
1029 dBinsPair[dim_val] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val]);
1030 axisTitlePair[dim_val] = "TrigCharge";
1032 dBinsPair[dim_val+1] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+1]);
1033 axisTitlePair[dim_val+1] = "AssoCharge";
1036 if(fcontainPIDtrig && !fcontainPIDasso){
1038 dBinsPair[dim_val] = GetBinning(fBinningString, "InvariantMass", iBinPair[dim_val]);
1039 axisTitlePair[dim_val] = "InvariantMass";
1042 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
1043 axisTitlePair[dim_val] = "PIDTrig";
1045 if(SetChargeAxis==2){
1046 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
1047 axisTitlePair[dim_val+1] = "TrigCharge";
1049 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
1050 axisTitlePair[dim_val+2] = "AssoCharge";
1054 if(!fcontainPIDtrig && fcontainPIDasso){
1055 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val]);
1056 axisTitlePair[dim_val] = "PIDAsso";
1058 if(SetChargeAxis==2){
1059 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
1060 axisTitlePair[dim_val+1] = "TrigCharge";
1062 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
1063 axisTitlePair[dim_val+2] = "AssoCharge";
1067 if(fcontainPIDtrig && fcontainPIDasso){
1070 dBinsPair[dim_val] = GetBinning(fBinningString, "InvariantMass", iBinPair[dim_val]);
1071 axisTitlePair[dim_val] = "InvariantMass";
1074 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
1075 axisTitlePair[dim_val] = "PIDTrig";
1078 dBinsPair[dim_val+1] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val+1]);
1079 axisTitlePair[dim_val+1] = "PIDAsso";
1081 if(SetChargeAxis==2){
1082 dBinsPair[dim_val+2] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+2]);
1083 axisTitlePair[dim_val+2] = "TrigCharge";
1085 dBinsPair[dim_val+3] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+3]);
1086 axisTitlePair[dim_val+3] = "AssoCharge";
1091 Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
1093 Int_t nPteffbin = -1;
1094 Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
1096 Int_t multmixbin = -1;
1097 Double_t* multmix = GetBinning(fBinningString, "multiplicity_mixing", multmixbin);
1100 //Set the limits from custom binning
1101 fminPtTrig=dBinsPair[2][0];
1102 fmaxPtTrig=dBinsPair[2][iBinPair[2]];
1103 fminPtAsso=dBinsPair[3][0];
1104 fmaxPtAsso=dBinsPair[3][iBinPair[3]];
1105 fmincentmult=dBinsPair[0][0];
1106 fmaxcentmult=dBinsPair[0][iBinPair[0]];
1108 //event pool manager
1109 Int_t MaxNofEvents=1000;
1110 const Int_t NofVrtxBins=10+(1+10)*2;
1111 Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
1112 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
1113 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210};
1116 if(fRequestEventPlanemixing){
1117 // Event plane angle (Psi) bins for event mixing
1120 Double_t* psibins = GetBinning(fBinningString, "eventPlanemixing", nPsiBins);
1121 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,multmixbin,multmix,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1122 if(psibins) delete [] psibins;
1126 const Int_t nPsiBinsd=1;
1127 Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1128 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,multmixbin,multmix,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1131 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1134 AliError("Event Mixing required, but Pool Manager not initialized...");
1138 //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
1139 //fmaxPtComboeff=fmaxPtTrig;
1140 //THnSparses for calculation of efficiency
1142 if((fAnalysisType =="MCAOD") && ffillefficiency) {
1145 effbin[0]=iBinPair[0];
1146 effbin[1]=iBinPair[1];
1147 effbin[2]=nPteffbin;
1149 Int_t effsteps=5;//for each species type::primMCParticles(0),primRecoTracksMatched(1),allRecoTracksMatched(2),primRecoTracksMatchedPID(3),allRecoTracksMatchedPID(4)
1150 for(Int_t jj=0;jj<6;jj++)//PID type binning
1152 if(jj==5) effsteps=3;//for unidentified particles
1153 Histrename="fTrackHistEfficiency";Histrename+=jj;
1154 fTrackHistEfficiency[jj] = new AliTHn(Histrename.Data(), "Tracking efficiency", effsteps, 4, effbin);
1155 fTrackHistEfficiency[jj]->SetBinLimits(0, dBinsPair[0]);
1156 fTrackHistEfficiency[jj]->SetVarTitle(0, "Centrality");
1157 fTrackHistEfficiency[jj]->SetBinLimits(1, dBinsPair[1]);
1158 fTrackHistEfficiency[jj]->SetVarTitle(1, "zvtx");
1159 fTrackHistEfficiency[jj]->SetBinLimits(2, Pteff);
1160 fTrackHistEfficiency[jj]->SetVarTitle(2, "p_{T} (GeV/c)");
1161 fTrackHistEfficiency[jj]->SetBinLimits(3, EtaBin);
1162 fTrackHistEfficiency[jj]->SetVarTitle(3, "#eta");
1163 fOutput->Add(fTrackHistEfficiency[jj]);
1167 //AliThns for Correlation plots(data & MC)
1169 if(ffilltrigassoUNID)
1171 fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1172 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1173 fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
1174 fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
1176 fOutput->Add(fTHnCorrUNID);
1178 fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1179 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1180 fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
1181 fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1183 fOutput->Add(fTHnCorrUNIDmix);
1186 if(ffilltrigIDassoID)
1188 fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
1189 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1190 fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
1191 fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
1193 fOutput->Add(fTHnCorrID);
1195 fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1196 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1197 fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
1198 fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
1200 fOutput->Add(fTHnCorrIDmix);
1203 if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
1205 fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1206 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1207 fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
1208 fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
1210 fOutput->Add(fTHnCorrIDUNID);
1213 fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1214 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1215 fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
1216 fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1218 fOutput->Add(fTHnCorrIDUNIDmix);
1223 //ThnSparse for Correlation plots(truth MC)
1224 if(ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
1226 fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
1227 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1228 fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
1229 fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
1231 fOutput->Add(fCorrelatonTruthPrimary);
1234 fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
1235 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1236 fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
1237 fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
1239 fOutput->Add(fCorrelatonTruthPrimarymix);
1242 //binning for trigger no. counting
1245 if(SetChargeAxis==2) ChargeAxis=1;
1248 Int_t dims=3+ChargeAxis+eventplaneaxis;
1249 if(fcontainPIDtrig) dims=4+ChargeAxis+eventplaneaxis;
1250 fBinst= new Int_t[dims];
1251 Double_t* dBinsTrig[dims]; // bins for track variables
1252 TString* axisTitleTrig; // axis titles for track variables
1253 axisTitleTrig=new TString[dims];
1255 for(Int_t i=0; i<3;i++)
1257 fBinst[i]=iBinPair[i];
1258 dBinsTrig[i]=dBinsPair[i];
1259 axisTitleTrig[i]=axisTitlePair[i];
1261 Int_t dim_val_trig=3;
1262 if(fRequestEventPlane){
1263 fBinst[dim_val_trig]=iBinPair[6];//if fRequestEventPlane=TRUE, dim_val already becomes 7.
1264 dBinsTrig[dim_val_trig]=dBinsPair[6];
1265 axisTitleTrig[dim_val_trig]=axisTitlePair[6];
1269 if(!fcontainPIDtrig && !fcontainPIDasso && ChargeAxis==1){
1270 fBinst[dim_val_trig]=iBinPair[dim_val];
1271 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1272 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1275 if(fcontainPIDtrig && !fcontainPIDasso){
1276 fBinst[dim_val_trig]=iBinPair[dim_val];
1277 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1278 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1280 fBinst[dim_val_trig+1]=iBinPair[dim_val+1];
1281 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+1];
1282 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+1];
1286 if(!fcontainPIDtrig && fcontainPIDasso){
1288 fBinst[dim_val_trig]=iBinPair[dim_val+1];
1289 dBinsTrig[dim_val_trig]=dBinsPair[dim_val+1];
1290 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val+1];
1294 if(fcontainPIDtrig && fcontainPIDasso){
1295 fBinst[dim_val_trig]=iBinPair[dim_val];
1296 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1297 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1299 fBinst[dim_val_trig+1]=iBinPair[dim_val+2];
1300 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+2];
1301 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+2];
1305 //ThSparse for trigger counting(data & reco MC)
1306 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
1308 fTHnTrigcount = new AliTHn("fTHnTrigcount", "fTHnTrigcount", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1309 for(Int_t i=0; i<dims;i++){
1310 fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1311 fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1313 fOutput->Add(fTHnTrigcount);
1316 if(ffilltrigIDassoIDMCTRUTH) {
1317 //AliTHns for trigger counting(truth MC)
1318 fTHnTrigcountMCTruthPrim = new AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1319 for(Int_t i=0; i<dims;i++){
1320 fTHnTrigcountMCTruthPrim->SetBinLimits(i, dBinsTrig[i]);
1321 fTHnTrigcountMCTruthPrim->SetVarTitle(i, axisTitleTrig[i]);
1323 fOutput->Add(fTHnTrigcountMCTruthPrim);
1326 if(fAnalysisType=="MCAOD" || fAnalysisType=="MC"){
1327 if(ffillhistQATruth)
1329 MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
1330 fOutputList->Add(MCtruthpt);
1332 MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
1333 fOutputList->Add(MCtrutheta);
1335 MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
1336 fOutputList->Add(MCtruthphi);
1338 MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
1339 fOutputList->Add(MCtruthpionpt);
1341 MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
1342 fOutputList->Add(MCtruthpioneta);
1344 MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
1345 fOutputList->Add(MCtruthpionphi);
1347 MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
1348 fOutputList->Add(MCtruthkaonpt);
1350 MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
1351 fOutputList->Add(MCtruthkaoneta);
1353 MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
1354 fOutputList->Add(MCtruthkaonphi);
1356 MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
1357 fOutputList->Add(MCtruthprotonpt);
1359 MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
1360 fOutputList->Add(MCtruthprotoneta);
1362 MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
1363 fOutputList->Add(MCtruthprotonphi);
1365 fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
1366 fOutputList->Add(fPioncont);
1368 fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
1369 fOutputList->Add(fKaoncont);
1371 fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
1372 fOutputList->Add(fProtoncont);
1374 fUNIDcont=new TH2F("fUNIDcont","fUNIDcont",10,-0.5,9.5,100,0.,10.);
1375 fOutputList->Add(fUNIDcont);
1378 fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1379 fEventno->GetXaxis()->SetTitle("Centrality");
1380 fEventno->GetYaxis()->SetTitle("Z_Vtx");
1381 fOutput->Add(fEventno);
1382 fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1383 fEventnobaryon->GetXaxis()->SetTitle("Centrality");
1384 fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
1385 fOutput->Add(fEventnobaryon);
1386 fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1387 fEventnomeson->GetXaxis()->SetTitle("Centrality");
1388 fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
1389 fOutput->Add(fEventnomeson);
1391 fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
1392 fOutput->Add(fhistJetTrigestimate);
1394 fTwoTrackDistancePtdip = new TH3F("fTwoTrackDistancePtdip", ";#Delta#eta;#Delta#varphi;#Delta p_{T}", 36, -1.8, 1.8, 72,-TMath::Pi()/2, 3*TMath::Pi()/2, 40, 0, 10);
1395 fOutput->Add(fTwoTrackDistancePtdip);
1397 fTwoTrackDistancePtdipmix = new TH3F("fTwoTrackDistancePtdipmix", ";#Delta#eta;#Delta#varphi;#Delta p_{T}", 36, -1.8, 1.8, 72,-TMath::Pi()/2, 3*TMath::Pi()/2, 40, 0, 10);
1398 fOutput->Add(fTwoTrackDistancePtdipmix);
1400 TString Histttrname;
1401 for(Int_t jj=0;jj<2;jj++)// PID type binning
1403 Histttrname="fTwoTrackDistancePt";Histttrname+=jj;
1404 fTwoTrackDistancePt[jj] = new TH3F(Histttrname.Data(), ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
1405 fOutput->Add(fTwoTrackDistancePt[jj]);
1407 Histttrname="fTwoTrackDistancePtmix";Histttrname+=jj;
1408 fTwoTrackDistancePtmix[jj] = new TH3F(Histttrname.Data(), ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
1409 fOutput->Add(fTwoTrackDistancePtmix[jj]);
1412 //DefineEventPool();
1414 if(fapplyTrigefficiency || fapplyAssoefficiency)
1416 const Int_t nDimt = 4;// cent zvtx pt eta
1417 Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
1418 Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
1419 Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
1421 TString Histrexname;
1422 for(Int_t jj=0;jj<6;jj++)// PID type binning
1424 Histrexname="effcorection";Histrexname+=jj;
1425 effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
1426 effcorection[jj]->Sumw2();
1427 effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
1428 effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
1429 effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
1430 effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
1431 effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
1432 effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1433 effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
1434 effcorection[jj]->GetAxis(3)->SetTitle("#eta");
1435 fOutput->Add(effcorection[jj]);
1437 // TFile *fsifile = new TFile(fefffilename,"READ");
1439 if (TString(fefffilename).BeginsWith("alien:"))
1440 TGrid::Connect("alien:");
1441 TFile *fileT=TFile::Open(fefffilename);
1443 for(Int_t jj=0;jj<6;jj++)//type binning
1445 Nameg="effmap";Nameg+=jj;
1446 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
1447 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
1449 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
1460 //******************************************************************V0 plots*********************************************//
1463 //Double_t BinsV0[]={1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.5,8.0};
1464 //const Int_t nbinsV0 =sizeof(BinsV0)/sizeof(Double_t)-1;
1468 fHistRawPtCentInvK0s= new TH3F("fHistRawPtCentInvK0s", "K^{0}_{s}: mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,0.398,0.598,100,0.0,10.0,100,0.0,100.);
1469 fOutput->Add(fHistRawPtCentInvK0s);
1472 fHistRawPtCentInvLambda= new TH3F("fHistRawPtCentInvLambda", "#Lambda: mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,1.065,1.165,100,0.0,10.0,100,0.0,100.);
1473 fOutput->Add(fHistRawPtCentInvLambda);
1476 fHistRawPtCentInvAntiLambda= new TH3F("fHistRawPtCentInvAntiLambda", "#bar{#Lambda} : mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,1.065,1.165,100,0.0,10.0,100,0.0,100.);
1477 fOutput->Add(fHistRawPtCentInvAntiLambda);
1480 fHistFinalPtCentInvK0s= new TH3F("fHistFinalPtCentInvK0s", "K^{0}_{s}: mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,0.398,0.598,100,0.0,10.0,100,0.0,100.);
1481 fOutput->Add(fHistFinalPtCentInvK0s);
1484 fHistFinalPtCentInvLambda= new TH3F("fHistFinalPtCentInvLambda", "#Lambda: mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,1.065,1.165,100,0.0,10.0,100,0.0,100.);
1485 fOutput->Add(fHistFinalPtCentInvLambda);
1488 fHistFinalPtCentInvAntiLambda= new TH3F("fHistFinalPtCentInvAntiLambda", "#bar{#Lambda} : mass vs #it{p}_{T};Mass (GeV/#it{c}^2);#it{p}_{T} (GeV/#it{c});centrality",100,1.065,1.165,100,0.0,10.0,100,0.0,100.);
1489 fOutput->Add(fHistFinalPtCentInvAntiLambda);
1493 //*************************************************************EP plots***********************************************//
1494 if(fRequestEventPlane){
1495 // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
1497 fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
1498 fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
1499 fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
1501 fList->Add(fHResTPCv0A2);
1502 fList->Add(fHResTPCv0C2);
1503 fList->Add(fHResv0Cv0A2);
1506 fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
1507 fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
1508 fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
1510 fList->Add(fHResTPCv0A3);
1511 fList->Add(fHResTPCv0C3);
1512 fList->Add(fHResv0Cv0A3);
1514 // MC as in the dataEP resolution (but using MC tracks)
1515 if(fAnalysisType == "MCAOD" && fV2){
1516 fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
1517 fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
1518 fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
1519 fList->Add(fHResMA2);
1520 fList->Add(fHResMC2);
1521 fList->Add(fHResAC2);
1523 if(fAnalysisType == "MCAOD" && fV3){
1524 fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
1525 fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
1526 fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
1527 fList->Add(fHResMA3);
1528 fList->Add(fHResMC3);
1529 fList->Add(fHResAC3);
1533 // V0A and V0C event plane distributions
1535 fPhiRPTPC = new TH2F("fPhiRPTPCv2","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1536 fPhiRPTPCv3 = new TH2F("fPhiRPTPCv3","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1537 fList->Add(fPhiRPTPC);
1538 fList->Add(fPhiRPTPCv3);
1540 fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1541 fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1542 fList->Add(fPhiRPv0A);
1543 fList->Add(fPhiRPv0C);
1546 fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1547 fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1548 fList->Add(fPhiRPv0Av3);
1549 fList->Add(fPhiRPv0Cv3);
1551 fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth","#phi distribution of EP MCTRUTHheader;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1552 fList->Add(fHistEventPlaneTruth);
1556 //*****************************************************PIDQA histos*****************************************************//
1560 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1561 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1564 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1565 TH2F *fHistoNSigma=new TH2F(Form("NSigma_%d_%d",ipart,ipid),Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1566 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1567 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1568 fOutputList->Add(fHistoNSigma);
1573 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1574 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1577 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1578 TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
1579 Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1580 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1581 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1582 fOutputList->Add(fHistoNSigma);
1587 if(fPIDType==Bayes){//use bayesianPID
1588 fPIDCombined = new AliPIDCombined();
1589 fPIDCombined->SetDefaultTPCPriors();//****************************************Need to know about it
1591 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1594 TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
1595 Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1596 fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1597 fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
1598 fOutputList->Add(fHistoBayes);
1601 TH2F *fHistoBayesTPC=new TH2F(Form("probBayes_TPC_%d",ipart),
1602 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1603 fHistoBayesTPC->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1604 fHistoBayesTPC->GetYaxis()->SetTitle(Form("Bayes prob TPC %s",kParticleSpeciesName[ipart]));
1605 fOutputList->Add(fHistoBayesTPC);
1607 TH2F *fHistoBayesTOF=new TH2F(Form("probBayes_TOF_%d",ipart),
1608 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1609 fHistoBayesTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1610 fHistoBayesTOF->GetYaxis()->SetTitle(Form("Bayes prob TOF %s",kParticleSpeciesName[ipart]));
1611 fOutputList->Add(fHistoBayesTOF);
1613 TH2F *fHistoBayesTPCTOF=new TH2F(Form("probBayes_TPCTOF_%d",ipart),
1614 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1615 fHistoBayesTPCTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1616 fHistoBayesTPCTOF->GetYaxis()->SetTitle(Form("Bayes prob TPCTOF %s",kParticleSpeciesName[ipart]));
1617 fOutputList->Add(fHistoBayesTPCTOF);
1621 //nsigma separation power plot
1622 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1625 TH2F *Pi_Ka_sep=new TH2F(Form("Pi_Ka_sep_%d",ipid),
1626 Form("Pi_Ka separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1627 Pi_Ka_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1628 Pi_Ka_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1629 fOutputList->Add(Pi_Ka_sep);
1631 TH2F *Pi_Pr_sep=new TH2F(Form("Pi_Pr_sep_%d",ipid),
1632 Form("Pi_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1633 Pi_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1634 Pi_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1635 fOutputList->Add(Pi_Pr_sep);
1637 TH2F *Ka_Pr_sep=new TH2F(Form("Ka_Pr_sep_%d",ipid),
1638 Form("Ka_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1639 Ka_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1640 Ka_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1641 fOutputList->Add(Ka_Pr_sep);
1645 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1646 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1649 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1650 TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
1651 Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1652 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1653 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1654 fOutputList->Add(fHistoNSigma);
1659 if (fAnalysisType == "MCAOD"){
1660 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1661 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1664 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1665 TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
1666 Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1667 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1668 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1669 fOutputList->Add(fHistoNSigma);
1674 for(Int_t idet=0;idet<fNDetectors;idet++){
1675 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1677 if(idet==fTOF)maxy=1.1;
1678 TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
1679 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1680 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1681 fOutputList->Add(fHistoPID);
1684 //PID signal plot, before PID cut
1685 for(Int_t idet=0;idet<fNDetectors;idet++){
1687 if(idet==fTOF)maxy=1.1;
1688 TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
1689 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1690 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1691 fOutputList->Add(fHistoPID);
1694 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
1695 PostData(2, fOutputList);
1696 if(fRequestEventPlane) PostData(3, fList);
1697 AliInfo("Finished setting up the Output");
1699 TH1::AddDirectory(oldStatus);
1701 //-------------------------------------------------------------------------------
1702 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1704 if(fAnalysisType == "AOD") {
1708 }//AOD--analysis-----
1710 else if(fAnalysisType == "MCAOD" || fAnalysisType == "MC") {
1719 //-------------------------------------------------------------------------
1720 void AliTwoParticlePIDCorr::doMCAODevent()
1723 // get the event (for generator level: MCEvent())
1724 AliVEvent* event = NULL;
1725 if(fAnalysisType == "MC") {
1726 event = dynamic_cast<AliVEvent*>(MCEvent());
1729 event = dynamic_cast<AliVEvent*>(InputEvent());
1732 AliError("eventMain not available");
1736 Double_t Inv_mass=0.0;//has no meaning for pions, kaons and protons(just set 0.0) to fill the LRCParticlePID position
1744 gReactionPlane = 999.;
1746 // get centrality object and check quality(valid for p-Pb and Pb-Pb; coming soon for pp 7 TeV)
1747 Double_t cent_v0=-1.0;
1748 Double_t effcent=1.0;
1749 Double_t refmultReco =0.0;
1750 Double_t nooftrackstruth=0.0;//in case of pp this will give the multiplicity(for truth case) after the track loop(only for unidentified particles that pass kinematic cuts)
1753 if(fAnalysisType == "MC"){
1755 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1758 AliError("mcEvent not available");
1761 // count all events(physics triggered)
1762 fEventCounter->Fill(1);
1764 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(gMCEvent->GenEventHeader());
1766 TArrayF gVertexArray;
1767 header->PrimaryVertex(gVertexArray);
1768 Float_t zVtxmc =gVertexArray.At(2);
1769 //cout<<"*****************************************************************************************************hi I am here"<<endl;
1772 cent_v0=GetAcceptedEventMultiplicity((AliVEvent*)gMCEvent,kFALSE); //b value; 2nd argument has no meaning
1774 if(cent_v0<0.) return;//mainly returns impact parameter
1776 //get the event plane in case of PbPb
1777 if(fRequestEventPlane){
1778 gReactionPlane=GetEventPlane((AliVEvent*)gMCEvent,kTRUE,cent_v0);//get the truth event plane,middle argument has no meaning in this case
1779 if(gReactionPlane==999.) return;
1782 TObjArray* tracksMCtruth=new TObjArray;//for truth MC particles with PID,here unidentified means any particle other than pion, kaon or proton(Basicaly Spundefined of AliHelperPID)******WARNING::different from data and reco MC
1783 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
1785 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
1786 AliMCParticle* partMC = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
1788 AliError(Form("Could not receive particle %d", iTracks));
1791 //exclude non stable particles
1792 if(fselectprimaryTruth && !(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
1794 //consider only charged particles
1795 if(partMC->Charge() == 0) continue;
1798 //give only kinematic cuts at the generator level
1799 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1800 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
1802 if(!partMC) continue;//for safety
1804 TParticle *particle = partMC->Particle();
1805 if(!particle) continue;
1806 Int_t particletypeTruth=-999;
1808 Int_t pdgtruth = particle->GetPdgCode();
1810 //To determine multiplicity in case of PP
1812 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1813 //only physical primary(all/unidentified)
1814 if(ffillhistQATruth)
1816 MCtruthpt->Fill(partMC->Pt());
1817 MCtrutheta->Fill(partMC->Eta());
1818 MCtruthphi->Fill(partMC->Phi());
1820 if (TMath::Abs(pdgtruth)==211)
1822 particletypeTruth=SpPion;
1823 if(ffillhistQATruth)
1825 MCtruthpionpt->Fill(partMC->Pt());
1826 MCtruthpioneta->Fill(partMC->Eta());
1827 MCtruthpionphi->Fill(partMC->Phi());
1830 if (TMath::Abs(pdgtruth)==321)
1832 particletypeTruth=SpKaon;
1833 if(ffillhistQATruth)
1835 MCtruthkaonpt->Fill(partMC->Pt());
1836 MCtruthkaoneta->Fill(partMC->Eta());
1837 MCtruthkaonphi->Fill(partMC->Phi());
1840 if(TMath::Abs(pdgtruth)==2212)
1842 particletypeTruth=SpProton;
1843 if(ffillhistQATruth)
1845 MCtruthprotonpt->Fill(partMC->Pt());
1846 MCtruthprotoneta->Fill(partMC->Eta());
1847 MCtruthprotonphi->Fill(partMC->Phi());
1850 if(TMath::Abs(pdgtruth)!=211 && TMath::Abs(pdgtruth)!=321 && TMath::Abs(pdgtruth)!=2212) particletypeTruth=unidentified;//*********************WARNING:: situation is different from reco MC and data case(here we don't have SpUndefined particles,because here unidentified=SpUndefined)
1852 if(fRequestEventPlane){
1853 FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
1856 //Exclude resonances
1857 if(fExcludeResonancesInMC) {
1858 TParticle *particle = track->Particle();
1859 if(!particle) continue;
1861 Bool_t kExcludeParticle = kFALSE;
1862 Int_t gMotherIndex = particle->GetFirstMother();
1863 if(gMotherIndex != -1) {
1864 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
1866 TParticle *motherParticle = motherTrack->Particle();
1867 if(motherParticle) {
1868 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
1869 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1871 if(pdgCodeOfMother == 113 // rho0
1872 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1873 // || pdgCodeOfMother == 221 // eta
1874 // || pdgCodeOfMother == 331 // eta'
1875 // || pdgCodeOfMother == 223 // omega
1876 // || pdgCodeOfMother == 333 // phi
1877 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1878 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1879 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1880 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1881 || pdgCodeOfMother == 111 // pi0 Dalitz
1883 kExcludeParticle = kTRUE;
1889 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1890 if(kExcludeParticle) continue;
1893 //Exclude electrons with PDG
1894 if(fExcludeElectronsInMC) {
1896 TParticle *particle = track->Particle();
1899 if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
1904 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1905 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1907 Short_t chargeval=0;
1908 if(partMC->Charge()>0) chargeval=1;
1909 if(partMC->Charge()<0) chargeval=-1;
1910 if(chargeval==0) continue;
1911 const TBits *clustermap=0;
1912 const TBits *sharemap=0;
1913 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,Inv_mass,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth,clustermap,sharemap);
1914 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1915 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1916 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1920 if (fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b") fCentralityCorrelation->Fill(cent_v0, nooftrackstruth);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
1922 if (fRandomizeReactionPlane)//only for TRuth MC??
1924 Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
1925 Double_t angle = TMath::TwoPi() * centralityDigits;
1926 AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
1927 ShiftTracks(tracksMCtruth, angle);
1931 Float_t weghtval=1.0;
1934 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1936 //Fill Correlations for MC truth particles(same event)
1937 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1938 Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1941 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
1942 if (pool2 && pool2->IsReady())
1943 {//start mixing only when pool->IsReady
1944 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1945 {//proceed only when no. of trigger particles >0 in current event
1946 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
1947 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1948 { //pool event loop start
1949 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1950 if(!bgTracks6) continue;
1951 Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1953 }// pool event loop ends mixing case
1954 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1955 } //if pool->IsReady() condition ends mixing case
1957 //still in main event loop
1960 if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1964 //still in main event loop
1966 if(tracksMCtruth) delete tracksMCtruth;
1971 //"MC" type analysis is finished but still in event loop
1973 else{//if(fAnalysisType=="MCAOD")
1975 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1977 AliError("Cannot get the AOD event");
1981 // count all events(physics triggered)
1982 fEventCounter->Fill(1);
1986 //check the PIDResponse handler
1987 fPID = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
1988 if (!fPID) AliFatal("This Task needs the PID response attached to the inputHandler");
1989 //if (!fPID) return;
1990 // get mag. field required for twotrack efficiency cut
1992 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1994 //check for TClonesArray(truth track MC information)
1995 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1997 AliFatal("Error: MC particles branch not found!\n");
2001 //check for AliAODMCHeader(truth event MC information)
2002 AliAODMCHeader *header=NULL;
2003 header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
2005 printf("MC header branch not found!\n");
2009 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
2010 Float_t zVtxmc =header->GetVtxZ();
2011 if(TMath::Abs(zVtxmc)>fzvtxcut) return;
2013 // For productions with injected signals, figure out above which label to skip particles/tracks
2015 if (fInjectedSignals)
2017 AliGenEventHeader* eventHeader = 0;
2022 AliFatal("fInjectedSignals set but no MC header found");
2024 headers = header->GetNCocktailHeaders();
2025 eventHeader = header->GetCocktailHeader(0);
2029 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
2030 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
2031 AliError("First event header not found. Skipping this event.");
2032 //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
2035 skipParticlesAbove = eventHeader->NProduced();
2036 AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
2039 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
2041 //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
2042 Double_t refmultTruth = GetAcceptedEventMultiplicity((AliVEvent*)aod,kTRUE); //incase of ref multiplicity it will return the truth MC ref mullt value; need to determine the ref mult value separately for reco Mc case; in case of centrality this is final and fine
2043 refmultReco = GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE);
2044 if(refmultTruth<=0 || refmultReco<=0) return;
2045 cent_v0=refmultTruth;
2048 cent_v0=GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE); //centrality value; 2nd argument has no meaning
2051 if(cent_v0<0.) return;
2052 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
2054 //get the event plane in case of PbPb
2055 if(fRequestEventPlane){
2056 gReactionPlane=GetEventPlane((AliVEvent*)aod,kTRUE,cent_v0);//get the truth event plane
2057 if(gReactionPlane==999.) return;
2060 //TObjArray* tracksMCtruth=0;
2061 TObjArray* tracksMCtruth=new TObjArray;//for truth MC particles with PID,here unidentified means any particle other than pion, kaon or proton(Basicaly Spundefined of AliHelperPID)******WARNING::different from data and reco MC
2062 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
2066 //There is a small difference between zvtx and zVtxmc??????
2067 //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
2068 //cout<<"***********************************************zvtx="<<zvtx<<endl;
2070 //now process the truth particles(for both efficiency & correlation function)
2071 Int_t nMCTrack = fArrayMC->GetEntriesFast();
2073 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
2074 { //MC truth track loop starts
2076 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
2079 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
2083 //consider only charged particles
2084 if(partMC->Charge() == 0) continue;
2086 //consider only primary particles; neglect all secondary particles including from weak decays
2087 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
2090 //remove injected signals(primaries above <maxLabel>)
2091 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
2094 Bool_t isduplicate=kFALSE;
2095 if (fRemoveDuplicates)
2097 for (Int_t j=iMC+1; j<nMCTrack; ++j)
2098 {//2nd trutuh loop starts
2099 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
2101 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
2104 if (partMC->GetLabel() == partMC2->GetLabel())
2109 }//2nd truth loop ends
2111 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
2113 //give only kinematic cuts at the truth level
2114 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
2115 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
2117 if(!partMC) continue;//for safety
2119 //To determine multiplicity in case of PP
2121 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
2122 //only physical primary(all/unidentified)
2123 if(ffillhistQATruth)
2125 MCtruthpt->Fill(partMC->Pt());
2126 MCtrutheta->Fill(partMC->Eta());
2127 MCtruthphi->Fill(partMC->Phi());
2130 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
2131 Int_t particletypeTruth=-999;
2132 if (TMath::Abs(pdgtruth)==211)
2134 particletypeTruth=SpPion;
2135 if(ffillhistQATruth)
2138 MCtruthpionpt->Fill(partMC->Pt());
2139 MCtruthpioneta->Fill(partMC->Eta());
2140 MCtruthpionphi->Fill(partMC->Phi());
2143 if (TMath::Abs(pdgtruth)==321)
2145 particletypeTruth=SpKaon;
2146 if(ffillhistQATruth)
2148 MCtruthkaonpt->Fill(partMC->Pt());
2149 MCtruthkaoneta->Fill(partMC->Eta());
2150 MCtruthkaonphi->Fill(partMC->Phi());
2153 if(TMath::Abs(pdgtruth)==2212)
2155 particletypeTruth=SpProton;
2156 if(ffillhistQATruth)
2158 MCtruthprotonpt->Fill(partMC->Pt());
2159 MCtruthprotoneta->Fill(partMC->Eta());
2160 MCtruthprotonphi->Fill(partMC->Phi());
2163 if(TMath::Abs(pdgtruth)!=211 && TMath::Abs(pdgtruth)!=321 && TMath::Abs(pdgtruth)!=2212) particletypeTruth=unidentified;//*********************WARNING:: situation is different from reco MC and data case(here we don't have SpUndefined particles,because here unidentified=SpUndefined)
2165 if(fRequestEventPlane){
2166 FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
2169 // -- Fill THnSparse for efficiency and contamination calculation
2170 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicity i.e. multiplicity intregated out )**************Important
2172 Double_t primmctruth[4] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
2175 fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
2176 if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for primary truth mesons(3)
2177 if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for primary truth kaons+protons(4)
2178 if (TMath::Abs(pdgtruth)==211) fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
2179 if (TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
2180 if(TMath::Abs(pdgtruth)==2212) fTrackHistEfficiency[2]->Fill(primmctruth,0);//for protons
2183 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
2184 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2186 Short_t chargeval=0;
2187 if(partMC->Charge()>0) chargeval=1;
2188 if(partMC->Charge()<0) chargeval=-1;
2189 if(chargeval==0) continue;
2190 const TBits *clustermap=0;
2191 const TBits *sharemap=0;
2192 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,Inv_mass,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth,clustermap,sharemap);
2193 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
2194 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
2195 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
2197 }//MC truth track loop ends
2199 //*********************still in event loop
2201 if (fRandomizeReactionPlane)//only for TRuth MC??
2203 Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
2204 Double_t angle = TMath::TwoPi() * centralityDigits;
2205 AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
2206 ShiftTracks(tracksMCtruth, angle);
2210 Float_t weghtval=1.0;
2212 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
2215 if (fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b") fCentralityCorrelation->Fill(cent_v0, nooftrackstruth);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2217 //Fill Correlations for MC truth particles(same event)
2218 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
2219 Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
2222 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
2223 if (pool2 && pool2->IsReady())
2224 {//start mixing only when pool->IsReady
2225 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
2226 {//proceed only when no. of trigger particles >0 in current event
2227 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
2228 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
2229 { //pool event loop start
2230 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
2231 if(!bgTracks6) continue;
2232 Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
2234 }// pool event loop ends mixing case
2235 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
2236 } //if pool->IsReady() condition ends mixing case
2238 //still in main event loop
2241 if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
2245 //still in main event loop
2247 if(tracksMCtruth) delete tracksMCtruth;
2249 //now deal with reco tracks
2252 //Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
2253 Float_t bSign1=aod->GetMagneticField() ;//used for reconstructed track dca cut
2256 //detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
2257 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) cent_v0=refmultReco;
2258 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
2260 if(fRequestEventPlane){
2261 gReactionPlane = GetEventPlane((AliVEvent*)aod,kFALSE,cent_v0);//get the reconstructed event plane
2262 if(gReactionPlane==999.) return;
2266 TExMap *trackMap = new TExMap();
2268 // --- track loop for mapping matrix
2271 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2272 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2273 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2274 if (!track) continue;
2275 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
2276 if(tracktype!=1) continue;
2278 if(!track) continue;//for safety
2280 Int_t gid = track->GetID();
2281 trackMap->Add(gid,itrk);
2287 //TObjArray* tracksUNID=0;
2288 TObjArray* tracksUNID = new TObjArray;
2289 tracksUNID->SetOwner(kTRUE);
2291 //TObjArray* tracksID=0;
2292 TObjArray* tracksID = new TObjArray;
2293 tracksID->SetOwner(kTRUE);
2296 //get the selected v0 particle TObjArray
2297 TObjArray* tracksIDV0 = new TObjArray;
2298 tracksIDV0->SetOwner(kTRUE);
2300 Double_t trackscount=0.0;
2301 // loop over reconstructed tracks
2302 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2303 { // reconstructed track loop starts
2304 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2305 if (!track) continue;
2306 //get the corresponding MC track at the truth level (doing reco matching)
2307 AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
2308 if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
2310 //remove injected signals
2311 if(fInjectedSignals)
2313 AliAODMCParticle* mother = recomatched;
2315 while (!mother->IsPhysicalPrimary())
2316 {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
2317 if (mother->GetMother() < 0)
2323 mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
2329 Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
2332 if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
2333 }//remove injected signals
2335 if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
2337 Bool_t isduplicate2=kFALSE;
2338 if (fRemoveDuplicates)
2340 for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
2342 AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
2343 if (!track2) continue;
2344 AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
2345 if(!recomatched2) continue;
2347 if (track->GetLabel() == track2->GetLabel())
2354 if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
2356 fHistQA[11]->Fill(track->GetTPCNcls());
2357 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2359 if(tracktype==0) continue;
2360 if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
2362 if(!track) continue;//for safety
2363 //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
2366 if(IsTrackFromV0(aod,track)) continue;// remove auto correlation
2369 AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2371 if(fFilterBit==128){
2372 Int_t gid1 = track->GetID();
2373 //if(gid1>=0) PIDtrack = track;
2374 PIDtrack =(AliAODTrack*) aod->GetTrack(trackMap->GetValue(-1-gid1));
2375 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2380 //check for eta , phi holes
2381 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2382 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2384 Int_t particletypeMC=-9999;
2386 //tag all particles as unidentified
2387 particletypeMC=unidentified;
2389 Float_t effmatrix=1.;
2391 // -- Fill THnSparse for efficiency calculation
2392 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicity i.e. multiplicity intregated out )**************Important
2393 //NOTE:: this will be used for fillinfg THnSparse of efficiency & also to get the the track by track eff. factor on the fly(only in case of pp)
2395 //Clone & Reduce track list(TObjArray) for unidentified particles
2396 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2398 Short_t chargeval=0;
2399 if(track->Charge()>0) chargeval=1;
2400 if(track->Charge()<0) chargeval=-1;
2401 if(chargeval==0) continue;
2402 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2403 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
2404 LRCParticlePID* copy = new LRCParticlePID(particletypeMC,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2405 copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
2406 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2409 //get the pdg code of the corresponding truth particle
2410 Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
2412 Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
2413 if(ffillefficiency) {
2414 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
2415 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
2416 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
2417 if(TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions
2418 if(TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
2419 if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
2421 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
2422 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
2423 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
2424 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
2425 if( TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions
2426 if( TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
2427 if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
2431 //now start the particle identification process:)
2433 Float_t dEdx = PIDtrack->GetTPCsignal();
2434 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2436 if(HasTOFPID(PIDtrack))
2438 Double_t beta = GetBeta(PIDtrack);
2439 fHistoTOFbeta->Fill(track->Pt(), beta);
2442 //do track identification(nsigma method)
2443 particletypeMC=GetParticle(PIDtrack,fFIllPIDQAHistos);//******************************problem is here
2444 switch(TMath::Abs(pdgCode)){
2446 if(fFIllPIDQAHistos){
2447 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2448 if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2449 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
2450 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
2455 if(fFIllPIDQAHistos){
2456 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2457 if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2458 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpKaon,ipid));
2459 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
2464 if(fFIllPIDQAHistos){
2465 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2466 if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2467 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpPion,ipid));
2468 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
2475 //2-d TPCTOF map(for each Pt interval)
2476 if(HasTOFPID(PIDtrack)){
2477 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2478 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2479 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2482 //Pt, Eta , Phi distribution of the reconstructed identified particles
2485 if (particletypeMC==SpPion)
2487 fPionPt->Fill(track->Pt());
2488 fPionEta->Fill(track->Eta());
2489 fPionPhi->Fill(track->Phi());
2491 if (particletypeMC==SpKaon)
2493 fKaonPt->Fill(track->Pt());
2494 fKaonEta->Fill(track->Eta());
2495 fKaonPhi->Fill(track->Phi());
2497 if (particletypeMC==SpProton)
2499 fProtonPt->Fill(track->Pt());
2500 fProtonEta->Fill(track->Eta());
2501 fProtonPhi->Fill(track->Phi());
2505 //for misidentification fraction calculation(do it with tuneonPID)
2506 if(particletypeMC==SpPion )
2508 if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
2509 if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
2510 if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
2511 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
2513 if(particletypeMC==SpKaon )
2515 if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
2516 if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
2517 if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
2518 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
2520 if(particletypeMC==SpProton )
2522 if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
2523 if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
2524 if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
2525 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
2527 if(particletypeMC==SpUndefined )
2529 if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
2530 if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
2531 if(TMath::Abs(pdgCode)==2212) fUNIDcont->Fill(5.,track->Pt());
2532 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fUNIDcont->Fill(7.,track->Pt());
2535 if(particletypeMC==SpUndefined) continue;
2537 if(fRequestEventPlane){
2538 FillPIDEventPlane(cent_v0,particletypeMC,track->Phi(),gReactionPlane);
2541 //fill tracking efficiency
2544 if(particletypeMC==SpPion || particletypeMC==SpKaon)
2546 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) {
2547 fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
2548 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
2551 if(particletypeMC==SpKaon || particletypeMC==SpProton)
2553 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) {
2554 fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
2555 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
2558 if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) {
2559 fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
2560 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
2562 if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
2563 fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
2564 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
2566 if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
2567 fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
2568 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
2572 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2574 Short_t chargeval=0;
2575 if(track->Charge()>0) chargeval=1;
2576 if(track->Charge()<0) chargeval=-1;
2577 if(chargeval==0) continue;
2578 if (fapplyTrigefficiency || fapplyAssoefficiency)
2579 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
2580 LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2581 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2582 tracksID->Add(copy1);
2584 }// if(tracktype==1) condition structure ands
2586 }//reco track loop ends
2588 //*************************************************************still in event loop
2593 //fill the centrality/multiplicity distribution of the selected events
2594 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2596 if (fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2598 //count selected events having centrality betn 0-100%
2599 fEventCounter->Fill(13);
2601 //***************************************event no. counting
2602 Bool_t isbaryontrig=kFALSE;
2603 Bool_t ismesontrig=kFALSE;
2604 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2606 if(tracksID && tracksID->GetEntriesFast()>0)
2608 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2609 { //trigger loop starts
2610 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2612 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2613 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2614 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2615 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2617 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2618 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2623 tracksIDV0=GetV0Particles((AliVEvent*) aod, cent_v0);
2624 if(tracksIDV0->GetEntriesFast()<=0) return;
2626 //same event delte-eta, delta-phi plot
2627 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2628 {//same event calculation starts
2629 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2630 if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(gReactionPlane,tracksUNID,tracksID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
2633 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2634 {//same event calculation starts
2635 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) {
2636 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2638 else Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2640 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2643 //still in main event loop
2645 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2646 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2647 if (pool && pool->IsReady())
2648 {//start mixing only when pool->IsReady
2649 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2650 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2651 { //pool event loop start
2652 TObjArray* bgTracks = pool->GetEvent(jMix);
2653 if(!bgTracks) continue;
2654 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2655 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2656 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2658 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2660 else Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2662 }// pool event loop ends mixing case
2664 } //if pool->IsReady() condition ends mixing case
2667 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2669 }//mixing with unidentified particles
2671 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2672 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2673 if (pool1 && pool1->IsReady())
2674 {//start mixing only when pool->IsReady
2675 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2676 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2677 { //pool event loop start
2678 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2679 if(!bgTracks2) continue;
2680 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2681 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2682 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2683 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2685 }// pool event loop ends mixing case
2686 } //if pool1->IsReady() condition ends mixing case
2690 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2692 }//mixing with identified particles
2694 //no. of events analyzed
2695 fEventCounter->Fill(15);
2698 if(tracksUNID) delete tracksUNID;
2700 if(tracksID) delete tracksID;
2702 if(tracksIDV0) delete tracksIDV0;
2706 }//AOD || MCAOD condition
2708 //still in the main event loop
2711 //________________________________________________________________________
2712 void AliTwoParticlePIDCorr::doAODevent()
2716 AliVEvent *event = InputEvent();
2717 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
2718 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
2720 AliError("Cannot get the AOD event");
2723 Double_t Inv_mass=0.0;
2726 fEventCounter->Fill(1);
2728 fPID = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
2729 if (!fPID) AliFatal("This Task needs the PID response attached to the inputHandler");
2730 //if (!fPID) return;//this should be available with each event even if we don't do PID selection
2738 gReactionPlane = 999.;
2740 Double_t cent_v0= -999.;
2741 Double_t effcent=1.0;
2743 Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
2746 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
2747 Float_t bSign1=aod->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
2750 // check event cuts and fill event histograms and return the centrality or reference multiplicity value
2751 if((cent_v0 = GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE)) < 0){
2754 effcent=cent_v0;//required for efficiency correction case********Extremely Important
2755 //get the event plane in case of PbPb
2756 if(fRequestEventPlane){
2757 gReactionPlane = GetEventPlane((AliVEvent*)aod,kFALSE,cent_v0);
2758 if(gReactionPlane==999.) return;
2762 TExMap *trackMap = new TExMap();
2763 // --- track loop for mapping matrix
2766 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2767 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2768 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2769 if (!track) continue;
2770 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
2771 if(tracktype!=1) continue;
2773 if(!track) continue;//for safety
2775 Int_t gid = track->GetID();
2776 trackMap->Add(gid,itrk);
2780 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
2781 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
2783 TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
2784 tracksID->SetOwner(kTRUE); // IMPORTANT!
2786 //get the selected v0 particle TObjArray
2787 TObjArray* tracksIDV0= new TObjArray;
2788 tracksIDV0->SetOwner(kTRUE); // IMPORTANT!
2792 Bool_t fTrigPtmin1=kFALSE;
2793 Bool_t fTrigPtmin2=kFALSE;
2794 Bool_t fTrigPtJet=kFALSE;
2796 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2797 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2798 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2799 if (!track) continue;
2800 fHistQA[11]->Fill(track->GetTPCNcls());
2801 Int_t particletype=-9999;//required for PID filling
2802 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2803 if(tracktype!=1) continue;
2805 if(!track) continue;//for safety
2808 if(IsTrackFromV0(aod,track)) continue;// remove auto correlation
2810 AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2812 if(fFilterBit==128){
2813 Int_t gid1 = track->GetID();
2814 //if(gid1>=0) PIDtrack = track;
2815 PIDtrack =(AliAODTrack*) aod->GetTrack(trackMap->GetValue(-1-gid1));
2816 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2819 //check for eta , phi holes
2820 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2821 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2825 //if no applyefficiency , set the eff factor=1.0
2826 Float_t effmatrix=1.0;
2828 //tag all particles as unidentified that passed filterbit & kinematic cuts
2829 particletype=unidentified;
2831 //To count the no. of tracks having an accepted track in a certain PT(e.g. Jet Pt) range
2832 if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
2833 if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
2834 if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
2837 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) effcent=15.0;//integrated over multiplicity [i.e each track has multiplicity 15.0](so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important for efficiency related issues
2840 //to reduce memory consumption in pool
2841 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2843 //Clone & Reduce track list(TObjArray) for unidentified particles
2844 Short_t chargeval=0;
2845 if(track->Charge()>0) chargeval=1;
2846 if(track->Charge()<0) chargeval=-1;
2847 if(chargeval==0) continue;
2848 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2849 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
2850 LRCParticlePID* copy = new LRCParticlePID(particletype,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2851 copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2852 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2855 //now start the particle identificaion process:)
2857 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
2859 Float_t dEdx = PIDtrack->GetTPCsignal();
2860 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2862 //fill beta vs Pt plots only for tracks having proper TOF response(much less tracks compared to the no. that pass the filterbit & kinematic cuts)
2863 if(HasTOFPID(PIDtrack))
2865 Double_t beta = GetBeta(PIDtrack);
2866 fHistoTOFbeta->Fill(track->Pt(), beta);
2870 //track identification(using nsigma method)
2871 particletype=GetParticle(PIDtrack,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
2873 //2-d TPCTOF map(for each Pt interval)
2874 if(HasTOFPID(PIDtrack)){
2875 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2876 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2877 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2880 //ignore the Spundefined particles as they also contain pion, kaon, proton outside the nsigma cut(also if tracks don't have proper TOF PID in a certain Pt interval) & these tracks are actually counted when we do the the efficiency correction, so considering them as unidentified particles & doing the efficiency correction(i.e defining unidentified=pion+Kaon+proton+SpUndefined is right only without efficiency correction) for them will be two times wrong!!!!!
2881 if (particletype==SpUndefined) continue;//this condition creating a modulated structure in delphi projection in mixed event case(only when we are dealing with identified particles i.e. tracksID)!!!!!!!!!!!
2883 if(fRequestEventPlane){
2884 FillPIDEventPlane(cent_v0,particletype,track->Phi(),gReactionPlane);
2888 //Pt, Eta , Phi distribution of the reconstructed identified particles
2891 if (particletype==SpPion)
2893 fPionPt->Fill(track->Pt());
2894 fPionEta->Fill(track->Eta());
2895 fPionPhi->Fill(track->Phi());
2897 if (particletype==SpKaon)
2899 fKaonPt->Fill(track->Pt());
2900 fKaonEta->Fill(track->Eta());
2901 fKaonPhi->Fill(track->Phi());
2903 if (particletype==SpProton)
2905 fProtonPt->Fill(track->Pt());
2906 fProtonEta->Fill(track->Eta());
2907 fProtonPhi->Fill(track->Phi());
2911 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2913 Short_t chargeval=0;
2914 if(track->Charge()>0) chargeval=1;
2915 if(track->Charge()<0) chargeval=-1;
2916 if(chargeval==0) continue;
2917 if (fapplyTrigefficiency || fapplyAssoefficiency)
2918 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles; Bool_t mesoneffrequired=kFALSE
2919 LRCParticlePID* copy1 = new LRCParticlePID(particletype,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2920 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2921 tracksID->Add(copy1);
2923 } //track loop ends but still in event loop
2925 if(trackscount<1.0){
2926 if(tracksUNID) delete tracksUNID;
2927 if(tracksID) delete tracksID;
2931 if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
2932 if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
2933 if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
2935 Float_t weightval=1.0;
2939 //fill the centrality/multiplicity distribution of the selected events
2940 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2942 if(fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE) fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2944 //count selected events having centrality betn 0-100%
2945 fEventCounter->Fill(13);
2947 //***************************************event no. counting
2948 Bool_t isbaryontrig=kFALSE;
2949 Bool_t ismesontrig=kFALSE;
2950 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2952 if(tracksID && tracksID->GetEntriesFast()>0)
2954 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2955 { //trigger loop starts
2956 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2958 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2959 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2960 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2961 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2963 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2964 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2967 //Get the TObjArray of V0 particles
2969 tracksIDV0=GetV0Particles((AliVEvent*) aod,cent_v0);
2970 if(tracksIDV0->GetEntriesFast()<=0) return;
2973 //same event delta-eta-deltaphi plot
2974 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2975 {//same event calculation starts
2976 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2977 if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(gReactionPlane,tracksUNID,tracksID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
2981 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2982 {//same event calculation starts
2983 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID){
2984 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2985 else Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2987 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2990 //still in main event loop
2994 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2995 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2996 if (pool && pool->IsReady())
2997 {//start mixing only when pool->IsReady
2998 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2999 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
3000 { //pool event loop start
3001 TObjArray* bgTracks = pool->GetEvent(jMix);
3002 if(!bgTracks) continue;
3003 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
3004 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
3005 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
3007 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
3008 else Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
3010 }// pool event loop ends mixing case
3011 } //if pool->IsReady() condition ends mixing case
3014 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
3016 }//mixing with unidentified particles
3019 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
3020 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
3021 if (pool1 && pool1->IsReady())
3022 {//start mixing only when pool->IsReady
3023 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
3024 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
3025 { //pool event loop start
3026 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
3027 if(!bgTracks2) continue;
3028 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
3029 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
3030 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
3031 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
3033 }// pool event loop ends mixing case
3034 } //if pool1->IsReady() condition ends mixing case
3038 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
3040 }//mixing with identified particles
3043 //no. of events analyzed
3044 fEventCounter->Fill(15);
3046 //still in main event loop
3049 if(tracksUNID) delete tracksUNID;
3051 if(tracksID) delete tracksID;
3054 if(tracksIDV0) delete tracksIDV0;
3057 } // *************************event loop ends******************************************
3058 //_______________________________________________________________________
3059 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
3061 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
3063 TObjArray* tracksClone = new TObjArray;
3064 tracksClone->SetOwner(kTRUE);
3066 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
3068 LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
3069 LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->GetInvMass(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval(),particle->GetTPCPadMap(),particle->GetTPCSharedMap());
3070 copy100->SetUniqueID(particle->GetUniqueID());
3071 tracksClone->Add(copy100);
3077 //--------------------------------------------------------------------------------
3078 void AliTwoParticlePIDCorr::Fillcorrelation(Float_t ReactionPlane,TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t weight,Bool_t firstTime,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)
3081 //before calling this function check that either trackstrig & tracksasso are available
3083 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
3084 TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
3085 TArrayF eta(input->GetEntriesFast());
3086 for (Int_t i=0; i<input->GetEntriesFast(); i++)
3087 eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
3090 Int_t jmax=trackstrig->GetEntriesFast();
3092 jmax=tracksasso->GetEntriesFast();
3094 // identify K, Lambda candidates and flag those particles
3095 // a TObject bit is used for this
3096 const UInt_t kResonanceDaughterFlag = 1 << 14;
3097 if (fRejectResonanceDaughters > 0)
3099 Double_t resonanceMass = -1;
3100 Double_t massDaughter1 = -1;
3101 Double_t massDaughter2 = -1;
3102 const Double_t interval = 0.02;
3103 switch (fRejectResonanceDaughters)
3105 case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
3106 case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
3107 case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
3108 default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
3111 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
3112 trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
3113 for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
3114 tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
3116 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
3118 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3119 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
3121 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
3123 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
3124 if (trig->IsEqual(asso)) continue;
3126 if (trig->Charge() * asso->Charge() > 0) continue;
3128 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
3130 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
3132 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
3134 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
3136 trig->SetBit(kResonanceDaughterFlag);
3137 asso->SetBit(kResonanceDaughterFlag);
3139 // Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
3146 //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
3148 Float_t TriggerPtMin=fminPtTrig;
3149 Float_t TriggerPtMax=fmaxPtTrig;
3150 Int_t HighestPtTriggerIndx=-99999;
3151 TH1* triggerWeighting = 0;
3153 if(fSelectHighestPtTrig || fWeightPerEvent)//**************add this data member to the constructor
3155 if (fWeightPerEvent)
3158 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID) axis = fTHnTrigcount->GetGrid(0)->GetGrid()->GetAxis(2);
3159 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) axis = fTHnTrigcountMCTruthPrim->GetGrid(0)->GetGrid()->GetAxis(2);
3160 triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
3162 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
3163 { //trigger loop starts(highest Pt trigger selection)
3164 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3166 Float_t trigpt=trig->Pt();
3167 //to avoid overflow qnd underflow
3168 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
3169 Int_t particlepidtrig=trig->getparticle();
3170 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
3172 Float_t trigeta=trig->Eta();
3174 // some optimization
3175 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
3178 if (fOnlyOneEtaSide != 0)
3180 if (fOnlyOneEtaSide * trigeta < 0)
3183 if (fTriggerSelectCharge != 0)
3184 if (trig->Charge() * fTriggerSelectCharge < 0)
3187 if (fRejectResonanceDaughters > 0)
3188 if (trig->TestBit(kResonanceDaughterFlag)) continue;
3190 if(fSelectHighestPtTrig){
3191 if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
3193 HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
3194 TriggerPtMin=trigpt;
3198 if (fWeightPerEvent) triggerWeighting->Fill(trigpt);
3200 }//trigger loop ends(highest Pt trigger selection)
3202 }//******************(fSelectHighestPtTrig || fWeightPerEvent) condition ends
3205 //two particle correlation filling
3206 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
3207 { //trigger loop starts
3208 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3210 Float_t trigpt=trig->Pt();
3211 //to avoid overflow qnd underflow
3212 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
3213 Double_t ParticlePID_InvMass=0.0;
3214 if(fV0TrigCorr) ParticlePID_InvMass=trig->GetInvMass();
3216 Int_t particlepidtrig=trig->getparticle();
3217 ParticlePID_InvMass=(Double_t) particlepidtrig;
3218 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}//***********************************forks,lam.Alam their PID numbers have no meaning, only their Inv_mass will be stored
3221 Float_t trigeta=trig->Eta();
3223 // some optimization
3224 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
3227 if (fOnlyOneEtaSide != 0)
3229 if (fOnlyOneEtaSide * trigeta < 0)
3232 if (fTriggerSelectCharge != 0)
3233 if (trig->Charge() * fTriggerSelectCharge < 0)
3236 if (fRejectResonanceDaughters > 0)
3237 if (trig->TestBit(kResonanceDaughterFlag)) continue;
3239 if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
3240 if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
3243 Float_t trigphi=trig->Phi();
3244 Float_t trackefftrig=1.0;
3245 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
3247 // Event plane (determine psi bin)
3248 Double_t gPsiMinusPhi = 0.;
3249 Double_t gPsiMinusPhiBin = -10.;
3250 if(fRequestEventPlane){
3251 gPsiMinusPhi = TMath::Abs(trigphi - ReactionPlane);
3252 //in-plane(Note thet event plane angle has to be defined within 0 to 180 degree(do not use this if event ), otherwise the definition of in plane and out plane particles is wrong)
3253 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
3254 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
3255 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
3256 gPsiMinusPhiBin = 0.0;
3258 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
3259 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
3260 gPsiMinusPhiBin = 0.0;
3263 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
3264 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
3265 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
3266 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
3267 gPsiMinusPhiBin = 1.0;
3269 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
3270 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
3271 gPsiMinusPhiBin = 2.0;
3274 gPsiMinusPhiBin = 3.0;
3276 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
3279 //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
3282 Int_t eventplaneAxis=0;
3283 if(fRequestEventPlane) eventplaneAxis=1;
3284 if(fcontainPIDtrig && SetChargeAxis==0) dim=4+eventplaneAxis;
3285 if(!fcontainPIDtrig && SetChargeAxis==2) dim=4+eventplaneAxis;
3286 if(fcontainPIDtrig && SetChargeAxis==2) dim=5+eventplaneAxis;
3287 trigval= new Double_t[dim];
3290 trigval[2] = trigpt;
3292 if(fRequestEventPlane){
3293 trigval[3] = gPsiMinusPhiBin;
3294 if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = ParticlePID_InvMass;
3295 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
3296 if(fcontainPIDtrig && SetChargeAxis==2) {
3297 trigval[4] = ParticlePID_InvMass;
3298 trigval[5] = trig->Charge();
3302 if(!fRequestEventPlane){
3303 if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = ParticlePID_InvMass;
3304 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
3305 if(fcontainPIDtrig && SetChargeAxis==2) {
3306 trigval[3] = ParticlePID_InvMass;
3307 trigval[4] = trig->Charge();
3313 if (fWeightPerEvent)
3315 // leads effectively to a filling of one entry per filled trigger particle pT bin
3316 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(trigval[2]);
3317 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3318 trackefftrig *= triggerWeighting->GetBinContent(weightBin);
3322 //trigger particle counting for both same and mixed event case;;;;;step=0->same event case;;;;;step=1->mixed event case;;;;;;
3323 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
3324 if(fillup=="trigassoUNID" ) {
3325 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3326 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3329 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
3330 if(fillup=="trigassoUNID" )
3332 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3333 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3336 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
3337 if(fillup=="trigUNIDassoID")
3339 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3340 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3343 //ensure that trigIDassoID , trigassoUNID, trigIDassoUNID & trigUNIDassoID case FillCorrelation called only once in the event loop for same event correlation function calculation, otherwise there will be multiple counting of pion, kaon,proton,unidentified
3344 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
3345 if(fillup=="trigIDassoID")
3347 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3348 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3351 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
3352 if(fillup=="trigIDassoUNID" )
3354 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3355 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3358 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
3359 if(fillup=="trigIDassoID")
3361 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3362 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3366 if(fillup=="trigIDassoIDMCTRUTH") { //In truth MC case "Unidentified" means any particle other than pion,kaon or proton and no efficiency correction(default value 1.0)************************be careful!!!!
3367 if(mixcase==kFALSE) fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig);
3368 if(mixcase==kTRUE && firstTime) fTHnTrigcountMCTruthPrim->Fill(trigval,1,1.0/trackefftrig);
3371 //asso loop starts within trigger loop
3372 for(Int_t j=0;j<jmax;j++)
3374 LRCParticlePID *asso=0;
3376 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
3378 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
3382 //to avoid overflow and underflow
3383 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
3385 //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
3387 if(!tracksasso && j==i) continue;
3389 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event,i.e. both Trig and asso TObjArray belongs to the same Pt range but say Trig is Unidentified but asso is identified then the serial no. wise particles are not same and and j==i doesn't aplly)
3390 // if (tracksasso && trig->IsEqual(asso)) continue;
3392 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
3395 if (asso->Pt() >= trig->Pt()) continue;
3397 Int_t particlepidasso=asso->getparticle();
3398 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
3401 if (fAssociatedSelectCharge != 0)
3402 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
3404 if (fSelectCharge > 0)
3407 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
3411 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
3417 if (trigeta < 0 && asso->Eta() < trigeta)
3419 if (trigeta > 0 && asso->Eta() > trigeta)
3423 if (fRejectResonanceDaughters > 0)
3424 if (asso->TestBit(kResonanceDaughterFlag))
3426 // Printf("Skipped j=%d", j);
3431 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
3433 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3437 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3439 fControlConvResoncances->Fill(0.0, mass);
3441 if (mass < 0.04*0.04)
3447 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3449 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
3451 const Float_t kK0smass = 0.4976;
3453 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
3455 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
3457 fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
3459 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
3465 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3467 Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
3468 Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3470 const Float_t kLambdaMass = 1.115;
3472 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
3474 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
3476 fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
3478 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3481 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
3483 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3485 fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
3487 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3492 if (twoTrackEfficiencyCut)
3494 // the variables & cuthave been developed by the HBT group
3495 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
3496 Float_t phi1 = trig->Phi();
3497 Float_t pt1 = trig->Pt();
3498 Float_t charge1 = trig->Charge();
3499 Float_t phi2 = asso->Phi();
3500 Float_t pt2 = asso->Pt();
3501 Float_t charge2 = asso->Charge();
3503 Float_t deta= trigeta - eta[j];
3506 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
3509 // check first boundaries to see if is worth to loop and find the minimum
3510 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMinRadius, bSign);
3511 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMaxRadius, bSign);
3513 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
3515 Float_t dphistarminabs = 1e5;
3516 Float_t dphistarmin = 1e5;
3518 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
3520 for (Double_t rad=fTwoTrackCutMinRadius; rad<=fTwoTrackCutMaxRadius; rad+=0.01)
3522 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
3524 Float_t dphistarabs = TMath::Abs(dphistar);
3526 if (dphistarabs < dphistarminabs)
3528 dphistarmin = dphistar;
3529 dphistarminabs = dphistarabs;
3532 if(mixcase==kFALSE) fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3533 if(mixcase==kTRUE) fTwoTrackDistancePtmix[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3535 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
3537 // Printf("Removed track pair %d %d with %f %f %f %f %f %f %f %f %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
3540 if(mixcase==kFALSE) fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3541 if(mixcase==kTRUE) fTwoTrackDistancePtmix[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3547 //pair sharedfraction cut(only between the trig and asso track)
3548 if(fillup!="trigIDassoIDMCTRUTH")//******************************************NOT for TRUTH MC particles
3550 if(fSharedfraction_Pair_cut>=0){
3551 Bool_t passsharedfractionpaircut=CalculateSharedFraction(trig->GetTPCPadMap(),asso->GetTPCPadMap(),trig->GetTPCSharedMap(),asso->GetTPCSharedMap());
3552 if(!passsharedfractionpaircut) continue;
3555 Float_t weightperevent=weight;
3556 Float_t trackeffasso=1.0;
3557 if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
3558 //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
3559 Float_t deleta=trigeta-eta[j];
3560 Float_t delphi=PhiRange(trigphi-asso->Phi());
3562 Float_t delpt=trigpt-asso->Pt();
3563 //fill it with/without two track efficiency cut
3564 if(mixcase==kFALSE) fTwoTrackDistancePtdip->Fill(deleta, delphi, TMath::Abs(delpt));//for same event
3565 if(mixcase==kTRUE) fTwoTrackDistancePtdipmix->Fill(deleta, delphi, TMath::Abs(delpt));//for mixed event
3567 //here get the two particle efficiency correction factor
3568 Float_t effweight=trackefftrig*trackeffasso*weightperevent;
3569 // if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
3571 Int_t dimused=kTrackVariablesPair+eventplaneAxis;
3572 vars= new Double_t[dimused];
3581 if(fRequestEventPlane)
3583 vars[6]=gPsiMinusPhiBin;
3587 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
3588 vars[dimension]=trig->Charge();
3589 vars[dimension+1]=asso->Charge();
3591 if(fcontainPIDtrig && !fcontainPIDasso){
3592 vars[dimension]=ParticlePID_InvMass;
3593 if(SetChargeAxis==2){
3594 vars[dimension+1]=trig->Charge();
3595 vars[dimension+2]=asso->Charge();
3598 if(!fcontainPIDtrig && fcontainPIDasso){
3599 vars[dimension]=particlepidasso;
3600 if(SetChargeAxis==2){
3601 vars[dimension+1]=trig->Charge();
3602 vars[dimension+2]=asso->Charge();
3605 if(fcontainPIDtrig && fcontainPIDasso){
3606 vars[dimension]=ParticlePID_InvMass;
3607 vars[dimension+1]=particlepidasso;
3608 if(SetChargeAxis==2){
3609 vars[dimension+2]=trig->Charge();
3610 vars[dimension+3]=asso->Charge();
3614 if (fWeightPerEvent)
3616 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
3617 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3618 effweight *= triggerWeighting->GetBinContent(weightBin);
3622 //Fill Correlation ThnSparses
3623 if(fillup=="trigassoUNID")
3625 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,0,1.0/effweight);
3626 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
3628 if(fillup=="trigIDassoID")
3630 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,0,1.0/effweight);
3631 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
3633 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
3634 {//in this case effweight should be 1 always
3635 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
3636 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
3638 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
3640 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
3641 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
3647 }//trigger loop ends
3649 if (triggerWeighting)
3651 delete triggerWeighting;
3652 triggerWeighting = 0;
3656 //------------------------------------------------------------------------------------------------
3657 Bool_t AliTwoParticlePIDCorr:: CalculateSharedFraction(const TBits *triggerPadMap,const TBits *assocPadMap,const TBits *triggerShareMap,const TBits *assocShareMap)
3658 {//source code-AliFemtoShareQualityPairCut.cxx
3660 Double_t nofsharedhits=0;
3662 for(UInt_t imap=0;imap< (triggerPadMap->GetNbits() );imap++)
3664 //if they are in same pad
3665 //cout<<triggerPadMap->TestBitNumber(imap)<<" "<< assocPadMap->TestBitNumber(imap)<<endl;
3666 if (triggerPadMap->TestBitNumber(imap) &&
3667 assocPadMap->TestBitNumber(imap))
3670 //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
3671 if (triggerShareMap->TestBitNumber(imap) &&
3672 assocShareMap->TestBitNumber(imap))
3674 //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
3691 //cout<< (triggerPadMap->TestBitNumber(imap) || assocPadMap->TestBitNumber(imap))<<endl;
3692 else if (triggerPadMap->TestBitNumber(imap) ||
3693 assocPadMap->TestBitNumber(imap)) {
3694 // One track has a hit, the other does not
3697 //cout<<"No hits :"<<nofhits<<endl;
3705 Double_t SharedFraction=0.0;
3706 if(nofhits>0) SharedFraction=(nofsharedhits/nofhits);
3708 //cout<<"Fraction shared hits :"<<SharedFraction<<endl;
3710 if(SharedFraction>fSharedfraction_Pair_cut) return kFALSE;
3716 //________________________________________________________________________________________________
3717 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
3719 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
3721 Float_t effvalue=1.;
3723 if(parpid==unidentified)
3725 effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
3726 effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
3727 effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
3728 effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
3729 effvalue=effcorection[5]->GetBinContent(effVars);
3731 if(parpid==SpPion || parpid==SpKaon)
3733 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3735 effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
3736 effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
3737 effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
3738 effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
3739 effvalue=effcorection[3]->GetBinContent(effVars);
3744 effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
3745 effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
3746 effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
3747 effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
3748 effvalue=effcorection[0]->GetBinContent(effVars);
3753 effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
3754 effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
3755 effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
3756 effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
3757 effvalue=effcorection[1]->GetBinContent(effVars);
3762 if(parpid==SpProton)
3764 effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
3765 effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
3766 effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
3767 effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
3768 effvalue=effcorection[2]->GetBinContent(effVars);
3771 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3773 if(parpid==SpProton || parpid==SpKaon)
3775 effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
3776 effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
3777 effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
3778 effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
3779 effvalue=effcorection[4]->GetBinContent(effVars);
3782 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
3783 if(effvalue==0.) effvalue=1.;
3788 //---------------------------------------------------------------------------------
3790 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield, Bool_t fill)
3793 if(!track) return 0;
3794 Bool_t trackOK = track->TestFilterBit(fFilterBit);
3795 if(!trackOK) return 0;
3796 if (fTrackStatus != 0 && !CheckTrack(track)) return 0;
3797 //select only primary traks(for data & reco MC tracks)
3798 if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
3799 if(track->Charge()==0) return 0;
3800 if (fill) fHistQA[12]->Fill(track->GetTPCNcls());
3803 dz = track->ZAtDCA();
3804 if (fill) fHistQA[6]->Fill(dxy);
3805 if (fill) fHistQA[7]->Fill(dz);
3806 Float_t chi2ndf = track->Chi2perNDF();
3807 if (fill) fHistQA[13]->Fill(chi2ndf);
3808 // Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
3809 Float_t nCrossedRowsTPC = track->GetTPCNCrossedRows();
3810 if (fill) fHistQA[14]->Fill(nCrossedRowsTPC);
3811 //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
3812 if (track->GetTPCNclsF()>0) {
3813 Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
3814 if (fill) fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
3817 Float_t pt=track->Pt();
3818 if(pt< fminPt || pt> fmaxPt) return 0;
3819 if(TMath::Abs(track->Eta())> fmaxeta) return 0;
3820 if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
3821 //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required???
3823 if (fdcacut && fDCAXYCut)
3830 AliAODTrack* clone =(AliAODTrack*) track->Clone();
3831 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
3836 // Printf("%f", ((AliAODTrack*)part)->DCA());
3837 // Printf("%f", pos[0]);
3838 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
3842 if (fSharedClusterCut >= 0)
3844 Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
3845 if (frac > fSharedClusterCut)
3849 // Rejects tracks with shared clusters after filling a control histogram
3850 // This overload is used for primaries
3852 // Get the shared maps
3853 const TBits sharedMap = track->GetTPCSharedMap();
3854 // Fill a control histogram
3855 fPriHistShare->Fill(sharedMap.CountBits());
3857 // Reject shared clusters
3858 if (fSharedTPCmapCut >= 0)
3860 if((sharedMap.CountBits()) >= 1) return 0;// Bad track, has too many shared clusters!
3863 if (fill) fHistQA[8]->Fill(pt);
3864 if (fill) fHistQA[9]->Fill(track->Eta());
3865 if (fill) fHistQA[10]->Fill(track->Phi());
3868 //________________________________________________________________________________
3869 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
3871 //This function is called within the func GetParticle() for accepted tracks only i.e.after call of Classifytrack() & for those tracks which have proper TPC PID response . combined nsigma(circular) cut only for particles having pt upto 4.0 Gev/c and beyond that use the asymmetric nsigma cut around pion's mean position in TPC ( while filling the TObjArray for trig & asso )
3872 Float_t pt=track->Pt();
3874 //plot the separation power
3876 Double_t bethe[AliPID::kSPECIES]={0.};
3877 Double_t sigma_TPC[AliPID::kSPECIES]={0.};
3879 Double_t Pi_Ka_sep[NSigmaPIDType+1]={0.};
3880 Double_t Pi_Pr_sep[NSigmaPIDType+1]={0.};
3881 Double_t Ka_Pr_sep[NSigmaPIDType+1]={0.};
3884 Double_t ptpc = track->GetTPCmomentum();
3885 Int_t dEdxN = track->GetTPCsignalN();
3886 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
3887 bethe[ipart] = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
3888 //Double_t diff = dEdx - bethe;
3889 sigma_TPC[ipart] = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
3890 //nSigma[ipart] = diff / sigma;
3892 Pi_Ka_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kKaon])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kKaon])/2.0);
3893 Pi_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kProton])/2.0);
3894 Ka_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kKaon]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kKaon]+sigma_TPC[AliPID::kProton])/2.0);
3897 Double_t sigma_TOF[AliPID::kSPECIES]={0.};
3899 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3901 Double_t timei[AliPID::kSPECIES];
3902 track->GetIntegratedTimes(timei);
3903 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { sigma_TOF[ipart]= fPID->GetTOFResponse().GetExpectedSigma(track->P(), timei[ipart], AliPID::ParticleMass(ipart));}
3904 Pi_Ka_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kKaon])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kKaon])/2.0);
3905 Pi_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kProton])/2.0);
3906 Ka_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kKaon]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kKaon]+sigma_TOF[AliPID::kProton])/2.0);
3908 Pi_Ka_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Ka_sep[NSigmaTPC]*Pi_Ka_sep[NSigmaTPC]+Pi_Ka_sep[NSigmaTOF]*Pi_Ka_sep[NSigmaTOF]);
3909 Pi_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Pr_sep[NSigmaTPC]*Pi_Pr_sep[NSigmaTPC]+Pi_Pr_sep[NSigmaTOF]*Pi_Pr_sep[NSigmaTOF]);
3910 Ka_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Ka_Pr_sep[NSigmaTPC]*Ka_Pr_sep[NSigmaTPC]+Ka_Pr_sep[NSigmaTOF]*Ka_Pr_sep[NSigmaTOF]);
3914 //fill separation power histograms
3915 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3917 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3918 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3919 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3920 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3921 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3922 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3924 if(HasTOFPID(track) && pt>fPtTOFPIDmin && ipid!=0){
3925 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3926 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3927 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3928 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3929 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3930 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3935 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
3936 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
3937 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
3938 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
3940 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
3941 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
3943 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3946 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
3947 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
3948 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
3950 nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
3951 nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
3952 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
3958 // if TOF is missing and below fPtTOFPID only the TPC information is used
3959 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
3960 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
3961 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
3965 //set data member fnsigmas
3966 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
3967 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
3968 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
3970 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
3971 fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
3972 fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
3973 fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
3975 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TPCTOF nsigma values will be TMath::Abs(TPC only nsigma)
3976 fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
3977 fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
3978 fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
3981 //Fill NSigma SeparationPlot
3982 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3983 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3984 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3985 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
3986 h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
3992 //----------------------------------------------------------------------------
3993 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
3995 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3996 if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && (!HasTOFPID(track)) )return SpUndefined;//so any track having Pt>0.6 withot having proper TOF response will be defined as SpUndefined
3997 //get the identity of the particle with the minimum Nsigma
3998 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
4001 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
4002 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
4003 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
4006 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
4007 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
4008 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
4010 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
4011 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4012 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4013 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4015 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
4016 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4017 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4018 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4023 if(fdiffPIDcutvalues){
4024 if(track->Pt()<=4) fNSigmaPID=fPIDCutval1;
4025 if(track->Pt()>4 && track->Pt()<=6) fNSigmaPID=fPIDCutval2;
4026 if(track->Pt()>6 && track->Pt()<=8) fNSigmaPID=fPIDCutval3;
4027 if(track->Pt()>8) fNSigmaPID=fPIDCutval4;
4030 // guess the particle based on the smaller nsigma (within fNSigmaPID)
4031 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
4033 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
4034 if((fHighPtKaonNSigmaPID>0) && (track->Pt()>fHighPtKaonSigma) && (nsigmaKaon > fHighPtKaonNSigmaPID)) return SpUndefined;//different nsigma cut for kaons above a particular Pt range(within the TPC-TOF PID range)
4036 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4037 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4038 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
4039 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
4044 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
4046 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4047 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4048 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
4049 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
4054 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
4056 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4057 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4058 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
4059 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
4065 // else, return undefined
4071 //------------------------------------------------------------------------------------------
4072 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
4073 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
4075 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
4077 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
4079 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
4082 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
4084 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
4087 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
4088 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
4089 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
4092 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
4093 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
4094 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
4096 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
4097 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4098 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4099 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4101 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
4102 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4103 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4104 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4108 // Actually the tracks in the overlapping region(in TPC-TOF nSigma plane) will be ignored
4110 if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
4111 if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
4112 if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
4117 //fill NSigma distr for double counting
4118 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4119 if(fHasDoubleCounting[ipart]){//this may be kTRUE only for particles having Pt<=4.0 GeV/C, so this histo contains all the particles having Pt<=4.0 GeV/C in the nsigma overlapping region in TPC/TPC-TOF plane
4120 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4121 if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)))continue;//not filling TOF and combined if no TOF PID
4122 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
4123 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
4130 return fHasDoubleCounting;
4133 //////////////////////////////////////////////////////////////////////////////////////////////////
4135 Bool_t* AliTwoParticlePIDCorr::GetAllCompatibleIdentitiesNSigma(AliAODTrack * trk,Bool_t FIllQAHistos){
4136 //mainly intended to check the probability of the PID of the tracks which are in the overlapping nSigma regions and near about the middle position from the mean position of two ID particle
4137 Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
4138 IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
4142 //////////////////////////////////////////////////////////////////////////////////////////////////
4144 UInt_t AliTwoParticlePIDCorr::CalcPIDCombined(AliAODTrack *track, Int_t detMask, Double_t* prob) const{
4146 // Bayesian PID calculation
4148 for(Int_t i=0;i<AliPID::kSPECIES;i++)
4152 fPIDCombined->SetDetectorMask(detMask);
4154 return fPIDCombined->ComputeProbabilities((AliAODTrack*)track, fPID, prob);
4157 //////////////////////////////////////////////////////////////////////////////////////////////////
4159 Int_t AliTwoParticlePIDCorr::GetIDBayes(AliAODTrack * trk, Bool_t FIllQAHistos){
4161 Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
4164 //Filling of Probability histos
4165 Double_t probTPC[AliPID::kSPECIES]={0.};
4166 Double_t probTOF[AliPID::kSPECIES]={0.};
4167 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
4169 UInt_t detUsedTPC = 0;
4170 UInt_t detUsedTOF = 0;
4171 UInt_t detUsedTPCTOF = 0;
4173 //get the TPC probability
4174 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
4175 detUsedTPC = fPIDCombined->ComputeProbabilities(trk, fPID, probTPC);
4176 if(detUsedTPC == AliPIDResponse::kDetTPC)
4178 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4180 TH2F *h=GetHistogram2D(Form("probBayes_TPC_%d",ipart));
4181 if(ipart==0) h->Fill(trk->Pt(),probTPC[AliPID::kPion]);
4182 if(ipart==1) h->Fill(trk->Pt(),probTPC[AliPID::kKaon]);
4183 if(ipart==2) h->Fill(trk->Pt(),probTPC[AliPID::kProton]);
4187 //get the TOF probability
4188 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
4189 detUsedTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTOF);
4190 if(detUsedTOF == AliPIDResponse::kDetTOF)
4192 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4193 TH2F *h=GetHistogram2D(Form("probBayes_TOF_%d",ipart));
4194 if(ipart==0) h->Fill(trk->Pt(),probTOF[AliPID::kPion]);
4195 if(ipart==1) h->Fill(trk->Pt(),probTOF[AliPID::kKaon]);
4196 if(ipart==2) h->Fill(trk->Pt(),probTOF[AliPID::kProton]);
4200 //get the TPC-TOF probability
4201 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
4202 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTPCTOF);
4203 if(detUsedTPCTOF == (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))
4205 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4206 TH2F *h=GetHistogram2D(Form("probBayes_TPCTOF_%d",ipart));
4207 if(ipart==0) h->Fill(trk->Pt(),probTPCTOF[AliPID::kPion]);
4208 if(ipart==1) h->Fill(trk->Pt(),probTPCTOF[AliPID::kKaon]);
4209 if(ipart==2) h->Fill(trk->Pt(),probTPCTOF[AliPID::kProton]);
4214 Double_t probBayes[AliPID::kSPECIES];
4217 if(HasTOFPID(trk) && trk->Pt()>fPtTOFPIDmin){//use TOF information
4218 detUsed = CalcPIDCombined(trk, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
4219 if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return SpUndefined;//check that TPC and TOF are used
4221 detUsed = CalcPIDCombined(trk,AliPIDResponse::kDetTPC, probBayes);
4222 if(detUsed != AliPIDResponse::kDetTPC)return SpUndefined;//check that TPC is used
4225 //the probability has to be normalized to one, we check it
4227 for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
4228 if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
4229 AliFatal("Bayesian probability not normalized to one");
4232 if(fdiffPIDcutvalues){
4233 if(trk->Pt()<=4) fBayesCut=fPIDCutval1;
4234 if(trk->Pt()>4 && trk->Pt()<=6) fBayesCut=fPIDCutval2;
4235 if(trk->Pt()>6 && trk->Pt()<=8) fBayesCut=fPIDCutval3;
4236 if(trk->Pt()>8) fBayesCut=fPIDCutval4;
4240 //probabilities are normalized to one, if the cut is above .5 there is no problem
4241 if(probBayes[AliPID::kPion]>fBayesCut && IDs[SpPion]==1){
4242 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpPion));
4243 h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
4246 else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[SpKaon]==1){
4247 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpKaon));
4248 h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
4251 else if(probBayes[AliPID::kProton]>fBayesCut && IDs[SpProton]==1){
4252 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpProton));
4253 h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
4262 //////////////////////////////////////////////////////////////////////////////////////////////////
4263 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
4264 //return the specie according to the minimum nsigma value
4265 //no double counting, this has to be evaluated using CheckDoubleCounting()
4267 Int_t ID=SpUndefined;
4269 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
4273 if(fPIDType==Bayes){//use bayesianPID
4276 AliFatal("PIDCombined object has to be set in the steering macro");
4279 ID = GetIDBayes(trk,FIllQAHistos);
4281 }else{ //use nsigma PID
4283 ID=FindMinNSigma(trk,FIllQAHistos);
4284 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
4286 HasDC=GetDoubleCounting(trk,FIllQAHistos);
4287 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4288 if(HasDC[ipart]==kTRUE) ID = SpUndefined;
4292 //Fill PID signal plot
4293 if(ID != SpUndefined){
4294 for(Int_t idet=0;idet<fNDetectors;idet++){
4295 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
4296 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
4297 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
4298 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
4301 //Fill PID signal plot without cuts
4302 for(Int_t idet=0;idet<fNDetectors;idet++){
4303 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
4304 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
4305 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
4306 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
4312 //-------------------------------------------------------------------------------------
4314 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
4317 AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
4318 if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
4319 //ULong_t status=track->GetStatus();
4320 //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
4321 //if (track->GetTPCsignal() <= 0.) return kFALSE;
4322 // if(track->GetTPCsignalN() < 60) return kFALSE;//tracks with TPCsignalN< 60 have questionable dEdx,cutting on TPCsignalN > 70 or > 60 shouldn't make too much difference in statistics,also it is IMO safe to use TPC also for MIPs.
4326 //___________________________________________________________
4329 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
4331 // check TOF matched track
4332 //ULong_t status=track->GetStatus();
4333 //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
4334 AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
4335 if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
4336 if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
4337 //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
4338 //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
4339 // if (probMis > 0.01) return kFALSE;
4340 if(fRemoveTracksT0Fill)
4342 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
4343 if (startTimeMask < 0)return kFALSE;
4348 //________________________________________________________________________
4349 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
4351 //it is called only when TOF PID is available
4352 //TOF beta calculation
4353 Double_t tofTime=track->GetTOFsignal();
4355 Double_t c=TMath::C()*1.E-9;// m/ns
4356 Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
4357 Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
4358 tofTime -= startTime; // subtract startTime to the signal
4359 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
4365 Double_t p = track->P();
4366 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
4368 track->GetIntegratedTimes(timei);
4369 return timei[0]/time;
4372 //------------------------------------------------------------------------------------------------------
4374 Float_t AliTwoParticlePIDCorr::GetInvMassSquared(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
4376 // calculate inv mass squared
4377 // same can be achieved, but with more computing time with
4378 /*TLorentzVector photon, p1, p2;
4379 p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
4380 p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
4384 Float_t tantheta1 = 1e10;
4386 if (eta1 < -1e-10 || eta1 > 1e-10)
4387 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
4389 Float_t tantheta2 = 1e10;
4390 if (eta2 < -1e-10 || eta2 > 1e-10)
4391 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
4393 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4394 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4396 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );
4400 //---------------------------------------------------------------------------------
4402 Float_t AliTwoParticlePIDCorr::GetInvMassSquaredCheap(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
4404 // calculate inv mass squared approximately
4406 Float_t tantheta1 = 1e10;
4408 if (eta1 < -1e-10 || eta1 > 1e-10)
4410 Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
4411 tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4414 Float_t tantheta2 = 1e10;
4415 if (eta2 < -1e-10 || eta2 > 1e-10)
4417 Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
4418 tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4421 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4422 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4425 Float_t deltaPhi = TMath::Abs(phi1 - phi2);
4426 while (deltaPhi > TMath::TwoPi())
4427 deltaPhi -= TMath::TwoPi();
4428 if (deltaPhi > TMath::Pi())
4429 deltaPhi = TMath::TwoPi() - deltaPhi;
4431 Float_t cosDeltaPhi = 0;
4432 if (deltaPhi < TMath::Pi()/3)
4433 cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
4434 else if (deltaPhi < 2*TMath::Pi()/3)
4435 cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
4437 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
4439 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
4441 // Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
4445 //--------------------------------------------------------------------------------
4446 Float_t AliTwoParticlePIDCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign)
4449 // calculates dphistar
4452 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
4454 static const Double_t kPi = TMath::Pi();
4457 // if (dphistar > 2 * kPi)
4458 // dphistar -= 2 * kPi;
4459 // if (dphistar < -2 * kPi)
4460 // dphistar += 2 * kPi;
4463 dphistar = kPi * 2 - dphistar;
4464 if (dphistar < -kPi)
4465 dphistar = -kPi * 2 - dphistar;
4466 if (dphistar > kPi) // might look funny but is needed
4467 dphistar = kPi * 2 - dphistar;
4472 //------------------------------------------------------------------------
4473 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
4475 // This method is a copy from AliUEHist::GetBinning
4476 // takes the binning from <configuration> identified by <tag>
4477 // configuration syntax example:
4478 // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
4481 // returns bin edges which have to be deleted by the caller
4483 TString config(configuration);
4484 TObjArray* lines = config.Tokenize("\n");
4485 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
4487 TString line(lines->At(i)->GetName());
4488 if (line.BeginsWith(TString(tag) + ":"))
4490 line.Remove(0, strlen(tag) + 1);
4491 line.ReplaceAll(" ", "");
4492 TObjArray* binning = line.Tokenize(",");
4493 Double_t* bins = new Double_t[binning->GetEntriesFast()];
4494 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
4495 bins[j] = TString(binning->At(j)->GetName()).Atof();
4497 nBins = binning->GetEntriesFast() - 1;
4506 AliFatal(Form("Tag %s not found in %s", tag, configuration));
4510 //____________________________________________________________________
4512 Bool_t AliTwoParticlePIDCorr::CheckTrack(AliAODTrack * part)
4514 // check if the track status flags are set
4516 UInt_t status=((AliAODTrack*)part)->GetStatus();
4517 if ((status & fTrackStatus) == fTrackStatus)
4521 //________________________________________________________________________
4523 Bool_t AliTwoParticlePIDCorr::AcceptEventCentralityWeight(Double_t centrality)
4525 // rejects "randomly" events such that the centrality gets flat
4526 // uses fCentralityWeights histogram
4528 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
4530 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
4531 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
4533 Bool_t result = kFALSE;
4534 if (centralityDigits < weight)
4537 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
4542 //____________________________________________________________________
4543 void AliTwoParticlePIDCorr::ShiftTracks(TObjArray* tracks, Double_t angle)
4545 // shifts the phi angle of all tracks by angle
4546 // 0 <= angle <= 2pi
4548 for (Int_t i=0; i<tracks->GetEntriesFast(); ++i)
4550 LRCParticlePID *part=(LRCParticlePID*)(tracks->UncheckedAt(i));
4552 Double_t newAngle = part->Phi() + angle;
4553 if (newAngle >= TMath::TwoPi())
4554 newAngle -= TMath::TwoPi();
4556 part->SetPhi(newAngle);
4561 //________________________________________________________________________
4562 void AliTwoParticlePIDCorr::SetVZEROCalibrationFile(const char* filename,const char* lhcPeriod) {
4563 //Function to setup the VZERO gain equalization
4564 //============Get the equilization map============//
4565 TFile *calibrationFile = TFile::Open(filename);
4566 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
4567 Printf("No calibration file found!!!");
4571 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
4573 Printf("Calibration TList not found!!!");
4577 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
4578 if(!fHistVZEROAGainEqualizationMap) {
4579 Printf("VZERO-A calibration object not found!!!");
4582 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
4583 if(!fHistVZEROCGainEqualizationMap) {
4584 Printf("VZERO-C calibration object not found!!!");
4588 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
4589 if(!fHistVZEROChannelGainEqualizationMap) {
4590 Printf("VZERO channel calibration object not found!!!");
4595 //________________________________________________________________________
4596 Double_t AliTwoParticlePIDCorr::GetChannelEqualizationFactor(Int_t run,Int_t channel) {
4598 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4600 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
4601 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4602 if(gRunNumber == run)
4603 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
4609 //________________________________________________________________________
4610 Double_t AliTwoParticlePIDCorr::GetEqualizationFactor(Int_t run, const char* side) {
4612 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4614 TString gVZEROSide = side;
4615 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
4616 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4617 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
4618 if(gRunNumber == run) {
4619 if(gVZEROSide == "A")
4620 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
4621 else if(gVZEROSide == "C")
4622 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
4628 //________________________________________________________________________
4629 Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliVEvent *mainevent){
4630 //Function that returns the reference multiplicity from AODs (data or reco MC, Not for Truth)
4631 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
4632 if(!mainevent) return -1;
4634 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
4636 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
4637 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
4639 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
4641 Printf("ERROR: AOD header not available");
4644 Int_t gRunNumber = header->GetRunNumber();
4645 Float_t bSign1=header->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
4648 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++)
4649 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
4650 AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
4651 if (!track) continue;
4652 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
4653 if(tracktype!=1) continue;
4655 if(!track) continue;//for safety
4657 gRefMultiplicityTPC += 1.0;
4661 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
4662 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
4663 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
4664 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
4667 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
4668 else if(iChannel >= 32)
4669 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
4672 //Equalization of gain
4673 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
4675 gRefMultiplicityVZEROA /= gFactorA;
4676 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
4678 gRefMultiplicityVZEROC /= gFactorC;
4679 if((gFactorA != 0)&&(gFactorC != 0))
4680 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
4683 //EQVZERO vs TPC multiplicity
4684 fHistEQVZEROvsTPCmultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
4685 fHistEQVZEROAvsTPCmultiplicity->Fill(gRefMultiplicityVZEROA,gRefMultiplicityTPC);
4686 fHistEQVZEROCvsTPCmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityTPC);
4688 //EQVZERO vs VZERO multiplicity
4689 fHistVZEROCvsEQVZEROCmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),gRefMultiplicityVZEROC);
4690 fHistVZEROAvsEQVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0A(),gRefMultiplicityVZEROA);
4692 //VZEROC vs VZEROA multiplicity
4693 fHistVZEROCvsVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),event->GetVZEROData()->GetMTotV0A());
4695 //EQVZEROC vs EQVZEROA multiplicity
4696 fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
4698 fHistRefmult->Fill(3.,gRefMultiplicityTPC);
4699 fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
4700 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
4701 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
4704 if(fCentralityMethod == "TRACKS_MANUAL") gRefMultiplicity = gRefMultiplicityTPC;
4706 else if(fCentralityMethod == "V0M_MANUAL") gRefMultiplicity = gRefMultiplicityVZERO;
4708 else if(fCentralityMethod == "V0A_MANUAL") gRefMultiplicity = gRefMultiplicityVZEROA;
4710 else if(fCentralityMethod == "V0C_MANUAL") gRefMultiplicity = gRefMultiplicityVZEROC;
4712 else gRefMultiplicity = gRefMultiplicityTPC;
4714 return gRefMultiplicity;
4717 //-------------------------------------------------------------------------------------------------------
4718 Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliVEvent *mainevent, Bool_t truth){
4720 if(!mainevent) return -1;
4721 // get centrality object and check quality
4722 Double_t cent_v0=-1;
4723 Bool_t shift_to_TRACKS_MANUAL=kFALSE;//in case of wrong setting automatic shift to Tracks_Manual method
4725 Double_t gRefMultiplicityTPC_Truth = 0.;
4726 Double_t gRefMultiplicityVZERO_Truth = 0., gRefMultiplicityVZEROA_Truth = 0., gRefMultiplicityVZEROC_Truth = 0.;
4728 if(fAnalysisType == "AOD"|| fAnalysisType == "MCAOD") { //centrality in AOD header //++++++++++++++
4729 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
4731 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")//for PbPb, pPb, pp7TeV(still to be introduced)//data or RecoMC and also for TRUTH
4734 if(fSampleType=="pp_7" && fPPVsMultUtils)
4735 {//for pp 7 TeV case only using Alianalysisutils class
4736 if(fAnalysisUtils) cent_v0 = fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,fCentralityMethod);
4738 fHistCentStats->Fill(0.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0A"));
4739 fHistCentStats->Fill(1.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0C"));
4740 fHistCentStats->Fill(2.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0M"));
4741 fHistCentStats->Fill(3.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0AEq"));//only available for LHC10d at present (Quantile info)
4742 fHistCentStats->Fill(4.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0CEq"));//only available for LHC10d at present (Quantile info)
4743 fHistCentStats->Fill(5.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0MEq"));//only available for LHC10d at present (Quantile info)
4746 else if(fSampleType=="pPb" || fSampleType=="PbPb")
4748 AliCentrality *centralityObj=0;
4749 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
4750 if(!header) return -1;
4751 centralityObj = header->GetCentralityP();
4752 // if (centrality->GetQuality() != 0) return ;
4754 fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
4755 fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
4756 fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
4757 fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
4758 fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
4759 fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
4761 fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
4762 fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
4764 cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
4768 else shift_to_TRACKS_MANUAL=kTRUE;
4770 }//centralitymethod condition
4772 else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL" || shift_to_TRACKS_MANUAL)//data or RecoMc and also for TRUTH
4774 if(!truth){//for data or RecoMC
4775 cent_v0 = GetReferenceMultiplicityVZEROFromAOD((AliVEvent*)event);
4776 }//for data or RecoMC
4778 if(truth && (fAnalysisType == "MCAOD")){//condition for TRUTH case
4779 //check for TClonesArray(truth track MC information)
4780 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
4782 //AliFatal("Error: MC particles branch not found!\n");
4785 //now process the truth particles(for both efficiency & correlation function)
4786 Int_t nMCTrack = fArrayMC->GetEntriesFast();
4788 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
4789 {//MC truth track loop starts
4791 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
4794 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
4798 //consider only charged particles
4799 if(partMC->Charge() == 0) continue;
4801 //consider only primary particles; neglect all secondary particles including from weak decays
4802 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
4805 //remove injected signals(primaries above <maxLabel>)
4806 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
4809 Bool_t isduplicate=kFALSE;
4810 if (fRemoveDuplicates)
4812 for (Int_t j=iMC+1; j<nMCTrack; ++j)
4813 {//2nd trutuh loop starts
4814 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
4816 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
4819 if (partMC->GetLabel() == partMC2->GetLabel())
4824 }//2nd truth loop ends
4826 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
4829 // if (fCentralityMethod=="V0M_MANUAL")
4830 if((partMC->Eta() < 5.1 && partMC->Eta() > 2.8) || (partMC->Eta() > -3.7 && partMC->Eta() < -1.7)) gRefMultiplicityVZERO_Truth+=1;
4831 // else if (fCentralityMethod=="V0A_MANUAL") {
4832 if(partMC->Eta() < 5.1 && partMC->Eta() > 2.8) gRefMultiplicityVZEROA_Truth+=1;
4833 // else if (fCentralityMethod=="V0C_MANUAL") {
4834 if(partMC->Eta() < -1.7 && partMC->Eta() > -3.7) gRefMultiplicityVZEROC_Truth+=1;
4835 //else if (fCentralityMethod=="TRACKS_MANUAL") {
4836 if (partMC->Eta() > fmineta && partMC->Eta() < fmaxeta) {
4837 if (partMC->Pt() > fminPt && partMC->Pt() < fmaxPt) gRefMultiplicityTPC_Truth+=1;
4840 }//truth track loop ends
4842 fHistRefmult->Fill(3.,gRefMultiplicityTPC_Truth);
4843 fHistRefmult->Fill(2.,gRefMultiplicityVZERO_Truth);
4844 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA_Truth);
4845 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC_Truth);
4847 if(fCentralityMethod == "TRACKS_MANUAL") cent_v0=gRefMultiplicityTPC_Truth;
4849 else if(fCentralityMethod == "V0M_MANUAL") cent_v0=gRefMultiplicityVZERO_Truth;
4851 else if(fCentralityMethod == "V0A_MANUAL") cent_v0=gRefMultiplicityVZEROA_Truth;
4853 else if(fCentralityMethod == "V0C_MANUAL") cent_v0=gRefMultiplicityVZEROC_Truth;
4855 else cent_v0=gRefMultiplicityTPC_Truth;
4857 }//condition for TRUTH case
4859 }//end of MANUAL method
4861 else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC in AOD production
4863 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4867 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4870 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
4871 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
4872 AliError("Event header not found. Skipping this event.");
4876 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4880 cent_v0 = collGeometry->ImpactParameter();
4881 fhistImpactParm->Fill(cent_v0);
4884 }//end of Impact parameter method
4888 }//AOD OR MCAOD condition
4891 else if(fAnalysisType == "MC"){
4892 Double_t gImpactParameter = -1.;
4893 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(mainevent);
4895 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
4897 gImpactParameter = headerH->ImpactParameter();
4899 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
4900 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
4902 AliError(Form("Could not receive particle %d", iParticle));
4906 //exclude non stable particles
4907 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
4909 if(track->Charge() == 0) continue;
4911 // if (fCentralityMethod=="V0M_MANUAL")
4912 if((track->Eta() < 5.1 && track->Eta() > 2.8) || (track->Eta() > -3.7 && track->Eta() < -1.7)) gRefMultiplicityVZERO_Truth+=1;
4913 // else if (fCentralityMethod=="V0A_MANUAL") {
4914 if(track->Eta() < 5.1 && track->Eta() > 2.8) gRefMultiplicityVZEROA_Truth+=1;
4915 // else if (fCentralityMethod=="V0C_MANUAL") {
4916 if(track->Eta() < -1.7 && track->Eta() > -3.7) gRefMultiplicityVZEROC_Truth+=1;
4917 //else if (fCentralityMethod=="TRACKS_MANUAL") {
4918 if (track->Eta() > fmineta && track->Eta() < fmaxeta) {
4919 if (track->Pt() > fminPt && track->Pt() < fmaxPt) gRefMultiplicityTPC_Truth+=1;}
4921 }//loop over primaries
4923 fHistRefmult->Fill(3.,gRefMultiplicityTPC_Truth);
4924 fHistRefmult->Fill(2.,gRefMultiplicityVZERO_Truth);
4925 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA_Truth);
4926 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC_Truth);
4927 if (fCentralityMethod == "MC_b"){
4928 cent_v0=gImpactParameter;
4929 fhistImpactParm->Fill(gImpactParameter);
4932 else if(fCentralityMethod == "TRACKS_MANUAL") cent_v0=gRefMultiplicityTPC_Truth;
4934 else if(fCentralityMethod == "V0M_MANUAL") cent_v0=gRefMultiplicityVZERO_Truth;
4936 else if(fCentralityMethod == "V0A_MANUAL") cent_v0=gRefMultiplicityVZEROA_Truth;
4938 else if(fCentralityMethod == "V0C_MANUAL") cent_v0=gRefMultiplicityVZEROC_Truth;
4940 else cent_v0=gImpactParameter;//default value is the impact parameter
4951 //-----------------------------------------------------------------------------------------
4952 Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliVEvent *event,Bool_t truth){
4953 //do the event selection(zvtx, pileup, centrality/multiplicity cut) and then return the value of the centrality of that event
4954 if(!event) return -1;
4956 Float_t gRefMultiplicity = -1.;
4958 //***********************************SOURCE CODE-TASKBFPsi
4961 if(fAnalysisType == "MC") {
4962 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
4964 AliError("mcEvent not available");
4969 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
4971 TArrayF gVertexArray;
4972 header->PrimaryVertex(gVertexArray);
4973 //count events having a proper vertex
4974 fEventCounter->Fill(5);
4976 fHistQA[0]->Fill((gVertexArray.At(0)));fHistQA[1]->Fill((gVertexArray.At(1)));fHistQA[2]->Fill((gVertexArray.At(2))); //for trkVtx only before vertex cut |zvtx|<10 cm
4978 if(TMath::Abs(gVertexArray.At(0)) < fVxMax_MC) {
4979 if(TMath::Abs(gVertexArray.At(1)) < fVyMax_MC) {
4980 if(TMath::Abs(gVertexArray.At(2)) < fVzMax_MC) {
4981 //count events after vertex cut
4982 fEventCounter->Fill(7);
4983 fHistQA[3]->Fill((gVertexArray.At(0)));fHistQA[4]->Fill((gVertexArray.At(1)));fHistQA[5]->Fill((gVertexArray.At(2)));//after vertex cut,for trkVtx only
4985 // get the reference multiplicty or centrality
4986 gRefMultiplicity = GetRefMultiOrCentrality((AliVEvent*)mcevent,kFALSE);//2nd argument has no meaning
4988 if(gRefMultiplicity<0) return -1;
4990 // take events only within the multiplicity class mentioned in the custom binning
4991 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
4993 //count events having proper centrality/ref multiplicity
4994 fEventCounter->Fill(9);
5003 else if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD"){// if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD"
5004 //vertex selection(is it fine for PP?)
5005 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
5008 // check first event in chunk (is not needed for new reconstructions)
5009 if(fCheckFirstEventInChunk){
5010 AliAnalysisUtils ut;
5011 if(ut.IsFirstEventInChunk(aod))
5016 AliAnalysisUtils ut;
5017 ut.SetUseMVPlpSelection(kTRUE);
5018 ut.SetUseOutOfBunchPileUp(kTRUE);
5019 if(ut.IsPileUpEvent(aod))
5023 //count events after pileup selection
5024 fEventCounter->Fill(3);
5026 if (fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
5027 trkVtx = aod->GetPrimaryVertex();
5028 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
5029 TString vtxTtl = trkVtx->GetTitle();
5030 if (!vtxTtl.Contains("VertexerTracks")) return -1;
5031 zvtx = trkVtx->GetZ();
5032 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
5033 if (!spdVtx || spdVtx->GetNContributors()<=0) return -1;
5034 TString vtxTyp = spdVtx->GetTitle();
5035 Double_t cov[6]={0};
5036 spdVtx->GetCovarianceMatrix(cov);
5037 Double_t zRes = TMath::Sqrt(cov[5]);
5038 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return -1;
5039 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return -1;
5041 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
5042 Int_t nVertex = aod->GetNumberOfVertices();
5044 trkVtx = aod->GetPrimaryVertex();
5045 Int_t nTracksPrim = trkVtx->GetNContributors();
5046 zvtx = trkVtx->GetZ();
5047 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
5048 // Reject TPC only vertex
5049 TString name(trkVtx->GetName());
5050 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return -1;
5052 // Select a quality vertex by number of tracks?
5053 if( nTracksPrim < fnTracksVertex ) {
5054 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
5057 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
5058 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
5060 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
5065 else if(fVertextype==0){//default case
5066 trkVtx =(AliAODVertex*) aod->GetPrimaryVertex();
5067 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
5068 zvtx = trkVtx->GetZ();
5070 trkVtx->GetCovarianceMatrix(fCov);
5071 if(fCov[5] == 0) return -1;//proper vertex resolution
5074 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
5075 return -1;//as there is no proper sample type
5078 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
5080 //count events having a proper vertex
5081 fEventCounter->Fill(5);
5083 if (TMath::Abs(zvtx) > fzvtxcut) return -1;
5085 //count events after vertex cut
5086 fEventCounter->Fill(7);
5089 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
5091 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
5093 //get the centrality or multiplicity
5094 if(truth) {gRefMultiplicity = GetRefMultiOrCentrality((AliVEvent*)aod,kTRUE);}//kTRUE-->for Truth case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
5096 else {gRefMultiplicity = GetRefMultiOrCentrality((AliVEvent*)aod,kFALSE);}//kFALSE-->for data and RecoMc case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
5098 if(gRefMultiplicity<0) return -1;
5100 // take events only within the multiplicity class mentioned in the custom binning
5101 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
5103 //count events having proper centrality/ref multiplicity
5104 fEventCounter->Fill(9);
5107 // centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
5108 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
5110 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
5114 //count events after rejection due to centrality weighting
5115 fEventCounter->Fill(11);
5118 else gRefMultiplicity=-1;
5120 return gRefMultiplicity;
5123 //--------------------------------------------------------------------------------------------------------
5124 Float_t AliTwoParticlePIDCorr::GetEventPlane(AliVEvent *mainevent,Bool_t truth, Double_t v0Centr)
5126 Float_t eventplane=999.;
5127 // Get the event plane
5128 if(!mainevent) return 999.;
5131 //MC: from reaction plane
5132 if(fAnalysisType == "MC"){
5134 AliError("mcEvent not available");
5138 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(mainevent);
5140 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
5143 // Impact parameter bins(it is only for Pb-Pb)
5144 if(v0Centr < 3.50) iC = 0;
5145 else if(v0Centr < 4.94) iC = 1;
5146 else if(v0Centr < 6.98) iC = 2;
5147 else if(v0Centr < 8.55) iC = 3;
5148 else if(v0Centr < 9.88) iC = 4;
5149 else if(v0Centr < 11.04) iC = 5;
5150 else if(v0Centr < 12.09) iC = 6;
5151 else if(v0Centr < 13.05) iC = 7;
5154 eventplane = headerH->ReactionPlaneAngle();
5155 if(eventplane > TMath::Pi()/2 && eventplane <= TMath::Pi()*3/2) eventplane-=TMath::Pi();
5156 if(eventplane > TMath::Pi()*3/2) eventplane-=2*TMath::Pi();
5157 fHistEventPlaneTruth->Fill(iC,eventplane);
5158 //gReactionPlane *= TMath::RadToDeg();
5163 else if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD") {
5164 //reset Q vector info
5166 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
5169 Int_t run = event->GetRunNumber();
5172 // Load the calibrations run dependent
5173 if(! fIsAfter2011) OpenInfoCalbration(run);
5179 if (v0Centr > 80) return 999.; // analysis only for 0-80% centrality classes
5181 if(v0Centr < 5) iC = 0;
5182 else if(v0Centr < 10) iC = 1;
5183 else if(v0Centr < 20) iC = 2;
5184 else if(v0Centr < 30) iC = 3;
5185 else if(v0Centr < 40) iC = 4;
5186 else if(v0Centr < 50) iC = 5;
5187 else if(v0Centr < 60) iC = 6;
5188 else if(v0Centr < 70) iC = 7;
5194 //reset Q vector info
5195 Double_t Qxa2 = 0, Qya2 = 0;
5196 Double_t Qxc2 = 0, Qyc2 = 0;
5197 Double_t Qxa3 = 0, Qya3 = 0;
5198 Double_t Qxc3 = 0, Qyc3 = 0;
5201 //MC: from reaction plane
5204 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
5206 evplaneMC = header->GetReactionPlaneAngle();//[0, 360]
5207 //make it within [-pi/2,pi/2] to make it general
5208 if(evplaneMC > TMath::Pi()/2 && evplaneMC <= TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
5209 if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
5210 fHistEventPlaneTruth->Fill(iC,evplaneMC);
5212 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
5216 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
5218 if (collGeometry){//get the reaction plane from MC header
5219 gReactionPlane = collGeometry->ReactionPlaneAngle();//[0,180]
5223 //taken from vnv0 code(get the TPC, V0A, V0C event plane using truth tracks)
5224 TClonesArray *mcArray = NULL;
5225 mcArray = (TClonesArray*)event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
5227 Float_t QxMCv2[3] = {0,0,0};
5228 Float_t QyMCv2[3] = {0,0,0};
5229 Float_t QxMCv3[3] = {0,0,0};
5230 Float_t QyMCv3[3] = {0,0,0};
5231 Float_t EvPlaneMCV2[3] = {0,0,0};
5232 Float_t EvPlaneMCV3[3] = {0,0,0};
5233 Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
5234 Float_t etaMax[3] = {4.88,-1.8,0.8};
5236 // analysis on MC tracks
5237 Int_t nMCtrack = mcArray->GetEntries() ;
5239 // EP computation with MC tracks
5240 for(Int_t iT=0;iT < nMCtrack;iT++){
5241 AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
5242 if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
5244 Float_t eta = mctr->Eta();
5245 for(Int_t iD=0;iD<3;iD++){
5246 if(eta > etaMin[iD] && eta < etaMax[iD]){
5247 Float_t phi = mctr->Phi();
5248 QxMCv2[iD] += TMath::Cos(2*phi);
5249 QyMCv2[iD] += TMath::Sin(2*phi);
5250 QxMCv3[iD] += TMath::Cos(3*phi);
5251 QyMCv3[iD] += TMath::Sin(3*phi);
5256 EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
5257 EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
5258 EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
5259 fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
5260 fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
5261 fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
5262 fgPsi2v0aMC = EvPlaneMCV2[0];
5263 fgPsi2v0cMC = EvPlaneMCV2[1];
5264 fgPsi2tpcMC = EvPlaneMCV2[2];
5267 EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
5268 EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
5269 EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
5270 fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
5271 fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
5272 fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
5273 fgPsi3v0aMC = EvPlaneMCV3[0];
5274 fgPsi3v0cMC = EvPlaneMCV3[1];
5275 fgPsi3tpcMC = EvPlaneMCV3[2];
5282 Int_t nAODTracks = event->GetNumberOfTracks();
5284 // TPC EP needed for resolution studies (TPC subevent)
5285 //AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
5286 //Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
5287 Double_t Qx2 = 0, Qy2 = 0;
5288 Double_t Qx3 = 0, Qy3 = 0;
5290 for(Int_t iT = 0; iT < nAODTracks; iT++) {
5292 AliAODTrack* aodTrack =(AliAODTrack*) event->GetTrack(iT);
5298 Bool_t trkFlag = aodTrack->TestFilterBit(1);
5300 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
5303 Double_t b[2] = {-99., -99.};
5304 Double_t bCov[3] = {-99., -99., -99.};
5307 AliAODTrack param(*aodTrack);
5308 if (!param.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov)){
5312 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
5315 Qx2 += TMath::Cos(2*aodTrack->Phi());
5316 Qy2 += TMath::Sin(2*aodTrack->Phi());
5317 Qx3 += TMath::Cos(3*aodTrack->Phi());
5318 Qy3 += TMath::Sin(3*aodTrack->Phi());
5322 Float_t evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
5323 Float_t evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
5325 fgPsi2tpc = evPlAng2;
5326 fgPsi3tpc = evPlAng3;
5328 fPhiRPTPC->Fill(iC,evPlAng2);
5329 fPhiRPTPCv3->Fill(iC,evPlAng3);
5334 AliAODVZERO* aodV0 = event->GetVZEROData();
5336 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
5337 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
5338 Float_t multv0 = aodV0->GetMultiplicity(iv0);
5341 if(fAnalysisType == "AOD"){//not for reco MC tracks, only for real data
5342 if (iv0 < 32){ // V0C
5343 Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5344 Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5345 Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5346 Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5348 Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5349 Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5350 Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5351 Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5355 if (iv0 < 32){ // V0C
5356 Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5357 Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5358 Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5359 Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5361 Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5362 Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5363 Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5364 Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5369 //grab for each centrality the proper histo with the Qx and Qy to do the recentering
5370 Double_t Qxamean2 = fMeanQ[iCcal][1][0];
5371 Double_t Qxarms2 = fWidthQ[iCcal][1][0];
5372 Double_t Qyamean2 = fMeanQ[iCcal][1][1];
5373 Double_t Qyarms2 = fWidthQ[iCcal][1][1];
5374 Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
5375 Double_t Qxarms3 = fWidthQv3[iCcal][1][0];
5376 Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
5377 Double_t Qyarms3 = fWidthQv3[iCcal][1][1];
5379 Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
5380 Double_t Qxcrms2 = fWidthQ[iCcal][0][0];
5381 Double_t Qycmean2 = fMeanQ[iCcal][0][1];
5382 Double_t Qycrms2 = fWidthQ[iCcal][0][1];
5383 Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
5384 Double_t Qxcrms3 = fWidthQv3[iCcal][0][0];
5385 Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
5386 Double_t Qycrms3 = fWidthQv3[iCcal][0][1];
5388 Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
5389 Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
5390 Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
5391 Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
5392 Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
5393 Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
5394 Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
5395 Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
5397 //to calculate 2nd order event plane with v0M
5398 Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
5399 /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
5400 Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
5401 /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
5403 //here the calculated event plane is within -Pi to +Pi(delete it , no use here , only for definition)
5404 Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
5405 Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
5406 Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
5410 Float_t evPlAngV0ACor2=999.;
5411 Float_t evPlAngV0CCor2=999.;
5412 Float_t evPlAngV0ACor3=999.;
5413 Float_t evPlAngV0CCor3=999.;
5416 if(fAnalysisType == "AOD"){
5417 evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
5418 evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
5419 evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
5420 evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
5423 evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
5424 evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
5425 evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
5426 evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
5430 AliEventplane *ep = event->GetEventplane();
5431 evPlAngV0ACor2 = ep->GetEventplane("V0A", event, 2);
5432 evPlAngV0CCor2 = ep->GetEventplane("V0C", event, 2);
5433 evPlAngV0ACor3 = ep->GetEventplane("V0A", event, 3);
5434 evPlAngV0CCor3 = ep->GetEventplane("V0C", event, 3);
5437 fgPsi2v0a = evPlAngV0ACor2;
5438 fgPsi2v0c = evPlAngV0CCor2;
5439 fgPsi3v0a = evPlAngV0ACor3;
5440 fgPsi3v0c = evPlAngV0CCor3;
5442 // Fill EP distribution histograms evPlAng
5444 fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
5445 fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
5447 fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
5448 fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
5450 // Fill histograms needed for resolution evaluation
5451 fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
5452 fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
5453 fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
5455 fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
5456 fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
5457 fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
5461 Float_t gVZEROEventPlane = -10.;
5462 Float_t gReactionPlane = -10.;
5463 Double_t qxTot = 0.0, qyTot = 0.0;
5465 AliEventplane *ep = event->GetEventplane();
5467 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
5468 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
5469 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
5470 gReactionPlane = gVZEROEventPlane;
5474 //return gReactionPlane;
5476 //make the final 2nd order event plane within 0 to Pi
5477 //using data and reco tracks only
5478 if(fgPsi2v0a!=999. && fgPsi2v0a < 0.) fgPsi2v0a += TMath::Pi();
5479 if(fgPsi2v0c!=999. && fgPsi2v0c < 0.) fgPsi2v0c += TMath::Pi();
5480 if(fgPsi2tpc!=999. && fgPsi2tpc < 0.) fgPsi2tpc += TMath::Pi();
5481 //using truth tracks only
5482 if(evplaneMC!=999. && evplaneMC < 0.) evplaneMC += TMath::Pi();
5483 if(fgPsi2v0aMC!=999. && fgPsi2v0aMC < 0.) fgPsi2v0aMC += TMath::Pi();
5484 if(fgPsi2v0cMC!=999. && fgPsi2v0cMC < 0.) fgPsi2v0cMC += TMath::Pi();
5485 if(fgPsi2tpcMC!=999. && fgPsi2tpcMC < 0.) fgPsi2tpcMC += TMath::Pi();
5486 //for the time being leave the 3rd order event planes within -pi/3 t0 +pi/3
5488 if(truth){//for truth MC
5489 if(fV2 && fEPdet=="header")eventplane=evplaneMC;
5490 if(fV2 && fEPdet=="V0A")eventplane=fgPsi2v0aMC;
5491 if(fV2 && fEPdet=="V0C")eventplane=fgPsi2v0cMC;
5492 if(fV2 && fEPdet=="TPC")eventplane=fgPsi2tpcMC;
5494 if(fV3 && fEPdet=="V0A")eventplane=fgPsi3v0aMC;
5495 if(fV3 && fEPdet=="V0C")eventplane=fgPsi3v0cMC;
5496 if(fV3 && fEPdet=="TPC")eventplane=fgPsi3tpcMC;
5498 else{//for data and recoMC
5499 if(fV2 && fEPdet=="V0A")eventplane=fgPsi2v0a;
5500 if(fV2 && fEPdet=="V0C")eventplane=fgPsi2v0c;
5501 if(fV2 && fEPdet=="TPC")eventplane=fgPsi2tpc;
5503 if(fV3 && fEPdet=="V0A")eventplane=fgPsi3v0a;
5504 if(fV3 && fEPdet=="V0C")eventplane=fgPsi3v0c;
5505 if(fV3 && fEPdet=="TPC")eventplane=fgPsi3tpc;
5510 }//AOD/MCAOD condition
5512 else eventplane=999.;
5517 //------------------------------------------------------------------------------------------------------------------
5518 void AliTwoParticlePIDCorr::OpenInfoCalbration(Int_t run){
5519 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
5520 TFile *foadb = TFile::Open(oadbfilename.Data());
5523 printf("OADB file %s cannot be opened\n",oadbfilename.Data());
5527 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
5529 printf("OADB object hMultV0BefCorr is not available in the file\n");
5533 if(!(cont->GetObject(run))){
5534 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
5537 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
5539 TF1 *fpol0 = new TF1("fpol0","pol0");
5540 fMultV0->Fit(fpol0,"","",0,31);
5541 fV0Cpol = fpol0->GetParameter(0);
5542 fMultV0->Fit(fpol0,"","",32,64);
5543 fV0Apol = fpol0->GetParameter(0);
5545 for(Int_t iside=0;iside<2;iside++){
5546 for(Int_t icoord=0;icoord<2;icoord++){
5547 for(Int_t i=0;i < 9;i++){
5549 if(iside==0 && icoord==0)
5550 snprintf(namecont,100,"hQxc2_%i",i);
5551 else if(iside==1 && icoord==0)
5552 snprintf(namecont,100,"hQxa2_%i",i);
5553 else if(iside==0 && icoord==1)
5554 snprintf(namecont,100,"hQyc2_%i",i);
5555 else if(iside==1 && icoord==1)
5556 snprintf(namecont,100,"hQya2_%i",i);
5558 cont = (AliOADBContainer*) foadb->Get(namecont);
5560 printf("OADB object %s is not available in the file\n",namecont);
5564 if(!(cont->GetObject(run))){
5565 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5568 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5569 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5572 if(iside==0 && icoord==0)
5573 snprintf(namecont,100,"hQxc3_%i",i);
5574 else if(iside==1 && icoord==0)
5575 snprintf(namecont,100,"hQxa3_%i",i);
5576 else if(iside==0 && icoord==1)
5577 snprintf(namecont,100,"hQyc3_%i",i);
5578 else if(iside==1 && icoord==1)
5579 snprintf(namecont,100,"hQya3_%i",i);
5581 cont = (AliOADBContainer*) foadb->Get(namecont);
5583 printf("OADB object %s is not available in the file\n",namecont);
5587 if(!(cont->GetObject(run))){
5588 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5591 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5592 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5598 //____________________________________________________________________
5599 void AliTwoParticlePIDCorr::FillPIDEventPlane(Double_t centrality,Int_t par,Float_t trigphi,Float_t fReactionPlane)
5602 // Event plane (determine psi bin)
5603 Double_t gPsiMinusPhi = 0.;
5604 Double_t gPsiMinusPhiBin = -10.;
5605 if(fRequestEventPlane){
5606 gPsiMinusPhi = TMath::Abs(trigphi - fReactionPlane);
5608 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
5609 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
5610 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
5611 gPsiMinusPhiBin = 0.0;
5613 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
5614 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
5615 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
5616 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
5617 gPsiMinusPhiBin = 1.0;
5619 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
5620 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
5621 gPsiMinusPhiBin = 2.0;
5624 gPsiMinusPhiBin = 3.0;
5626 fEventPlanePID->Fill(centrality,gPsiMinusPhiBin,(Float_t)par);
5630 //____________________________________________________________________________________________________
5632 TObjArray* AliTwoParticlePIDCorr::GetV0Particles(AliVEvent* event,Double_t Centrality)
5635 AliAODEvent* fAOD = dynamic_cast<AliAODEvent*>(event);
5637 //function to select v0's from AODs
5638 trkVtx=fAOD->GetPrimaryVertex();
5639 Float_t xv=trkVtx->GetX(), yv=trkVtx->GetY(), zv=trkVtx->GetZ();
5640 Int_t nV0sTot = fAOD->GetNumberOfV0s();
5642 TObjArray * selectedV0s = new TObjArray;
5643 selectedV0s->SetOwner(kTRUE);
5645 for (Int_t iV0 = 0; iV0 < nV0sTot; iV0++)
5648 AliAODv0 *v0=fAOD->GetV0(iV0);
5650 if(!CheckStatusv0(v0)) continue;
5652 AliAODTrack *ptrack=(AliAODTrack*)v0->GetDaughter(0);
5653 AliAODTrack *ntrack=(AliAODTrack*)v0->GetDaughter(1);
5655 Bool_t cutK0sPID=kFALSE;
5656 Bool_t cutLambdaPID=kFALSE;
5657 Bool_t cutAntiLambdaPID=kFALSE;
5659 if(fUsev0DaughterPID)
5661 //use fHelperPID check PID of the daughter tracks
5662 //v0 daughter PID may be helpful in distangling k0S and (Anti)Lamda
5664 Int_t PIDptrack = GetParticle(ptrack,kFALSE);
5665 Int_t PIDntrack = GetParticle(ntrack ,kFALSE);
5667 if(PIDptrack ==0 && PIDntrack == 0) cutK0sPID=kTRUE;
5669 if(PIDptrack==2 && PIDntrack ==0) cutLambdaPID=kTRUE;
5671 if (PIDptrack ==0 && PIDntrack == 2) cutAntiLambdaPID=kTRUE;
5675 // effective mass calculations for each hypothesis(without daughter PID)
5676 Double_t InvMassK0s = v0->MassK0Short();
5677 Double_t InvMassAntiLambda = v0->MassAntiLambda();
5678 Double_t InvMassLambda = v0->MassLambda();
5680 Float_t v0Pt=TMath::Sqrt(v0->Pt2V0());
5681 Float_t v0Eta=v0->Eta();
5682 Float_t v0Phi=v0->Phi();
5684 //This is simply raw v0 without any specialised cut
5685 fHistRawPtCentInvK0s->Fill(InvMassK0s,v0Pt,Centrality);
5686 fHistRawPtCentInvLambda->Fill(InvMassLambda,v0Pt,Centrality);
5687 fHistRawPtCentInvAntiLambda->Fill(InvMassAntiLambda,v0Pt,Centrality);
5691 v0->GetSecondaryVtx(xyz);
5693 dx=xyz[0]-xv, dy=xyz[1]-yv, dz=xyz[2]-zv;
5695 // Float_t v0DecayRadius=TMath::Sqrt(dx*dx + dy*dy);
5696 Float_t v0DecayLength=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
5697 // VO's main characteristics to check the reconstruction cuts
5698 // Float_t DcaV0Daughters = v0->DcaV0Daughters();
5699 Float_t V0cosPointAngle = v0->CosPointingAngle(trkVtx);
5700 // Float_t DcaPosToPrimVertex = v0->DcaPosToPrimVertex();
5701 //Float_t DcaNegToPrimVertex = v0->DcaNegToPrimVertex();
5702 //Float_t Dcav0PVz = v0->DcaV0ToPrimVertex();
5703 Float_t v0Pz=v0->Pz();
5704 Float_t v0P= TMath::Sqrt(v0Pt*v0Pt + v0Pz*v0Pz);
5706 Float_t ctauLambda =999.;
5707 Float_t ctauAntiLambda = 999.;
5708 Float_t ctauK0s = 999.;
5711 ctauLambda = (v0DecayLength*InvMassLambda)/v0P;
5712 ctauAntiLambda = (v0DecayLength*InvMassAntiLambda)/v0P;
5713 ctauK0s = (v0DecayLength*InvMassK0s)/v0P;
5716 Bool_t ctauCutK0s= ctauK0s < NCtau*fCutctauK0s ; //ctauK0s 2.68 cm, mean life time of K0s is 8.95 x10^(-11)
5717 Bool_t ctauCutLambda = ctauLambda < NCtau*fCutctauLambda; //ctauLambda 7.8 cm ,mean life is 2.6 x10 ^(-10) ***** 3xctau is the accepted limit
5718 Bool_t ctauCutAntiLambda= ctauAntiLambda < NCtau*fCutctauAntiLambda;
5720 Bool_t RapCutK0s = v0->RapK0Short() < fRapCutK0s;
5721 Bool_t RapCutLambda = v0->RapLambda() < fRapCutLambda;
5722 Bool_t RapCutAntiLambda = v0->Y(-3122) < fRapCutLambda;
5724 Bool_t CPACut= V0cosPointAngle > fMinCPA; //cosine of pointing angle with v0 should be greater than 0.998
5726 //Now we put a loose mass cut which will be tightened later
5727 Bool_t MassCutLooseK0s=(TMath::Abs(InvMassK0s - 0.497614) < 0.1);
5728 Bool_t MassCutLooseLambda=(TMath::Abs(InvMassLambda - 1.115683) < 0.1); // cut is same for Anti-Lambda
5729 Bool_t MassCutLooseAntiLambda=(TMath::Abs(InvMassAntiLambda - 1.115683) < 0.1); // cut is same for Anti-Lambda
5731 //Special Cut for Kshort arementeros podalanski plot
5732 Bool_t ArmenterosCut =kFALSE;
5733 if(ctauCutK0s && RapCutK0s && CPACut && MassCutLooseK0s)
5736 Float_t lAlphaV0 = v0->AlphaV0();
5737 Float_t lPtArmV0 = v0->PtArmV0();
5739 ArmenterosCut = lPtArmV0 > TMath::Abs(0.2*lAlphaV0);
5743 Bool_t IskShortOk=(ctauCutK0s && RapCutK0s && CPACut && MassCutLooseK0s && ArmenterosCut);
5745 Bool_t IsLambdaOk=(ctauCutLambda && RapCutLambda && CPACut && MassCutLooseLambda);
5747 Bool_t IsAntiLambdaOk=(ctauCutAntiLambda && RapCutAntiLambda && CPACut && MassCutLooseAntiLambda);
5749 //Difference on Lambda and Anti-Lambda can be made through daughter PID
5752 Int_t particletype=999;
5754 if( IskShortOk || cutK0sPID )
5756 fHistFinalPtCentInvK0s->Fill(InvMassK0s,v0Pt,Centrality);
5759 Short_t chargeval=0;
5760 Float_t effmatrix=1.0;
5761 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassK0s,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5762 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5763 selectedV0s->Add(copy1);
5768 if(IsLambdaOk || cutLambdaPID)
5770 fHistFinalPtCentInvLambda->Fill(InvMassLambda,v0Pt,Centrality);
5771 //Add in the LRCParticle and give Lambda a tag 5
5774 Short_t chargeval=0;
5775 Float_t effmatrix=1.0;
5776 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassLambda,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5777 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5778 selectedV0s->Add(copy1);
5781 if(IsAntiLambdaOk || cutAntiLambdaPID)
5783 fHistFinalPtCentInvLambda->Fill(InvMassAntiLambda,v0Pt,Centrality);
5784 //Add in the LRCParticle and give Lambda a tag 6
5785 particletype=SpALam;
5786 Short_t chargeval=0;
5787 Float_t effmatrix=1.0;
5788 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassAntiLambda,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5789 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5790 selectedV0s->Add(copy1);
5799 //___________________________________________________________________
5800 Bool_t AliTwoParticlePIDCorr :: CheckStatusv0Daughter(AliAODTrack *t1 ,AliAODTrack *t2)
5802 if (!t1->IsOn(AliAODTrack::kTPCrefit) || !t2->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
5803 // Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
5804 if(t1->GetTPCClusterInfo(2,1)<fDaugNClsTPC || t2->GetTPCClusterInfo(2,1)<fDaugNClsTPC) return kFALSE ;
5806 // ---------------- Fraction of TPC Shared Cluster
5807 Float_t fracPosDaugTPCSharedMap = GetFractionTPCSharedCls(t1);
5808 Float_t fracNegDaugTPCSharedMap = GetFractionTPCSharedCls(t2);
5810 if( (fracPosDaugTPCSharedMap > fFracTPCcls) || (fracNegDaugTPCSharedMap > fFracTPCcls) )
5816 //___________________________________________________________________________________________
5818 Float_t AliTwoParticlePIDCorr :: GetFractionTPCSharedCls( AliAODTrack *track)
5820 // Rejects tracks with shared clusters after filling a control histogram
5821 // This overload is used for primaries
5823 // Get the shared maps
5824 const TBits sharedMap = track->GetTPCSharedMap();
5826 return 1.*sharedMap.CountBits()/track->GetTPCNclsF();
5829 //______________________________________________________________________
5830 Bool_t AliTwoParticlePIDCorr :: CheckStatusv0(AliAODv0 *v1)
5833 // Offline reconstructed V0 only
5834 if (v1->GetOnFlyStatus()) return kFALSE;
5836 AliAODTrack *ptrack=(AliAODTrack *)v1->GetDaughter(0);
5837 AliAODTrack *ntrack=(AliAODTrack *)v1->GetDaughter(1);
5839 if(!ptrack || !ntrack) return kFALSE;
5841 if(ptrack->Charge()==-1 || ntrack->Charge()==1) return kFALSE; //remove wrongly identified charge pairs
5843 if(ptrack->Charge()==0 || ntrack->Charge()==0) return kFALSE; //remove uncharged pairs
5845 if(ptrack->Charge() == ntrack->Charge()) return kFALSE; //remove like sign pairs
5847 if(!CheckStatusv0Daughter(ptrack,ntrack)) return kFALSE;//daughters need to pass some basic cuts
5849 if(TMath::Abs(ptrack->Eta()) > fmaxeta || TMath::Abs(ntrack->Eta()) > fmaxeta) return kFALSE; // remove daughters beyond eta bound |0.8|
5851 if(ptrack->Pt() < fMinPtDaughter || ntrack->Pt() < fMinPtDaughter) return kFALSE; // remove daughter tracks below minmum p |1.0 GeV/c|
5853 if(ptrack->Pt() > fMaxPtDaughter || ntrack->Pt() > fMaxPtDaughter) return kFALSE; // remove daughter tracks above maximum p ** to make it compatiable with AliHelperPID**|4.0 GeV/C|
5855 // Daughters: Impact parameter of daughter to prim vtx
5856 Float_t xy = v1->DcaPosToPrimVertex();
5857 if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE; //0.1 cm
5858 xy = v1->DcaNegToPrimVertex();
5859 if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE; //0.1 cm
5862 Float_t dca = v1->DcaV0Daughters();
5863 if (dca>fMaxDCADaughter) return kFALSE; //1.0 cm
5865 // V0: Cosine of the pointing angle
5866 Float_t cpa=v1->CosPointingAngle(trkVtx); //0.997
5867 if (cpa<fMinCPA) return kFALSE;
5869 // V0: Fiducial volume
5870 Double_t xyz[3]; v1->GetSecondaryVtx(xyz);
5871 Float_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
5872 if (r2<5.*5.) return kFALSE;
5873 if (r2>lMax*lMax) return kFALSE; //lmax=100 cm
5879 //__________________________________________________________________________________________
5880 Bool_t AliTwoParticlePIDCorr::IsTrackFromV0(AliAODEvent* fAOD,AliAODTrack* track)
5882 //to check whether a daughter being taken as associated
5883 Int_t assoID = track->GetID();
5885 for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
5886 AliAODv0* aodV0 = fAOD->GetV0(i);
5888 AliAODTrack *trackPos=(AliAODTrack *)(aodV0->GetDaughter(0));
5889 AliAODTrack *trackNeg=(AliAODTrack *)(aodV0->GetDaughter(1));
5892 Int_t negID = trackNeg->GetID();
5893 Int_t posID = trackPos->GetID();
5895 if ((TMath::Abs(negID)+1)==(TMath::Abs(assoID))){ return kTRUE;}
5896 if ((TMath::Abs(posID)+1)==(TMath::Abs(assoID))){ return kTRUE;}
5897 //----------------------------------
5904 //----------------------------------------------------------
5905 void AliTwoParticlePIDCorr::Terminate(Option_t *)
5907 // Draw result to screen, or perform fitting, normalizations
5908 // Called once at the end of the query
5909 fOutput = dynamic_cast<TList*> (GetOutputData(1));
5910 if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
5914 //------------------------------------------------------------------