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 //fill tracking efficiency
2508 if(particletypeMC==SpPion || particletypeMC==SpKaon)
2510 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) {
2511 fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
2512 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
2515 if(particletypeMC==SpKaon || particletypeMC==SpProton)
2517 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) {
2518 fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
2519 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
2522 if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) {
2523 fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
2524 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
2526 if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
2527 fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
2528 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
2530 if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
2531 fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
2532 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
2536 //remove the tracks which don't have proper TOF response-otherwise the misIDentification rate values will be wrong(do it after the efficiency filling is complete)
2537 if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && (!HasTOFPID(PIDtrack)) ) continue;
2539 //for misidentification fraction calculation(do it with tuneonPID)
2540 if(particletypeMC==SpPion )
2542 if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
2543 if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
2544 if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
2545 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
2547 if(particletypeMC==SpKaon )
2549 if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
2550 if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
2551 if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
2552 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
2554 if(particletypeMC==SpProton )
2556 if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
2557 if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
2558 if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
2559 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
2561 if(particletypeMC==SpUndefined )//these undefined are not due to absence of proper TOF response, rather due to the PID method only
2563 if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
2564 if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
2565 if(TMath::Abs(pdgCode)==2212) fUNIDcont->Fill(5.,track->Pt());
2566 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fUNIDcont->Fill(7.,track->Pt());
2569 if(particletypeMC==SpUndefined) continue;
2571 if(fRequestEventPlane){
2572 FillPIDEventPlane(cent_v0,particletypeMC,track->Phi(),gReactionPlane);
2575 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2577 Short_t chargeval=0;
2578 if(track->Charge()>0) chargeval=1;
2579 if(track->Charge()<0) chargeval=-1;
2580 if(chargeval==0) continue;
2581 if (fapplyTrigefficiency || fapplyAssoefficiency)
2582 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
2583 LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2584 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2585 tracksID->Add(copy1);
2587 }// if(tracktype==1) condition structure ands
2589 }//reco track loop ends
2591 //*************************************************************still in event loop
2596 //fill the centrality/multiplicity distribution of the selected events
2597 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2599 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??????
2601 //count selected events having centrality betn 0-100%
2602 fEventCounter->Fill(13);
2604 //***************************************event no. counting
2605 Bool_t isbaryontrig=kFALSE;
2606 Bool_t ismesontrig=kFALSE;
2607 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2609 if(tracksID && tracksID->GetEntriesFast()>0)
2611 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2612 { //trigger loop starts
2613 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2615 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2616 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2617 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2618 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2620 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2621 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2626 tracksIDV0=GetV0Particles((AliVEvent*) aod, cent_v0);
2627 if(tracksIDV0->GetEntriesFast()<=0) return;
2629 //same event delte-eta, delta-phi plot
2630 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2631 {//same event calculation starts
2632 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2633 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)
2636 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2637 {//same event calculation starts
2638 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) {
2639 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2641 else Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2643 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2646 //still in main event loop
2648 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2649 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2650 if (pool && pool->IsReady())
2651 {//start mixing only when pool->IsReady
2652 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2653 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2654 { //pool event loop start
2655 TObjArray* bgTracks = pool->GetEvent(jMix);
2656 if(!bgTracks) continue;
2657 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2658 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2659 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2661 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2663 else Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2665 }// pool event loop ends mixing case
2667 } //if pool->IsReady() condition ends mixing case
2670 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2672 }//mixing with unidentified particles
2674 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2675 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2676 if (pool1 && pool1->IsReady())
2677 {//start mixing only when pool->IsReady
2678 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2679 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2680 { //pool event loop start
2681 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2682 if(!bgTracks2) continue;
2683 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2684 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2685 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2686 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2688 }// pool event loop ends mixing case
2689 } //if pool1->IsReady() condition ends mixing case
2693 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2695 }//mixing with identified particles
2697 //no. of events analyzed
2698 fEventCounter->Fill(15);
2701 if(tracksUNID) delete tracksUNID;
2703 if(tracksID) delete tracksID;
2705 if(tracksIDV0) delete tracksIDV0;
2709 }//AOD || MCAOD condition
2711 //still in the main event loop
2714 //________________________________________________________________________
2715 void AliTwoParticlePIDCorr::doAODevent()
2719 AliVEvent *event = InputEvent();
2720 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
2721 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
2723 AliError("Cannot get the AOD event");
2726 Double_t Inv_mass=0.0;
2729 fEventCounter->Fill(1);
2731 fPID = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
2732 if (!fPID) AliFatal("This Task needs the PID response attached to the inputHandler");
2733 //if (!fPID) return;//this should be available with each event even if we don't do PID selection
2741 gReactionPlane = 999.;
2743 Double_t cent_v0= -999.;
2744 Double_t effcent=1.0;
2746 Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
2749 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
2750 Float_t bSign1=aod->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
2753 // check event cuts and fill event histograms and return the centrality or reference multiplicity value
2754 if((cent_v0 = GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE)) < 0){
2757 effcent=cent_v0;//required for efficiency correction case********Extremely Important
2758 //get the event plane in case of PbPb
2759 if(fRequestEventPlane){
2760 gReactionPlane = GetEventPlane((AliVEvent*)aod,kFALSE,cent_v0);
2761 if(gReactionPlane==999.) return;
2765 TExMap *trackMap = new TExMap();
2766 // --- track loop for mapping matrix
2769 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2770 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2771 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2772 if (!track) continue;
2773 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
2774 if(tracktype!=1) continue;
2776 if(!track) continue;//for safety
2778 Int_t gid = track->GetID();
2779 trackMap->Add(gid,itrk);
2783 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
2784 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
2786 TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
2787 tracksID->SetOwner(kTRUE); // IMPORTANT!
2789 //get the selected v0 particle TObjArray
2790 TObjArray* tracksIDV0= new TObjArray;
2791 tracksIDV0->SetOwner(kTRUE); // IMPORTANT!
2795 Bool_t fTrigPtmin1=kFALSE;
2796 Bool_t fTrigPtmin2=kFALSE;
2797 Bool_t fTrigPtJet=kFALSE;
2799 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2800 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2801 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2802 if (!track) continue;
2803 fHistQA[11]->Fill(track->GetTPCNcls());
2804 Int_t particletype=-9999;//required for PID filling
2805 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2806 if(tracktype!=1) continue;
2808 if(!track) continue;//for safety
2811 if(IsTrackFromV0(aod,track)) continue;// remove auto correlation
2813 AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2815 if(fFilterBit==128){
2816 Int_t gid1 = track->GetID();
2817 //if(gid1>=0) PIDtrack = track;
2818 PIDtrack =(AliAODTrack*) aod->GetTrack(trackMap->GetValue(-1-gid1));
2819 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2822 //check for eta , phi holes
2823 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2824 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2828 //if no applyefficiency , set the eff factor=1.0
2829 Float_t effmatrix=1.0;
2831 //tag all particles as unidentified that passed filterbit & kinematic cuts
2832 particletype=unidentified;
2834 //To count the no. of tracks having an accepted track in a certain PT(e.g. Jet Pt) range
2835 if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
2836 if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
2837 if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
2840 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
2843 //to reduce memory consumption in pool
2844 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2846 //Clone & Reduce track list(TObjArray) for unidentified particles
2847 Short_t chargeval=0;
2848 if(track->Charge()>0) chargeval=1;
2849 if(track->Charge()<0) chargeval=-1;
2850 if(chargeval==0) continue;
2851 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2852 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
2853 LRCParticlePID* copy = new LRCParticlePID(particletype,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2854 copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2855 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2858 //now start the particle identificaion process:)
2860 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
2862 Float_t dEdx = PIDtrack->GetTPCsignal();
2863 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2865 //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)
2866 if(HasTOFPID(PIDtrack))
2868 Double_t beta = GetBeta(PIDtrack);
2869 fHistoTOFbeta->Fill(track->Pt(), beta);
2873 //track identification(using nsigma method)
2874 particletype=GetParticle(PIDtrack,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
2876 //2-d TPCTOF map(for each Pt interval)
2877 if(HasTOFPID(PIDtrack)){
2878 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2879 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2880 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2883 //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!!!!!
2884 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)!!!!!!!!!!!
2886 if(fRequestEventPlane){
2887 FillPIDEventPlane(cent_v0,particletype,track->Phi(),gReactionPlane);
2891 //Pt, Eta , Phi distribution of the reconstructed identified particles
2894 if (particletype==SpPion)
2896 fPionPt->Fill(track->Pt());
2897 fPionEta->Fill(track->Eta());
2898 fPionPhi->Fill(track->Phi());
2900 if (particletype==SpKaon)
2902 fKaonPt->Fill(track->Pt());
2903 fKaonEta->Fill(track->Eta());
2904 fKaonPhi->Fill(track->Phi());
2906 if (particletype==SpProton)
2908 fProtonPt->Fill(track->Pt());
2909 fProtonEta->Fill(track->Eta());
2910 fProtonPhi->Fill(track->Phi());
2914 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2916 Short_t chargeval=0;
2917 if(track->Charge()>0) chargeval=1;
2918 if(track->Charge()<0) chargeval=-1;
2919 if(chargeval==0) continue;
2920 if (fapplyTrigefficiency || fapplyAssoefficiency)
2921 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
2922 LRCParticlePID* copy1 = new LRCParticlePID(particletype,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2923 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2924 tracksID->Add(copy1);
2926 } //track loop ends but still in event loop
2928 if(trackscount<1.0){
2929 if(tracksUNID) delete tracksUNID;
2930 if(tracksID) delete tracksID;
2934 if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
2935 if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
2936 if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
2938 Float_t weightval=1.0;
2942 //fill the centrality/multiplicity distribution of the selected events
2943 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2945 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??????
2947 //count selected events having centrality betn 0-100%
2948 fEventCounter->Fill(13);
2950 //***************************************event no. counting
2951 Bool_t isbaryontrig=kFALSE;
2952 Bool_t ismesontrig=kFALSE;
2953 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2955 if(tracksID && tracksID->GetEntriesFast()>0)
2957 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2958 { //trigger loop starts
2959 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2961 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2962 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2963 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2964 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2966 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2967 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2970 //Get the TObjArray of V0 particles
2972 tracksIDV0=GetV0Particles((AliVEvent*) aod,cent_v0);
2973 if(tracksIDV0->GetEntriesFast()<=0) return;
2976 //same event delta-eta-deltaphi plot
2977 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2978 {//same event calculation starts
2979 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2980 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)
2984 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2985 {//same event calculation starts
2986 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID){
2987 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2988 else Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2990 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2993 //still in main event loop
2997 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2998 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2999 if (pool && pool->IsReady())
3000 {//start mixing only when pool->IsReady
3001 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
3002 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
3003 { //pool event loop start
3004 TObjArray* bgTracks = pool->GetEvent(jMix);
3005 if(!bgTracks) continue;
3006 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
3007 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
3008 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
3010 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
3011 else Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
3013 }// pool event loop ends mixing case
3014 } //if pool->IsReady() condition ends mixing case
3017 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
3019 }//mixing with unidentified particles
3022 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
3023 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
3024 if (pool1 && pool1->IsReady())
3025 {//start mixing only when pool->IsReady
3026 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
3027 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
3028 { //pool event loop start
3029 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
3030 if(!bgTracks2) continue;
3031 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
3032 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
3033 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
3034 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
3036 }// pool event loop ends mixing case
3037 } //if pool1->IsReady() condition ends mixing case
3041 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
3043 }//mixing with identified particles
3046 //no. of events analyzed
3047 fEventCounter->Fill(15);
3049 //still in main event loop
3052 if(tracksUNID) delete tracksUNID;
3054 if(tracksID) delete tracksID;
3057 if(tracksIDV0) delete tracksIDV0;
3060 } // *************************event loop ends******************************************
3061 //_______________________________________________________________________
3062 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
3064 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
3066 TObjArray* tracksClone = new TObjArray;
3067 tracksClone->SetOwner(kTRUE);
3069 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
3071 LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
3072 LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->GetInvMass(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval(),particle->GetTPCPadMap(),particle->GetTPCSharedMap());
3073 copy100->SetUniqueID(particle->GetUniqueID());
3074 tracksClone->Add(copy100);
3080 //--------------------------------------------------------------------------------
3081 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)
3084 //before calling this function check that either trackstrig & tracksasso are available
3086 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
3087 TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
3088 TArrayF eta(input->GetEntriesFast());
3089 for (Int_t i=0; i<input->GetEntriesFast(); i++)
3090 eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
3093 Int_t jmax=trackstrig->GetEntriesFast();
3095 jmax=tracksasso->GetEntriesFast();
3097 // identify K, Lambda candidates and flag those particles
3098 // a TObject bit is used for this
3099 const UInt_t kResonanceDaughterFlag = 1 << 14;
3100 if (fRejectResonanceDaughters > 0)
3102 Double_t resonanceMass = -1;
3103 Double_t massDaughter1 = -1;
3104 Double_t massDaughter2 = -1;
3105 const Double_t interval = 0.02;
3106 switch (fRejectResonanceDaughters)
3108 case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
3109 case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
3110 case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
3111 default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
3114 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
3115 trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
3116 for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
3117 tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
3119 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
3121 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3122 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
3124 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
3126 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
3127 if (trig->IsEqual(asso)) continue;
3129 if (trig->Charge() * asso->Charge() > 0) continue;
3131 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
3133 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
3135 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
3137 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
3139 trig->SetBit(kResonanceDaughterFlag);
3140 asso->SetBit(kResonanceDaughterFlag);
3142 // Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
3149 //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
3151 Float_t TriggerPtMin=fminPtTrig;
3152 Float_t TriggerPtMax=fmaxPtTrig;
3153 Int_t HighestPtTriggerIndx=-99999;
3154 TH1* triggerWeighting = 0;
3156 if(fSelectHighestPtTrig || fWeightPerEvent)//**************add this data member to the constructor
3158 if (fWeightPerEvent)
3161 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID) axis = fTHnTrigcount->GetGrid(0)->GetGrid()->GetAxis(2);
3162 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) axis = fTHnTrigcountMCTruthPrim->GetGrid(0)->GetGrid()->GetAxis(2);
3163 triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
3165 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
3166 { //trigger loop starts(highest Pt trigger selection)
3167 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3169 Float_t trigpt=trig->Pt();
3170 //to avoid overflow qnd underflow
3171 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
3172 Int_t particlepidtrig=trig->getparticle();
3173 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
3175 Float_t trigeta=trig->Eta();
3177 // some optimization
3178 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
3181 if (fOnlyOneEtaSide != 0)
3183 if (fOnlyOneEtaSide * trigeta < 0)
3186 if (fTriggerSelectCharge != 0)
3187 if (trig->Charge() * fTriggerSelectCharge < 0)
3190 if (fRejectResonanceDaughters > 0)
3191 if (trig->TestBit(kResonanceDaughterFlag)) continue;
3193 if(fSelectHighestPtTrig){
3194 if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
3196 HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
3197 TriggerPtMin=trigpt;
3201 if (fWeightPerEvent) triggerWeighting->Fill(trigpt);
3203 }//trigger loop ends(highest Pt trigger selection)
3205 }//******************(fSelectHighestPtTrig || fWeightPerEvent) condition ends
3208 //two particle correlation filling
3209 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
3210 { //trigger loop starts
3211 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3213 Float_t trigpt=trig->Pt();
3214 //to avoid overflow qnd underflow
3215 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
3216 Double_t ParticlePID_InvMass=0.0;
3217 if(fV0TrigCorr) ParticlePID_InvMass=trig->GetInvMass();
3219 Int_t particlepidtrig=trig->getparticle();
3220 ParticlePID_InvMass=(Double_t) particlepidtrig;
3221 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}//***********************************forks,lam.Alam their PID numbers have no meaning, only their Inv_mass will be stored
3224 Float_t trigeta=trig->Eta();
3226 // some optimization
3227 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
3230 if (fOnlyOneEtaSide != 0)
3232 if (fOnlyOneEtaSide * trigeta < 0)
3235 if (fTriggerSelectCharge != 0)
3236 if (trig->Charge() * fTriggerSelectCharge < 0)
3239 if (fRejectResonanceDaughters > 0)
3240 if (trig->TestBit(kResonanceDaughterFlag)) continue;
3242 if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
3243 if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
3246 Float_t trigphi=trig->Phi();
3247 Float_t trackefftrig=1.0;
3248 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
3250 // Event plane (determine psi bin)
3251 Double_t gPsiMinusPhi = 0.;
3252 Double_t gPsiMinusPhiBin = -10.;
3253 if(fRequestEventPlane){
3254 gPsiMinusPhi = TMath::Abs(trigphi - ReactionPlane);
3255 //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)
3256 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
3257 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
3258 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
3259 gPsiMinusPhiBin = 0.0;
3261 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
3262 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
3263 gPsiMinusPhiBin = 0.0;
3266 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
3267 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
3268 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
3269 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
3270 gPsiMinusPhiBin = 1.0;
3272 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
3273 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
3274 gPsiMinusPhiBin = 2.0;
3277 gPsiMinusPhiBin = 3.0;
3279 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
3282 //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
3285 Int_t eventplaneAxis=0;
3286 if(fRequestEventPlane) eventplaneAxis=1;
3287 if(fcontainPIDtrig && SetChargeAxis==0) dim=4+eventplaneAxis;
3288 if(!fcontainPIDtrig && SetChargeAxis==2) dim=4+eventplaneAxis;
3289 if(fcontainPIDtrig && SetChargeAxis==2) dim=5+eventplaneAxis;
3290 trigval= new Double_t[dim];
3293 trigval[2] = trigpt;
3295 if(fRequestEventPlane){
3296 trigval[3] = gPsiMinusPhiBin;
3297 if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = ParticlePID_InvMass;
3298 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
3299 if(fcontainPIDtrig && SetChargeAxis==2) {
3300 trigval[4] = ParticlePID_InvMass;
3301 trigval[5] = trig->Charge();
3305 if(!fRequestEventPlane){
3306 if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = ParticlePID_InvMass;
3307 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
3308 if(fcontainPIDtrig && SetChargeAxis==2) {
3309 trigval[3] = ParticlePID_InvMass;
3310 trigval[4] = trig->Charge();
3316 if (fWeightPerEvent)
3318 // leads effectively to a filling of one entry per filled trigger particle pT bin
3319 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(trigval[2]);
3320 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3321 trackefftrig *= triggerWeighting->GetBinContent(weightBin);
3325 //trigger particle counting for both same and mixed event case;;;;;step=0->same event case;;;;;step=1->mixed event case;;;;;;
3326 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
3327 if(fillup=="trigassoUNID" ) {
3328 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3329 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3332 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
3333 if(fillup=="trigassoUNID" )
3335 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3336 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3339 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
3340 if(fillup=="trigUNIDassoID")
3342 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3343 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3346 //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
3347 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
3348 if(fillup=="trigIDassoID")
3350 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3351 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3354 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
3355 if(fillup=="trigIDassoUNID" )
3357 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3358 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3361 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
3362 if(fillup=="trigIDassoID")
3364 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3365 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3369 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!!!!
3370 if(mixcase==kFALSE) fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig);
3371 if(mixcase==kTRUE && firstTime) fTHnTrigcountMCTruthPrim->Fill(trigval,1,1.0/trackefftrig);
3374 //asso loop starts within trigger loop
3375 for(Int_t j=0;j<jmax;j++)
3377 LRCParticlePID *asso=0;
3379 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
3381 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
3385 //to avoid overflow and underflow
3386 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
3388 //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
3390 if(!tracksasso && j==i) continue;
3392 // 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)
3393 // if (tracksasso && trig->IsEqual(asso)) continue;
3395 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
3398 if (asso->Pt() >= trig->Pt()) continue;
3400 Int_t particlepidasso=asso->getparticle();
3401 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
3404 if (fAssociatedSelectCharge != 0)
3405 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
3407 if (fSelectCharge > 0)
3410 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
3414 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
3420 if (trigeta < 0 && asso->Eta() < trigeta)
3422 if (trigeta > 0 && asso->Eta() > trigeta)
3426 if (fRejectResonanceDaughters > 0)
3427 if (asso->TestBit(kResonanceDaughterFlag))
3429 // Printf("Skipped j=%d", j);
3434 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
3436 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3440 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3442 fControlConvResoncances->Fill(0.0, mass);
3444 if (mass < 0.04*0.04)
3450 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3452 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
3454 const Float_t kK0smass = 0.4976;
3456 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
3458 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
3460 fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
3462 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
3468 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3470 Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
3471 Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3473 const Float_t kLambdaMass = 1.115;
3475 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
3477 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
3479 fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
3481 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3484 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
3486 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3488 fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
3490 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3495 if (twoTrackEfficiencyCut)
3497 // the variables & cuthave been developed by the HBT group
3498 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
3499 Float_t phi1 = trig->Phi();
3500 Float_t pt1 = trig->Pt();
3501 Float_t charge1 = trig->Charge();
3502 Float_t phi2 = asso->Phi();
3503 Float_t pt2 = asso->Pt();
3504 Float_t charge2 = asso->Charge();
3506 Float_t deta= trigeta - eta[j];
3509 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
3512 // check first boundaries to see if is worth to loop and find the minimum
3513 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMinRadius, bSign);
3514 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMaxRadius, bSign);
3516 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
3518 Float_t dphistarminabs = 1e5;
3519 Float_t dphistarmin = 1e5;
3521 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
3523 for (Double_t rad=fTwoTrackCutMinRadius; rad<=fTwoTrackCutMaxRadius; rad+=0.01)
3525 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
3527 Float_t dphistarabs = TMath::Abs(dphistar);
3529 if (dphistarabs < dphistarminabs)
3531 dphistarmin = dphistar;
3532 dphistarminabs = dphistarabs;
3535 if(mixcase==kFALSE) fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3536 if(mixcase==kTRUE) fTwoTrackDistancePtmix[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3538 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
3540 // 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);
3543 if(mixcase==kFALSE) fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3544 if(mixcase==kTRUE) fTwoTrackDistancePtmix[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3550 //pair sharedfraction cut(only between the trig and asso track)
3551 if(fillup!="trigIDassoIDMCTRUTH")//******************************************NOT for TRUTH MC particles
3553 if(fSharedfraction_Pair_cut>=0){
3554 Bool_t passsharedfractionpaircut=CalculateSharedFraction(trig->GetTPCPadMap(),asso->GetTPCPadMap(),trig->GetTPCSharedMap(),asso->GetTPCSharedMap());
3555 if(!passsharedfractionpaircut) continue;
3558 Float_t weightperevent=weight;
3559 Float_t trackeffasso=1.0;
3560 if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
3561 //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
3562 Float_t deleta=trigeta-eta[j];
3563 Float_t delphi=PhiRange(trigphi-asso->Phi());
3565 Float_t delpt=trigpt-asso->Pt();
3566 //fill it with/without two track efficiency cut
3567 if(mixcase==kFALSE) fTwoTrackDistancePtdip->Fill(deleta, delphi, TMath::Abs(delpt));//for same event
3568 if(mixcase==kTRUE) fTwoTrackDistancePtdipmix->Fill(deleta, delphi, TMath::Abs(delpt));//for mixed event
3570 //here get the two particle efficiency correction factor
3571 Float_t effweight=trackefftrig*trackeffasso*weightperevent;
3572 // if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
3574 Int_t dimused=kTrackVariablesPair+eventplaneAxis;
3575 vars= new Double_t[dimused];
3584 if(fRequestEventPlane)
3586 vars[6]=gPsiMinusPhiBin;
3590 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
3591 vars[dimension]=trig->Charge();
3592 vars[dimension+1]=asso->Charge();
3594 if(fcontainPIDtrig && !fcontainPIDasso){
3595 vars[dimension]=ParticlePID_InvMass;
3596 if(SetChargeAxis==2){
3597 vars[dimension+1]=trig->Charge();
3598 vars[dimension+2]=asso->Charge();
3601 if(!fcontainPIDtrig && fcontainPIDasso){
3602 vars[dimension]=particlepidasso;
3603 if(SetChargeAxis==2){
3604 vars[dimension+1]=trig->Charge();
3605 vars[dimension+2]=asso->Charge();
3608 if(fcontainPIDtrig && fcontainPIDasso){
3609 vars[dimension]=ParticlePID_InvMass;
3610 vars[dimension+1]=particlepidasso;
3611 if(SetChargeAxis==2){
3612 vars[dimension+2]=trig->Charge();
3613 vars[dimension+3]=asso->Charge();
3617 if (fWeightPerEvent)
3619 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
3620 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3621 effweight *= triggerWeighting->GetBinContent(weightBin);
3625 //Fill Correlation ThnSparses
3626 if(fillup=="trigassoUNID")
3628 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,0,1.0/effweight);
3629 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
3631 if(fillup=="trigIDassoID")
3633 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,0,1.0/effweight);
3634 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
3636 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
3637 {//in this case effweight should be 1 always
3638 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
3639 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
3641 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
3643 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
3644 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
3650 }//trigger loop ends
3652 if (triggerWeighting)
3654 delete triggerWeighting;
3655 triggerWeighting = 0;
3659 //------------------------------------------------------------------------------------------------
3660 Bool_t AliTwoParticlePIDCorr:: CalculateSharedFraction(const TBits *triggerPadMap,const TBits *assocPadMap,const TBits *triggerShareMap,const TBits *assocShareMap)
3661 {//source code-AliFemtoShareQualityPairCut.cxx
3663 Double_t nofsharedhits=0;
3665 for(UInt_t imap=0;imap< (triggerPadMap->GetNbits() );imap++)
3667 //if they are in same pad
3668 //cout<<triggerPadMap->TestBitNumber(imap)<<" "<< assocPadMap->TestBitNumber(imap)<<endl;
3669 if (triggerPadMap->TestBitNumber(imap) &&
3670 assocPadMap->TestBitNumber(imap))
3673 //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
3674 if (triggerShareMap->TestBitNumber(imap) &&
3675 assocShareMap->TestBitNumber(imap))
3677 //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
3694 //cout<< (triggerPadMap->TestBitNumber(imap) || assocPadMap->TestBitNumber(imap))<<endl;
3695 else if (triggerPadMap->TestBitNumber(imap) ||
3696 assocPadMap->TestBitNumber(imap)) {
3697 // One track has a hit, the other does not
3700 //cout<<"No hits :"<<nofhits<<endl;
3708 Double_t SharedFraction=0.0;
3709 if(nofhits>0) SharedFraction=(nofsharedhits/nofhits);
3711 //cout<<"Fraction shared hits :"<<SharedFraction<<endl;
3713 if(SharedFraction>fSharedfraction_Pair_cut) return kFALSE;
3719 //________________________________________________________________________________________________
3720 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
3722 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
3724 Float_t effvalue=1.;
3726 if(parpid==unidentified)
3728 effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
3729 effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
3730 effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
3731 effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
3732 effvalue=effcorection[5]->GetBinContent(effVars);
3734 if(parpid==SpPion || parpid==SpKaon)
3736 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3738 effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
3739 effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
3740 effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
3741 effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
3742 effvalue=effcorection[3]->GetBinContent(effVars);
3747 effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
3748 effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
3749 effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
3750 effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
3751 effvalue=effcorection[0]->GetBinContent(effVars);
3756 effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
3757 effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
3758 effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
3759 effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
3760 effvalue=effcorection[1]->GetBinContent(effVars);
3765 if(parpid==SpProton)
3767 effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
3768 effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
3769 effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
3770 effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
3771 effvalue=effcorection[2]->GetBinContent(effVars);
3774 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3776 if(parpid==SpProton || parpid==SpKaon)
3778 effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
3779 effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
3780 effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
3781 effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
3782 effvalue=effcorection[4]->GetBinContent(effVars);
3785 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
3786 if(effvalue==0.) effvalue=1.;
3791 //---------------------------------------------------------------------------------
3793 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield, Bool_t fill)
3796 if(!track) return 0;
3797 Bool_t trackOK = track->TestFilterBit(fFilterBit);
3798 if(!trackOK) return 0;
3799 if (fTrackStatus != 0 && !CheckTrack(track)) return 0;
3800 //select only primary traks(for data & reco MC tracks)
3801 if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
3802 if(track->Charge()==0) return 0;
3803 if (fill) fHistQA[12]->Fill(track->GetTPCNcls());
3806 dz = track->ZAtDCA();
3807 if (fill) fHistQA[6]->Fill(dxy);
3808 if (fill) fHistQA[7]->Fill(dz);
3809 Float_t chi2ndf = track->Chi2perNDF();
3810 if (fill) fHistQA[13]->Fill(chi2ndf);
3811 // Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
3812 Float_t nCrossedRowsTPC = track->GetTPCNCrossedRows();
3813 if (fill) fHistQA[14]->Fill(nCrossedRowsTPC);
3814 //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
3815 if (track->GetTPCNclsF()>0) {
3816 Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
3817 if (fill) fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
3820 Float_t pt=track->Pt();
3821 if(pt< fminPt || pt> fmaxPt) return 0;
3822 if(TMath::Abs(track->Eta())> fmaxeta) return 0;
3823 if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
3824 //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required???
3826 if (fdcacut && fDCAXYCut)
3833 AliAODTrack* clone =(AliAODTrack*) track->Clone();
3834 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
3839 // Printf("%f", ((AliAODTrack*)part)->DCA());
3840 // Printf("%f", pos[0]);
3841 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
3845 if (fSharedClusterCut >= 0)
3847 Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
3848 if (frac > fSharedClusterCut)
3852 // Rejects tracks with shared clusters after filling a control histogram
3853 // This overload is used for primaries
3855 // Get the shared maps
3856 const TBits sharedMap = track->GetTPCSharedMap();
3857 // Fill a control histogram
3858 fPriHistShare->Fill(sharedMap.CountBits());
3860 // Reject shared clusters
3861 if (fSharedTPCmapCut >= 0)
3863 if((sharedMap.CountBits()) >= 1) return 0;// Bad track, has too many shared clusters!
3866 if (fill) fHistQA[8]->Fill(pt);
3867 if (fill) fHistQA[9]->Fill(track->Eta());
3868 if (fill) fHistQA[10]->Fill(track->Phi());
3871 //________________________________________________________________________________
3872 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
3874 //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 )
3875 Float_t pt=track->Pt();
3877 //plot the separation power
3879 Double_t bethe[AliPID::kSPECIES]={0.};
3880 Double_t sigma_TPC[AliPID::kSPECIES]={0.};
3882 Double_t Pi_Ka_sep[NSigmaPIDType+1]={0.};
3883 Double_t Pi_Pr_sep[NSigmaPIDType+1]={0.};
3884 Double_t Ka_Pr_sep[NSigmaPIDType+1]={0.};
3887 Double_t ptpc = track->GetTPCmomentum();
3888 Int_t dEdxN = track->GetTPCsignalN();
3889 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
3890 bethe[ipart] = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
3891 //Double_t diff = dEdx - bethe;
3892 sigma_TPC[ipart] = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
3893 //nSigma[ipart] = diff / sigma;
3895 Pi_Ka_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kKaon])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kKaon])/2.0);
3896 Pi_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kProton])/2.0);
3897 Ka_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kKaon]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kKaon]+sigma_TPC[AliPID::kProton])/2.0);
3900 Double_t sigma_TOF[AliPID::kSPECIES]={0.};
3902 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3904 Double_t timei[AliPID::kSPECIES];
3905 track->GetIntegratedTimes(timei);
3906 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { sigma_TOF[ipart]= fPID->GetTOFResponse().GetExpectedSigma(track->P(), timei[ipart], AliPID::ParticleMass(ipart));}
3907 Pi_Ka_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kKaon])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kKaon])/2.0);
3908 Pi_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kProton])/2.0);
3909 Ka_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kKaon]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kKaon]+sigma_TOF[AliPID::kProton])/2.0);
3911 Pi_Ka_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Ka_sep[NSigmaTPC]*Pi_Ka_sep[NSigmaTPC]+Pi_Ka_sep[NSigmaTOF]*Pi_Ka_sep[NSigmaTOF]);
3912 Pi_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Pr_sep[NSigmaTPC]*Pi_Pr_sep[NSigmaTPC]+Pi_Pr_sep[NSigmaTOF]*Pi_Pr_sep[NSigmaTOF]);
3913 Ka_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Ka_Pr_sep[NSigmaTPC]*Ka_Pr_sep[NSigmaTPC]+Ka_Pr_sep[NSigmaTOF]*Ka_Pr_sep[NSigmaTOF]);
3917 //fill separation power histograms
3918 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3920 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3921 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3922 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3923 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3924 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3925 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3927 if(HasTOFPID(track) && pt>fPtTOFPIDmin && ipid!=0){
3928 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3929 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3930 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3931 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3932 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3933 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3938 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
3939 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
3940 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
3941 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
3943 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
3944 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
3946 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3949 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
3950 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
3951 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
3953 nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
3954 nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
3955 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
3961 // if TOF is missing and below fPtTOFPID only the TPC information is used
3962 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
3963 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
3964 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
3968 //set data member fnsigmas
3969 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
3970 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
3971 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
3973 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
3974 fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
3975 fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
3976 fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
3978 //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)
3979 fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
3980 fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
3981 fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
3984 //Fill NSigma SeparationPlot
3985 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3986 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3987 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3988 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
3989 h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
3995 //----------------------------------------------------------------------------
3996 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
3998 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3999 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
4000 //get the identity of the particle with the minimum Nsigma
4001 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
4004 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
4005 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
4006 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
4009 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
4010 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
4011 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
4013 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
4014 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4015 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4016 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4018 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
4019 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4020 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4021 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4026 if(fdiffPIDcutvalues){
4027 if(track->Pt()<=4) fNSigmaPID=fPIDCutval1;
4028 if(track->Pt()>4 && track->Pt()<=6) fNSigmaPID=fPIDCutval2;
4029 if(track->Pt()>6 && track->Pt()<=8) fNSigmaPID=fPIDCutval3;
4030 if(track->Pt()>8) fNSigmaPID=fPIDCutval4;
4033 // guess the particle based on the smaller nsigma (within fNSigmaPID)
4034 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
4036 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
4037 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)
4039 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4040 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4041 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
4042 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
4047 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
4049 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4050 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4051 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
4052 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
4057 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
4059 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4060 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4061 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
4062 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
4068 // else, return undefined
4074 //------------------------------------------------------------------------------------------
4075 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
4076 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
4078 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
4080 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
4082 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
4085 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
4087 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
4090 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
4091 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
4092 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
4095 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
4096 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
4097 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
4099 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
4100 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4101 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4102 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4104 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
4105 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4106 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4107 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4111 // Actually the tracks in the overlapping region(in TPC-TOF nSigma plane) will be ignored
4113 if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
4114 if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
4115 if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
4120 //fill NSigma distr for double counting
4121 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4122 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
4123 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4124 if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)))continue;//not filling TOF and combined if no TOF PID
4125 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
4126 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
4133 return fHasDoubleCounting;
4136 //////////////////////////////////////////////////////////////////////////////////////////////////
4138 Bool_t* AliTwoParticlePIDCorr::GetAllCompatibleIdentitiesNSigma(AliAODTrack * trk,Bool_t FIllQAHistos){
4139 //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
4140 Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
4141 IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
4145 //////////////////////////////////////////////////////////////////////////////////////////////////
4147 UInt_t AliTwoParticlePIDCorr::CalcPIDCombined(AliAODTrack *track, Int_t detMask, Double_t* prob) const{
4149 // Bayesian PID calculation
4151 for(Int_t i=0;i<AliPID::kSPECIES;i++)
4155 fPIDCombined->SetDetectorMask(detMask);
4157 return fPIDCombined->ComputeProbabilities((AliAODTrack*)track, fPID, prob);
4160 //////////////////////////////////////////////////////////////////////////////////////////////////
4162 Int_t AliTwoParticlePIDCorr::GetIDBayes(AliAODTrack * trk, Bool_t FIllQAHistos){
4164 Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
4167 //Filling of Probability histos
4168 Double_t probTPC[AliPID::kSPECIES]={0.};
4169 Double_t probTOF[AliPID::kSPECIES]={0.};
4170 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
4172 UInt_t detUsedTPC = 0;
4173 UInt_t detUsedTOF = 0;
4174 UInt_t detUsedTPCTOF = 0;
4176 //get the TPC probability
4177 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
4178 detUsedTPC = fPIDCombined->ComputeProbabilities(trk, fPID, probTPC);
4179 if(detUsedTPC == AliPIDResponse::kDetTPC)
4181 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4183 TH2F *h=GetHistogram2D(Form("probBayes_TPC_%d",ipart));
4184 if(ipart==0) h->Fill(trk->Pt(),probTPC[AliPID::kPion]);
4185 if(ipart==1) h->Fill(trk->Pt(),probTPC[AliPID::kKaon]);
4186 if(ipart==2) h->Fill(trk->Pt(),probTPC[AliPID::kProton]);
4190 //get the TOF probability
4191 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
4192 detUsedTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTOF);
4193 if(detUsedTOF == AliPIDResponse::kDetTOF)
4195 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4196 TH2F *h=GetHistogram2D(Form("probBayes_TOF_%d",ipart));
4197 if(ipart==0) h->Fill(trk->Pt(),probTOF[AliPID::kPion]);
4198 if(ipart==1) h->Fill(trk->Pt(),probTOF[AliPID::kKaon]);
4199 if(ipart==2) h->Fill(trk->Pt(),probTOF[AliPID::kProton]);
4203 //get the TPC-TOF probability
4204 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
4205 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTPCTOF);
4206 if(detUsedTPCTOF == (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))
4208 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4209 TH2F *h=GetHistogram2D(Form("probBayes_TPCTOF_%d",ipart));
4210 if(ipart==0) h->Fill(trk->Pt(),probTPCTOF[AliPID::kPion]);
4211 if(ipart==1) h->Fill(trk->Pt(),probTPCTOF[AliPID::kKaon]);
4212 if(ipart==2) h->Fill(trk->Pt(),probTPCTOF[AliPID::kProton]);
4217 Double_t probBayes[AliPID::kSPECIES];
4220 if(HasTOFPID(trk) && trk->Pt()>fPtTOFPIDmin){//use TOF information
4221 detUsed = CalcPIDCombined(trk, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
4222 if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return SpUndefined;//check that TPC and TOF are used
4224 detUsed = CalcPIDCombined(trk,AliPIDResponse::kDetTPC, probBayes);
4225 if(detUsed != AliPIDResponse::kDetTPC)return SpUndefined;//check that TPC is used
4228 //the probability has to be normalized to one, we check it
4230 for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
4231 if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
4232 AliFatal("Bayesian probability not normalized to one");
4235 if(fdiffPIDcutvalues){
4236 if(trk->Pt()<=4) fBayesCut=fPIDCutval1;
4237 if(trk->Pt()>4 && trk->Pt()<=6) fBayesCut=fPIDCutval2;
4238 if(trk->Pt()>6 && trk->Pt()<=8) fBayesCut=fPIDCutval3;
4239 if(trk->Pt()>8) fBayesCut=fPIDCutval4;
4243 //probabilities are normalized to one, if the cut is above .5 there is no problem
4244 if(probBayes[AliPID::kPion]>fBayesCut && IDs[SpPion]==1){
4245 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpPion));
4246 h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
4249 else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[SpKaon]==1){
4250 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpKaon));
4251 h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
4254 else if(probBayes[AliPID::kProton]>fBayesCut && IDs[SpProton]==1){
4255 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpProton));
4256 h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
4265 //////////////////////////////////////////////////////////////////////////////////////////////////
4266 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
4267 //return the specie according to the minimum nsigma value
4268 //no double counting, this has to be evaluated using CheckDoubleCounting()
4270 Int_t ID=SpUndefined;
4272 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
4276 if(fPIDType==Bayes){//use bayesianPID
4279 AliFatal("PIDCombined object has to be set in the steering macro");
4282 ID = GetIDBayes(trk,FIllQAHistos);
4284 }else{ //use nsigma PID
4286 ID=FindMinNSigma(trk,FIllQAHistos);
4287 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
4289 HasDC=GetDoubleCounting(trk,FIllQAHistos);
4290 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4291 if(HasDC[ipart]==kTRUE) ID = SpUndefined;
4295 //Fill PID signal plot
4296 if(ID != SpUndefined){
4297 for(Int_t idet=0;idet<fNDetectors;idet++){
4298 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
4299 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
4300 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
4301 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
4304 //Fill PID signal plot without cuts
4305 for(Int_t idet=0;idet<fNDetectors;idet++){
4306 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
4307 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
4308 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
4309 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
4315 //-------------------------------------------------------------------------------------
4317 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
4320 AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
4321 if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
4322 //ULong_t status=track->GetStatus();
4323 //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
4324 //if (track->GetTPCsignal() <= 0.) return kFALSE;
4325 // 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.
4329 //___________________________________________________________
4332 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
4334 // check TOF matched track
4335 //ULong_t status=track->GetStatus();
4336 //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
4337 AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
4338 if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
4339 if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
4340 //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
4341 //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
4342 // if (probMis > 0.01) return kFALSE;
4343 if(fRemoveTracksT0Fill)
4345 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
4346 if (startTimeMask < 0)return kFALSE;
4351 //________________________________________________________________________
4352 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
4354 //it is called only when TOF PID is available
4355 //TOF beta calculation
4356 Double_t tofTime=track->GetTOFsignal();
4358 Double_t c=TMath::C()*1.E-9;// m/ns
4359 Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
4360 Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
4361 tofTime -= startTime; // subtract startTime to the signal
4362 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
4368 Double_t p = track->P();
4369 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
4371 track->GetIntegratedTimes(timei);
4372 return timei[0]/time;
4375 //------------------------------------------------------------------------------------------------------
4377 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)
4379 // calculate inv mass squared
4380 // same can be achieved, but with more computing time with
4381 /*TLorentzVector photon, p1, p2;
4382 p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
4383 p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
4387 Float_t tantheta1 = 1e10;
4389 if (eta1 < -1e-10 || eta1 > 1e-10)
4390 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
4392 Float_t tantheta2 = 1e10;
4393 if (eta2 < -1e-10 || eta2 > 1e-10)
4394 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
4396 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4397 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4399 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 ) ) );
4403 //---------------------------------------------------------------------------------
4405 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)
4407 // calculate inv mass squared approximately
4409 Float_t tantheta1 = 1e10;
4411 if (eta1 < -1e-10 || eta1 > 1e-10)
4413 Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
4414 tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4417 Float_t tantheta2 = 1e10;
4418 if (eta2 < -1e-10 || eta2 > 1e-10)
4420 Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
4421 tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4424 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4425 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4428 Float_t deltaPhi = TMath::Abs(phi1 - phi2);
4429 while (deltaPhi > TMath::TwoPi())
4430 deltaPhi -= TMath::TwoPi();
4431 if (deltaPhi > TMath::Pi())
4432 deltaPhi = TMath::TwoPi() - deltaPhi;
4434 Float_t cosDeltaPhi = 0;
4435 if (deltaPhi < TMath::Pi()/3)
4436 cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
4437 else if (deltaPhi < 2*TMath::Pi()/3)
4438 cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
4440 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
4442 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
4444 // Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
4448 //--------------------------------------------------------------------------------
4449 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)
4452 // calculates dphistar
4455 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
4457 static const Double_t kPi = TMath::Pi();
4460 // if (dphistar > 2 * kPi)
4461 // dphistar -= 2 * kPi;
4462 // if (dphistar < -2 * kPi)
4463 // dphistar += 2 * kPi;
4466 dphistar = kPi * 2 - dphistar;
4467 if (dphistar < -kPi)
4468 dphistar = -kPi * 2 - dphistar;
4469 if (dphistar > kPi) // might look funny but is needed
4470 dphistar = kPi * 2 - dphistar;
4475 //------------------------------------------------------------------------
4476 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
4478 // This method is a copy from AliUEHist::GetBinning
4479 // takes the binning from <configuration> identified by <tag>
4480 // configuration syntax example:
4481 // 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
4484 // returns bin edges which have to be deleted by the caller
4486 TString config(configuration);
4487 TObjArray* lines = config.Tokenize("\n");
4488 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
4490 TString line(lines->At(i)->GetName());
4491 if (line.BeginsWith(TString(tag) + ":"))
4493 line.Remove(0, strlen(tag) + 1);
4494 line.ReplaceAll(" ", "");
4495 TObjArray* binning = line.Tokenize(",");
4496 Double_t* bins = new Double_t[binning->GetEntriesFast()];
4497 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
4498 bins[j] = TString(binning->At(j)->GetName()).Atof();
4500 nBins = binning->GetEntriesFast() - 1;
4509 AliFatal(Form("Tag %s not found in %s", tag, configuration));
4513 //____________________________________________________________________
4515 Bool_t AliTwoParticlePIDCorr::CheckTrack(AliAODTrack * part)
4517 // check if the track status flags are set
4519 UInt_t status=((AliAODTrack*)part)->GetStatus();
4520 if ((status & fTrackStatus) == fTrackStatus)
4524 //________________________________________________________________________
4526 Bool_t AliTwoParticlePIDCorr::AcceptEventCentralityWeight(Double_t centrality)
4528 // rejects "randomly" events such that the centrality gets flat
4529 // uses fCentralityWeights histogram
4531 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
4533 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
4534 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
4536 Bool_t result = kFALSE;
4537 if (centralityDigits < weight)
4540 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
4545 //____________________________________________________________________
4546 void AliTwoParticlePIDCorr::ShiftTracks(TObjArray* tracks, Double_t angle)
4548 // shifts the phi angle of all tracks by angle
4549 // 0 <= angle <= 2pi
4551 for (Int_t i=0; i<tracks->GetEntriesFast(); ++i)
4553 LRCParticlePID *part=(LRCParticlePID*)(tracks->UncheckedAt(i));
4555 Double_t newAngle = part->Phi() + angle;
4556 if (newAngle >= TMath::TwoPi())
4557 newAngle -= TMath::TwoPi();
4559 part->SetPhi(newAngle);
4564 //________________________________________________________________________
4565 void AliTwoParticlePIDCorr::SetVZEROCalibrationFile(const char* filename,const char* lhcPeriod) {
4566 //Function to setup the VZERO gain equalization
4567 //============Get the equilization map============//
4568 TFile *calibrationFile = TFile::Open(filename);
4569 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
4570 Printf("No calibration file found!!!");
4574 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
4576 Printf("Calibration TList not found!!!");
4580 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
4581 if(!fHistVZEROAGainEqualizationMap) {
4582 Printf("VZERO-A calibration object not found!!!");
4585 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
4586 if(!fHistVZEROCGainEqualizationMap) {
4587 Printf("VZERO-C calibration object not found!!!");
4591 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
4592 if(!fHistVZEROChannelGainEqualizationMap) {
4593 Printf("VZERO channel calibration object not found!!!");
4598 //________________________________________________________________________
4599 Double_t AliTwoParticlePIDCorr::GetChannelEqualizationFactor(Int_t run,Int_t channel) {
4601 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4603 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
4604 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4605 if(gRunNumber == run)
4606 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
4612 //________________________________________________________________________
4613 Double_t AliTwoParticlePIDCorr::GetEqualizationFactor(Int_t run, const char* side) {
4615 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4617 TString gVZEROSide = side;
4618 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
4619 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4620 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
4621 if(gRunNumber == run) {
4622 if(gVZEROSide == "A")
4623 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
4624 else if(gVZEROSide == "C")
4625 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
4631 //________________________________________________________________________
4632 Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliVEvent *mainevent){
4633 //Function that returns the reference multiplicity from AODs (data or reco MC, Not for Truth)
4634 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
4635 if(!mainevent) return -1;
4637 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
4639 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
4640 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
4642 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
4644 Printf("ERROR: AOD header not available");
4647 Int_t gRunNumber = header->GetRunNumber();
4648 Float_t bSign1=header->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
4651 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++)
4652 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
4653 AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
4654 if (!track) continue;
4655 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
4656 if(tracktype!=1) continue;
4658 if(!track) continue;//for safety
4660 gRefMultiplicityTPC += 1.0;
4664 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
4665 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
4666 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
4667 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
4670 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
4671 else if(iChannel >= 32)
4672 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
4675 //Equalization of gain
4676 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
4678 gRefMultiplicityVZEROA /= gFactorA;
4679 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
4681 gRefMultiplicityVZEROC /= gFactorC;
4682 if((gFactorA != 0)&&(gFactorC != 0))
4683 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
4686 //EQVZERO vs TPC multiplicity
4687 fHistEQVZEROvsTPCmultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
4688 fHistEQVZEROAvsTPCmultiplicity->Fill(gRefMultiplicityVZEROA,gRefMultiplicityTPC);
4689 fHistEQVZEROCvsTPCmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityTPC);
4691 //EQVZERO vs VZERO multiplicity
4692 fHistVZEROCvsEQVZEROCmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),gRefMultiplicityVZEROC);
4693 fHistVZEROAvsEQVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0A(),gRefMultiplicityVZEROA);
4695 //VZEROC vs VZEROA multiplicity
4696 fHistVZEROCvsVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),event->GetVZEROData()->GetMTotV0A());
4698 //EQVZEROC vs EQVZEROA multiplicity
4699 fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
4701 fHistRefmult->Fill(3.,gRefMultiplicityTPC);
4702 fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
4703 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
4704 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
4707 if(fCentralityMethod == "TRACKS_MANUAL") gRefMultiplicity = gRefMultiplicityTPC;
4709 else if(fCentralityMethod == "V0M_MANUAL") gRefMultiplicity = gRefMultiplicityVZERO;
4711 else if(fCentralityMethod == "V0A_MANUAL") gRefMultiplicity = gRefMultiplicityVZEROA;
4713 else if(fCentralityMethod == "V0C_MANUAL") gRefMultiplicity = gRefMultiplicityVZEROC;
4715 else gRefMultiplicity = gRefMultiplicityTPC;
4717 return gRefMultiplicity;
4720 //-------------------------------------------------------------------------------------------------------
4721 Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliVEvent *mainevent, Bool_t truth){
4723 if(!mainevent) return -1;
4724 // get centrality object and check quality
4725 Double_t cent_v0=-1;
4726 Bool_t shift_to_TRACKS_MANUAL=kFALSE;//in case of wrong setting automatic shift to Tracks_Manual method
4728 Double_t gRefMultiplicityTPC_Truth = 0.;
4729 Double_t gRefMultiplicityVZERO_Truth = 0., gRefMultiplicityVZEROA_Truth = 0., gRefMultiplicityVZEROC_Truth = 0.;
4731 if(fAnalysisType == "AOD"|| fAnalysisType == "MCAOD") { //centrality in AOD header //++++++++++++++
4732 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
4734 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
4737 if(fSampleType=="pp_7" && fPPVsMultUtils)
4738 {//for pp 7 TeV case only using Alianalysisutils class
4739 if(fAnalysisUtils) cent_v0 = fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,fCentralityMethod);
4741 fHistCentStats->Fill(0.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0A"));
4742 fHistCentStats->Fill(1.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0C"));
4743 fHistCentStats->Fill(2.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0M"));
4744 fHistCentStats->Fill(3.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0AEq"));//only available for LHC10d at present (Quantile info)
4745 fHistCentStats->Fill(4.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0CEq"));//only available for LHC10d at present (Quantile info)
4746 fHistCentStats->Fill(5.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0MEq"));//only available for LHC10d at present (Quantile info)
4749 else if(fSampleType=="pPb" || fSampleType=="PbPb")
4751 AliCentrality *centralityObj=0;
4752 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
4753 if(!header) return -1;
4754 centralityObj = header->GetCentralityP();
4755 // if (centrality->GetQuality() != 0) return ;
4757 fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
4758 fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
4759 fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
4760 fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
4761 fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
4762 fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
4764 fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
4765 fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
4767 cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
4771 else shift_to_TRACKS_MANUAL=kTRUE;
4773 }//centralitymethod condition
4775 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
4777 if(!truth){//for data or RecoMC
4778 cent_v0 = GetReferenceMultiplicityVZEROFromAOD((AliVEvent*)event);
4779 }//for data or RecoMC
4781 if(truth && (fAnalysisType == "MCAOD")){//condition for TRUTH case
4782 //check for TClonesArray(truth track MC information)
4783 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
4785 //AliFatal("Error: MC particles branch not found!\n");
4788 //now process the truth particles(for both efficiency & correlation function)
4789 Int_t nMCTrack = fArrayMC->GetEntriesFast();
4791 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
4792 {//MC truth track loop starts
4794 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
4797 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
4801 //consider only charged particles
4802 if(partMC->Charge() == 0) continue;
4804 //consider only primary particles; neglect all secondary particles including from weak decays
4805 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
4808 //remove injected signals(primaries above <maxLabel>)
4809 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
4812 Bool_t isduplicate=kFALSE;
4813 if (fRemoveDuplicates)
4815 for (Int_t j=iMC+1; j<nMCTrack; ++j)
4816 {//2nd trutuh loop starts
4817 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
4819 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
4822 if (partMC->GetLabel() == partMC2->GetLabel())
4827 }//2nd truth loop ends
4829 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
4832 // if (fCentralityMethod=="V0M_MANUAL")
4833 if((partMC->Eta() < 5.1 && partMC->Eta() > 2.8) || (partMC->Eta() > -3.7 && partMC->Eta() < -1.7)) gRefMultiplicityVZERO_Truth+=1;
4834 // else if (fCentralityMethod=="V0A_MANUAL") {
4835 if(partMC->Eta() < 5.1 && partMC->Eta() > 2.8) gRefMultiplicityVZEROA_Truth+=1;
4836 // else if (fCentralityMethod=="V0C_MANUAL") {
4837 if(partMC->Eta() < -1.7 && partMC->Eta() > -3.7) gRefMultiplicityVZEROC_Truth+=1;
4838 //else if (fCentralityMethod=="TRACKS_MANUAL") {
4839 if (partMC->Eta() > fmineta && partMC->Eta() < fmaxeta) {
4840 if (partMC->Pt() > fminPt && partMC->Pt() < fmaxPt) gRefMultiplicityTPC_Truth+=1;
4843 }//truth track loop ends
4845 fHistRefmult->Fill(3.,gRefMultiplicityTPC_Truth);
4846 fHistRefmult->Fill(2.,gRefMultiplicityVZERO_Truth);
4847 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA_Truth);
4848 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC_Truth);
4850 if(fCentralityMethod == "TRACKS_MANUAL") cent_v0=gRefMultiplicityTPC_Truth;
4852 else if(fCentralityMethod == "V0M_MANUAL") cent_v0=gRefMultiplicityVZERO_Truth;
4854 else if(fCentralityMethod == "V0A_MANUAL") cent_v0=gRefMultiplicityVZEROA_Truth;
4856 else if(fCentralityMethod == "V0C_MANUAL") cent_v0=gRefMultiplicityVZEROC_Truth;
4858 else cent_v0=gRefMultiplicityTPC_Truth;
4860 }//condition for TRUTH case
4862 }//end of MANUAL method
4864 else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC in AOD production
4866 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4870 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4873 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
4874 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
4875 AliError("Event header not found. Skipping this event.");
4879 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4883 cent_v0 = collGeometry->ImpactParameter();
4884 fhistImpactParm->Fill(cent_v0);
4887 }//end of Impact parameter method
4891 }//AOD OR MCAOD condition
4894 else if(fAnalysisType == "MC"){
4895 Double_t gImpactParameter = -1.;
4896 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(mainevent);
4898 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
4900 gImpactParameter = headerH->ImpactParameter();
4902 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
4903 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
4905 AliError(Form("Could not receive particle %d", iParticle));
4909 //exclude non stable particles
4910 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
4912 if(track->Charge() == 0) continue;
4914 // if (fCentralityMethod=="V0M_MANUAL")
4915 if((track->Eta() < 5.1 && track->Eta() > 2.8) || (track->Eta() > -3.7 && track->Eta() < -1.7)) gRefMultiplicityVZERO_Truth+=1;
4916 // else if (fCentralityMethod=="V0A_MANUAL") {
4917 if(track->Eta() < 5.1 && track->Eta() > 2.8) gRefMultiplicityVZEROA_Truth+=1;
4918 // else if (fCentralityMethod=="V0C_MANUAL") {
4919 if(track->Eta() < -1.7 && track->Eta() > -3.7) gRefMultiplicityVZEROC_Truth+=1;
4920 //else if (fCentralityMethod=="TRACKS_MANUAL") {
4921 if (track->Eta() > fmineta && track->Eta() < fmaxeta) {
4922 if (track->Pt() > fminPt && track->Pt() < fmaxPt) gRefMultiplicityTPC_Truth+=1;}
4924 }//loop over primaries
4926 fHistRefmult->Fill(3.,gRefMultiplicityTPC_Truth);
4927 fHistRefmult->Fill(2.,gRefMultiplicityVZERO_Truth);
4928 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA_Truth);
4929 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC_Truth);
4930 if (fCentralityMethod == "MC_b"){
4931 cent_v0=gImpactParameter;
4932 fhistImpactParm->Fill(gImpactParameter);
4935 else if(fCentralityMethod == "TRACKS_MANUAL") cent_v0=gRefMultiplicityTPC_Truth;
4937 else if(fCentralityMethod == "V0M_MANUAL") cent_v0=gRefMultiplicityVZERO_Truth;
4939 else if(fCentralityMethod == "V0A_MANUAL") cent_v0=gRefMultiplicityVZEROA_Truth;
4941 else if(fCentralityMethod == "V0C_MANUAL") cent_v0=gRefMultiplicityVZEROC_Truth;
4943 else cent_v0=gImpactParameter;//default value is the impact parameter
4954 //-----------------------------------------------------------------------------------------
4955 Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliVEvent *event,Bool_t truth){
4956 //do the event selection(zvtx, pileup, centrality/multiplicity cut) and then return the value of the centrality of that event
4957 if(!event) return -1;
4959 Float_t gRefMultiplicity = -1.;
4961 //***********************************SOURCE CODE-TASKBFPsi
4964 if(fAnalysisType == "MC") {
4965 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
4967 AliError("mcEvent not available");
4972 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
4974 TArrayF gVertexArray;
4975 header->PrimaryVertex(gVertexArray);
4976 //count events having a proper vertex
4977 fEventCounter->Fill(5);
4979 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
4981 if(TMath::Abs(gVertexArray.At(0)) < fVxMax_MC) {
4982 if(TMath::Abs(gVertexArray.At(1)) < fVyMax_MC) {
4983 if(TMath::Abs(gVertexArray.At(2)) < fVzMax_MC) {
4984 //count events after vertex cut
4985 fEventCounter->Fill(7);
4986 fHistQA[3]->Fill((gVertexArray.At(0)));fHistQA[4]->Fill((gVertexArray.At(1)));fHistQA[5]->Fill((gVertexArray.At(2)));//after vertex cut,for trkVtx only
4988 // get the reference multiplicty or centrality
4989 gRefMultiplicity = GetRefMultiOrCentrality((AliVEvent*)mcevent,kFALSE);//2nd argument has no meaning
4991 if(gRefMultiplicity<0) return -1;
4993 // take events only within the multiplicity class mentioned in the custom binning
4994 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
4996 //count events having proper centrality/ref multiplicity
4997 fEventCounter->Fill(9);
5006 else if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD"){// if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD"
5007 //vertex selection(is it fine for PP?)
5008 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
5011 // check first event in chunk (is not needed for new reconstructions)
5012 if(fCheckFirstEventInChunk){
5013 AliAnalysisUtils ut;
5014 if(ut.IsFirstEventInChunk(aod))
5019 AliAnalysisUtils ut;
5020 ut.SetUseMVPlpSelection(kTRUE);
5021 ut.SetUseOutOfBunchPileUp(kTRUE);
5022 if(ut.IsPileUpEvent(aod))
5026 //count events after pileup selection
5027 fEventCounter->Fill(3);
5029 if (fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
5030 trkVtx = aod->GetPrimaryVertex();
5031 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
5032 TString vtxTtl = trkVtx->GetTitle();
5033 if (!vtxTtl.Contains("VertexerTracks")) return -1;
5034 zvtx = trkVtx->GetZ();
5035 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
5036 if (!spdVtx || spdVtx->GetNContributors()<=0) return -1;
5037 TString vtxTyp = spdVtx->GetTitle();
5038 Double_t cov[6]={0};
5039 spdVtx->GetCovarianceMatrix(cov);
5040 Double_t zRes = TMath::Sqrt(cov[5]);
5041 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return -1;
5042 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return -1;
5044 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
5045 Int_t nVertex = aod->GetNumberOfVertices();
5047 trkVtx = aod->GetPrimaryVertex();
5048 Int_t nTracksPrim = trkVtx->GetNContributors();
5049 zvtx = trkVtx->GetZ();
5050 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
5051 // Reject TPC only vertex
5052 TString name(trkVtx->GetName());
5053 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return -1;
5055 // Select a quality vertex by number of tracks?
5056 if( nTracksPrim < fnTracksVertex ) {
5057 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
5060 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
5061 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
5063 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
5068 else if(fVertextype==0){//default case
5069 trkVtx =(AliAODVertex*) aod->GetPrimaryVertex();
5070 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
5071 zvtx = trkVtx->GetZ();
5073 trkVtx->GetCovarianceMatrix(fCov);
5074 if(fCov[5] == 0) return -1;//proper vertex resolution
5077 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
5078 return -1;//as there is no proper sample type
5081 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
5083 //count events having a proper vertex
5084 fEventCounter->Fill(5);
5086 if (TMath::Abs(zvtx) > fzvtxcut) return -1;
5088 //count events after vertex cut
5089 fEventCounter->Fill(7);
5092 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
5094 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
5096 //get the centrality or multiplicity
5097 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)
5099 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)
5101 if(gRefMultiplicity<0) return -1;
5103 // take events only within the multiplicity class mentioned in the custom binning
5104 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
5106 //count events having proper centrality/ref multiplicity
5107 fEventCounter->Fill(9);
5110 // centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
5111 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
5113 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
5117 //count events after rejection due to centrality weighting
5118 fEventCounter->Fill(11);
5121 else gRefMultiplicity=-1;
5123 return gRefMultiplicity;
5126 //--------------------------------------------------------------------------------------------------------
5127 Float_t AliTwoParticlePIDCorr::GetEventPlane(AliVEvent *mainevent,Bool_t truth, Double_t v0Centr)
5129 Float_t eventplane=999.;
5130 // Get the event plane
5131 if(!mainevent) return 999.;
5134 //MC: from reaction plane
5135 if(fAnalysisType == "MC"){
5137 AliError("mcEvent not available");
5141 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(mainevent);
5143 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
5146 // Impact parameter bins(it is only for Pb-Pb)
5147 if(v0Centr < 3.50) iC = 0;
5148 else if(v0Centr < 4.94) iC = 1;
5149 else if(v0Centr < 6.98) iC = 2;
5150 else if(v0Centr < 8.55) iC = 3;
5151 else if(v0Centr < 9.88) iC = 4;
5152 else if(v0Centr < 11.04) iC = 5;
5153 else if(v0Centr < 12.09) iC = 6;
5154 else if(v0Centr < 13.05) iC = 7;
5157 eventplane = headerH->ReactionPlaneAngle();
5158 if(eventplane > TMath::Pi()/2 && eventplane <= TMath::Pi()*3/2) eventplane-=TMath::Pi();
5159 if(eventplane > TMath::Pi()*3/2) eventplane-=2*TMath::Pi();
5160 fHistEventPlaneTruth->Fill(iC,eventplane);
5161 //gReactionPlane *= TMath::RadToDeg();
5166 else if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD") {
5167 //reset Q vector info
5169 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
5172 Int_t run = event->GetRunNumber();
5175 // Load the calibrations run dependent
5176 if(! fIsAfter2011) OpenInfoCalbration(run);
5182 if (v0Centr > 80) return 999.; // analysis only for 0-80% centrality classes
5184 if(v0Centr < 5) iC = 0;
5185 else if(v0Centr < 10) iC = 1;
5186 else if(v0Centr < 20) iC = 2;
5187 else if(v0Centr < 30) iC = 3;
5188 else if(v0Centr < 40) iC = 4;
5189 else if(v0Centr < 50) iC = 5;
5190 else if(v0Centr < 60) iC = 6;
5191 else if(v0Centr < 70) iC = 7;
5197 //reset Q vector info
5198 Double_t Qxa2 = 0, Qya2 = 0;
5199 Double_t Qxc2 = 0, Qyc2 = 0;
5200 Double_t Qxa3 = 0, Qya3 = 0;
5201 Double_t Qxc3 = 0, Qyc3 = 0;
5204 //MC: from reaction plane
5207 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
5209 evplaneMC = header->GetReactionPlaneAngle();//[0, 360]
5210 //make it within [-pi/2,pi/2] to make it general
5211 if(evplaneMC > TMath::Pi()/2 && evplaneMC <= TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
5212 if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
5213 fHistEventPlaneTruth->Fill(iC,evplaneMC);
5215 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
5219 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
5221 if (collGeometry){//get the reaction plane from MC header
5222 gReactionPlane = collGeometry->ReactionPlaneAngle();//[0,180]
5226 //taken from vnv0 code(get the TPC, V0A, V0C event plane using truth tracks)
5227 TClonesArray *mcArray = NULL;
5228 mcArray = (TClonesArray*)event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
5230 Float_t QxMCv2[3] = {0,0,0};
5231 Float_t QyMCv2[3] = {0,0,0};
5232 Float_t QxMCv3[3] = {0,0,0};
5233 Float_t QyMCv3[3] = {0,0,0};
5234 Float_t EvPlaneMCV2[3] = {0,0,0};
5235 Float_t EvPlaneMCV3[3] = {0,0,0};
5236 Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
5237 Float_t etaMax[3] = {4.88,-1.8,0.8};
5239 // analysis on MC tracks
5240 Int_t nMCtrack = mcArray->GetEntries() ;
5242 // EP computation with MC tracks
5243 for(Int_t iT=0;iT < nMCtrack;iT++){
5244 AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
5245 if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
5247 Float_t eta = mctr->Eta();
5248 for(Int_t iD=0;iD<3;iD++){
5249 if(eta > etaMin[iD] && eta < etaMax[iD]){
5250 Float_t phi = mctr->Phi();
5251 QxMCv2[iD] += TMath::Cos(2*phi);
5252 QyMCv2[iD] += TMath::Sin(2*phi);
5253 QxMCv3[iD] += TMath::Cos(3*phi);
5254 QyMCv3[iD] += TMath::Sin(3*phi);
5259 EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
5260 EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
5261 EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
5262 fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
5263 fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
5264 fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
5265 fgPsi2v0aMC = EvPlaneMCV2[0];
5266 fgPsi2v0cMC = EvPlaneMCV2[1];
5267 fgPsi2tpcMC = EvPlaneMCV2[2];
5270 EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
5271 EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
5272 EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
5273 fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
5274 fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
5275 fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
5276 fgPsi3v0aMC = EvPlaneMCV3[0];
5277 fgPsi3v0cMC = EvPlaneMCV3[1];
5278 fgPsi3tpcMC = EvPlaneMCV3[2];
5285 Int_t nAODTracks = event->GetNumberOfTracks();
5287 // TPC EP needed for resolution studies (TPC subevent)
5288 //AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
5289 //Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
5290 Double_t Qx2 = 0, Qy2 = 0;
5291 Double_t Qx3 = 0, Qy3 = 0;
5293 for(Int_t iT = 0; iT < nAODTracks; iT++) {
5295 AliAODTrack* aodTrack =(AliAODTrack*) event->GetTrack(iT);
5301 Bool_t trkFlag = aodTrack->TestFilterBit(1);
5303 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
5306 Double_t b[2] = {-99., -99.};
5307 Double_t bCov[3] = {-99., -99., -99.};
5310 AliAODTrack param(*aodTrack);
5311 if (!param.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov)){
5315 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
5318 Qx2 += TMath::Cos(2*aodTrack->Phi());
5319 Qy2 += TMath::Sin(2*aodTrack->Phi());
5320 Qx3 += TMath::Cos(3*aodTrack->Phi());
5321 Qy3 += TMath::Sin(3*aodTrack->Phi());
5325 Float_t evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
5326 Float_t evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
5328 fgPsi2tpc = evPlAng2;
5329 fgPsi3tpc = evPlAng3;
5331 fPhiRPTPC->Fill(iC,evPlAng2);
5332 fPhiRPTPCv3->Fill(iC,evPlAng3);
5337 AliAODVZERO* aodV0 = event->GetVZEROData();
5339 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
5340 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
5341 Float_t multv0 = aodV0->GetMultiplicity(iv0);
5344 if(fAnalysisType == "AOD"){//not for reco MC tracks, only for real data
5345 if (iv0 < 32){ // V0C
5346 Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5347 Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5348 Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5349 Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5351 Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5352 Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5353 Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5354 Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5358 if (iv0 < 32){ // V0C
5359 Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5360 Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5361 Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5362 Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5364 Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5365 Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5366 Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5367 Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5372 //grab for each centrality the proper histo with the Qx and Qy to do the recentering
5373 Double_t Qxamean2 = fMeanQ[iCcal][1][0];
5374 Double_t Qxarms2 = fWidthQ[iCcal][1][0];
5375 Double_t Qyamean2 = fMeanQ[iCcal][1][1];
5376 Double_t Qyarms2 = fWidthQ[iCcal][1][1];
5377 Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
5378 Double_t Qxarms3 = fWidthQv3[iCcal][1][0];
5379 Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
5380 Double_t Qyarms3 = fWidthQv3[iCcal][1][1];
5382 Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
5383 Double_t Qxcrms2 = fWidthQ[iCcal][0][0];
5384 Double_t Qycmean2 = fMeanQ[iCcal][0][1];
5385 Double_t Qycrms2 = fWidthQ[iCcal][0][1];
5386 Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
5387 Double_t Qxcrms3 = fWidthQv3[iCcal][0][0];
5388 Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
5389 Double_t Qycrms3 = fWidthQv3[iCcal][0][1];
5391 Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
5392 Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
5393 Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
5394 Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
5395 Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
5396 Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
5397 Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
5398 Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
5400 //to calculate 2nd order event plane with v0M
5401 Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
5402 /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
5403 Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
5404 /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
5406 //here the calculated event plane is within -Pi to +Pi(delete it , no use here , only for definition)
5407 Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
5408 Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
5409 Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
5413 Float_t evPlAngV0ACor2=999.;
5414 Float_t evPlAngV0CCor2=999.;
5415 Float_t evPlAngV0ACor3=999.;
5416 Float_t evPlAngV0CCor3=999.;
5419 if(fAnalysisType == "AOD"){
5420 evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
5421 evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
5422 evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
5423 evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
5426 evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
5427 evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
5428 evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
5429 evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
5433 AliEventplane *ep = event->GetEventplane();
5434 evPlAngV0ACor2 = ep->GetEventplane("V0A", event, 2);
5435 evPlAngV0CCor2 = ep->GetEventplane("V0C", event, 2);
5436 evPlAngV0ACor3 = ep->GetEventplane("V0A", event, 3);
5437 evPlAngV0CCor3 = ep->GetEventplane("V0C", event, 3);
5440 fgPsi2v0a = evPlAngV0ACor2;
5441 fgPsi2v0c = evPlAngV0CCor2;
5442 fgPsi3v0a = evPlAngV0ACor3;
5443 fgPsi3v0c = evPlAngV0CCor3;
5445 // Fill EP distribution histograms evPlAng
5447 fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
5448 fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
5450 fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
5451 fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
5453 // Fill histograms needed for resolution evaluation
5454 fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
5455 fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
5456 fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
5458 fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
5459 fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
5460 fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
5464 Float_t gVZEROEventPlane = -10.;
5465 Float_t gReactionPlane = -10.;
5466 Double_t qxTot = 0.0, qyTot = 0.0;
5468 AliEventplane *ep = event->GetEventplane();
5470 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
5471 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
5472 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
5473 gReactionPlane = gVZEROEventPlane;
5477 //return gReactionPlane;
5479 //make the final 2nd order event plane within 0 to Pi
5480 //using data and reco tracks only
5481 if(fgPsi2v0a!=999. && fgPsi2v0a < 0.) fgPsi2v0a += TMath::Pi();
5482 if(fgPsi2v0c!=999. && fgPsi2v0c < 0.) fgPsi2v0c += TMath::Pi();
5483 if(fgPsi2tpc!=999. && fgPsi2tpc < 0.) fgPsi2tpc += TMath::Pi();
5484 //using truth tracks only
5485 if(evplaneMC!=999. && evplaneMC < 0.) evplaneMC += TMath::Pi();
5486 if(fgPsi2v0aMC!=999. && fgPsi2v0aMC < 0.) fgPsi2v0aMC += TMath::Pi();
5487 if(fgPsi2v0cMC!=999. && fgPsi2v0cMC < 0.) fgPsi2v0cMC += TMath::Pi();
5488 if(fgPsi2tpcMC!=999. && fgPsi2tpcMC < 0.) fgPsi2tpcMC += TMath::Pi();
5489 //for the time being leave the 3rd order event planes within -pi/3 t0 +pi/3
5491 if(truth){//for truth MC
5492 if(fV2 && fEPdet=="header")eventplane=evplaneMC;
5493 if(fV2 && fEPdet=="V0A")eventplane=fgPsi2v0aMC;
5494 if(fV2 && fEPdet=="V0C")eventplane=fgPsi2v0cMC;
5495 if(fV2 && fEPdet=="TPC")eventplane=fgPsi2tpcMC;
5497 if(fV3 && fEPdet=="V0A")eventplane=fgPsi3v0aMC;
5498 if(fV3 && fEPdet=="V0C")eventplane=fgPsi3v0cMC;
5499 if(fV3 && fEPdet=="TPC")eventplane=fgPsi3tpcMC;
5501 else{//for data and recoMC
5502 if(fV2 && fEPdet=="V0A")eventplane=fgPsi2v0a;
5503 if(fV2 && fEPdet=="V0C")eventplane=fgPsi2v0c;
5504 if(fV2 && fEPdet=="TPC")eventplane=fgPsi2tpc;
5506 if(fV3 && fEPdet=="V0A")eventplane=fgPsi3v0a;
5507 if(fV3 && fEPdet=="V0C")eventplane=fgPsi3v0c;
5508 if(fV3 && fEPdet=="TPC")eventplane=fgPsi3tpc;
5513 }//AOD/MCAOD condition
5515 else eventplane=999.;
5520 //------------------------------------------------------------------------------------------------------------------
5521 void AliTwoParticlePIDCorr::OpenInfoCalbration(Int_t run){
5522 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
5523 TFile *foadb = TFile::Open(oadbfilename.Data());
5526 printf("OADB file %s cannot be opened\n",oadbfilename.Data());
5530 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
5532 printf("OADB object hMultV0BefCorr is not available in the file\n");
5536 if(!(cont->GetObject(run))){
5537 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
5540 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
5542 TF1 *fpol0 = new TF1("fpol0","pol0");
5543 fMultV0->Fit(fpol0,"","",0,31);
5544 fV0Cpol = fpol0->GetParameter(0);
5545 fMultV0->Fit(fpol0,"","",32,64);
5546 fV0Apol = fpol0->GetParameter(0);
5548 for(Int_t iside=0;iside<2;iside++){
5549 for(Int_t icoord=0;icoord<2;icoord++){
5550 for(Int_t i=0;i < 9;i++){
5552 if(iside==0 && icoord==0)
5553 snprintf(namecont,100,"hQxc2_%i",i);
5554 else if(iside==1 && icoord==0)
5555 snprintf(namecont,100,"hQxa2_%i",i);
5556 else if(iside==0 && icoord==1)
5557 snprintf(namecont,100,"hQyc2_%i",i);
5558 else if(iside==1 && icoord==1)
5559 snprintf(namecont,100,"hQya2_%i",i);
5561 cont = (AliOADBContainer*) foadb->Get(namecont);
5563 printf("OADB object %s is not available in the file\n",namecont);
5567 if(!(cont->GetObject(run))){
5568 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5571 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5572 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5575 if(iside==0 && icoord==0)
5576 snprintf(namecont,100,"hQxc3_%i",i);
5577 else if(iside==1 && icoord==0)
5578 snprintf(namecont,100,"hQxa3_%i",i);
5579 else if(iside==0 && icoord==1)
5580 snprintf(namecont,100,"hQyc3_%i",i);
5581 else if(iside==1 && icoord==1)
5582 snprintf(namecont,100,"hQya3_%i",i);
5584 cont = (AliOADBContainer*) foadb->Get(namecont);
5586 printf("OADB object %s is not available in the file\n",namecont);
5590 if(!(cont->GetObject(run))){
5591 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5594 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5595 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5601 //____________________________________________________________________
5602 void AliTwoParticlePIDCorr::FillPIDEventPlane(Double_t centrality,Int_t par,Float_t trigphi,Float_t fReactionPlane)
5605 // Event plane (determine psi bin)
5606 Double_t gPsiMinusPhi = 0.;
5607 Double_t gPsiMinusPhiBin = -10.;
5608 if(fRequestEventPlane){
5609 gPsiMinusPhi = TMath::Abs(trigphi - fReactionPlane);
5611 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
5612 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
5613 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
5614 gPsiMinusPhiBin = 0.0;
5616 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
5617 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
5618 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
5619 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
5620 gPsiMinusPhiBin = 1.0;
5622 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
5623 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
5624 gPsiMinusPhiBin = 2.0;
5627 gPsiMinusPhiBin = 3.0;
5629 fEventPlanePID->Fill(centrality,gPsiMinusPhiBin,(Float_t)par);
5633 //____________________________________________________________________________________________________
5635 TObjArray* AliTwoParticlePIDCorr::GetV0Particles(AliVEvent* event,Double_t Centrality)
5638 AliAODEvent* fAOD = dynamic_cast<AliAODEvent*>(event);
5640 //function to select v0's from AODs
5641 trkVtx=fAOD->GetPrimaryVertex();
5642 Float_t xv=trkVtx->GetX(), yv=trkVtx->GetY(), zv=trkVtx->GetZ();
5643 Int_t nV0sTot = fAOD->GetNumberOfV0s();
5645 TObjArray * selectedV0s = new TObjArray;
5646 selectedV0s->SetOwner(kTRUE);
5648 for (Int_t iV0 = 0; iV0 < nV0sTot; iV0++)
5651 AliAODv0 *v0=fAOD->GetV0(iV0);
5653 if(!CheckStatusv0(v0)) continue;
5655 AliAODTrack *ptrack=(AliAODTrack*)v0->GetDaughter(0);
5656 AliAODTrack *ntrack=(AliAODTrack*)v0->GetDaughter(1);
5658 Bool_t cutK0sPID=kFALSE;
5659 Bool_t cutLambdaPID=kFALSE;
5660 Bool_t cutAntiLambdaPID=kFALSE;
5662 if(fUsev0DaughterPID)
5664 //use fHelperPID check PID of the daughter tracks
5665 //v0 daughter PID may be helpful in distangling k0S and (Anti)Lamda
5667 Int_t PIDptrack = GetParticle(ptrack,kFALSE);
5668 Int_t PIDntrack = GetParticle(ntrack ,kFALSE);
5670 if(PIDptrack ==0 && PIDntrack == 0) cutK0sPID=kTRUE;
5672 if(PIDptrack==2 && PIDntrack ==0) cutLambdaPID=kTRUE;
5674 if (PIDptrack ==0 && PIDntrack == 2) cutAntiLambdaPID=kTRUE;
5678 // effective mass calculations for each hypothesis(without daughter PID)
5679 Double_t InvMassK0s = v0->MassK0Short();
5680 Double_t InvMassAntiLambda = v0->MassAntiLambda();
5681 Double_t InvMassLambda = v0->MassLambda();
5683 Float_t v0Pt=TMath::Sqrt(v0->Pt2V0());
5684 Float_t v0Eta=v0->Eta();
5685 Float_t v0Phi=v0->Phi();
5687 //This is simply raw v0 without any specialised cut
5688 fHistRawPtCentInvK0s->Fill(InvMassK0s,v0Pt,Centrality);
5689 fHistRawPtCentInvLambda->Fill(InvMassLambda,v0Pt,Centrality);
5690 fHistRawPtCentInvAntiLambda->Fill(InvMassAntiLambda,v0Pt,Centrality);
5694 v0->GetSecondaryVtx(xyz);
5696 dx=xyz[0]-xv, dy=xyz[1]-yv, dz=xyz[2]-zv;
5698 // Float_t v0DecayRadius=TMath::Sqrt(dx*dx + dy*dy);
5699 Float_t v0DecayLength=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
5700 // VO's main characteristics to check the reconstruction cuts
5701 // Float_t DcaV0Daughters = v0->DcaV0Daughters();
5702 Float_t V0cosPointAngle = v0->CosPointingAngle(trkVtx);
5703 // Float_t DcaPosToPrimVertex = v0->DcaPosToPrimVertex();
5704 //Float_t DcaNegToPrimVertex = v0->DcaNegToPrimVertex();
5705 //Float_t Dcav0PVz = v0->DcaV0ToPrimVertex();
5706 Float_t v0Pz=v0->Pz();
5707 Float_t v0P= TMath::Sqrt(v0Pt*v0Pt + v0Pz*v0Pz);
5709 Float_t ctauLambda =999.;
5710 Float_t ctauAntiLambda = 999.;
5711 Float_t ctauK0s = 999.;
5714 ctauLambda = (v0DecayLength*InvMassLambda)/v0P;
5715 ctauAntiLambda = (v0DecayLength*InvMassAntiLambda)/v0P;
5716 ctauK0s = (v0DecayLength*InvMassK0s)/v0P;
5719 Bool_t ctauCutK0s= ctauK0s < NCtau*fCutctauK0s ; //ctauK0s 2.68 cm, mean life time of K0s is 8.95 x10^(-11)
5720 Bool_t ctauCutLambda = ctauLambda < NCtau*fCutctauLambda; //ctauLambda 7.8 cm ,mean life is 2.6 x10 ^(-10) ***** 3xctau is the accepted limit
5721 Bool_t ctauCutAntiLambda= ctauAntiLambda < NCtau*fCutctauAntiLambda;
5723 Bool_t RapCutK0s = v0->RapK0Short() < fRapCutK0s;
5724 Bool_t RapCutLambda = v0->RapLambda() < fRapCutLambda;
5725 Bool_t RapCutAntiLambda = v0->Y(-3122) < fRapCutLambda;
5727 Bool_t CPACut= V0cosPointAngle > fMinCPA; //cosine of pointing angle with v0 should be greater than 0.998
5729 //Now we put a loose mass cut which will be tightened later
5730 Bool_t MassCutLooseK0s=(TMath::Abs(InvMassK0s - 0.497614) < 0.1);
5731 Bool_t MassCutLooseLambda=(TMath::Abs(InvMassLambda - 1.115683) < 0.1); // cut is same for Anti-Lambda
5732 Bool_t MassCutLooseAntiLambda=(TMath::Abs(InvMassAntiLambda - 1.115683) < 0.1); // cut is same for Anti-Lambda
5734 //Special Cut for Kshort arementeros podalanski plot
5735 Bool_t ArmenterosCut =kFALSE;
5736 if(ctauCutK0s && RapCutK0s && CPACut && MassCutLooseK0s)
5739 Float_t lAlphaV0 = v0->AlphaV0();
5740 Float_t lPtArmV0 = v0->PtArmV0();
5742 ArmenterosCut = lPtArmV0 > TMath::Abs(0.2*lAlphaV0);
5746 Bool_t IskShortOk=(ctauCutK0s && RapCutK0s && CPACut && MassCutLooseK0s && ArmenterosCut);
5748 Bool_t IsLambdaOk=(ctauCutLambda && RapCutLambda && CPACut && MassCutLooseLambda);
5750 Bool_t IsAntiLambdaOk=(ctauCutAntiLambda && RapCutAntiLambda && CPACut && MassCutLooseAntiLambda);
5752 //Difference on Lambda and Anti-Lambda can be made through daughter PID
5755 Int_t particletype=999;
5757 if( IskShortOk || cutK0sPID )
5759 fHistFinalPtCentInvK0s->Fill(InvMassK0s,v0Pt,Centrality);
5762 Short_t chargeval=0;
5763 Float_t effmatrix=1.0;
5764 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassK0s,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5765 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5766 selectedV0s->Add(copy1);
5771 if(IsLambdaOk || cutLambdaPID)
5773 fHistFinalPtCentInvLambda->Fill(InvMassLambda,v0Pt,Centrality);
5774 //Add in the LRCParticle and give Lambda a tag 5
5777 Short_t chargeval=0;
5778 Float_t effmatrix=1.0;
5779 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassLambda,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5780 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5781 selectedV0s->Add(copy1);
5784 if(IsAntiLambdaOk || cutAntiLambdaPID)
5786 fHistFinalPtCentInvLambda->Fill(InvMassAntiLambda,v0Pt,Centrality);
5787 //Add in the LRCParticle and give Lambda a tag 6
5788 particletype=SpALam;
5789 Short_t chargeval=0;
5790 Float_t effmatrix=1.0;
5791 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassAntiLambda,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5792 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5793 selectedV0s->Add(copy1);
5802 //___________________________________________________________________
5803 Bool_t AliTwoParticlePIDCorr :: CheckStatusv0Daughter(AliAODTrack *t1 ,AliAODTrack *t2)
5805 if (!t1->IsOn(AliAODTrack::kTPCrefit) || !t2->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
5806 // Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
5807 if(t1->GetTPCClusterInfo(2,1)<fDaugNClsTPC || t2->GetTPCClusterInfo(2,1)<fDaugNClsTPC) return kFALSE ;
5809 // ---------------- Fraction of TPC Shared Cluster
5810 Float_t fracPosDaugTPCSharedMap = GetFractionTPCSharedCls(t1);
5811 Float_t fracNegDaugTPCSharedMap = GetFractionTPCSharedCls(t2);
5813 if( (fracPosDaugTPCSharedMap > fFracTPCcls) || (fracNegDaugTPCSharedMap > fFracTPCcls) )
5819 //___________________________________________________________________________________________
5821 Float_t AliTwoParticlePIDCorr :: GetFractionTPCSharedCls( AliAODTrack *track)
5823 // Rejects tracks with shared clusters after filling a control histogram
5824 // This overload is used for primaries
5826 // Get the shared maps
5827 const TBits sharedMap = track->GetTPCSharedMap();
5829 return 1.*sharedMap.CountBits()/track->GetTPCNclsF();
5832 //______________________________________________________________________
5833 Bool_t AliTwoParticlePIDCorr :: CheckStatusv0(AliAODv0 *v1)
5836 // Offline reconstructed V0 only
5837 if (v1->GetOnFlyStatus()) return kFALSE;
5839 AliAODTrack *ptrack=(AliAODTrack *)v1->GetDaughter(0);
5840 AliAODTrack *ntrack=(AliAODTrack *)v1->GetDaughter(1);
5842 if(!ptrack || !ntrack) return kFALSE;
5844 if(ptrack->Charge()==-1 || ntrack->Charge()==1) return kFALSE; //remove wrongly identified charge pairs
5846 if(ptrack->Charge()==0 || ntrack->Charge()==0) return kFALSE; //remove uncharged pairs
5848 if(ptrack->Charge() == ntrack->Charge()) return kFALSE; //remove like sign pairs
5850 if(!CheckStatusv0Daughter(ptrack,ntrack)) return kFALSE;//daughters need to pass some basic cuts
5852 if(TMath::Abs(ptrack->Eta()) > fmaxeta || TMath::Abs(ntrack->Eta()) > fmaxeta) return kFALSE; // remove daughters beyond eta bound |0.8|
5854 if(ptrack->Pt() < fMinPtDaughter || ntrack->Pt() < fMinPtDaughter) return kFALSE; // remove daughter tracks below minmum p |1.0 GeV/c|
5856 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|
5858 // Daughters: Impact parameter of daughter to prim vtx
5859 Float_t xy = v1->DcaPosToPrimVertex();
5860 if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE; //0.1 cm
5861 xy = v1->DcaNegToPrimVertex();
5862 if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE; //0.1 cm
5865 Float_t dca = v1->DcaV0Daughters();
5866 if (dca>fMaxDCADaughter) return kFALSE; //1.0 cm
5868 // V0: Cosine of the pointing angle
5869 Float_t cpa=v1->CosPointingAngle(trkVtx); //0.997
5870 if (cpa<fMinCPA) return kFALSE;
5872 // V0: Fiducial volume
5873 Double_t xyz[3]; v1->GetSecondaryVtx(xyz);
5874 Float_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
5875 if (r2<5.*5.) return kFALSE;
5876 if (r2>lMax*lMax) return kFALSE; //lmax=100 cm
5882 //__________________________________________________________________________________________
5883 Bool_t AliTwoParticlePIDCorr::IsTrackFromV0(AliAODEvent* fAOD,AliAODTrack* track)
5885 //to check whether a daughter being taken as associated
5886 Int_t assoID = track->GetID();
5888 for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
5889 AliAODv0* aodV0 = fAOD->GetV0(i);
5891 AliAODTrack *trackPos=(AliAODTrack *)(aodV0->GetDaughter(0));
5892 AliAODTrack *trackNeg=(AliAODTrack *)(aodV0->GetDaughter(1));
5895 Int_t negID = trackNeg->GetID();
5896 Int_t posID = trackPos->GetID();
5898 if ((TMath::Abs(negID)+1)==(TMath::Abs(assoID))){ return kTRUE;}
5899 if ((TMath::Abs(posID)+1)==(TMath::Abs(assoID))){ return kTRUE;}
5900 //----------------------------------
5907 //----------------------------------------------------------
5908 void AliTwoParticlePIDCorr::Terminate(Option_t *)
5910 // Draw result to screen, or perform fitting, normalizations
5911 // Called once at the end of the query
5912 fOutput = dynamic_cast<TList*> (GetOutputData(1));
5913 if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
5917 //------------------------------------------------------------------