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),
159 fhistImpactParmvsMult(0),
182 fhistJetTrigestimate(0),
183 fTwoTrackDistancePtdip(0x0),
184 fTwoTrackDistancePtdipmix(0x0),
185 fCentralityCorrelation(0x0),
186 fHistVZEROAGainEqualizationMap(0),
187 fHistVZEROCGainEqualizationMap(0),
188 fHistVZEROChannelGainEqualizationMap(0),
189 fCentralityWeights(0),
192 fHistEQVZEROvsTPCmultiplicity(0x0),
193 fHistEQVZEROAvsTPCmultiplicity(0x0),
194 fHistEQVZEROCvsTPCmultiplicity(0x0),
195 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
196 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
197 fHistVZEROCvsVZEROAmultiplicity(0x0),
198 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
199 fHistVZEROSignal(0x0),
200 fHistEventPlaneTruth(0x0),
201 fHistPsiMinusPhi(0x0),
216 gReactionPlane(999.),
244 fControlConvResoncances(0),
259 fCorrelatonTruthPrimary(0),
260 fCorrelatonTruthPrimarymix(0),
266 fTHnCorrIDUNIDmix(0),
268 fTHnTrigcountMCTruthPrim(0),
271 fAnalysisType("AOD"),
273 ftwoTrackEfficiencyCutDataReco(kTRUE),
274 fTwoTrackCutMinRadius(0.8),
275 fTwoTrackCutMaxRadius(2.5),
276 twoTrackEfficiencyCutValue(0.02),
282 fRequestTOFPID(kTRUE),
283 fPIDType(NSigmaTPCTOF),
284 fFIllPIDQAHistos(kTRUE),
287 fdiffPIDcutvalues(kFALSE),
292 fHighPtKaonNSigmaPID(-1),
293 fHighPtKaonSigma(3.5),
294 fUseExclusiveNSigma(kFALSE),
295 fRemoveTracksT0Fill(kFALSE),
297 fTriggerSelectCharge(0),
298 fAssociatedSelectCharge(0),
299 fTriggerRestrictEta(-1),
300 fEtaOrdering(kFALSE),
301 fCutConversions(kFALSE),
302 fCutResonances(kFALSE),
303 fRejectResonanceDaughters(-1),
305 fInjectedSignals(kFALSE),
306 fRemoveWeakDecays(kFALSE),
307 fRemoveDuplicates(kFALSE),
308 fapplyTrigefficiency(kFALSE),
309 fapplyAssoefficiency(kFALSE),
310 ffillefficiency(kFALSE),
311 fmesoneffrequired(kFALSE),
312 fkaonprotoneffrequired(kFALSE),
316 fUsev0DaughterPID(kFALSE),
317 fMinPtDaughter(1.0),// v0 related cut starts here
320 fMaxDCADaughter(1.0),
323 fHistRawPtCentInvK0s(0x0),
324 fHistRawPtCentInvLambda(0x0),
325 fHistRawPtCentInvAntiLambda(0x0),
326 fHistFinalPtCentInvK0s(0x0),
327 fHistFinalPtCentInvLambda(0x0),
328 fHistFinalPtCentInvAntiLambda(0x0),
332 fCutctauAntiLambda(7.8),
339 for ( Int_t i = 0; i < 16; i++) {
343 for ( Int_t i = 0; i < 6; i++ ){
344 fTrackHistEfficiency[i] = NULL;
345 effcorection[i]=NULL;
348 for ( Int_t i = 0; i < 2; i++ ){
349 fTwoTrackDistancePt[i]=NULL;
350 fTwoTrackDistancePtmix[i]=NULL;
353 for(Int_t ipart=0;ipart<NSpecies;ipart++)
354 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
355 fnsigmas[ipart][ipid]=999.;
357 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
359 for(Int_t i = 0; i != 2; ++i)
360 for(Int_t j = 0; j != 2; ++j)
361 for(Int_t iC = 0; iC < 9; iC++){
362 fMeanQ[iC][i][j] = 0.;
363 fWidthQ[iC][i][j] = 1.;
364 fMeanQv3[iC][i][j] = 0.;
365 fWidthQv3[iC][i][j] = 1.;
369 //________________________________________________________________________
370 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
371 :AliAnalysisTaskSE(name),
375 fCentralityMethod("V0A"),
376 fPPVsMultUtils(kFALSE),
378 fRequestEventPlane(kFALSE),
379 fRequestEventPlanemixing(kFALSE),
380 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
385 fSharedClusterCut(-1),
386 fSharedTPCmapCut(-1),
387 fSharedfraction_Pair_cut(-1),
389 skipParticlesAbove(0),
394 ffilltrigassoUNID(kFALSE),
395 ffilltrigUNIDassoID(kFALSE),
396 ffilltrigIDassoUNID(kTRUE),
397 ffilltrigIDassoID(kFALSE),
398 ffilltrigIDassoIDMCTRUTH(kFALSE),
399 fMaxNofMixingTracks(50000),
400 fPtOrderMCTruth(kTRUE),
401 fPtOrderDataReco(kTRUE),
402 fWeightPerEvent(kFALSE),
403 fTriggerSpeciesSelection(kFALSE),
404 fAssociatedSpeciesSelection(kFALSE),
405 fRandomizeReactionPlane(kFALSE),
406 fTriggerSpecies(SpPion),
407 fAssociatedSpecies(SpPion),
410 fSelectHighestPtTrig(kFALSE),
411 fcontainPIDtrig(kTRUE),
412 fcontainPIDasso(kFALSE),
414 frejectPileUp(kFALSE),
419 fselectprimaryTruth(kTRUE),
420 fonlyprimarydatareco(kFALSE),
423 ffillhistQAReco(kFALSE),
424 ffillhistQATruth(kFALSE),
425 kTrackVariablesPair(0),
437 fhistImpactParmvsMult(0),
460 fhistJetTrigestimate(0),
461 fTwoTrackDistancePtdip(0x0),
462 fTwoTrackDistancePtdipmix(0x0),
463 fCentralityCorrelation(0x0),
464 fHistVZEROAGainEqualizationMap(0),
465 fHistVZEROCGainEqualizationMap(0),
466 fHistVZEROChannelGainEqualizationMap(0),
467 fCentralityWeights(0),
470 fHistEQVZEROvsTPCmultiplicity(0x0),
471 fHistEQVZEROAvsTPCmultiplicity(0x0),
472 fHistEQVZEROCvsTPCmultiplicity(0x0),
473 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
474 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
475 fHistVZEROCvsVZEROAmultiplicity(0x0),
476 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
477 fHistVZEROSignal(0x0),
478 fHistEventPlaneTruth(0x0),
479 fHistPsiMinusPhi(0x0),
494 gReactionPlane(999.),
522 fControlConvResoncances(0),
537 fCorrelatonTruthPrimary(0),
538 fCorrelatonTruthPrimarymix(0),
544 fTHnCorrIDUNIDmix(0),
546 fTHnTrigcountMCTruthPrim(0),
549 fAnalysisType("AOD"),
551 ftwoTrackEfficiencyCutDataReco(kTRUE),
552 fTwoTrackCutMinRadius(0.8),
553 fTwoTrackCutMaxRadius(2.5),
554 twoTrackEfficiencyCutValue(0.02),
560 fRequestTOFPID(kTRUE),
561 fPIDType(NSigmaTPCTOF),
562 fFIllPIDQAHistos(kTRUE),
565 fdiffPIDcutvalues(kFALSE),
570 fHighPtKaonNSigmaPID(-1),
571 fHighPtKaonSigma(3.5),
572 fUseExclusiveNSigma(kFALSE),
573 fRemoveTracksT0Fill(kFALSE),
575 fTriggerSelectCharge(0),
576 fAssociatedSelectCharge(0),
577 fTriggerRestrictEta(-1),
578 fEtaOrdering(kFALSE),
579 fCutConversions(kFALSE),
580 fCutResonances(kFALSE),
581 fRejectResonanceDaughters(-1),
583 fInjectedSignals(kFALSE),
584 fRemoveWeakDecays(kFALSE),
585 fRemoveDuplicates(kFALSE),
586 fapplyTrigefficiency(kFALSE),
587 fapplyAssoefficiency(kFALSE),
588 ffillefficiency(kFALSE),
589 fmesoneffrequired(kFALSE),
590 fkaonprotoneffrequired(kFALSE),
594 fUsev0DaughterPID(kFALSE),
595 fMinPtDaughter(1.0),// v0 related cut starts here
598 fMaxDCADaughter(1.0),
601 fHistRawPtCentInvK0s(0x0),
602 fHistRawPtCentInvLambda(0x0),
603 fHistRawPtCentInvAntiLambda(0x0),
604 fHistFinalPtCentInvK0s(0x0),
605 fHistFinalPtCentInvLambda(0x0),
606 fHistFinalPtCentInvAntiLambda(0x0),
610 fCutctauAntiLambda(7.8),
616 // The last in the above list should not have a comma after it
620 for ( Int_t i = 0; i < 16; i++) {
624 for ( Int_t i = 0; i < 6; i++ ){
625 fTrackHistEfficiency[i] = NULL;
626 effcorection[i]=NULL;
630 for ( Int_t i = 0; i < 2; i++ ){
631 fTwoTrackDistancePt[i]=NULL;
632 fTwoTrackDistancePtmix[i]=NULL;
635 for(Int_t ipart=0;ipart<NSpecies;ipart++)
636 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
637 fnsigmas[ipart][ipid]=999.;
639 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
641 for(Int_t i = 0; i != 2; ++i)
642 for(Int_t j = 0; j != 2; ++j)
643 for(Int_t iC = 0; iC < 9; iC++){
644 fMeanQ[iC][i][j] = 0.;
645 fWidthQ[iC][i][j] = 1.;
646 fMeanQv3[iC][i][j] = 0.;
647 fWidthQv3[iC][i][j] = 1.;
651 // Define input and output slots here (never in the dummy constructor)
652 // Input slot #0 works with a TChain - it is connected to the default input container
653 // Output slot #1 writes into a TH1 container
654 DefineInput(0, TChain::Class());
656 DefineOutput(1, TList::Class()); // for output list
657 DefineOutput(2, TList::Class());
658 DefineOutput(3, TList::Class());
662 //________________________________________________________________________
663 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
665 // Destructor. Clean-up the output list, but not the histograms that are put inside
666 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
667 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
672 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
677 if(fRequestEventPlane){
678 if (fList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
683 if (fPID) delete fPID;
684 if (fPIDCombined) delete fPIDCombined;
687 //________________________________________________________________________
689 //////////////////////////////////////////////////////////////////////////////////////////////////
691 TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
692 // returns histo named name
693 return (TH2F*) fOutputList->FindObject(name);
696 //////////////////////////////////////////////////////////////////////////////////////////////////
698 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
702 // Puts the argument in the range [-pi/2,3 pi/2].
705 if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
706 if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
711 //________________________________________________________________________
712 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
715 // Called once (on the worker node)
717 // global switch disabling the reference
718 // (to avoid "Replacing existing TH1" if several wagons are created in train)
719 Bool_t oldStatus = TH1::AddDirectoryStatus();
720 TH1::AddDirectory(kFALSE);
722 const Int_t nPsiTOF = 10;
723 const Int_t nCentrBin = 9;
726 fOutput = new TList();
727 fOutput->SetOwner(); // IMPORTANT!
729 fOutputList = new TList;
730 fOutputList->SetOwner();
731 fOutputList->SetName("PIDQAList");
733 if(fRequestEventPlane){
736 fList->SetName("EPQAList");
738 fEventCounter = new TH1F("fEventCounter","EventCounter", 19, 0.5,19.5);
739 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
740 fEventCounter->GetXaxis()->SetBinLabel(3,"After PileUP Cut");//only for Data
741 fEventCounter->GetXaxis()->SetBinLabel(5,"Have A Vertex");
742 fEventCounter->GetXaxis()->SetBinLabel(7,"After vertex Cut");
743 fEventCounter->GetXaxis()->SetBinLabel(9,"Getting centrality");
744 fEventCounter->GetXaxis()->SetBinLabel(11,"After centrality flattening");
745 fEventCounter->GetXaxis()->SetBinLabel(13,"Within 0-100% centrality");
746 fEventCounter->GetXaxis()->SetBinLabel(15,"Event Analyzed");
747 //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
748 fOutput->Add(fEventCounter);
750 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
751 fOutput->Add(fEtaSpectrasso);
753 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
754 fOutput->Add(fphiSpectraasso);
756 if(fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE || fCentralityMethod == "MC_b"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality_ImpactParam;multiplicity", 101, 0, 101, 20000, 0,40000);
757 fOutput->Add(fCentralityCorrelation);
760 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")
762 TString gCentName[8] = {"V0A","V0C","V0M","V0AEq","V0CEq","V0MEq","CL1","ZNA"};
763 fHistCentStats = new TH2F("fHistCentStats",
764 "Centrality statistics;;Cent percentile",
765 8,-0.5,7.5,220,-5,105);
766 for(Int_t i = 1; i <= 8; i++){
767 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
768 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
770 fOutput->Add(fHistCentStats);
773 if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
775 fhistcentrality=new TH1F("fhistcentrality","referencemultiplicity",30001,-0.5,30000.5);
776 fOutput->Add(fhistcentrality);
779 fhistcentrality=new TH1F("fhistcentrality","centrality",220,-5,105);
780 fOutput->Add(fhistcentrality);
782 if(fCentralityMethod=="MC_b"){
783 fhistImpactParm=new TH1F("fhistImpactParm","Impact_Parameter",300,0,300);
784 fOutput->Add(fhistImpactParm);
785 fhistImpactParmvsMult=new TH2F("fhistImpactParmvsMult","Impact_Parameter_vs_Multiplicity",300,0,300,50001,-0.5,50000.5);
786 fOutput->Add(fhistImpactParmvsMult);
789 TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
790 fHistRefmult = new TH2F("fHistRefmult",
791 "Reference multiplicity",
792 4,-0.5,3.5,10000,0,20000);
793 for(Int_t i = 1; i <= 4; i++){
794 fHistRefmult->GetXaxis()->SetBinLabel(i,gmultName[i-1].Data());
796 fOutput->Add(fHistRefmult);
798 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
799 //TPC vs EQVZERO multiplicity
800 fHistEQVZEROvsTPCmultiplicity = new TH2F("fHistEQVZEROvsTPCmultiplicity","EqVZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
801 fHistEQVZEROvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO multiplicity (a.u.)");
802 fHistEQVZEROvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
803 fOutput->Add(fHistEQVZEROvsTPCmultiplicity);
806 fHistEQVZEROAvsTPCmultiplicity = new TH2F("fHistEQVZEROAvsTPCmultiplicity","EqVZERO_A vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
807 fHistEQVZEROAvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
808 fHistEQVZEROAvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
809 fOutput->Add(fHistEQVZEROAvsTPCmultiplicity);
812 fHistEQVZEROCvsTPCmultiplicity = new TH2F("fHistEQVZEROCvsTPCmultiplicity","EqVZERO_C vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
813 fHistEQVZEROCvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
814 fHistEQVZEROCvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
815 fOutput->Add(fHistEQVZEROCvsTPCmultiplicity);
817 //EQVZERO vs VZERO multiplicity
818 fHistVZEROCvsEQVZEROCmultiplicity = new TH2F("fHistVZEROCvsEQVZEROCmultiplicity","EqVZERO_C vs VZERO_C multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
819 fHistVZEROCvsEQVZEROCmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
820 fHistVZEROCvsEQVZEROCmultiplicity->GetYaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
821 fOutput->Add(fHistVZEROCvsEQVZEROCmultiplicity);
824 fHistVZEROAvsEQVZEROAmultiplicity = new TH2F("fHistVZEROAvsEQVZEROAmultiplicity","EqVZERO_A vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
825 fHistVZEROAvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
826 fHistVZEROAvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
827 fOutput->Add(fHistVZEROAvsEQVZEROAmultiplicity);
830 //VZEROC vs VZEROA multiplicity
831 fHistVZEROCvsVZEROAmultiplicity = new TH2F("fHistVZEROCvsVZEROAmultiplicity","VZERO_C vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
832 fHistVZEROCvsVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
833 fHistVZEROCvsVZEROAmultiplicity->GetYaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
834 fOutput->Add(fHistVZEROCvsVZEROAmultiplicity);
838 //EQVZEROC vs EQVZEROA multiplicity
839 fHistEQVZEROCvsEQVZEROAmultiplicity = new TH2F("fHistEQVZEROCvsEQVZEROAmultiplicity","EqVZERO_C vs EqVZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
840 fHistEQVZEROCvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
841 fHistEQVZEROCvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
842 fOutput->Add(fHistEQVZEROCvsEQVZEROAmultiplicity);
844 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
845 fOutput->Add(fHistVZEROSignal);
849 if(fRequestEventPlane){
852 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
853 fList->Add(fHistPsiMinusPhi);
855 fEventPlanePID = new TH3F("fEventPlanePID",";centrality;eventplane;PID",20,0.0,100.0,4,-0.5,3.5,4,-0.5,3.5);
856 fList->Add(fEventPlanePID);
860 if(fCutConversions || fCutResonances)
862 fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
863 fOutput->Add(fControlConvResoncances);
866 fHistoTPCdEdx = new TH2F("fHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
867 fOutputList->Add(fHistoTPCdEdx);
868 fHistoTOFbeta = new TH2F(Form("fHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
869 fOutputList->Add(fHistoTOFbeta);
871 fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
872 fOutputList->Add(fTPCTOFPion3d);
874 fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
875 fOutputList->Add(fTPCTOFKaon3d);
877 fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
878 fOutputList->Add(fTPCTOFProton3d);
882 fPionPt = new TH1F("fPionPt","p_{T} distribution",200,0.,10.);
883 fOutputList->Add(fPionPt);
884 fPionEta= new TH1F("fPionEta","#eta distribution",360,-1.8,1.8);
885 fOutputList->Add(fPionEta);
886 fPionPhi = new TH1F("fPionPhi","#phi distribution",340,0,6.8);
887 fOutputList->Add(fPionPhi);
889 fKaonPt = new TH1F("fKaonPt","p_{T} distribution",200,0.,10.);
890 fOutputList->Add(fKaonPt);
891 fKaonEta= new TH1F("fKaonEta","#eta distribution",360,-1.8,1.8);
892 fOutputList->Add(fKaonEta);
893 fKaonPhi = new TH1F("fKaonPhi","#phi distribution",340,0,6.8);
894 fOutputList->Add(fKaonPhi);
896 fProtonPt = new TH1F("fProtonPt","p_{T} distribution",200,0.,10.);
897 fOutputList->Add(fProtonPt);
898 fProtonEta= new TH1F("fProtonEta","#eta distribution",360,-1.8,1.8);
899 fOutputList->Add(fProtonEta);
900 fProtonPhi= new TH1F("fProtonPhi","#phi distribution",340,0,6.8);
901 fOutputList->Add(fProtonPhi);
904 fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
905 fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
906 fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
907 fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);
908 fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
909 fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
910 fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
911 fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
912 fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
913 fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
914 fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
915 fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
916 fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
917 fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
918 fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
919 fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
921 for(Int_t i = 0; i < 16; i++)
923 fOutput->Add(fHistQA[i]);
926 fPriHistShare = new TH1F ("fPriHistShare","Shared clusters, primaries;#shared clusters;counts",160,0,160);
927 fOutput->Add(fPriHistShare);
929 Int_t eventplaneaxis=0;
931 if (fRequestEventPlane) eventplaneaxis=1;
933 kTrackVariablesPair=6+SetChargeAxis+eventplaneaxis;
935 if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
937 if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
939 if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8+SetChargeAxis+eventplaneaxis;
942 // two particle histograms
943 Int_t anaSteps = 1; // analysis steps
944 const char* title = "d^{2}N_{ch}/d#varphid#eta";
946 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
947 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
948 TString* axisTitlePair; // axis titles for track variables
949 axisTitlePair=new TString[kTrackVariablesPair];
951 TString defaultBinningStr;
952 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"
953 "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"
954 "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
955 "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"
956 "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
957 "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
958 "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"
959 "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n"
960 "multiplicity_mixing: 0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1\n";
964 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";
966 if(fRequestEventPlane){
967 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))
969 if(fRequestEventPlanemixing){
970 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";
973 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
976 defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
979 if(SetChargeAxis==2){
980 defaultBinningStr += "TrigCharge: -2.0,0.0,2.0\n"; // course
981 defaultBinningStr += "AssoCharge: -2.0,0.0,2.0\n"; // course
983 // =========================================================
984 // Customization (adopted from AliUEHistograms)
985 // =========================================================
987 TObjArray* lines = defaultBinningStr.Tokenize("\n");
988 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
990 TString line(lines->At(i)->GetName());
991 TString tag = line(0, line.Index(":")+1);
992 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
993 fBinningString += line + "\n";
995 AliInfo(Form("Using custom binning for %s", tag.Data()));
998 fBinningString += fCustomBinning;
1000 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
1002 // =========================================================
1004 // =========================================================
1006 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
1007 axisTitlePair[0] = " multiplicity";
1009 dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
1010 axisTitlePair[1] = "v_{Z} (cm)";
1012 dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
1013 axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
1015 dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
1016 axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
1018 dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
1019 axisTitlePair[4] = "#Delta#eta";
1021 dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
1022 axisTitlePair[5] = "#Delta#varphi (rad)";
1026 if(fRequestEventPlane){
1027 dBinsPair[dim_val] = GetBinning(fBinningString, "eventPlane", iBinPair[dim_val]);
1028 axisTitlePair[dim_val] = "#varphi - #Psi_{2} (a.u.)";
1032 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
1033 dBinsPair[dim_val] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val]);
1034 axisTitlePair[dim_val] = "TrigCharge";
1036 dBinsPair[dim_val+1] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+1]);
1037 axisTitlePair[dim_val+1] = "AssoCharge";
1040 if(fcontainPIDtrig && !fcontainPIDasso){
1042 dBinsPair[dim_val] = GetBinning(fBinningString, "InvariantMass", iBinPair[dim_val]);
1043 axisTitlePair[dim_val] = "InvariantMass";
1046 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
1047 axisTitlePair[dim_val] = "PIDTrig";
1049 if(SetChargeAxis==2){
1050 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
1051 axisTitlePair[dim_val+1] = "TrigCharge";
1053 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
1054 axisTitlePair[dim_val+2] = "AssoCharge";
1058 if(!fcontainPIDtrig && fcontainPIDasso){
1059 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val]);
1060 axisTitlePair[dim_val] = "PIDAsso";
1062 if(SetChargeAxis==2){
1063 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
1064 axisTitlePair[dim_val+1] = "TrigCharge";
1066 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
1067 axisTitlePair[dim_val+2] = "AssoCharge";
1071 if(fcontainPIDtrig && fcontainPIDasso){
1074 dBinsPair[dim_val] = GetBinning(fBinningString, "InvariantMass", iBinPair[dim_val]);
1075 axisTitlePair[dim_val] = "InvariantMass";
1078 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
1079 axisTitlePair[dim_val] = "PIDTrig";
1082 dBinsPair[dim_val+1] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val+1]);
1083 axisTitlePair[dim_val+1] = "PIDAsso";
1085 if(SetChargeAxis==2){
1086 dBinsPair[dim_val+2] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+2]);
1087 axisTitlePair[dim_val+2] = "TrigCharge";
1089 dBinsPair[dim_val+3] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+3]);
1090 axisTitlePair[dim_val+3] = "AssoCharge";
1095 Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
1097 Int_t nPteffbin = -1;
1098 Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
1100 Int_t multmixbin = -1;
1101 Double_t* multmix = GetBinning(fBinningString, "multiplicity_mixing", multmixbin);
1104 //Set the limits from custom binning
1105 fminPtTrig=dBinsPair[2][0];
1106 fmaxPtTrig=dBinsPair[2][iBinPair[2]];
1107 fminPtAsso=dBinsPair[3][0];
1108 fmaxPtAsso=dBinsPair[3][iBinPair[3]];
1109 fmincentmult=dBinsPair[0][0];
1110 fmaxcentmult=dBinsPair[0][iBinPair[0]];
1112 //event pool manager
1113 Int_t MaxNofEvents=1000;
1114 const Int_t NofVrtxBins=10+(1+10)*2;
1115 Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
1116 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
1117 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210};
1120 if(fRequestEventPlanemixing){
1121 // Event plane angle (Psi) bins for event mixing
1124 Double_t* psibins = GetBinning(fBinningString, "eventPlanemixing", nPsiBins);
1125 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,multmixbin,multmix,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1126 if(psibins) delete [] psibins;
1130 const Int_t nPsiBinsd=1;
1131 Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1132 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,multmixbin,multmix,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1135 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1138 AliError("Event Mixing required, but Pool Manager not initialized...");
1142 //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
1143 //fmaxPtComboeff=fmaxPtTrig;
1144 //THnSparses for calculation of efficiency
1146 if((fAnalysisType =="MCAOD") && ffillefficiency) {
1149 effbin[0]=iBinPair[0];
1150 effbin[1]=iBinPair[1];
1151 effbin[2]=nPteffbin;
1153 Int_t effsteps=5;//for each species type::primMCParticles(0),primRecoTracksMatched(1),allRecoTracksMatched(2),primRecoTracksMatchedPID(3),allRecoTracksMatchedPID(4)
1154 for(Int_t jj=0;jj<6;jj++)//PID type binning
1156 if(jj==5) effsteps=3;//for unidentified particles
1157 Histrename="fTrackHistEfficiency";Histrename+=jj;
1158 fTrackHistEfficiency[jj] = new AliTHn(Histrename.Data(), "Tracking efficiency", effsteps, 4, effbin);
1159 fTrackHistEfficiency[jj]->SetBinLimits(0, dBinsPair[0]);
1160 fTrackHistEfficiency[jj]->SetVarTitle(0, "Centrality");
1161 fTrackHistEfficiency[jj]->SetBinLimits(1, dBinsPair[1]);
1162 fTrackHistEfficiency[jj]->SetVarTitle(1, "zvtx");
1163 fTrackHistEfficiency[jj]->SetBinLimits(2, Pteff);
1164 fTrackHistEfficiency[jj]->SetVarTitle(2, "p_{T} (GeV/c)");
1165 fTrackHistEfficiency[jj]->SetBinLimits(3, EtaBin);
1166 fTrackHistEfficiency[jj]->SetVarTitle(3, "#eta");
1167 fOutput->Add(fTrackHistEfficiency[jj]);
1171 //AliThns for Correlation plots(data & MC)
1173 if(ffilltrigassoUNID)
1175 fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1176 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1177 fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
1178 fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
1180 fOutput->Add(fTHnCorrUNID);
1182 fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1183 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1184 fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
1185 fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1187 fOutput->Add(fTHnCorrUNIDmix);
1190 if(ffilltrigIDassoID)
1192 fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
1193 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1194 fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
1195 fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
1197 fOutput->Add(fTHnCorrID);
1199 fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1200 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1201 fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
1202 fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
1204 fOutput->Add(fTHnCorrIDmix);
1207 if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
1209 fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1210 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1211 fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
1212 fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
1214 fOutput->Add(fTHnCorrIDUNID);
1217 fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1218 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1219 fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
1220 fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1222 fOutput->Add(fTHnCorrIDUNIDmix);
1227 //ThnSparse for Correlation plots(truth MC)
1228 if(ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
1230 fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
1231 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1232 fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
1233 fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
1235 fOutput->Add(fCorrelatonTruthPrimary);
1238 fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
1239 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1240 fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
1241 fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
1243 fOutput->Add(fCorrelatonTruthPrimarymix);
1246 //binning for trigger no. counting
1249 if(SetChargeAxis==2) ChargeAxis=1;
1252 Int_t dims=3+ChargeAxis+eventplaneaxis;
1253 if(fcontainPIDtrig) dims=4+ChargeAxis+eventplaneaxis;
1254 fBinst= new Int_t[dims];
1255 Double_t* dBinsTrig[dims]; // bins for track variables
1256 TString* axisTitleTrig; // axis titles for track variables
1257 axisTitleTrig=new TString[dims];
1259 for(Int_t i=0; i<3;i++)
1261 fBinst[i]=iBinPair[i];
1262 dBinsTrig[i]=dBinsPair[i];
1263 axisTitleTrig[i]=axisTitlePair[i];
1265 Int_t dim_val_trig=3;
1266 if(fRequestEventPlane){
1267 fBinst[dim_val_trig]=iBinPair[6];//if fRequestEventPlane=TRUE, dim_val already becomes 7.
1268 dBinsTrig[dim_val_trig]=dBinsPair[6];
1269 axisTitleTrig[dim_val_trig]=axisTitlePair[6];
1273 if(!fcontainPIDtrig && !fcontainPIDasso && ChargeAxis==1){
1274 fBinst[dim_val_trig]=iBinPair[dim_val];
1275 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1276 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1279 if(fcontainPIDtrig && !fcontainPIDasso){
1280 fBinst[dim_val_trig]=iBinPair[dim_val];
1281 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1282 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1284 fBinst[dim_val_trig+1]=iBinPair[dim_val+1];
1285 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+1];
1286 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+1];
1290 if(!fcontainPIDtrig && fcontainPIDasso){
1292 fBinst[dim_val_trig]=iBinPair[dim_val+1];
1293 dBinsTrig[dim_val_trig]=dBinsPair[dim_val+1];
1294 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val+1];
1298 if(fcontainPIDtrig && fcontainPIDasso){
1299 fBinst[dim_val_trig]=iBinPair[dim_val];
1300 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1301 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1303 fBinst[dim_val_trig+1]=iBinPair[dim_val+2];
1304 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+2];
1305 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+2];
1309 //ThSparse for trigger counting(data & reco MC)
1310 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
1312 fTHnTrigcount = new AliTHn("fTHnTrigcount", "fTHnTrigcount", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1313 for(Int_t i=0; i<dims;i++){
1314 fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1315 fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1317 fOutput->Add(fTHnTrigcount);
1320 if(ffilltrigIDassoIDMCTRUTH) {
1321 //AliTHns for trigger counting(truth MC)
1322 fTHnTrigcountMCTruthPrim = new AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1323 for(Int_t i=0; i<dims;i++){
1324 fTHnTrigcountMCTruthPrim->SetBinLimits(i, dBinsTrig[i]);
1325 fTHnTrigcountMCTruthPrim->SetVarTitle(i, axisTitleTrig[i]);
1327 fOutput->Add(fTHnTrigcountMCTruthPrim);
1330 if(fAnalysisType=="MCAOD" || fAnalysisType=="MC"){
1331 if(ffillhistQATruth)
1333 MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
1334 fOutputList->Add(MCtruthpt);
1336 MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
1337 fOutputList->Add(MCtrutheta);
1339 MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
1340 fOutputList->Add(MCtruthphi);
1342 MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
1343 fOutputList->Add(MCtruthpionpt);
1345 MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
1346 fOutputList->Add(MCtruthpioneta);
1348 MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
1349 fOutputList->Add(MCtruthpionphi);
1351 MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
1352 fOutputList->Add(MCtruthkaonpt);
1354 MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
1355 fOutputList->Add(MCtruthkaoneta);
1357 MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
1358 fOutputList->Add(MCtruthkaonphi);
1360 MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
1361 fOutputList->Add(MCtruthprotonpt);
1363 MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
1364 fOutputList->Add(MCtruthprotoneta);
1366 MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
1367 fOutputList->Add(MCtruthprotonphi);
1369 fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
1370 fOutputList->Add(fPioncont);
1372 fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
1373 fOutputList->Add(fKaoncont);
1375 fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
1376 fOutputList->Add(fProtoncont);
1378 fUNIDcont=new TH2F("fUNIDcont","fUNIDcont",10,-0.5,9.5,100,0.,10.);
1379 fOutputList->Add(fUNIDcont);
1382 fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1383 fEventno->GetXaxis()->SetTitle("Centrality");
1384 fEventno->GetYaxis()->SetTitle("Z_Vtx");
1385 fOutput->Add(fEventno);
1386 fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1387 fEventnobaryon->GetXaxis()->SetTitle("Centrality");
1388 fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
1389 fOutput->Add(fEventnobaryon);
1390 fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1391 fEventnomeson->GetXaxis()->SetTitle("Centrality");
1392 fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
1393 fOutput->Add(fEventnomeson);
1395 fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
1396 fOutput->Add(fhistJetTrigestimate);
1398 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);
1399 fOutput->Add(fTwoTrackDistancePtdip);
1401 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);
1402 fOutput->Add(fTwoTrackDistancePtdipmix);
1404 TString Histttrname;
1405 for(Int_t jj=0;jj<2;jj++)// PID type binning
1407 Histttrname="fTwoTrackDistancePt";Histttrname+=jj;
1408 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);
1409 fOutput->Add(fTwoTrackDistancePt[jj]);
1411 Histttrname="fTwoTrackDistancePtmix";Histttrname+=jj;
1412 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);
1413 fOutput->Add(fTwoTrackDistancePtmix[jj]);
1416 //DefineEventPool();
1418 if(fapplyTrigefficiency || fapplyAssoefficiency)
1420 const Int_t nDimt = 4;// cent zvtx pt eta
1421 Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
1422 Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
1423 Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
1425 TString Histrexname;
1426 for(Int_t jj=0;jj<6;jj++)// PID type binning
1428 Histrexname="effcorection";Histrexname+=jj;
1429 effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
1430 effcorection[jj]->Sumw2();
1431 effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
1432 effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
1433 effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
1434 effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
1435 effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
1436 effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1437 effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
1438 effcorection[jj]->GetAxis(3)->SetTitle("#eta");
1439 fOutput->Add(effcorection[jj]);
1441 // TFile *fsifile = new TFile(fefffilename,"READ");
1443 if (TString(fefffilename).BeginsWith("alien:"))
1444 TGrid::Connect("alien:");
1445 TFile *fileT=TFile::Open(fefffilename);
1447 for(Int_t jj=0;jj<6;jj++)//type binning
1449 Nameg="effmap";Nameg+=jj;
1450 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
1451 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
1453 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
1464 //******************************************************************V0 plots*********************************************//
1467 //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};
1468 //const Int_t nbinsV0 =sizeof(BinsV0)/sizeof(Double_t)-1;
1472 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.);
1473 fOutput->Add(fHistRawPtCentInvK0s);
1476 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.);
1477 fOutput->Add(fHistRawPtCentInvLambda);
1480 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.);
1481 fOutput->Add(fHistRawPtCentInvAntiLambda);
1484 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.);
1485 fOutput->Add(fHistFinalPtCentInvK0s);
1488 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.);
1489 fOutput->Add(fHistFinalPtCentInvLambda);
1492 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.);
1493 fOutput->Add(fHistFinalPtCentInvAntiLambda);
1497 //*************************************************************EP plots***********************************************//
1498 if(fRequestEventPlane){
1499 // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
1501 fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
1502 fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
1503 fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
1505 fList->Add(fHResTPCv0A2);
1506 fList->Add(fHResTPCv0C2);
1507 fList->Add(fHResv0Cv0A2);
1510 fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
1511 fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
1512 fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
1514 fList->Add(fHResTPCv0A3);
1515 fList->Add(fHResTPCv0C3);
1516 fList->Add(fHResv0Cv0A3);
1518 // MC as in the dataEP resolution (but using MC tracks)
1519 if(fAnalysisType == "MCAOD" && fV2){
1520 fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
1521 fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
1522 fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
1523 fList->Add(fHResMA2);
1524 fList->Add(fHResMC2);
1525 fList->Add(fHResAC2);
1527 if(fAnalysisType == "MCAOD" && fV3){
1528 fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
1529 fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
1530 fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
1531 fList->Add(fHResMA3);
1532 fList->Add(fHResMC3);
1533 fList->Add(fHResAC3);
1537 // V0A and V0C event plane distributions
1539 fPhiRPTPC = new TH2F("fPhiRPTPCv2","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1540 fPhiRPTPCv3 = new TH2F("fPhiRPTPCv3","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1541 fList->Add(fPhiRPTPC);
1542 fList->Add(fPhiRPTPCv3);
1544 fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1545 fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1546 fList->Add(fPhiRPv0A);
1547 fList->Add(fPhiRPv0C);
1550 fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1551 fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1552 fList->Add(fPhiRPv0Av3);
1553 fList->Add(fPhiRPv0Cv3);
1555 fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth","#phi distribution of EP MCTRUTHheader;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1556 fList->Add(fHistEventPlaneTruth);
1560 //*****************************************************PIDQA histos*****************************************************//
1564 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1565 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1568 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1569 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);
1570 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1571 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1572 fOutputList->Add(fHistoNSigma);
1577 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1578 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1581 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1582 TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
1583 Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1584 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1585 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1586 fOutputList->Add(fHistoNSigma);
1591 if(fPIDType==Bayes){//use bayesianPID
1592 fPIDCombined = new AliPIDCombined();
1593 fPIDCombined->SetDefaultTPCPriors();//****************************************Need to know about it
1595 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1598 TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
1599 Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1600 fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1601 fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
1602 fOutputList->Add(fHistoBayes);
1605 TH2F *fHistoBayesTPC=new TH2F(Form("probBayes_TPC_%d",ipart),
1606 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1607 fHistoBayesTPC->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1608 fHistoBayesTPC->GetYaxis()->SetTitle(Form("Bayes prob TPC %s",kParticleSpeciesName[ipart]));
1609 fOutputList->Add(fHistoBayesTPC);
1611 TH2F *fHistoBayesTOF=new TH2F(Form("probBayes_TOF_%d",ipart),
1612 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1613 fHistoBayesTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1614 fHistoBayesTOF->GetYaxis()->SetTitle(Form("Bayes prob TOF %s",kParticleSpeciesName[ipart]));
1615 fOutputList->Add(fHistoBayesTOF);
1617 TH2F *fHistoBayesTPCTOF=new TH2F(Form("probBayes_TPCTOF_%d",ipart),
1618 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1619 fHistoBayesTPCTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1620 fHistoBayesTPCTOF->GetYaxis()->SetTitle(Form("Bayes prob TPCTOF %s",kParticleSpeciesName[ipart]));
1621 fOutputList->Add(fHistoBayesTPCTOF);
1625 //nsigma separation power plot
1626 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1629 TH2F *Pi_Ka_sep=new TH2F(Form("Pi_Ka_sep_%d",ipid),
1630 Form("Pi_Ka separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1631 Pi_Ka_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1632 Pi_Ka_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1633 fOutputList->Add(Pi_Ka_sep);
1635 TH2F *Pi_Pr_sep=new TH2F(Form("Pi_Pr_sep_%d",ipid),
1636 Form("Pi_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1637 Pi_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1638 Pi_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1639 fOutputList->Add(Pi_Pr_sep);
1641 TH2F *Ka_Pr_sep=new TH2F(Form("Ka_Pr_sep_%d",ipid),
1642 Form("Ka_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1643 Ka_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1644 Ka_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1645 fOutputList->Add(Ka_Pr_sep);
1649 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1650 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1653 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1654 TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
1655 Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1656 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1657 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1658 fOutputList->Add(fHistoNSigma);
1663 if (fAnalysisType == "MCAOD"){
1664 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1665 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1668 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1669 TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
1670 Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1671 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1672 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1673 fOutputList->Add(fHistoNSigma);
1678 for(Int_t idet=0;idet<fNDetectors;idet++){
1679 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1681 if(idet==fTOF)maxy=1.1;
1682 TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
1683 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1684 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1685 fOutputList->Add(fHistoPID);
1688 //PID signal plot, before PID cut
1689 for(Int_t idet=0;idet<fNDetectors;idet++){
1691 if(idet==fTOF)maxy=1.1;
1692 TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
1693 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1694 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1695 fOutputList->Add(fHistoPID);
1698 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
1699 PostData(2, fOutputList);
1700 if(fRequestEventPlane) PostData(3, fList);
1701 AliInfo("Finished setting up the Output");
1703 TH1::AddDirectory(oldStatus);
1705 //-------------------------------------------------------------------------------
1706 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1708 if(fAnalysisType == "AOD") {
1712 }//AOD--analysis-----
1714 else if(fAnalysisType == "MCAOD" || fAnalysisType == "MC") {
1723 //-------------------------------------------------------------------------
1724 void AliTwoParticlePIDCorr::doMCAODevent()
1727 // get the event (for generator level: MCEvent())
1728 AliVEvent* event = NULL;
1729 if(fAnalysisType == "MC") {
1730 event = dynamic_cast<AliVEvent*>(MCEvent());
1733 event = dynamic_cast<AliVEvent*>(InputEvent());
1736 AliError("eventMain not available");
1740 Double_t Inv_mass=0.0;//has no meaning for pions, kaons and protons(just set 0.0) to fill the LRCParticlePID position
1748 gReactionPlane = 999.;
1750 // get centrality object and check quality(valid for p-Pb and Pb-Pb; coming soon for pp 7 TeV)
1751 Double_t cent_v0=-1.0;
1752 Double_t effcent=1.0;
1753 Double_t refmultReco =0.0;
1754 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)
1757 if(fAnalysisType == "MC"){
1759 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1762 AliError("mcEvent not available");
1765 // count all events(physics triggered)
1766 fEventCounter->Fill(1);
1768 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(gMCEvent->GenEventHeader());
1770 TArrayF gVertexArray;
1771 header->PrimaryVertex(gVertexArray);
1772 Float_t zVtxmc =gVertexArray.At(2);
1773 //cout<<"*****************************************************************************************************hi I am here"<<endl;
1776 cent_v0=GetAcceptedEventMultiplicity((AliVEvent*)gMCEvent,kFALSE); //b value; 2nd argument has no meaning
1778 if(cent_v0<0.) return;//mainly returns impact parameter
1780 //get the event plane in case of PbPb
1781 if(fRequestEventPlane){
1782 gReactionPlane=GetEventPlane((AliVEvent*)gMCEvent,kTRUE,cent_v0);//get the truth event plane,middle argument has no meaning in this case
1783 if(gReactionPlane==999.) return;
1786 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
1787 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
1789 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
1790 AliMCParticle* partMC = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
1792 AliError(Form("Could not receive particle %d", iTracks));
1795 //exclude non stable particles
1796 if(fselectprimaryTruth && !(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
1798 //consider only charged particles
1799 if(partMC->Charge() == 0) continue;
1802 //give only kinematic cuts at the generator level
1803 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1804 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
1806 if(!partMC) continue;//for safety
1808 TParticle *particle = partMC->Particle();
1809 if(!particle) continue;
1810 Int_t particletypeTruth=-999;
1812 Int_t pdgtruth = particle->GetPdgCode();
1814 //To determine multiplicity in case of PP
1816 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1817 //only physical primary(all/unidentified)
1818 if(ffillhistQATruth)
1820 MCtruthpt->Fill(partMC->Pt());
1821 MCtrutheta->Fill(partMC->Eta());
1822 MCtruthphi->Fill(partMC->Phi());
1824 if (TMath::Abs(pdgtruth)==211)
1826 particletypeTruth=SpPion;
1827 if(ffillhistQATruth)
1829 MCtruthpionpt->Fill(partMC->Pt());
1830 MCtruthpioneta->Fill(partMC->Eta());
1831 MCtruthpionphi->Fill(partMC->Phi());
1834 if (TMath::Abs(pdgtruth)==321)
1836 particletypeTruth=SpKaon;
1837 if(ffillhistQATruth)
1839 MCtruthkaonpt->Fill(partMC->Pt());
1840 MCtruthkaoneta->Fill(partMC->Eta());
1841 MCtruthkaonphi->Fill(partMC->Phi());
1844 if(TMath::Abs(pdgtruth)==2212)
1846 particletypeTruth=SpProton;
1847 if(ffillhistQATruth)
1849 MCtruthprotonpt->Fill(partMC->Pt());
1850 MCtruthprotoneta->Fill(partMC->Eta());
1851 MCtruthprotonphi->Fill(partMC->Phi());
1854 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)
1856 if(fRequestEventPlane){
1857 FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
1860 //Exclude resonances
1861 if(fExcludeResonancesInMC) {
1862 TParticle *particle = track->Particle();
1863 if(!particle) continue;
1865 Bool_t kExcludeParticle = kFALSE;
1866 Int_t gMotherIndex = particle->GetFirstMother();
1867 if(gMotherIndex != -1) {
1868 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
1870 TParticle *motherParticle = motherTrack->Particle();
1871 if(motherParticle) {
1872 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
1873 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1875 if(pdgCodeOfMother == 113 // rho0
1876 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1877 // || pdgCodeOfMother == 221 // eta
1878 // || pdgCodeOfMother == 331 // eta'
1879 // || pdgCodeOfMother == 223 // omega
1880 // || pdgCodeOfMother == 333 // phi
1881 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1882 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1883 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1884 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1885 || pdgCodeOfMother == 111 // pi0 Dalitz
1887 kExcludeParticle = kTRUE;
1893 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1894 if(kExcludeParticle) continue;
1897 //Exclude electrons with PDG
1898 if(fExcludeElectronsInMC) {
1900 TParticle *particle = track->Particle();
1903 if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
1908 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1909 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1911 Short_t chargeval=0;
1912 if(partMC->Charge()>0) chargeval=1;
1913 if(partMC->Charge()<0) chargeval=-1;
1914 if(chargeval==0) continue;
1915 const TBits *clustermap=0;
1916 const TBits *sharemap=0;
1917 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,Inv_mass,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth,clustermap,sharemap);
1918 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1919 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1920 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1924 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??????
1926 if (fRandomizeReactionPlane)//only for TRuth MC??
1928 Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
1929 Double_t angle = TMath::TwoPi() * centralityDigits;
1930 AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
1931 ShiftTracks(tracksMCtruth, angle);
1935 Float_t weghtval=1.0;
1938 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1940 //Fill Correlations for MC truth particles(same event)
1941 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1942 Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1945 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
1946 if (pool2 && pool2->IsReady())
1947 {//start mixing only when pool->IsReady
1948 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1949 {//proceed only when no. of trigger particles >0 in current event
1950 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
1951 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1952 { //pool event loop start
1953 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1954 if(!bgTracks6) continue;
1955 Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1957 }// pool event loop ends mixing case
1958 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1959 } //if pool->IsReady() condition ends mixing case
1961 //still in main event loop
1964 if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1968 //still in main event loop
1970 if(tracksMCtruth) delete tracksMCtruth;
1975 //"MC" type analysis is finished but still in event loop
1977 else{//if(fAnalysisType=="MCAOD")
1979 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1981 AliError("Cannot get the AOD event");
1985 // count all events(physics triggered)
1986 fEventCounter->Fill(1);
1990 //check the PIDResponse handler
1991 fPID = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
1992 if (!fPID) AliFatal("This Task needs the PID response attached to the inputHandler");
1993 //if (!fPID) return;
1994 // get mag. field required for twotrack efficiency cut
1996 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1998 //check for TClonesArray(truth track MC information)
1999 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2001 AliFatal("Error: MC particles branch not found!\n");
2005 //check for AliAODMCHeader(truth event MC information)
2006 AliAODMCHeader *header=NULL;
2007 header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
2009 printf("MC header branch not found!\n");
2013 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
2014 Float_t zVtxmc =header->GetVtxZ();
2015 if(TMath::Abs(zVtxmc)>fzvtxcut) return;
2017 // For productions with injected signals, figure out above which label to skip particles/tracks
2019 if (fInjectedSignals)
2021 AliGenEventHeader* eventHeader = 0;
2026 AliFatal("fInjectedSignals set but no MC header found");
2028 headers = header->GetNCocktailHeaders();
2029 eventHeader = header->GetCocktailHeader(0);
2033 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
2034 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
2035 AliError("First event header not found. Skipping this event.");
2036 //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
2039 skipParticlesAbove = eventHeader->NProduced();
2040 AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
2043 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
2045 //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
2046 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
2047 refmultReco = GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE);
2048 if(refmultTruth<=0 || refmultReco<=0) return;
2049 cent_v0=refmultTruth;
2052 cent_v0=GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE); //centrality value; 2nd argument has no meaning
2055 if(cent_v0<0.) return;
2056 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
2058 //get the event plane in case of PbPb
2059 if(fRequestEventPlane){
2060 gReactionPlane=GetEventPlane((AliVEvent*)aod,kTRUE,cent_v0);//get the truth event plane
2061 if(gReactionPlane==999.) return;
2064 //TObjArray* tracksMCtruth=0;
2065 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
2066 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
2070 //There is a small difference between zvtx and zVtxmc??????
2071 //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
2072 //cout<<"***********************************************zvtx="<<zvtx<<endl;
2074 //now process the truth particles(for both efficiency & correlation function)
2075 Int_t nMCTrack = fArrayMC->GetEntriesFast();
2077 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
2078 { //MC truth track loop starts
2080 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
2083 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
2087 //consider only charged particles
2088 if(partMC->Charge() == 0) continue;
2090 //consider only primary particles; neglect all secondary particles including from weak decays
2091 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
2094 //remove injected signals(primaries above <maxLabel>)
2095 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
2098 Bool_t isduplicate=kFALSE;
2099 if (fRemoveDuplicates)
2101 for (Int_t j=iMC+1; j<nMCTrack; ++j)
2102 {//2nd trutuh loop starts
2103 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
2105 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
2108 if (partMC->GetLabel() == partMC2->GetLabel())
2113 }//2nd truth loop ends
2115 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
2117 //give only kinematic cuts at the truth level
2118 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
2119 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
2121 if(!partMC) continue;//for safety
2123 //To determine multiplicity in case of PP
2125 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
2126 //only physical primary(all/unidentified)
2127 if(ffillhistQATruth)
2129 MCtruthpt->Fill(partMC->Pt());
2130 MCtrutheta->Fill(partMC->Eta());
2131 MCtruthphi->Fill(partMC->Phi());
2134 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
2135 Int_t particletypeTruth=-999;
2136 if (TMath::Abs(pdgtruth)==211)
2138 particletypeTruth=SpPion;
2139 if(ffillhistQATruth)
2142 MCtruthpionpt->Fill(partMC->Pt());
2143 MCtruthpioneta->Fill(partMC->Eta());
2144 MCtruthpionphi->Fill(partMC->Phi());
2147 if (TMath::Abs(pdgtruth)==321)
2149 particletypeTruth=SpKaon;
2150 if(ffillhistQATruth)
2152 MCtruthkaonpt->Fill(partMC->Pt());
2153 MCtruthkaoneta->Fill(partMC->Eta());
2154 MCtruthkaonphi->Fill(partMC->Phi());
2157 if(TMath::Abs(pdgtruth)==2212)
2159 particletypeTruth=SpProton;
2160 if(ffillhistQATruth)
2162 MCtruthprotonpt->Fill(partMC->Pt());
2163 MCtruthprotoneta->Fill(partMC->Eta());
2164 MCtruthprotonphi->Fill(partMC->Phi());
2167 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)
2169 if(fRequestEventPlane){
2170 FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
2173 // -- Fill THnSparse for efficiency and contamination calculation
2174 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
2176 Double_t primmctruth[4] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
2179 fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
2180 if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for primary truth mesons(3)
2181 if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for primary truth kaons+protons(4)
2182 if (TMath::Abs(pdgtruth)==211) fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
2183 if (TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
2184 if(TMath::Abs(pdgtruth)==2212) fTrackHistEfficiency[2]->Fill(primmctruth,0);//for protons
2187 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
2188 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2190 Short_t chargeval=0;
2191 if(partMC->Charge()>0) chargeval=1;
2192 if(partMC->Charge()<0) chargeval=-1;
2193 if(chargeval==0) continue;
2194 const TBits *clustermap=0;
2195 const TBits *sharemap=0;
2196 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,Inv_mass,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth,clustermap,sharemap);
2197 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
2198 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
2199 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
2201 }//MC truth track loop ends
2203 //*********************still in event loop
2205 if (fRandomizeReactionPlane)//only for TRuth MC??
2207 Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
2208 Double_t angle = TMath::TwoPi() * centralityDigits;
2209 AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
2210 ShiftTracks(tracksMCtruth, angle);
2214 Float_t weghtval=1.0;
2216 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
2219 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??????
2221 //Fill Correlations for MC truth particles(same event)
2222 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
2223 Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
2226 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
2227 if (pool2 && pool2->IsReady())
2228 {//start mixing only when pool->IsReady
2229 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
2230 {//proceed only when no. of trigger particles >0 in current event
2231 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
2232 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
2233 { //pool event loop start
2234 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
2235 if(!bgTracks6) continue;
2236 Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
2238 }// pool event loop ends mixing case
2239 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
2240 } //if pool->IsReady() condition ends mixing case
2242 //still in main event loop
2245 if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
2249 //still in main event loop
2251 if(tracksMCtruth) delete tracksMCtruth;
2253 //now deal with reco tracks
2256 //Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
2257 Float_t bSign1=aod->GetMagneticField() ;//used for reconstructed track dca cut
2260 //detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
2261 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) cent_v0=refmultReco;
2262 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
2264 if(fRequestEventPlane){
2265 gReactionPlane = GetEventPlane((AliVEvent*)aod,kFALSE,cent_v0);//get the reconstructed event plane
2266 if(gReactionPlane==999.) return;
2270 TExMap *trackMap = new TExMap();
2272 // --- track loop for mapping matrix
2275 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2276 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2277 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2278 if (!track) continue;
2279 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
2280 if(tracktype!=1) continue;
2282 if(!track) continue;//for safety
2284 Int_t gid = track->GetID();
2285 trackMap->Add(gid,itrk);
2291 //TObjArray* tracksUNID=0;
2292 TObjArray* tracksUNID = new TObjArray;
2293 tracksUNID->SetOwner(kTRUE);
2295 //TObjArray* tracksID=0;
2296 TObjArray* tracksID = new TObjArray;
2297 tracksID->SetOwner(kTRUE);
2300 //get the selected v0 particle TObjArray
2301 TObjArray* tracksIDV0 = new TObjArray;
2302 tracksIDV0->SetOwner(kTRUE);
2304 Double_t trackscount=0.0;
2305 // loop over reconstructed tracks
2306 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2307 { // reconstructed track loop starts
2308 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2309 if (!track) continue;
2310 //get the corresponding MC track at the truth level (doing reco matching)
2311 AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
2312 if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
2314 //remove injected signals
2315 if(fInjectedSignals)
2317 AliAODMCParticle* mother = recomatched;
2319 while (!mother->IsPhysicalPrimary())
2320 {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
2321 if (mother->GetMother() < 0)
2327 mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
2333 Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
2336 if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
2337 }//remove injected signals
2339 if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
2341 Bool_t isduplicate2=kFALSE;
2342 if (fRemoveDuplicates)
2344 for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
2346 AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
2347 if (!track2) continue;
2348 AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
2349 if(!recomatched2) continue;
2351 if (track->GetLabel() == track2->GetLabel())
2358 if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
2360 fHistQA[11]->Fill(track->GetTPCNcls());
2361 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2363 if(tracktype==0) continue;
2364 if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
2366 if(!track) continue;//for safety
2367 //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
2370 if(IsTrackFromV0(aod,track)) continue;// remove auto correlation
2373 AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2375 if(fFilterBit==128){
2376 Int_t gid1 = track->GetID();
2377 //if(gid1>=0) PIDtrack = track;
2378 PIDtrack =(AliAODTrack*) aod->GetTrack(trackMap->GetValue(-1-gid1));
2379 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2384 //check for eta , phi holes
2385 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2386 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2388 Int_t particletypeMC=-9999;
2390 //tag all particles as unidentified
2391 particletypeMC=unidentified;
2393 Float_t effmatrix=1.;
2395 // -- Fill THnSparse for efficiency calculation
2396 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
2397 //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)
2399 //Clone & Reduce track list(TObjArray) for unidentified particles
2400 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2402 Short_t chargeval=0;
2403 if(track->Charge()>0) chargeval=1;
2404 if(track->Charge()<0) chargeval=-1;
2405 if(chargeval==0) continue;
2406 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2407 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
2408 LRCParticlePID* copy = new LRCParticlePID(particletypeMC,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2409 copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
2410 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2413 //get the pdg code of the corresponding truth particle
2414 Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
2416 Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
2417 if(ffillefficiency) {
2418 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
2419 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
2420 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
2421 if(TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions
2422 if(TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
2423 if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
2425 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
2426 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
2427 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
2428 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
2429 if( TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions
2430 if( TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
2431 if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
2435 switch(TMath::Abs(pdgCode)){
2437 if(fFIllPIDQAHistos){
2438 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2439 if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2440 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
2441 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
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",SpKaon,ipid));
2450 h->Fill(track->Pt(),fnsigmas[SpKaon][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",SpPion,ipid));
2459 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
2466 //now start the particle identification process:)
2468 Float_t dEdx = PIDtrack->GetTPCsignal();
2469 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2471 if(HasTOFPID(PIDtrack))
2473 Double_t beta = GetBeta(PIDtrack);
2474 fHistoTOFbeta->Fill(track->Pt(), beta);
2477 //remove the tracks which don't have proper TOF response-otherwise the misIDentification rate values will be wrong
2478 if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && (!HasTOFPID(PIDtrack)) ) continue;
2481 //do track identification(nsigma method)
2482 particletypeMC=GetParticle(PIDtrack,fFIllPIDQAHistos);//******************************problem is here
2484 //2-d TPCTOF map(for each Pt interval)
2485 if(HasTOFPID(PIDtrack)){
2486 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2487 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2488 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2491 //Pt, Eta , Phi distribution of the reconstructed identified particles
2494 if (particletypeMC==SpPion)
2496 fPionPt->Fill(track->Pt());
2497 fPionEta->Fill(track->Eta());
2498 fPionPhi->Fill(track->Phi());
2500 if (particletypeMC==SpKaon)
2502 fKaonPt->Fill(track->Pt());
2503 fKaonEta->Fill(track->Eta());
2504 fKaonPhi->Fill(track->Phi());
2506 if (particletypeMC==SpProton)
2508 fProtonPt->Fill(track->Pt());
2509 fProtonEta->Fill(track->Eta());
2510 fProtonPhi->Fill(track->Phi());
2514 //fill tracking efficiency
2517 if(particletypeMC==SpPion || particletypeMC==SpKaon)
2519 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) {
2520 fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
2521 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
2524 if(particletypeMC==SpKaon || particletypeMC==SpProton)
2526 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) {
2527 fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
2528 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
2531 if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) {
2532 fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
2533 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
2535 if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
2536 fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
2537 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
2539 if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
2540 fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
2541 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
2546 //for misidentification fraction calculation(do it with tuneonPID)
2547 if(particletypeMC==SpPion )
2549 if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
2550 if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
2551 if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
2552 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
2554 if(particletypeMC==SpKaon )
2556 if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
2557 if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
2558 if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
2559 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
2561 if(particletypeMC==SpProton )
2563 if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
2564 if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
2565 if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
2566 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
2568 if(particletypeMC==SpUndefined )//these undefined are not due to absence of proper TOF response, rather due to the PID method only
2570 if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
2571 if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
2572 if(TMath::Abs(pdgCode)==2212) fUNIDcont->Fill(5.,track->Pt());
2573 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fUNIDcont->Fill(7.,track->Pt());
2576 if(particletypeMC==SpUndefined) continue;
2578 if(fRequestEventPlane){
2579 FillPIDEventPlane(cent_v0,particletypeMC,track->Phi(),gReactionPlane);
2582 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2584 Short_t chargeval=0;
2585 if(track->Charge()>0) chargeval=1;
2586 if(track->Charge()<0) chargeval=-1;
2587 if(chargeval==0) continue;
2588 if (fapplyTrigefficiency || fapplyAssoefficiency)
2589 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
2590 LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2591 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2592 tracksID->Add(copy1);
2594 }// if(tracktype==1) condition structure ands
2596 }//reco track loop ends
2598 //*************************************************************still in event loop
2603 //fill the centrality/multiplicity distribution of the selected events
2604 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2606 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??????
2608 //count selected events having centrality betn 0-100%
2609 fEventCounter->Fill(13);
2611 //***************************************event no. counting
2612 Bool_t isbaryontrig=kFALSE;
2613 Bool_t ismesontrig=kFALSE;
2614 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2616 if(tracksID && tracksID->GetEntriesFast()>0)
2618 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2619 { //trigger loop starts
2620 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2622 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2623 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2624 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2625 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2627 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2628 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2633 tracksIDV0=GetV0Particles((AliVEvent*) aod, cent_v0);
2634 if(tracksIDV0->GetEntriesFast()<=0) return;
2636 //same event delte-eta, delta-phi plot
2637 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2638 {//same event calculation starts
2639 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2640 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)
2643 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2644 {//same event calculation starts
2645 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) {
2646 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2648 else Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2650 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2653 //still in main event loop
2655 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2656 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2657 if (pool && pool->IsReady())
2658 {//start mixing only when pool->IsReady
2659 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2660 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2661 { //pool event loop start
2662 TObjArray* bgTracks = pool->GetEvent(jMix);
2663 if(!bgTracks) continue;
2664 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2665 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2666 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2668 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2670 else Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2672 }// pool event loop ends mixing case
2674 } //if pool->IsReady() condition ends mixing case
2677 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2679 }//mixing with unidentified particles
2681 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2682 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2683 if (pool1 && pool1->IsReady())
2684 {//start mixing only when pool->IsReady
2685 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2686 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2687 { //pool event loop start
2688 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2689 if(!bgTracks2) continue;
2690 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2691 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2692 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2693 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2695 }// pool event loop ends mixing case
2696 } //if pool1->IsReady() condition ends mixing case
2700 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2702 }//mixing with identified particles
2704 //no. of events analyzed
2705 fEventCounter->Fill(15);
2708 if(tracksUNID) delete tracksUNID;
2710 if(tracksID) delete tracksID;
2712 if(tracksIDV0) delete tracksIDV0;
2716 }//AOD || MCAOD condition
2718 //still in the main event loop
2721 //________________________________________________________________________
2722 void AliTwoParticlePIDCorr::doAODevent()
2726 AliVEvent *event = InputEvent();
2727 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
2728 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
2730 AliError("Cannot get the AOD event");
2733 Double_t Inv_mass=0.0;
2736 fEventCounter->Fill(1);
2738 fPID = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
2739 if (!fPID) AliFatal("This Task needs the PID response attached to the inputHandler");
2740 //if (!fPID) return;//this should be available with each event even if we don't do PID selection
2748 gReactionPlane = 999.;
2750 Double_t cent_v0= -999.;
2751 Double_t effcent=1.0;
2753 Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
2756 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
2757 Float_t bSign1=aod->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
2760 // check event cuts and fill event histograms and return the centrality or reference multiplicity value
2761 if((cent_v0 = GetAcceptedEventMultiplicity((AliVEvent*)aod,kFALSE)) < 0){
2764 effcent=cent_v0;//required for efficiency correction case********Extremely Important
2765 //get the event plane in case of PbPb
2766 if(fRequestEventPlane){
2767 gReactionPlane = GetEventPlane((AliVEvent*)aod,kFALSE,cent_v0);
2768 if(gReactionPlane==999.) return;
2772 TExMap *trackMap = new TExMap();
2773 // --- track loop for mapping matrix
2776 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2777 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2778 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2779 if (!track) continue;
2780 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
2781 if(tracktype!=1) continue;
2783 if(!track) continue;//for safety
2785 Int_t gid = track->GetID();
2786 trackMap->Add(gid,itrk);
2790 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
2791 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
2793 TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
2794 tracksID->SetOwner(kTRUE); // IMPORTANT!
2796 //get the selected v0 particle TObjArray
2797 TObjArray* tracksIDV0= new TObjArray;
2798 tracksIDV0->SetOwner(kTRUE); // IMPORTANT!
2802 Bool_t fTrigPtmin1=kFALSE;
2803 Bool_t fTrigPtmin2=kFALSE;
2804 Bool_t fTrigPtJet=kFALSE;
2806 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2807 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2808 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2809 if (!track) continue;
2810 fHistQA[11]->Fill(track->GetTPCNcls());
2811 Int_t particletype=-9999;//required for PID filling
2812 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2813 if(tracktype!=1) continue;
2815 if(!track) continue;//for safety
2818 if(IsTrackFromV0(aod,track)) continue;// remove auto correlation
2820 AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2822 if(fFilterBit==128){
2823 Int_t gid1 = track->GetID();
2824 //if(gid1>=0) PIDtrack = track;
2825 PIDtrack =(AliAODTrack*) aod->GetTrack(trackMap->GetValue(-1-gid1));
2826 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2829 //check for eta , phi holes
2830 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2831 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2835 //if no applyefficiency , set the eff factor=1.0
2836 Float_t effmatrix=1.0;
2838 //tag all particles as unidentified that passed filterbit & kinematic cuts
2839 particletype=unidentified;
2841 //To count the no. of tracks having an accepted track in a certain PT(e.g. Jet Pt) range
2842 if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
2843 if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
2844 if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
2847 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
2850 //to reduce memory consumption in pool
2851 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2853 //Clone & Reduce track list(TObjArray) for unidentified particles
2854 Short_t chargeval=0;
2855 if(track->Charge()>0) chargeval=1;
2856 if(track->Charge()<0) chargeval=-1;
2857 if(chargeval==0) continue;
2858 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2859 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
2860 LRCParticlePID* copy = new LRCParticlePID(particletype,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2861 copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2862 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2865 //now start the particle identificaion process:)
2867 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
2869 Float_t dEdx = PIDtrack->GetTPCsignal();
2870 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2872 //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)
2873 if(HasTOFPID(PIDtrack))
2875 Double_t beta = GetBeta(PIDtrack);
2876 fHistoTOFbeta->Fill(track->Pt(), beta);
2879 //remove the tracks which don't have proper TOF response-otherwise the misIDentification rate values will be wrong(in MC)
2880 if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && (!HasTOFPID(PIDtrack)) ) continue;
2882 //track identification(using nsigma method)
2883 particletype=GetParticle(PIDtrack,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
2885 //2-d TPCTOF map(for each Pt interval)
2886 if(HasTOFPID(PIDtrack)){
2887 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2888 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2889 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2892 //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!!!!!
2893 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)!!!!!!!!!!!
2895 if(fRequestEventPlane){
2896 FillPIDEventPlane(cent_v0,particletype,track->Phi(),gReactionPlane);
2900 //Pt, Eta , Phi distribution of the reconstructed identified particles
2903 if (particletype==SpPion)
2905 fPionPt->Fill(track->Pt());
2906 fPionEta->Fill(track->Eta());
2907 fPionPhi->Fill(track->Phi());
2909 if (particletype==SpKaon)
2911 fKaonPt->Fill(track->Pt());
2912 fKaonEta->Fill(track->Eta());
2913 fKaonPhi->Fill(track->Phi());
2915 if (particletype==SpProton)
2917 fProtonPt->Fill(track->Pt());
2918 fProtonEta->Fill(track->Eta());
2919 fProtonPhi->Fill(track->Phi());
2923 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2925 Short_t chargeval=0;
2926 if(track->Charge()>0) chargeval=1;
2927 if(track->Charge()<0) chargeval=-1;
2928 if(chargeval==0) continue;
2929 if (fapplyTrigefficiency || fapplyAssoefficiency)
2930 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
2931 LRCParticlePID* copy1 = new LRCParticlePID(particletype,Inv_mass,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2932 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2933 tracksID->Add(copy1);
2935 } //track loop ends but still in event loop
2937 if(trackscount<1.0){
2938 if(tracksUNID) delete tracksUNID;
2939 if(tracksID) delete tracksID;
2943 if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
2944 if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
2945 if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
2947 Float_t weightval=1.0;
2951 //fill the centrality/multiplicity distribution of the selected events
2952 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2954 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??????
2956 //count selected events having centrality betn 0-100%
2957 fEventCounter->Fill(13);
2959 //***************************************event no. counting
2960 Bool_t isbaryontrig=kFALSE;
2961 Bool_t ismesontrig=kFALSE;
2962 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2964 if(tracksID && tracksID->GetEntriesFast()>0)
2966 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2967 { //trigger loop starts
2968 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2970 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2971 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2972 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2973 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2975 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2976 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2979 //Get the TObjArray of V0 particles
2981 tracksIDV0=GetV0Particles((AliVEvent*) aod,cent_v0);
2982 if(tracksIDV0->GetEntriesFast()<=0) return;
2985 //same event delta-eta-deltaphi plot
2986 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2987 {//same event calculation starts
2988 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2989 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)
2993 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2994 {//same event calculation starts
2995 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID){
2996 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2997 else Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2999 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
3002 //still in main event loop
3006 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
3007 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
3008 if (pool && pool->IsReady())
3009 {//start mixing only when pool->IsReady
3010 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
3011 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
3012 { //pool event loop start
3013 TObjArray* bgTracks = pool->GetEvent(jMix);
3014 if(!bgTracks) continue;
3015 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
3016 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
3017 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
3019 if(fV0TrigCorr) Fillcorrelation(gReactionPlane,tracksIDV0,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
3020 else Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
3022 }// pool event loop ends mixing case
3023 } //if pool->IsReady() condition ends mixing case
3026 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
3028 }//mixing with unidentified particles
3031 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
3032 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
3033 if (pool1 && pool1->IsReady())
3034 {//start mixing only when pool->IsReady
3035 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
3036 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
3037 { //pool event loop start
3038 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
3039 if(!bgTracks2) continue;
3040 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
3041 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
3042 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
3043 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
3045 }// pool event loop ends mixing case
3046 } //if pool1->IsReady() condition ends mixing case
3050 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
3052 }//mixing with identified particles
3055 //no. of events analyzed
3056 fEventCounter->Fill(15);
3058 //still in main event loop
3061 if(tracksUNID) delete tracksUNID;
3063 if(tracksID) delete tracksID;
3066 if(tracksIDV0) delete tracksIDV0;
3069 } // *************************event loop ends******************************************
3070 //_______________________________________________________________________
3071 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
3073 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
3075 TObjArray* tracksClone = new TObjArray;
3076 tracksClone->SetOwner(kTRUE);
3078 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
3080 LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
3081 LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->GetInvMass(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval(),particle->GetTPCPadMap(),particle->GetTPCSharedMap());
3082 copy100->SetUniqueID(particle->GetUniqueID());
3083 tracksClone->Add(copy100);
3089 //--------------------------------------------------------------------------------
3090 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)
3093 //before calling this function check that either trackstrig & tracksasso are available
3095 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
3096 TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
3097 TArrayF eta(input->GetEntriesFast());
3098 for (Int_t i=0; i<input->GetEntriesFast(); i++)
3099 eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
3102 Int_t jmax=trackstrig->GetEntriesFast();
3104 jmax=tracksasso->GetEntriesFast();
3106 // identify K, Lambda candidates and flag those particles
3107 // a TObject bit is used for this
3108 const UInt_t kResonanceDaughterFlag = 1 << 14;
3109 if (fRejectResonanceDaughters > 0)
3111 Double_t resonanceMass = -1;
3112 Double_t massDaughter1 = -1;
3113 Double_t massDaughter2 = -1;
3114 const Double_t interval = 0.02;
3115 switch (fRejectResonanceDaughters)
3117 case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
3118 case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
3119 case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
3120 default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
3123 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
3124 trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
3125 for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
3126 tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
3128 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
3130 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3131 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
3133 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
3135 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
3136 if (trig->IsEqual(asso)) continue;
3138 if (trig->Charge() * asso->Charge() > 0) continue;
3140 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
3142 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
3144 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
3146 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
3148 trig->SetBit(kResonanceDaughterFlag);
3149 asso->SetBit(kResonanceDaughterFlag);
3151 // Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
3158 //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
3160 Float_t TriggerPtMin=fminPtTrig;
3161 Float_t TriggerPtMax=fmaxPtTrig;
3162 Int_t HighestPtTriggerIndx=-99999;
3163 TH1* triggerWeighting = 0;
3165 if(fSelectHighestPtTrig || fWeightPerEvent)//**************add this data member to the constructor
3167 if (fWeightPerEvent)
3170 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID) axis = fTHnTrigcount->GetGrid(0)->GetGrid()->GetAxis(2);
3171 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) axis = fTHnTrigcountMCTruthPrim->GetGrid(0)->GetGrid()->GetAxis(2);
3172 triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
3174 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
3175 { //trigger loop starts(highest Pt trigger selection)
3176 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3178 Float_t trigpt=trig->Pt();
3179 //to avoid overflow qnd underflow
3180 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
3181 Int_t particlepidtrig=trig->getparticle();
3182 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
3184 Float_t trigeta=trig->Eta();
3186 // some optimization
3187 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
3190 if (fOnlyOneEtaSide != 0)
3192 if (fOnlyOneEtaSide * trigeta < 0)
3195 if (fTriggerSelectCharge != 0)
3196 if (trig->Charge() * fTriggerSelectCharge < 0)
3199 if (fRejectResonanceDaughters > 0)
3200 if (trig->TestBit(kResonanceDaughterFlag)) continue;
3202 if(fSelectHighestPtTrig){
3203 if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
3205 HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
3206 TriggerPtMin=trigpt;
3210 if (fWeightPerEvent) triggerWeighting->Fill(trigpt);
3212 }//trigger loop ends(highest Pt trigger selection)
3214 }//******************(fSelectHighestPtTrig || fWeightPerEvent) condition ends
3217 //two particle correlation filling
3218 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
3219 { //trigger loop starts
3220 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
3222 Float_t trigpt=trig->Pt();
3223 //to avoid overflow qnd underflow
3224 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
3225 Double_t ParticlePID_InvMass=0.0;
3226 if(fV0TrigCorr) ParticlePID_InvMass=trig->GetInvMass();
3228 Int_t particlepidtrig=trig->getparticle();
3229 ParticlePID_InvMass=(Double_t) particlepidtrig;
3230 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}//***********************************forks,lam.Alam their PID numbers have no meaning, only their Inv_mass will be stored
3233 Float_t trigeta=trig->Eta();
3235 // some optimization
3236 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
3239 if (fOnlyOneEtaSide != 0)
3241 if (fOnlyOneEtaSide * trigeta < 0)
3244 if (fTriggerSelectCharge != 0)
3245 if (trig->Charge() * fTriggerSelectCharge < 0)
3248 if (fRejectResonanceDaughters > 0)
3249 if (trig->TestBit(kResonanceDaughterFlag)) continue;
3251 if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
3252 if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
3255 Float_t trigphi=trig->Phi();
3256 Float_t trackefftrig=1.0;
3257 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
3259 // Event plane (determine psi bin)
3260 Double_t gPsiMinusPhi = 0.;
3261 Double_t gPsiMinusPhiBin = -10.;
3262 if(fRequestEventPlane){
3263 gPsiMinusPhi = TMath::Abs(trigphi - ReactionPlane);
3264 //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)
3265 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
3266 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
3267 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
3268 gPsiMinusPhiBin = 0.0;
3270 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
3271 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
3272 gPsiMinusPhiBin = 0.0;
3275 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
3276 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
3277 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
3278 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
3279 gPsiMinusPhiBin = 1.0;
3281 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
3282 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
3283 gPsiMinusPhiBin = 2.0;
3286 gPsiMinusPhiBin = 3.0;
3288 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
3291 //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
3294 Int_t eventplaneAxis=0;
3295 if(fRequestEventPlane) eventplaneAxis=1;
3296 if(fcontainPIDtrig && SetChargeAxis==0) dim=4+eventplaneAxis;
3297 if(!fcontainPIDtrig && SetChargeAxis==2) dim=4+eventplaneAxis;
3298 if(fcontainPIDtrig && SetChargeAxis==2) dim=5+eventplaneAxis;
3299 trigval= new Double_t[dim];
3302 trigval[2] = trigpt;
3304 if(fRequestEventPlane){
3305 trigval[3] = gPsiMinusPhiBin;
3306 if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = ParticlePID_InvMass;
3307 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
3308 if(fcontainPIDtrig && SetChargeAxis==2) {
3309 trigval[4] = ParticlePID_InvMass;
3310 trigval[5] = trig->Charge();
3314 if(!fRequestEventPlane){
3315 if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = ParticlePID_InvMass;
3316 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
3317 if(fcontainPIDtrig && SetChargeAxis==2) {
3318 trigval[3] = ParticlePID_InvMass;
3319 trigval[4] = trig->Charge();
3325 if (fWeightPerEvent)
3327 // leads effectively to a filling of one entry per filled trigger particle pT bin
3328 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(trigval[2]);
3329 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3330 trackefftrig *= triggerWeighting->GetBinContent(weightBin);
3334 //trigger particle counting for both same and mixed event case;;;;;step=0->same event case;;;;;step=1->mixed event case;;;;;;
3335 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
3336 if(fillup=="trigassoUNID" ) {
3337 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3338 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3341 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
3342 if(fillup=="trigassoUNID" )
3344 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3345 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3348 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
3349 if(fillup=="trigUNIDassoID")
3351 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3352 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3355 //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
3356 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
3357 if(fillup=="trigIDassoID")
3359 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3360 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3363 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
3364 if(fillup=="trigIDassoUNID" )
3366 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3367 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3370 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
3371 if(fillup=="trigIDassoID")
3373 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
3374 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
3378 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!!!!
3379 if(mixcase==kFALSE) fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig);
3380 if(mixcase==kTRUE && firstTime) fTHnTrigcountMCTruthPrim->Fill(trigval,1,1.0/trackefftrig);
3383 //asso loop starts within trigger loop
3384 for(Int_t j=0;j<jmax;j++)
3386 LRCParticlePID *asso=0;
3388 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
3390 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
3394 //to avoid overflow and underflow
3395 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
3397 //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
3399 if(!tracksasso && j==i) continue;
3401 // 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)
3402 // if (tracksasso && trig->IsEqual(asso)) continue;
3404 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
3407 if (asso->Pt() >= trig->Pt()) continue;
3409 Int_t particlepidasso=asso->getparticle();
3410 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
3413 if (fAssociatedSelectCharge != 0)
3414 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
3416 if (fSelectCharge > 0)
3419 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
3423 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
3429 if (trigeta < 0 && asso->Eta() < trigeta)
3431 if (trigeta > 0 && asso->Eta() > trigeta)
3435 if (fRejectResonanceDaughters > 0)
3436 if (asso->TestBit(kResonanceDaughterFlag))
3438 // Printf("Skipped j=%d", j);
3443 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
3445 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3449 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3451 fControlConvResoncances->Fill(0.0, mass);
3453 if (mass < 0.04*0.04)
3459 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3461 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
3463 const Float_t kK0smass = 0.4976;
3465 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
3467 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
3469 fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
3471 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
3477 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3479 Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
3480 Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3482 const Float_t kLambdaMass = 1.115;
3484 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
3486 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
3488 fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
3490 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3493 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
3495 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3497 fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
3499 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3504 if (twoTrackEfficiencyCut)
3506 // the variables & cuthave been developed by the HBT group
3507 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
3508 Float_t phi1 = trig->Phi();
3509 Float_t pt1 = trig->Pt();
3510 Float_t charge1 = trig->Charge();
3511 Float_t phi2 = asso->Phi();
3512 Float_t pt2 = asso->Pt();
3513 Float_t charge2 = asso->Charge();
3515 Float_t deta= trigeta - eta[j];
3518 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
3521 // check first boundaries to see if is worth to loop and find the minimum
3522 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMinRadius, bSign);
3523 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMaxRadius, bSign);
3525 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
3527 Float_t dphistarminabs = 1e5;
3528 Float_t dphistarmin = 1e5;
3530 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
3532 for (Double_t rad=fTwoTrackCutMinRadius; rad<=fTwoTrackCutMaxRadius; rad+=0.01)
3534 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
3536 Float_t dphistarabs = TMath::Abs(dphistar);
3538 if (dphistarabs < dphistarminabs)
3540 dphistarmin = dphistar;
3541 dphistarminabs = dphistarabs;
3544 if(mixcase==kFALSE) fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3545 if(mixcase==kTRUE) fTwoTrackDistancePtmix[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3547 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
3549 // 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);
3552 if(mixcase==kFALSE) fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3553 if(mixcase==kTRUE) fTwoTrackDistancePtmix[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3559 //pair sharedfraction cut(only between the trig and asso track)
3560 if(fillup!="trigIDassoIDMCTRUTH")//******************************************NOT for TRUTH MC particles
3562 if(fSharedfraction_Pair_cut>=0){
3563 Bool_t passsharedfractionpaircut=CalculateSharedFraction(trig->GetTPCPadMap(),asso->GetTPCPadMap(),trig->GetTPCSharedMap(),asso->GetTPCSharedMap());
3564 if(!passsharedfractionpaircut) continue;
3567 Float_t weightperevent=weight;
3568 Float_t trackeffasso=1.0;
3569 if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
3570 //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
3571 Float_t deleta=trigeta-eta[j];
3572 Float_t delphi=PhiRange(trigphi-asso->Phi());
3574 Float_t delpt=trigpt-asso->Pt();
3575 //fill it with/without two track efficiency cut
3576 if(mixcase==kFALSE) fTwoTrackDistancePtdip->Fill(deleta, delphi, TMath::Abs(delpt));//for same event
3577 if(mixcase==kTRUE) fTwoTrackDistancePtdipmix->Fill(deleta, delphi, TMath::Abs(delpt));//for mixed event
3579 //here get the two particle efficiency correction factor
3580 Float_t effweight=trackefftrig*trackeffasso*weightperevent;
3581 // if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
3583 Int_t dimused=kTrackVariablesPair+eventplaneAxis;
3584 vars= new Double_t[dimused];
3593 if(fRequestEventPlane)
3595 vars[6]=gPsiMinusPhiBin;
3599 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
3600 vars[dimension]=trig->Charge();
3601 vars[dimension+1]=asso->Charge();
3603 if(fcontainPIDtrig && !fcontainPIDasso){
3604 vars[dimension]=ParticlePID_InvMass;
3605 if(SetChargeAxis==2){
3606 vars[dimension+1]=trig->Charge();
3607 vars[dimension+2]=asso->Charge();
3610 if(!fcontainPIDtrig && fcontainPIDasso){
3611 vars[dimension]=particlepidasso;
3612 if(SetChargeAxis==2){
3613 vars[dimension+1]=trig->Charge();
3614 vars[dimension+2]=asso->Charge();
3617 if(fcontainPIDtrig && fcontainPIDasso){
3618 vars[dimension]=ParticlePID_InvMass;
3619 vars[dimension+1]=particlepidasso;
3620 if(SetChargeAxis==2){
3621 vars[dimension+2]=trig->Charge();
3622 vars[dimension+3]=asso->Charge();
3626 if (fWeightPerEvent)
3628 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
3629 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3630 effweight *= triggerWeighting->GetBinContent(weightBin);
3634 //Fill Correlation ThnSparses
3635 if(fillup=="trigassoUNID")
3637 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,0,1.0/effweight);
3638 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
3640 if(fillup=="trigIDassoID")
3642 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,0,1.0/effweight);
3643 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
3645 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
3646 {//in this case effweight should be 1 always
3647 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
3648 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
3650 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
3652 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
3653 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
3659 }//trigger loop ends
3661 if (triggerWeighting)
3663 delete triggerWeighting;
3664 triggerWeighting = 0;
3668 //------------------------------------------------------------------------------------------------
3669 Bool_t AliTwoParticlePIDCorr:: CalculateSharedFraction(const TBits *triggerPadMap,const TBits *assocPadMap,const TBits *triggerShareMap,const TBits *assocShareMap)
3670 {//source code-AliFemtoShareQualityPairCut.cxx
3672 Double_t nofsharedhits=0;
3674 for(UInt_t imap=0;imap< (triggerPadMap->GetNbits() );imap++)
3676 //if they are in same pad
3677 //cout<<triggerPadMap->TestBitNumber(imap)<<" "<< assocPadMap->TestBitNumber(imap)<<endl;
3678 if (triggerPadMap->TestBitNumber(imap) &&
3679 assocPadMap->TestBitNumber(imap))
3682 //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
3683 if (triggerShareMap->TestBitNumber(imap) &&
3684 assocShareMap->TestBitNumber(imap))
3686 //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
3703 //cout<< (triggerPadMap->TestBitNumber(imap) || assocPadMap->TestBitNumber(imap))<<endl;
3704 else if (triggerPadMap->TestBitNumber(imap) ||
3705 assocPadMap->TestBitNumber(imap)) {
3706 // One track has a hit, the other does not
3709 //cout<<"No hits :"<<nofhits<<endl;
3717 Double_t SharedFraction=0.0;
3718 if(nofhits>0) SharedFraction=(nofsharedhits/nofhits);
3720 //cout<<"Fraction shared hits :"<<SharedFraction<<endl;
3722 if(SharedFraction>fSharedfraction_Pair_cut) return kFALSE;
3728 //________________________________________________________________________________________________
3729 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
3731 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
3733 Float_t effvalue=1.;
3735 if(parpid==unidentified)
3737 effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
3738 effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
3739 effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
3740 effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
3741 effvalue=effcorection[5]->GetBinContent(effVars);
3743 if(parpid==SpPion || parpid==SpKaon)
3745 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3747 effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
3748 effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
3749 effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
3750 effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
3751 effvalue=effcorection[3]->GetBinContent(effVars);
3756 effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
3757 effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
3758 effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
3759 effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
3760 effvalue=effcorection[0]->GetBinContent(effVars);
3765 effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
3766 effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
3767 effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
3768 effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
3769 effvalue=effcorection[1]->GetBinContent(effVars);
3774 if(parpid==SpProton)
3776 effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
3777 effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
3778 effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
3779 effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
3780 effvalue=effcorection[2]->GetBinContent(effVars);
3783 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3785 if(parpid==SpProton || parpid==SpKaon)
3787 effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
3788 effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
3789 effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
3790 effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
3791 effvalue=effcorection[4]->GetBinContent(effVars);
3794 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
3795 if(effvalue==0.) effvalue=1.;
3800 //---------------------------------------------------------------------------------
3802 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield, Bool_t fill)
3805 if(!track) return 0;
3806 Bool_t trackOK = track->TestFilterBit(fFilterBit);
3807 if(!trackOK) return 0;
3808 if (fTrackStatus != 0 && !CheckTrack(track)) return 0;
3809 //select only primary traks(for data & reco MC tracks)
3810 if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
3811 if(track->Charge()==0) return 0;
3812 if (fill) fHistQA[12]->Fill(track->GetTPCNcls());
3815 dz = track->ZAtDCA();
3816 if (fill) fHistQA[6]->Fill(dxy);
3817 if (fill) fHistQA[7]->Fill(dz);
3818 Float_t chi2ndf = track->Chi2perNDF();
3819 if (fill) fHistQA[13]->Fill(chi2ndf);
3820 // Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
3821 Float_t nCrossedRowsTPC = track->GetTPCNCrossedRows();
3822 if (fill) fHistQA[14]->Fill(nCrossedRowsTPC);
3823 //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
3824 if (track->GetTPCNclsF()>0) {
3825 Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
3826 if (fill) fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
3829 Float_t pt=track->Pt();
3830 if(pt< fminPt || pt> fmaxPt) return 0;
3831 if(TMath::Abs(track->Eta())> fmaxeta) return 0;
3832 if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
3833 //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required???
3835 if (fdcacut && fDCAXYCut)
3842 AliAODTrack* clone =(AliAODTrack*) track->Clone();
3843 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
3848 // Printf("%f", ((AliAODTrack*)part)->DCA());
3849 // Printf("%f", pos[0]);
3850 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
3854 if (fSharedClusterCut >= 0)
3856 Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
3857 if (frac > fSharedClusterCut)
3861 // Rejects tracks with shared clusters after filling a control histogram
3862 // This overload is used for primaries
3864 // Get the shared maps
3865 const TBits sharedMap = track->GetTPCSharedMap();
3866 // Fill a control histogram
3867 fPriHistShare->Fill(sharedMap.CountBits());
3869 // Reject shared clusters
3870 if (fSharedTPCmapCut >= 0)
3872 if((sharedMap.CountBits()) >= 1) return 0;// Bad track, has too many shared clusters!
3875 if (fill) fHistQA[8]->Fill(pt);
3876 if (fill) fHistQA[9]->Fill(track->Eta());
3877 if (fill) fHistQA[10]->Fill(track->Phi());
3880 //________________________________________________________________________________
3881 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
3883 //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 )
3884 Float_t pt=track->Pt();
3886 //plot the separation power
3888 Double_t bethe[AliPID::kSPECIES]={0.};
3889 Double_t sigma_TPC[AliPID::kSPECIES]={0.};
3891 Double_t Pi_Ka_sep[NSigmaPIDType+1]={0.};
3892 Double_t Pi_Pr_sep[NSigmaPIDType+1]={0.};
3893 Double_t Ka_Pr_sep[NSigmaPIDType+1]={0.};
3896 Double_t ptpc = track->GetTPCmomentum();
3897 Int_t dEdxN = track->GetTPCsignalN();
3898 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
3899 bethe[ipart] = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
3900 //Double_t diff = dEdx - bethe;
3901 sigma_TPC[ipart] = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
3902 //nSigma[ipart] = diff / sigma;
3904 Pi_Ka_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kKaon])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kKaon])/2.0);
3905 Pi_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kProton])/2.0);
3906 Ka_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kKaon]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kKaon]+sigma_TPC[AliPID::kProton])/2.0);
3909 Double_t sigma_TOF[AliPID::kSPECIES]={0.};
3911 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3913 Double_t timei[AliPID::kSPECIES];
3914 track->GetIntegratedTimes(timei);
3915 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { sigma_TOF[ipart]= fPID->GetTOFResponse().GetExpectedSigma(track->P(), timei[ipart], AliPID::ParticleMass(ipart));}
3916 Pi_Ka_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kKaon])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kKaon])/2.0);
3917 Pi_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kProton])/2.0);
3918 Ka_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kKaon]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kKaon]+sigma_TOF[AliPID::kProton])/2.0);
3920 Pi_Ka_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Ka_sep[NSigmaTPC]*Pi_Ka_sep[NSigmaTPC]+Pi_Ka_sep[NSigmaTOF]*Pi_Ka_sep[NSigmaTOF]);
3921 Pi_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Pr_sep[NSigmaTPC]*Pi_Pr_sep[NSigmaTPC]+Pi_Pr_sep[NSigmaTOF]*Pi_Pr_sep[NSigmaTOF]);
3922 Ka_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Ka_Pr_sep[NSigmaTPC]*Ka_Pr_sep[NSigmaTPC]+Ka_Pr_sep[NSigmaTOF]*Ka_Pr_sep[NSigmaTOF]);
3926 //fill separation power histograms
3927 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3929 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3930 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3931 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3932 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3933 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3934 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3936 if(HasTOFPID(track) && pt>fPtTOFPIDmin && ipid!=0){
3937 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3938 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3939 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3940 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3941 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3942 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3947 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
3948 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
3949 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
3950 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
3952 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
3953 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
3955 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3958 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
3959 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
3960 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
3962 nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
3963 nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
3964 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
3970 // if TOF is missing and below fPtTOFPID only the TPC information is used
3971 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
3972 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
3973 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
3977 //set data member fnsigmas
3978 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
3979 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
3980 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
3982 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
3983 fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
3984 fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
3985 fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
3987 //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)
3988 fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
3989 fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
3990 fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
3993 //Fill NSigma SeparationPlot
3994 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3995 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3996 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3997 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
3998 h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
4004 //----------------------------------------------------------------------------
4005 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
4007 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
4008 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
4009 //get the identity of the particle with the minimum Nsigma
4010 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
4013 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
4014 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
4015 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
4018 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
4019 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
4020 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
4022 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
4023 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4024 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4025 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4027 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
4028 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4029 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4030 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4035 if(fdiffPIDcutvalues){
4036 if(track->Pt()<=4) fNSigmaPID=fPIDCutval1;
4037 if(track->Pt()>4 && track->Pt()<=6) fNSigmaPID=fPIDCutval2;
4038 if(track->Pt()>6 && track->Pt()<=8) fNSigmaPID=fPIDCutval3;
4039 if(track->Pt()>8) fNSigmaPID=fPIDCutval4;
4042 // guess the particle based on the smaller nsigma (within fNSigmaPID)
4043 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
4045 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
4046 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)
4048 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4049 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4050 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
4051 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
4056 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
4058 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4059 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4060 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
4061 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
4066 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
4068 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4069 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
4070 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
4071 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
4077 // else, return undefined
4083 //------------------------------------------------------------------------------------------
4084 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
4085 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
4087 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
4089 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
4091 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
4094 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
4096 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
4099 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
4100 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
4101 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
4104 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
4105 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
4106 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
4108 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
4109 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4110 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4111 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4113 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
4114 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
4115 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
4116 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
4120 // Actually the tracks in the overlapping region(in TPC-TOF nSigma plane) will be ignored
4122 if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
4123 if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
4124 if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
4129 //fill NSigma distr for double counting
4130 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4131 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
4132 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
4133 if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)))continue;//not filling TOF and combined if no TOF PID
4134 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
4135 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
4142 return fHasDoubleCounting;
4145 //////////////////////////////////////////////////////////////////////////////////////////////////
4147 Bool_t* AliTwoParticlePIDCorr::GetAllCompatibleIdentitiesNSigma(AliAODTrack * trk,Bool_t FIllQAHistos){
4148 //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
4149 Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
4150 IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
4154 //////////////////////////////////////////////////////////////////////////////////////////////////
4156 UInt_t AliTwoParticlePIDCorr::CalcPIDCombined(AliAODTrack *track, Int_t detMask, Double_t* prob) const{
4158 // Bayesian PID calculation
4160 for(Int_t i=0;i<AliPID::kSPECIES;i++)
4164 fPIDCombined->SetDetectorMask(detMask);
4166 return fPIDCombined->ComputeProbabilities((AliAODTrack*)track, fPID, prob);
4169 //////////////////////////////////////////////////////////////////////////////////////////////////
4171 Int_t AliTwoParticlePIDCorr::GetIDBayes(AliAODTrack * trk, Bool_t FIllQAHistos){
4173 Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
4176 //Filling of Probability histos
4177 Double_t probTPC[AliPID::kSPECIES]={0.};
4178 Double_t probTOF[AliPID::kSPECIES]={0.};
4179 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
4181 UInt_t detUsedTPC = 0;
4182 UInt_t detUsedTOF = 0;
4183 UInt_t detUsedTPCTOF = 0;
4185 //get the TPC probability
4186 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
4187 detUsedTPC = fPIDCombined->ComputeProbabilities(trk, fPID, probTPC);
4188 if(detUsedTPC == AliPIDResponse::kDetTPC)
4190 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4192 TH2F *h=GetHistogram2D(Form("probBayes_TPC_%d",ipart));
4193 if(ipart==0) h->Fill(trk->Pt(),probTPC[AliPID::kPion]);
4194 if(ipart==1) h->Fill(trk->Pt(),probTPC[AliPID::kKaon]);
4195 if(ipart==2) h->Fill(trk->Pt(),probTPC[AliPID::kProton]);
4199 //get the TOF probability
4200 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
4201 detUsedTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTOF);
4202 if(detUsedTOF == AliPIDResponse::kDetTOF)
4204 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4205 TH2F *h=GetHistogram2D(Form("probBayes_TOF_%d",ipart));
4206 if(ipart==0) h->Fill(trk->Pt(),probTOF[AliPID::kPion]);
4207 if(ipart==1) h->Fill(trk->Pt(),probTOF[AliPID::kKaon]);
4208 if(ipart==2) h->Fill(trk->Pt(),probTOF[AliPID::kProton]);
4212 //get the TPC-TOF probability
4213 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
4214 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTPCTOF);
4215 if(detUsedTPCTOF == (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))
4217 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4218 TH2F *h=GetHistogram2D(Form("probBayes_TPCTOF_%d",ipart));
4219 if(ipart==0) h->Fill(trk->Pt(),probTPCTOF[AliPID::kPion]);
4220 if(ipart==1) h->Fill(trk->Pt(),probTPCTOF[AliPID::kKaon]);
4221 if(ipart==2) h->Fill(trk->Pt(),probTPCTOF[AliPID::kProton]);
4226 Double_t probBayes[AliPID::kSPECIES];
4229 if(HasTOFPID(trk) && trk->Pt()>fPtTOFPIDmin){//use TOF information
4230 detUsed = CalcPIDCombined(trk, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
4231 if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return SpUndefined;//check that TPC and TOF are used
4233 detUsed = CalcPIDCombined(trk,AliPIDResponse::kDetTPC, probBayes);
4234 if(detUsed != AliPIDResponse::kDetTPC)return SpUndefined;//check that TPC is used
4237 //the probability has to be normalized to one, we check it
4239 for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
4240 if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
4241 AliFatal("Bayesian probability not normalized to one");
4244 if(fdiffPIDcutvalues){
4245 if(trk->Pt()<=4) fBayesCut=fPIDCutval1;
4246 if(trk->Pt()>4 && trk->Pt()<=6) fBayesCut=fPIDCutval2;
4247 if(trk->Pt()>6 && trk->Pt()<=8) fBayesCut=fPIDCutval3;
4248 if(trk->Pt()>8) fBayesCut=fPIDCutval4;
4252 //probabilities are normalized to one, if the cut is above .5 there is no problem
4253 if(probBayes[AliPID::kPion]>fBayesCut && IDs[SpPion]==1){
4254 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpPion));
4255 h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
4258 else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[SpKaon]==1){
4259 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpKaon));
4260 h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
4263 else if(probBayes[AliPID::kProton]>fBayesCut && IDs[SpProton]==1){
4264 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpProton));
4265 h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
4274 //////////////////////////////////////////////////////////////////////////////////////////////////
4275 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
4276 //return the specie according to the minimum nsigma value
4277 //no double counting, this has to be evaluated using CheckDoubleCounting()
4279 Int_t ID=SpUndefined;
4281 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
4285 if(fPIDType==Bayes){//use bayesianPID
4288 AliFatal("PIDCombined object has to be set in the steering macro");
4291 ID = GetIDBayes(trk,FIllQAHistos);
4293 }else{ //use nsigma PID
4295 ID=FindMinNSigma(trk,FIllQAHistos);
4296 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
4298 HasDC=GetDoubleCounting(trk,FIllQAHistos);
4299 for(Int_t ipart=0;ipart<NSpecies;ipart++){
4300 if(HasDC[ipart]==kTRUE) ID = SpUndefined;
4304 //Fill PID signal plot
4305 if(ID != SpUndefined){
4306 for(Int_t idet=0;idet<fNDetectors;idet++){
4307 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
4308 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
4309 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
4310 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
4313 //Fill PID signal plot without cuts
4314 for(Int_t idet=0;idet<fNDetectors;idet++){
4315 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
4316 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
4317 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
4318 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
4324 //-------------------------------------------------------------------------------------
4326 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
4329 AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
4330 if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
4331 //ULong_t status=track->GetStatus();
4332 //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
4333 //if (track->GetTPCsignal() <= 0.) return kFALSE;
4334 // 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.
4338 //___________________________________________________________
4341 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
4343 // check TOF matched track
4344 //ULong_t status=track->GetStatus();
4345 //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
4346 AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
4347 if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
4348 if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
4349 //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
4350 //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
4351 // if (probMis > 0.01) return kFALSE;
4352 if(fRemoveTracksT0Fill)
4354 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
4355 if (startTimeMask < 0)return kFALSE;
4360 //________________________________________________________________________
4361 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
4363 //it is called only when TOF PID is available
4364 //TOF beta calculation
4365 Double_t tofTime=track->GetTOFsignal();
4367 Double_t c=TMath::C()*1.E-9;// m/ns
4368 Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
4369 Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
4370 tofTime -= startTime; // subtract startTime to the signal
4371 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
4377 Double_t p = track->P();
4378 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
4380 track->GetIntegratedTimes(timei);
4381 return timei[0]/time;
4384 //------------------------------------------------------------------------------------------------------
4386 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)
4388 // calculate inv mass squared
4389 // same can be achieved, but with more computing time with
4390 /*TLorentzVector photon, p1, p2;
4391 p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
4392 p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
4396 Float_t tantheta1 = 1e10;
4398 if (eta1 < -1e-10 || eta1 > 1e-10)
4399 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
4401 Float_t tantheta2 = 1e10;
4402 if (eta2 < -1e-10 || eta2 > 1e-10)
4403 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
4405 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4406 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4408 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 ) ) );
4412 //---------------------------------------------------------------------------------
4414 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)
4416 // calculate inv mass squared approximately
4418 Float_t tantheta1 = 1e10;
4420 if (eta1 < -1e-10 || eta1 > 1e-10)
4422 Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
4423 tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4426 Float_t tantheta2 = 1e10;
4427 if (eta2 < -1e-10 || eta2 > 1e-10)
4429 Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
4430 tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4433 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4434 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4437 Float_t deltaPhi = TMath::Abs(phi1 - phi2);
4438 while (deltaPhi > TMath::TwoPi())
4439 deltaPhi -= TMath::TwoPi();
4440 if (deltaPhi > TMath::Pi())
4441 deltaPhi = TMath::TwoPi() - deltaPhi;
4443 Float_t cosDeltaPhi = 0;
4444 if (deltaPhi < TMath::Pi()/3)
4445 cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
4446 else if (deltaPhi < 2*TMath::Pi()/3)
4447 cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
4449 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
4451 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
4453 // Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
4457 //--------------------------------------------------------------------------------
4458 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)
4461 // calculates dphistar
4464 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
4466 static const Double_t kPi = TMath::Pi();
4469 // if (dphistar > 2 * kPi)
4470 // dphistar -= 2 * kPi;
4471 // if (dphistar < -2 * kPi)
4472 // dphistar += 2 * kPi;
4475 dphistar = kPi * 2 - dphistar;
4476 if (dphistar < -kPi)
4477 dphistar = -kPi * 2 - dphistar;
4478 if (dphistar > kPi) // might look funny but is needed
4479 dphistar = kPi * 2 - dphistar;
4484 //------------------------------------------------------------------------
4485 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
4487 // This method is a copy from AliUEHist::GetBinning
4488 // takes the binning from <configuration> identified by <tag>
4489 // configuration syntax example:
4490 // 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
4493 // returns bin edges which have to be deleted by the caller
4495 TString config(configuration);
4496 TObjArray* lines = config.Tokenize("\n");
4497 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
4499 TString line(lines->At(i)->GetName());
4500 if (line.BeginsWith(TString(tag) + ":"))
4502 line.Remove(0, strlen(tag) + 1);
4503 line.ReplaceAll(" ", "");
4504 TObjArray* binning = line.Tokenize(",");
4505 Double_t* bins = new Double_t[binning->GetEntriesFast()];
4506 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
4507 bins[j] = TString(binning->At(j)->GetName()).Atof();
4509 nBins = binning->GetEntriesFast() - 1;
4518 AliFatal(Form("Tag %s not found in %s", tag, configuration));
4522 //____________________________________________________________________
4524 Bool_t AliTwoParticlePIDCorr::CheckTrack(AliAODTrack * part)
4526 // check if the track status flags are set
4528 UInt_t status=((AliAODTrack*)part)->GetStatus();
4529 if ((status & fTrackStatus) == fTrackStatus)
4533 //________________________________________________________________________
4535 Bool_t AliTwoParticlePIDCorr::AcceptEventCentralityWeight(Double_t centrality)
4537 // rejects "randomly" events such that the centrality gets flat
4538 // uses fCentralityWeights histogram
4540 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
4542 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
4543 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
4545 Bool_t result = kFALSE;
4546 if (centralityDigits < weight)
4549 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
4554 //____________________________________________________________________
4555 void AliTwoParticlePIDCorr::ShiftTracks(TObjArray* tracks, Double_t angle)
4557 // shifts the phi angle of all tracks by angle
4558 // 0 <= angle <= 2pi
4560 for (Int_t i=0; i<tracks->GetEntriesFast(); ++i)
4562 LRCParticlePID *part=(LRCParticlePID*)(tracks->UncheckedAt(i));
4564 Double_t newAngle = part->Phi() + angle;
4565 if (newAngle >= TMath::TwoPi())
4566 newAngle -= TMath::TwoPi();
4568 part->SetPhi(newAngle);
4573 //________________________________________________________________________
4574 void AliTwoParticlePIDCorr::SetVZEROCalibrationFile(const char* filename,const char* lhcPeriod) {
4575 //Function to setup the VZERO gain equalization
4576 //============Get the equilization map============//
4577 TFile *calibrationFile = TFile::Open(filename);
4578 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
4579 Printf("No calibration file found!!!");
4583 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
4585 Printf("Calibration TList not found!!!");
4589 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
4590 if(!fHistVZEROAGainEqualizationMap) {
4591 Printf("VZERO-A calibration object not found!!!");
4594 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
4595 if(!fHistVZEROCGainEqualizationMap) {
4596 Printf("VZERO-C calibration object not found!!!");
4600 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
4601 if(!fHistVZEROChannelGainEqualizationMap) {
4602 Printf("VZERO channel calibration object not found!!!");
4607 //________________________________________________________________________
4608 Double_t AliTwoParticlePIDCorr::GetChannelEqualizationFactor(Int_t run,Int_t channel) {
4610 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4612 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
4613 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4614 if(gRunNumber == run)
4615 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
4621 //________________________________________________________________________
4622 Double_t AliTwoParticlePIDCorr::GetEqualizationFactor(Int_t run, const char* side) {
4624 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4626 TString gVZEROSide = side;
4627 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
4628 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4629 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
4630 if(gRunNumber == run) {
4631 if(gVZEROSide == "A")
4632 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
4633 else if(gVZEROSide == "C")
4634 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
4640 //________________________________________________________________________
4641 Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliVEvent *mainevent){
4642 //Function that returns the reference multiplicity from AODs (data or reco MC, Not for Truth)
4643 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
4644 if(!mainevent) return -1;
4646 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
4648 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
4649 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
4651 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
4653 Printf("ERROR: AOD header not available");
4656 Int_t gRunNumber = header->GetRunNumber();
4657 Float_t bSign1=header->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
4660 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++)
4661 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
4662 AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
4663 if (!track) continue;
4664 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
4665 if(tracktype!=1) continue;
4667 if(!track) continue;//for safety
4669 gRefMultiplicityTPC += 1.0;
4673 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
4674 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
4675 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
4676 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
4679 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
4680 else if(iChannel >= 32)
4681 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
4684 //Equalization of gain
4685 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
4687 gRefMultiplicityVZEROA /= gFactorA;
4688 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
4690 gRefMultiplicityVZEROC /= gFactorC;
4691 if((gFactorA != 0)&&(gFactorC != 0))
4692 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
4695 //EQVZERO vs TPC multiplicity
4696 fHistEQVZEROvsTPCmultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
4697 fHistEQVZEROAvsTPCmultiplicity->Fill(gRefMultiplicityVZEROA,gRefMultiplicityTPC);
4698 fHistEQVZEROCvsTPCmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityTPC);
4700 //EQVZERO vs VZERO multiplicity
4701 fHistVZEROCvsEQVZEROCmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),gRefMultiplicityVZEROC);
4702 fHistVZEROAvsEQVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0A(),gRefMultiplicityVZEROA);
4704 //VZEROC vs VZEROA multiplicity
4705 fHistVZEROCvsVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),event->GetVZEROData()->GetMTotV0A());
4707 //EQVZEROC vs EQVZEROA multiplicity
4708 fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
4710 fHistRefmult->Fill(3.,gRefMultiplicityTPC);
4711 fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
4712 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
4713 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
4716 if(fCentralityMethod == "TRACKS_MANUAL") gRefMultiplicity = gRefMultiplicityTPC;
4718 else if(fCentralityMethod == "V0M_MANUAL") gRefMultiplicity = gRefMultiplicityVZERO;
4720 else if(fCentralityMethod == "V0A_MANUAL") gRefMultiplicity = gRefMultiplicityVZEROA;
4722 else if(fCentralityMethod == "V0C_MANUAL") gRefMultiplicity = gRefMultiplicityVZEROC;
4724 else gRefMultiplicity = gRefMultiplicityTPC;
4726 return gRefMultiplicity;
4729 //-------------------------------------------------------------------------------------------------------
4730 Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliVEvent *mainevent, Bool_t truth){
4732 if(!mainevent) return -1;
4733 // get centrality object and check quality
4734 Double_t cent_v0=-1;
4735 Bool_t shift_to_TRACKS_MANUAL=kFALSE;//in case of wrong setting automatic shift to Tracks_Manual method
4737 Double_t gRefMultiplicityTPC_Truth = 0.;
4738 Double_t gRefMultiplicityVZERO_Truth = 0., gRefMultiplicityVZEROA_Truth = 0., gRefMultiplicityVZEROC_Truth = 0.;
4740 if(fAnalysisType == "AOD"|| fAnalysisType == "MCAOD") { //centrality in AOD header //++++++++++++++
4741 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
4743 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
4746 if(fSampleType=="pp_7" && fPPVsMultUtils)
4747 {//for pp 7 TeV case only using Alianalysisutils class
4748 if(fAnalysisUtils) cent_v0 = fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,fCentralityMethod);
4750 fHistCentStats->Fill(0.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0A"));
4751 fHistCentStats->Fill(1.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0C"));
4752 fHistCentStats->Fill(2.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0M"));
4753 fHistCentStats->Fill(3.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0AEq"));//only available for LHC10d at present (Quantile info)
4754 fHistCentStats->Fill(4.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0CEq"));//only available for LHC10d at present (Quantile info)
4755 fHistCentStats->Fill(5.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0MEq"));//only available for LHC10d at present (Quantile info)
4759 else if(fSampleType=="pPb" || fSampleType=="PbPb")
4761 AliCentrality *centralityObj=0;
4762 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
4763 if(!header) return -1;
4764 centralityObj = header->GetCentralityP();
4765 // if (centrality->GetQuality() != 0) return ;
4767 fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
4768 fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
4769 fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
4770 fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
4771 fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
4772 fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
4774 fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
4775 fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
4777 cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
4781 else shift_to_TRACKS_MANUAL=kTRUE;
4783 }//centralitymethod condition
4785 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
4787 if(!truth){//for data or RecoMC
4788 cent_v0 = GetReferenceMultiplicityVZEROFromAOD((AliVEvent*)event);
4789 }//for data or RecoMC
4791 if(truth && (fAnalysisType == "MCAOD")){//condition for TRUTH case
4792 //check for TClonesArray(truth track MC information)
4793 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
4795 //AliFatal("Error: MC particles branch not found!\n");
4798 //now process the truth particles(for both efficiency & correlation function)
4799 Int_t nMCTrack = fArrayMC->GetEntriesFast();
4801 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
4802 {//MC truth track loop starts
4804 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
4807 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
4811 //consider only charged particles
4812 if(partMC->Charge() == 0) continue;
4814 //consider only primary particles; neglect all secondary particles including from weak decays
4815 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
4818 //remove injected signals(primaries above <maxLabel>)
4819 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
4822 Bool_t isduplicate=kFALSE;
4823 if (fRemoveDuplicates)
4825 for (Int_t j=iMC+1; j<nMCTrack; ++j)
4826 {//2nd trutuh loop starts
4827 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
4829 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
4832 if (partMC->GetLabel() == partMC2->GetLabel())
4837 }//2nd truth loop ends
4839 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
4842 // if (fCentralityMethod=="V0M_MANUAL")
4843 if((partMC->Eta() < 5.1 && partMC->Eta() > 2.8) || (partMC->Eta() > -3.7 && partMC->Eta() < -1.7)) gRefMultiplicityVZERO_Truth+=1;
4844 // else if (fCentralityMethod=="V0A_MANUAL") {
4845 if(partMC->Eta() < 5.1 && partMC->Eta() > 2.8) gRefMultiplicityVZEROA_Truth+=1;
4846 // else if (fCentralityMethod=="V0C_MANUAL") {
4847 if(partMC->Eta() < -1.7 && partMC->Eta() > -3.7) gRefMultiplicityVZEROC_Truth+=1;
4848 //else if (fCentralityMethod=="TRACKS_MANUAL") {
4849 if (partMC->Eta() > fmineta && partMC->Eta() < fmaxeta) {
4850 if (partMC->Pt() > fminPt && partMC->Pt() < fmaxPt) gRefMultiplicityTPC_Truth+=1;
4853 }//truth track loop ends
4855 fHistRefmult->Fill(3.,gRefMultiplicityTPC_Truth);
4856 fHistRefmult->Fill(2.,gRefMultiplicityVZERO_Truth);
4857 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA_Truth);
4858 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC_Truth);
4860 if(fCentralityMethod == "TRACKS_MANUAL") cent_v0=gRefMultiplicityTPC_Truth;
4862 else if(fCentralityMethod == "V0M_MANUAL") cent_v0=gRefMultiplicityVZERO_Truth;
4864 else if(fCentralityMethod == "V0A_MANUAL") cent_v0=gRefMultiplicityVZEROA_Truth;
4866 else if(fCentralityMethod == "V0C_MANUAL") cent_v0=gRefMultiplicityVZEROC_Truth;
4868 else cent_v0=gRefMultiplicityTPC_Truth;
4870 }//condition for TRUTH case
4872 }//end of MANUAL method
4874 else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC in AOD production
4876 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4880 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4883 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
4884 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
4885 AliError("Event header not found. Skipping this event.");
4889 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4893 cent_v0 = collGeometry->ImpactParameter();
4894 fhistImpactParm->Fill(cent_v0);
4897 }//end of Impact parameter method
4901 }//AOD OR MCAOD condition
4904 else if(fAnalysisType == "MC"){
4905 Double_t gImpactParameter = -1.;
4906 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(mainevent);
4908 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
4910 gImpactParameter = headerH->ImpactParameter();
4912 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
4913 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
4915 AliError(Form("Could not receive particle %d", iParticle));
4919 //exclude non stable particles
4920 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
4922 if(track->Charge() == 0) continue;
4924 // if (fCentralityMethod=="V0M_MANUAL")
4925 if((track->Eta() < 5.1 && track->Eta() > 2.8) || (track->Eta() > -3.7 && track->Eta() < -1.7)) gRefMultiplicityVZERO_Truth+=1;
4926 // else if (fCentralityMethod=="V0A_MANUAL") {
4927 if(track->Eta() < 5.1 && track->Eta() > 2.8) gRefMultiplicityVZEROA_Truth+=1;
4928 // else if (fCentralityMethod=="V0C_MANUAL") {
4929 if(track->Eta() < -1.7 && track->Eta() > -3.7) gRefMultiplicityVZEROC_Truth+=1;
4930 //else if (fCentralityMethod=="TRACKS_MANUAL") {
4931 if (track->Eta() > fmineta && track->Eta() < fmaxeta) {
4932 if (track->Pt() > fminPt && track->Pt() < fmaxPt) gRefMultiplicityTPC_Truth+=1;}
4934 }//loop over primaries
4936 fHistRefmult->Fill(3.,gRefMultiplicityTPC_Truth);
4937 fHistRefmult->Fill(2.,gRefMultiplicityVZERO_Truth);
4938 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA_Truth);
4939 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC_Truth);
4940 if (fCentralityMethod == "MC_b"){
4941 cent_v0=gImpactParameter;
4942 fhistImpactParm->Fill(gImpactParameter);
4943 fhistImpactParmvsMult->Fill(gImpactParameter,gRefMultiplicityTPC_Truth);
4946 else if(fCentralityMethod == "TRACKS_MANUAL") cent_v0=gRefMultiplicityTPC_Truth;
4948 else if(fCentralityMethod == "V0M_MANUAL") cent_v0=gRefMultiplicityVZERO_Truth;
4950 else if(fCentralityMethod == "V0A_MANUAL") cent_v0=gRefMultiplicityVZEROA_Truth;
4952 else if(fCentralityMethod == "V0C_MANUAL") cent_v0=gRefMultiplicityVZEROC_Truth;
4954 else cent_v0=gImpactParameter;//default value is the impact parameter
4965 //-----------------------------------------------------------------------------------------
4966 Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliVEvent *event,Bool_t truth){
4967 //do the event selection(zvtx, pileup, centrality/multiplicity cut) and then return the value of the centrality of that event
4968 if(!event) return -1;
4970 Float_t gRefMultiplicity = -1.;
4972 //***********************************SOURCE CODE-TASKBFPsi
4975 if(fAnalysisType == "MC") {
4976 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
4978 AliError("mcEvent not available");
4983 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
4985 TArrayF gVertexArray;
4986 header->PrimaryVertex(gVertexArray);
4987 //count events having a proper vertex
4988 fEventCounter->Fill(5);
4990 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
4992 if(TMath::Abs(gVertexArray.At(0)) < fVxMax_MC) {
4993 if(TMath::Abs(gVertexArray.At(1)) < fVyMax_MC) {
4994 if(TMath::Abs(gVertexArray.At(2)) < fVzMax_MC) {
4995 //count events after vertex cut
4996 fEventCounter->Fill(7);
4997 fHistQA[3]->Fill((gVertexArray.At(0)));fHistQA[4]->Fill((gVertexArray.At(1)));fHistQA[5]->Fill((gVertexArray.At(2)));//after vertex cut,for trkVtx only
4999 // get the reference multiplicty or centrality
5000 gRefMultiplicity = GetRefMultiOrCentrality((AliVEvent*)mcevent,kFALSE);//2nd argument has no meaning
5002 if(gRefMultiplicity<0) return -1;
5004 // take events only within the multiplicity class mentioned in the custom binning
5005 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
5007 //count events having proper centrality/ref multiplicity
5008 fEventCounter->Fill(9);
5017 else if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD"){// if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD"
5018 //vertex selection(is it fine for PP?)
5019 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
5022 // check first event in chunk (is not needed for new reconstructions)
5023 if(fCheckFirstEventInChunk){
5024 AliAnalysisUtils ut;
5025 if(ut.IsFirstEventInChunk(aod))
5030 AliAnalysisUtils ut;
5031 ut.SetUseMVPlpSelection(kTRUE);
5032 ut.SetUseOutOfBunchPileUp(kTRUE);
5033 if(ut.IsPileUpEvent(aod))
5037 //count events after pileup selection
5038 fEventCounter->Fill(3);
5040 if (fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
5041 trkVtx = aod->GetPrimaryVertex();
5042 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
5043 TString vtxTtl = trkVtx->GetTitle();
5044 if (!vtxTtl.Contains("VertexerTracks")) return -1;
5045 zvtx = trkVtx->GetZ();
5046 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
5047 if (!spdVtx || spdVtx->GetNContributors()<=0) return -1;
5048 TString vtxTyp = spdVtx->GetTitle();
5049 Double_t cov[6]={0};
5050 spdVtx->GetCovarianceMatrix(cov);
5051 Double_t zRes = TMath::Sqrt(cov[5]);
5052 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return -1;
5053 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return -1;
5055 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
5056 Int_t nVertex = aod->GetNumberOfVertices();
5058 trkVtx = aod->GetPrimaryVertex();
5059 Int_t nTracksPrim = trkVtx->GetNContributors();
5060 zvtx = trkVtx->GetZ();
5061 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
5062 // Reject TPC only vertex
5063 TString name(trkVtx->GetName());
5064 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return -1;
5066 // Select a quality vertex by number of tracks?
5067 if( nTracksPrim < fnTracksVertex ) {
5068 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
5071 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
5072 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
5074 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
5079 else if(fVertextype==0){//default case
5080 trkVtx =(AliAODVertex*) aod->GetPrimaryVertex();
5081 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
5082 zvtx = trkVtx->GetZ();
5084 trkVtx->GetCovarianceMatrix(fCov);
5085 if(fCov[5] == 0) return -1;//proper vertex resolution
5088 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
5089 return -1;//as there is no proper sample type
5092 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
5094 //count events having a proper vertex
5095 fEventCounter->Fill(5);
5097 if (TMath::Abs(zvtx) > fzvtxcut) return -1;
5099 //count events after vertex cut
5100 fEventCounter->Fill(7);
5103 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
5105 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
5107 //get the centrality or multiplicity
5108 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)
5110 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)
5112 if(gRefMultiplicity<0) return -1;
5114 // take events only within the multiplicity class mentioned in the custom binning
5115 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
5117 //count events having proper centrality/ref multiplicity
5118 fEventCounter->Fill(9);
5121 // centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
5122 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
5124 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
5128 //count events after rejection due to centrality weighting
5129 fEventCounter->Fill(11);
5132 else gRefMultiplicity=-1;
5134 return gRefMultiplicity;
5137 //--------------------------------------------------------------------------------------------------------
5138 Float_t AliTwoParticlePIDCorr::GetEventPlane(AliVEvent *mainevent,Bool_t truth, Double_t v0Centr)
5140 Float_t eventplane=999.;
5141 // Get the event plane
5142 if(!mainevent) return 999.;
5145 //MC: from reaction plane
5146 if(fAnalysisType == "MC"){
5148 AliError("mcEvent not available");
5152 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(mainevent);
5154 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
5157 // Impact parameter bins(it is only for Pb-Pb)
5158 if(v0Centr < 3.50) iC = 0;
5159 else if(v0Centr < 4.94) iC = 1;
5160 else if(v0Centr < 6.98) iC = 2;
5161 else if(v0Centr < 8.55) iC = 3;
5162 else if(v0Centr < 9.88) iC = 4;
5163 else if(v0Centr < 11.04) iC = 5;
5164 else if(v0Centr < 12.09) iC = 6;
5165 else if(v0Centr < 13.05) iC = 7;
5168 eventplane = headerH->ReactionPlaneAngle();
5169 if(eventplane > TMath::Pi()/2 && eventplane <= TMath::Pi()*3/2) eventplane-=TMath::Pi();
5170 if(eventplane > TMath::Pi()*3/2) eventplane-=2*TMath::Pi();
5171 fHistEventPlaneTruth->Fill(iC,eventplane);
5172 //gReactionPlane *= TMath::RadToDeg();
5177 else if(fAnalysisType == "MCAOD" || fAnalysisType == "AOD") {
5178 //reset Q vector info
5180 AliAODEvent* event = dynamic_cast<AliAODEvent*>(mainevent);
5183 Int_t run = event->GetRunNumber();
5186 // Load the calibrations run dependent
5187 if(! fIsAfter2011) OpenInfoCalbration(run);
5193 if (v0Centr > 80) return 999.; // analysis only for 0-80% centrality classes
5195 if(v0Centr < 5) iC = 0;
5196 else if(v0Centr < 10) iC = 1;
5197 else if(v0Centr < 20) iC = 2;
5198 else if(v0Centr < 30) iC = 3;
5199 else if(v0Centr < 40) iC = 4;
5200 else if(v0Centr < 50) iC = 5;
5201 else if(v0Centr < 60) iC = 6;
5202 else if(v0Centr < 70) iC = 7;
5208 //reset Q vector info
5209 Double_t Qxa2 = 0, Qya2 = 0;
5210 Double_t Qxc2 = 0, Qyc2 = 0;
5211 Double_t Qxa3 = 0, Qya3 = 0;
5212 Double_t Qxc3 = 0, Qyc3 = 0;
5215 //MC: from reaction plane
5218 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
5220 evplaneMC = header->GetReactionPlaneAngle();//[0, 360]
5221 //make it within [-pi/2,pi/2] to make it general
5222 if(evplaneMC > TMath::Pi()/2 && evplaneMC <= TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
5223 if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
5224 fHistEventPlaneTruth->Fill(iC,evplaneMC);
5226 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
5230 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
5232 if (collGeometry){//get the reaction plane from MC header
5233 gReactionPlane = collGeometry->ReactionPlaneAngle();//[0,180]
5237 //taken from vnv0 code(get the TPC, V0A, V0C event plane using truth tracks)
5238 TClonesArray *mcArray = NULL;
5239 mcArray = (TClonesArray*)event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
5241 Float_t QxMCv2[3] = {0,0,0};
5242 Float_t QyMCv2[3] = {0,0,0};
5243 Float_t QxMCv3[3] = {0,0,0};
5244 Float_t QyMCv3[3] = {0,0,0};
5245 Float_t EvPlaneMCV2[3] = {0,0,0};
5246 Float_t EvPlaneMCV3[3] = {0,0,0};
5247 Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
5248 Float_t etaMax[3] = {4.88,-1.8,0.8};
5250 // analysis on MC tracks
5251 Int_t nMCtrack = mcArray->GetEntries() ;
5253 // EP computation with MC tracks
5254 for(Int_t iT=0;iT < nMCtrack;iT++){
5255 AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
5256 if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
5258 Float_t eta = mctr->Eta();
5259 for(Int_t iD=0;iD<3;iD++){
5260 if(eta > etaMin[iD] && eta < etaMax[iD]){
5261 Float_t phi = mctr->Phi();
5262 QxMCv2[iD] += TMath::Cos(2*phi);
5263 QyMCv2[iD] += TMath::Sin(2*phi);
5264 QxMCv3[iD] += TMath::Cos(3*phi);
5265 QyMCv3[iD] += TMath::Sin(3*phi);
5270 EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
5271 EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
5272 EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
5273 fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
5274 fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
5275 fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
5276 fgPsi2v0aMC = EvPlaneMCV2[0];
5277 fgPsi2v0cMC = EvPlaneMCV2[1];
5278 fgPsi2tpcMC = EvPlaneMCV2[2];
5281 EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
5282 EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
5283 EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
5284 fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
5285 fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
5286 fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
5287 fgPsi3v0aMC = EvPlaneMCV3[0];
5288 fgPsi3v0cMC = EvPlaneMCV3[1];
5289 fgPsi3tpcMC = EvPlaneMCV3[2];
5296 Int_t nAODTracks = event->GetNumberOfTracks();
5298 // TPC EP needed for resolution studies (TPC subevent)
5299 //AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
5300 //Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
5301 Double_t Qx2 = 0, Qy2 = 0;
5302 Double_t Qx3 = 0, Qy3 = 0;
5304 for(Int_t iT = 0; iT < nAODTracks; iT++) {
5306 AliAODTrack* aodTrack =(AliAODTrack*) event->GetTrack(iT);
5312 Bool_t trkFlag = aodTrack->TestFilterBit(1);
5314 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
5317 Double_t b[2] = {-99., -99.};
5318 Double_t bCov[3] = {-99., -99., -99.};
5321 AliAODTrack param(*aodTrack);
5322 if (!param.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov)){
5326 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
5329 Qx2 += TMath::Cos(2*aodTrack->Phi());
5330 Qy2 += TMath::Sin(2*aodTrack->Phi());
5331 Qx3 += TMath::Cos(3*aodTrack->Phi());
5332 Qy3 += TMath::Sin(3*aodTrack->Phi());
5336 Float_t evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
5337 Float_t evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
5339 fgPsi2tpc = evPlAng2;
5340 fgPsi3tpc = evPlAng3;
5342 fPhiRPTPC->Fill(iC,evPlAng2);
5343 fPhiRPTPCv3->Fill(iC,evPlAng3);
5348 AliAODVZERO* aodV0 = event->GetVZEROData();
5350 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
5351 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
5352 Float_t multv0 = aodV0->GetMultiplicity(iv0);
5355 if(fAnalysisType == "AOD"){//not for reco MC tracks, only for real data
5356 if (iv0 < 32){ // V0C
5357 Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5358 Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5359 Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5360 Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5362 Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5363 Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5364 Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5365 Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
5369 if (iv0 < 32){ // V0C
5370 Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5371 Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5372 Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5373 Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
5375 Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5376 Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5377 Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5378 Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
5383 //grab for each centrality the proper histo with the Qx and Qy to do the recentering
5384 Double_t Qxamean2 = fMeanQ[iCcal][1][0];
5385 Double_t Qxarms2 = fWidthQ[iCcal][1][0];
5386 Double_t Qyamean2 = fMeanQ[iCcal][1][1];
5387 Double_t Qyarms2 = fWidthQ[iCcal][1][1];
5388 Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
5389 Double_t Qxarms3 = fWidthQv3[iCcal][1][0];
5390 Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
5391 Double_t Qyarms3 = fWidthQv3[iCcal][1][1];
5393 Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
5394 Double_t Qxcrms2 = fWidthQ[iCcal][0][0];
5395 Double_t Qycmean2 = fMeanQ[iCcal][0][1];
5396 Double_t Qycrms2 = fWidthQ[iCcal][0][1];
5397 Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
5398 Double_t Qxcrms3 = fWidthQv3[iCcal][0][0];
5399 Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
5400 Double_t Qycrms3 = fWidthQv3[iCcal][0][1];
5402 Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
5403 Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
5404 Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
5405 Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
5406 Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
5407 Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
5408 Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
5409 Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
5411 //to calculate 2nd order event plane with v0M
5412 Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
5413 /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
5414 Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
5415 /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
5417 //here the calculated event plane is within -Pi to +Pi(delete it , no use here , only for definition)
5418 Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
5419 Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
5420 Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
5424 Float_t evPlAngV0ACor2=999.;
5425 Float_t evPlAngV0CCor2=999.;
5426 Float_t evPlAngV0ACor3=999.;
5427 Float_t evPlAngV0CCor3=999.;
5430 if(fAnalysisType == "AOD"){
5431 evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
5432 evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
5433 evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
5434 evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
5437 evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
5438 evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
5439 evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
5440 evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
5444 AliEventplane *ep = event->GetEventplane();
5445 evPlAngV0ACor2 = ep->GetEventplane("V0A", event, 2);
5446 evPlAngV0CCor2 = ep->GetEventplane("V0C", event, 2);
5447 evPlAngV0ACor3 = ep->GetEventplane("V0A", event, 3);
5448 evPlAngV0CCor3 = ep->GetEventplane("V0C", event, 3);
5451 fgPsi2v0a = evPlAngV0ACor2;
5452 fgPsi2v0c = evPlAngV0CCor2;
5453 fgPsi3v0a = evPlAngV0ACor3;
5454 fgPsi3v0c = evPlAngV0CCor3;
5456 // Fill EP distribution histograms evPlAng
5458 fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
5459 fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
5461 fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
5462 fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
5464 // Fill histograms needed for resolution evaluation
5465 fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
5466 fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
5467 fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
5469 fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
5470 fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
5471 fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
5475 Float_t gVZEROEventPlane = -10.;
5476 Float_t gReactionPlane = -10.;
5477 Double_t qxTot = 0.0, qyTot = 0.0;
5479 AliEventplane *ep = event->GetEventplane();
5481 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
5482 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
5483 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
5484 gReactionPlane = gVZEROEventPlane;
5488 //return gReactionPlane;
5490 //make the final 2nd order event plane within 0 to Pi
5491 //using data and reco tracks only
5492 if(fgPsi2v0a!=999. && fgPsi2v0a < 0.) fgPsi2v0a += TMath::Pi();
5493 if(fgPsi2v0c!=999. && fgPsi2v0c < 0.) fgPsi2v0c += TMath::Pi();
5494 if(fgPsi2tpc!=999. && fgPsi2tpc < 0.) fgPsi2tpc += TMath::Pi();
5495 //using truth tracks only
5496 if(evplaneMC!=999. && evplaneMC < 0.) evplaneMC += TMath::Pi();
5497 if(fgPsi2v0aMC!=999. && fgPsi2v0aMC < 0.) fgPsi2v0aMC += TMath::Pi();
5498 if(fgPsi2v0cMC!=999. && fgPsi2v0cMC < 0.) fgPsi2v0cMC += TMath::Pi();
5499 if(fgPsi2tpcMC!=999. && fgPsi2tpcMC < 0.) fgPsi2tpcMC += TMath::Pi();
5500 //for the time being leave the 3rd order event planes within -pi/3 t0 +pi/3
5502 if(truth){//for truth MC
5503 if(fV2 && fEPdet=="header")eventplane=evplaneMC;
5504 if(fV2 && fEPdet=="V0A")eventplane=fgPsi2v0aMC;
5505 if(fV2 && fEPdet=="V0C")eventplane=fgPsi2v0cMC;
5506 if(fV2 && fEPdet=="TPC")eventplane=fgPsi2tpcMC;
5508 if(fV3 && fEPdet=="V0A")eventplane=fgPsi3v0aMC;
5509 if(fV3 && fEPdet=="V0C")eventplane=fgPsi3v0cMC;
5510 if(fV3 && fEPdet=="TPC")eventplane=fgPsi3tpcMC;
5512 else{//for data and recoMC
5513 if(fV2 && fEPdet=="V0A")eventplane=fgPsi2v0a;
5514 if(fV2 && fEPdet=="V0C")eventplane=fgPsi2v0c;
5515 if(fV2 && fEPdet=="TPC")eventplane=fgPsi2tpc;
5517 if(fV3 && fEPdet=="V0A")eventplane=fgPsi3v0a;
5518 if(fV3 && fEPdet=="V0C")eventplane=fgPsi3v0c;
5519 if(fV3 && fEPdet=="TPC")eventplane=fgPsi3tpc;
5524 }//AOD/MCAOD condition
5526 else eventplane=999.;
5531 //------------------------------------------------------------------------------------------------------------------
5532 void AliTwoParticlePIDCorr::OpenInfoCalbration(Int_t run){
5533 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
5534 TFile *foadb = TFile::Open(oadbfilename.Data());
5537 printf("OADB file %s cannot be opened\n",oadbfilename.Data());
5541 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
5543 printf("OADB object hMultV0BefCorr is not available in the file\n");
5547 if(!(cont->GetObject(run))){
5548 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
5551 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
5553 TF1 *fpol0 = new TF1("fpol0","pol0");
5554 fMultV0->Fit(fpol0,"","",0,31);
5555 fV0Cpol = fpol0->GetParameter(0);
5556 fMultV0->Fit(fpol0,"","",32,64);
5557 fV0Apol = fpol0->GetParameter(0);
5559 for(Int_t iside=0;iside<2;iside++){
5560 for(Int_t icoord=0;icoord<2;icoord++){
5561 for(Int_t i=0;i < 9;i++){
5563 if(iside==0 && icoord==0)
5564 snprintf(namecont,100,"hQxc2_%i",i);
5565 else if(iside==1 && icoord==0)
5566 snprintf(namecont,100,"hQxa2_%i",i);
5567 else if(iside==0 && icoord==1)
5568 snprintf(namecont,100,"hQyc2_%i",i);
5569 else if(iside==1 && icoord==1)
5570 snprintf(namecont,100,"hQya2_%i",i);
5572 cont = (AliOADBContainer*) foadb->Get(namecont);
5574 printf("OADB object %s is not available in the file\n",namecont);
5578 if(!(cont->GetObject(run))){
5579 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5582 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5583 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5586 if(iside==0 && icoord==0)
5587 snprintf(namecont,100,"hQxc3_%i",i);
5588 else if(iside==1 && icoord==0)
5589 snprintf(namecont,100,"hQxa3_%i",i);
5590 else if(iside==0 && icoord==1)
5591 snprintf(namecont,100,"hQyc3_%i",i);
5592 else if(iside==1 && icoord==1)
5593 snprintf(namecont,100,"hQya3_%i",i);
5595 cont = (AliOADBContainer*) foadb->Get(namecont);
5597 printf("OADB object %s is not available in the file\n",namecont);
5601 if(!(cont->GetObject(run))){
5602 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5605 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5606 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5612 //____________________________________________________________________
5613 void AliTwoParticlePIDCorr::FillPIDEventPlane(Double_t centrality,Int_t par,Float_t trigphi,Float_t fReactionPlane)
5616 // Event plane (determine psi bin)
5617 Double_t gPsiMinusPhi = 0.;
5618 Double_t gPsiMinusPhiBin = -10.;
5619 if(fRequestEventPlane){
5620 gPsiMinusPhi = TMath::Abs(trigphi - fReactionPlane);
5622 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
5623 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
5624 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
5625 gPsiMinusPhiBin = 0.0;
5627 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
5628 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
5629 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
5630 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
5631 gPsiMinusPhiBin = 1.0;
5633 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
5634 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
5635 gPsiMinusPhiBin = 2.0;
5638 gPsiMinusPhiBin = 3.0;
5640 fEventPlanePID->Fill(centrality,gPsiMinusPhiBin,(Float_t)par);
5644 //____________________________________________________________________________________________________
5646 TObjArray* AliTwoParticlePIDCorr::GetV0Particles(AliVEvent* event,Double_t Centrality)
5649 AliAODEvent* fAOD = dynamic_cast<AliAODEvent*>(event);
5651 //function to select v0's from AODs
5652 trkVtx=fAOD->GetPrimaryVertex();
5653 Float_t xv=trkVtx->GetX(), yv=trkVtx->GetY(), zv=trkVtx->GetZ();
5654 Int_t nV0sTot = fAOD->GetNumberOfV0s();
5656 TObjArray * selectedV0s = new TObjArray;
5657 selectedV0s->SetOwner(kTRUE);
5659 for (Int_t iV0 = 0; iV0 < nV0sTot; iV0++)
5662 AliAODv0 *v0=fAOD->GetV0(iV0);
5664 if(!CheckStatusv0(v0)) continue;
5666 AliAODTrack *ptrack=(AliAODTrack*)v0->GetDaughter(0);
5667 AliAODTrack *ntrack=(AliAODTrack*)v0->GetDaughter(1);
5669 Bool_t cutK0sPID=kFALSE;
5670 Bool_t cutLambdaPID=kFALSE;
5671 Bool_t cutAntiLambdaPID=kFALSE;
5673 if(fUsev0DaughterPID)
5675 //use fHelperPID check PID of the daughter tracks
5676 //v0 daughter PID may be helpful in distangling k0S and (Anti)Lamda
5678 Int_t PIDptrack = GetParticle(ptrack,kFALSE);
5679 Int_t PIDntrack = GetParticle(ntrack ,kFALSE);
5681 if(PIDptrack ==0 && PIDntrack == 0) cutK0sPID=kTRUE;
5683 if(PIDptrack==2 && PIDntrack ==0) cutLambdaPID=kTRUE;
5685 if (PIDptrack ==0 && PIDntrack == 2) cutAntiLambdaPID=kTRUE;
5689 // effective mass calculations for each hypothesis(without daughter PID)
5690 Double_t InvMassK0s = v0->MassK0Short();
5691 Double_t InvMassAntiLambda = v0->MassAntiLambda();
5692 Double_t InvMassLambda = v0->MassLambda();
5694 Float_t v0Pt=TMath::Sqrt(v0->Pt2V0());
5695 Float_t v0Eta=v0->Eta();
5696 Float_t v0Phi=v0->Phi();
5698 //This is simply raw v0 without any specialised cut
5699 fHistRawPtCentInvK0s->Fill(InvMassK0s,v0Pt,Centrality);
5700 fHistRawPtCentInvLambda->Fill(InvMassLambda,v0Pt,Centrality);
5701 fHistRawPtCentInvAntiLambda->Fill(InvMassAntiLambda,v0Pt,Centrality);
5705 v0->GetSecondaryVtx(xyz);
5707 dx=xyz[0]-xv, dy=xyz[1]-yv, dz=xyz[2]-zv;
5709 // Float_t v0DecayRadius=TMath::Sqrt(dx*dx + dy*dy);
5710 Float_t v0DecayLength=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
5711 // VO's main characteristics to check the reconstruction cuts
5712 // Float_t DcaV0Daughters = v0->DcaV0Daughters();
5713 Float_t V0cosPointAngle = v0->CosPointingAngle(trkVtx);
5714 // Float_t DcaPosToPrimVertex = v0->DcaPosToPrimVertex();
5715 //Float_t DcaNegToPrimVertex = v0->DcaNegToPrimVertex();
5716 //Float_t Dcav0PVz = v0->DcaV0ToPrimVertex();
5717 Float_t v0Pz=v0->Pz();
5718 Float_t v0P= TMath::Sqrt(v0Pt*v0Pt + v0Pz*v0Pz);
5720 Float_t ctauLambda =999.;
5721 Float_t ctauAntiLambda = 999.;
5722 Float_t ctauK0s = 999.;
5725 ctauLambda = (v0DecayLength*InvMassLambda)/v0P;
5726 ctauAntiLambda = (v0DecayLength*InvMassAntiLambda)/v0P;
5727 ctauK0s = (v0DecayLength*InvMassK0s)/v0P;
5730 Bool_t ctauCutK0s= ctauK0s < NCtau*fCutctauK0s ; //ctauK0s 2.68 cm, mean life time of K0s is 8.95 x10^(-11)
5731 Bool_t ctauCutLambda = ctauLambda < NCtau*fCutctauLambda; //ctauLambda 7.8 cm ,mean life is 2.6 x10 ^(-10) ***** 3xctau is the accepted limit
5732 Bool_t ctauCutAntiLambda= ctauAntiLambda < NCtau*fCutctauAntiLambda;
5734 Bool_t RapCutK0s = v0->RapK0Short() < fRapCutK0s;
5735 Bool_t RapCutLambda = v0->RapLambda() < fRapCutLambda;
5736 Bool_t RapCutAntiLambda = v0->Y(-3122) < fRapCutLambda;
5738 Bool_t CPACut= V0cosPointAngle > fMinCPA; //cosine of pointing angle with v0 should be greater than 0.998
5740 //Now we put a loose mass cut which will be tightened later
5741 Bool_t MassCutLooseK0s=(TMath::Abs(InvMassK0s - 0.497614) < 0.1);
5742 Bool_t MassCutLooseLambda=(TMath::Abs(InvMassLambda - 1.115683) < 0.1); // cut is same for Anti-Lambda
5743 Bool_t MassCutLooseAntiLambda=(TMath::Abs(InvMassAntiLambda - 1.115683) < 0.1); // cut is same for Anti-Lambda
5745 //Special Cut for Kshort arementeros podalanski plot
5746 Bool_t ArmenterosCut =kFALSE;
5747 if(ctauCutK0s && RapCutK0s && CPACut && MassCutLooseK0s)
5750 Float_t lAlphaV0 = v0->AlphaV0();
5751 Float_t lPtArmV0 = v0->PtArmV0();
5753 ArmenterosCut = lPtArmV0 > TMath::Abs(0.2*lAlphaV0);
5757 Bool_t IskShortOk=(ctauCutK0s && RapCutK0s && CPACut && MassCutLooseK0s && ArmenterosCut);
5759 Bool_t IsLambdaOk=(ctauCutLambda && RapCutLambda && CPACut && MassCutLooseLambda);
5761 Bool_t IsAntiLambdaOk=(ctauCutAntiLambda && RapCutAntiLambda && CPACut && MassCutLooseAntiLambda);
5763 //Difference on Lambda and Anti-Lambda can be made through daughter PID
5766 Int_t particletype=999;
5768 if( IskShortOk || cutK0sPID )
5770 fHistFinalPtCentInvK0s->Fill(InvMassK0s,v0Pt,Centrality);
5773 Short_t chargeval=0;
5774 Float_t effmatrix=1.0;
5775 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassK0s,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5776 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5777 selectedV0s->Add(copy1);
5782 if(IsLambdaOk || cutLambdaPID)
5784 fHistFinalPtCentInvLambda->Fill(InvMassLambda,v0Pt,Centrality);
5785 //Add in the LRCParticle and give Lambda a tag 5
5788 Short_t chargeval=0;
5789 Float_t effmatrix=1.0;
5790 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassLambda,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5791 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5792 selectedV0s->Add(copy1);
5795 if(IsAntiLambdaOk || cutAntiLambdaPID)
5797 fHistFinalPtCentInvLambda->Fill(InvMassAntiLambda,v0Pt,Centrality);
5798 //Add in the LRCParticle and give Lambda a tag 6
5799 particletype=SpALam;
5800 Short_t chargeval=0;
5801 Float_t effmatrix=1.0;
5802 LRCParticlePID* copy1 = new LRCParticlePID(particletype,InvMassAntiLambda,chargeval,v0Pt,v0Eta, v0Phi,effmatrix,ptrack->GetTPCSharedMapPtr(),ntrack->GetTPCSharedMapPtr());
5803 copy1->SetUniqueID(eventno * 200000 + (Int_t)iV0);
5804 selectedV0s->Add(copy1);
5813 //___________________________________________________________________
5814 Bool_t AliTwoParticlePIDCorr :: CheckStatusv0Daughter(AliAODTrack *t1 ,AliAODTrack *t2)
5816 if (!t1->IsOn(AliAODTrack::kTPCrefit) || !t2->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
5817 // Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
5818 if(t1->GetTPCClusterInfo(2,1)<fDaugNClsTPC || t2->GetTPCClusterInfo(2,1)<fDaugNClsTPC) return kFALSE ;
5820 // ---------------- Fraction of TPC Shared Cluster
5821 Float_t fracPosDaugTPCSharedMap = GetFractionTPCSharedCls(t1);
5822 Float_t fracNegDaugTPCSharedMap = GetFractionTPCSharedCls(t2);
5824 if( (fracPosDaugTPCSharedMap > fFracTPCcls) || (fracNegDaugTPCSharedMap > fFracTPCcls) )
5830 //___________________________________________________________________________________________
5832 Float_t AliTwoParticlePIDCorr :: GetFractionTPCSharedCls( AliAODTrack *track)
5834 // Rejects tracks with shared clusters after filling a control histogram
5835 // This overload is used for primaries
5837 // Get the shared maps
5838 const TBits sharedMap = track->GetTPCSharedMap();
5840 return 1.*sharedMap.CountBits()/track->GetTPCNclsF();
5843 //______________________________________________________________________
5844 Bool_t AliTwoParticlePIDCorr :: CheckStatusv0(AliAODv0 *v1)
5847 // Offline reconstructed V0 only
5848 if (v1->GetOnFlyStatus()) return kFALSE;
5850 AliAODTrack *ptrack=(AliAODTrack *)v1->GetDaughter(0);
5851 AliAODTrack *ntrack=(AliAODTrack *)v1->GetDaughter(1);
5853 if(!ptrack || !ntrack) return kFALSE;
5855 if(ptrack->Charge()==-1 || ntrack->Charge()==1) return kFALSE; //remove wrongly identified charge pairs
5857 if(ptrack->Charge()==0 || ntrack->Charge()==0) return kFALSE; //remove uncharged pairs
5859 if(ptrack->Charge() == ntrack->Charge()) return kFALSE; //remove like sign pairs
5861 if(!CheckStatusv0Daughter(ptrack,ntrack)) return kFALSE;//daughters need to pass some basic cuts
5863 if(TMath::Abs(ptrack->Eta()) > fmaxeta || TMath::Abs(ntrack->Eta()) > fmaxeta) return kFALSE; // remove daughters beyond eta bound |0.8|
5865 if(ptrack->Pt() < fMinPtDaughter || ntrack->Pt() < fMinPtDaughter) return kFALSE; // remove daughter tracks below minmum p |1.0 GeV/c|
5867 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|
5869 // Daughters: Impact parameter of daughter to prim vtx
5870 Float_t xy = v1->DcaPosToPrimVertex();
5871 if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE; //0.1 cm
5872 xy = v1->DcaNegToPrimVertex();
5873 if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE; //0.1 cm
5876 Float_t dca = v1->DcaV0Daughters();
5877 if (dca>fMaxDCADaughter) return kFALSE; //1.0 cm
5879 // V0: Cosine of the pointing angle
5880 Float_t cpa=v1->CosPointingAngle(trkVtx); //0.997
5881 if (cpa<fMinCPA) return kFALSE;
5883 // V0: Fiducial volume
5884 Double_t xyz[3]; v1->GetSecondaryVtx(xyz);
5885 Float_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
5886 if (r2<5.*5.) return kFALSE;
5887 if (r2>lMax*lMax) return kFALSE; //lmax=100 cm
5893 //__________________________________________________________________________________________
5894 Bool_t AliTwoParticlePIDCorr::IsTrackFromV0(AliAODEvent* fAOD,AliAODTrack* track)
5896 //to check whether a daughter being taken as associated
5897 Int_t assoID = track->GetID();
5899 for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
5900 AliAODv0* aodV0 = fAOD->GetV0(i);
5902 AliAODTrack *trackPos=(AliAODTrack *)(aodV0->GetDaughter(0));
5903 AliAODTrack *trackNeg=(AliAODTrack *)(aodV0->GetDaughter(1));
5906 Int_t negID = trackNeg->GetID();
5907 Int_t posID = trackPos->GetID();
5909 if ((TMath::Abs(negID)+1)==(TMath::Abs(assoID))){ return kTRUE;}
5910 if ((TMath::Abs(posID)+1)==(TMath::Abs(assoID))){ return kTRUE;}
5911 //----------------------------------
5918 //----------------------------------------------------------
5919 void AliTwoParticlePIDCorr::Terminate(Option_t *)
5921 // Draw result to screen, or perform fitting, normalizations
5922 // Called once at the end of the query
5923 fOutput = dynamic_cast<TList*> (GetOutputData(1));
5924 if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
5928 //------------------------------------------------------------------