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"
34 #include "AliCFContainer.h"
36 #include "THnSparse.h"
40 #include "AliESDpid.h"
41 #include "AliAODpidUtil.h"
42 #include <AliPIDResponse.h>
43 #include "AliPIDCombined.h"
45 #include <AliAnalysisManager.h>
46 #include <AliInputEventHandler.h>
47 #include "AliAODInputHandler.h"
49 #include "AliAnalysisTaskSE.h"
50 #include "AliAnalysisManager.h"
51 #include "AliCentrality.h"
53 #include "AliVEvent.h"
54 #include "AliAODEvent.h"
55 #include "AliAODTrack.h"
56 #include "AliVTrack.h"
58 #include "THnSparse.h"
60 #include "AliAODMCHeader.h"
61 #include "AliAODMCParticle.h"
62 #include "AliMCEventHandler.h"
63 #include "AliMCEvent.h"
64 #include "AliMCParticle.h"
65 #include "TParticle.h"
66 #include <TDatabasePDG.h>
67 #include <TParticlePDG.h>
69 #include "AliGenCocktailEventHeader.h"
70 #include "AliGenEventHeader.h"
71 #include "AliCollisionGeometry.h"
72 #include "AliOADBContainer.h"
74 #include "AliEventPoolManager.h"
75 #include "AliAnalysisUtils.h"
76 using namespace AliPIDNameSpace;
79 ClassImp(AliTwoParticlePIDCorr)
80 ClassImp(LRCParticlePID)
82 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
83 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
84 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
85 //Source code::dphicorrelations,VnV0, TaskBFpsi, AliHelperPID,
87 //________________________________________________________________________
88 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
93 fCentralityMethod("V0A"),
94 fPPVsMultUtils(kFALSE),
96 fRequestEventPlane(kFALSE),
97 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
102 fSharedClusterCut(-1),
103 fSharedTPCmapCut(-1),
104 fSharedfraction_Pair_cut(-1),
106 skipParticlesAbove(0),
108 ffilltrigassoUNID(kFALSE),
109 ffilltrigUNIDassoID(kFALSE),
110 ffilltrigIDassoUNID(kTRUE),
111 ffilltrigIDassoID(kFALSE),
112 ffilltrigIDassoIDMCTRUTH(kFALSE),
113 fMaxNofMixingTracks(50000),
114 fPtOrderMCTruth(kTRUE),
115 fPtOrderDataReco(kTRUE),
116 fWeightPerEvent(kFALSE),
117 fTriggerSpeciesSelection(kFALSE),
118 fAssociatedSpeciesSelection(kFALSE),
119 fRandomizeReactionPlane(kFALSE),
120 fTriggerSpecies(SpPion),
121 fAssociatedSpecies(SpPion),
124 fSelectHighestPtTrig(kFALSE),
125 fcontainPIDtrig(kTRUE),
126 fcontainPIDasso(kFALSE),
128 frejectPileUp(kFALSE),
133 fselectprimaryTruth(kTRUE),
134 fonlyprimarydatareco(kFALSE),
137 ffillhistQAReco(kFALSE),
138 ffillhistQATruth(kFALSE),
139 kTrackVariablesPair(0),
172 fhistJetTrigestimate(0),
173 fTwoTrackDistancePtdip(0x0),
174 fTwoTrackDistancePtdipmix(0x0),
175 fCentralityCorrelation(0x0),
176 fHistVZEROAGainEqualizationMap(0),
177 fHistVZEROCGainEqualizationMap(0),
178 fHistVZEROChannelGainEqualizationMap(0),
179 fCentralityWeights(0),
182 fHistEQVZEROvsTPCmultiplicity(0x0),
183 fHistEQVZEROAvsTPCmultiplicity(0x0),
184 fHistEQVZEROCvsTPCmultiplicity(0x0),
185 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
186 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
187 fHistVZEROCvsVZEROAmultiplicity(0x0),
188 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
189 fHistVZEROSignal(0x0),
190 fHistEventPlaneTruth(0x0),
191 fHistPsiMinusPhi(0x0),
206 gReactionPlane(999.),
234 fControlConvResoncances(0),
249 fCorrelatonTruthPrimary(0),
250 fCorrelatonTruthPrimarymix(0),
256 fTHnCorrIDUNIDmix(0),
258 fTHnTrigcountMCTruthPrim(0),
261 fAnalysisType("AOD"),
263 ftwoTrackEfficiencyCutDataReco(kTRUE),
264 fTwoTrackCutMinRadius(0.8),
265 fTwoTrackCutMaxRadius(2.5),
266 twoTrackEfficiencyCutValue(0.02),
272 fRequestTOFPID(kTRUE),
273 fPIDType(NSigmaTPCTOF),
274 fFIllPIDQAHistos(kTRUE),
277 fdiffPIDcutvalues(kFALSE),
282 fHighPtKaonNSigmaPID(-1),
283 fHighPtKaonSigma(3.5),
284 fUseExclusiveNSigma(kFALSE),
285 fRemoveTracksT0Fill(kFALSE),
287 fTriggerSelectCharge(0),
288 fAssociatedSelectCharge(0),
289 fTriggerRestrictEta(-1),
290 fEtaOrdering(kFALSE),
291 fCutConversions(kFALSE),
292 fCutResonances(kFALSE),
293 fRejectResonanceDaughters(-1),
295 fInjectedSignals(kFALSE),
296 fRemoveWeakDecays(kFALSE),
297 fRemoveDuplicates(kFALSE),
298 fapplyTrigefficiency(kFALSE),
299 fapplyAssoefficiency(kFALSE),
300 ffillefficiency(kFALSE),
301 fmesoneffrequired(kFALSE),
302 fkaonprotoneffrequired(kFALSE),
307 for ( Int_t i = 0; i < 16; i++) {
311 for ( Int_t i = 0; i < 6; i++ ){
312 fTrackHistEfficiency[i] = NULL;
313 effcorection[i]=NULL;
316 for ( Int_t i = 0; i < 2; i++ ){
317 fTwoTrackDistancePt[i]=NULL;
318 fTwoTrackDistancePtmix[i]=NULL;
321 for(Int_t ipart=0;ipart<NSpecies;ipart++)
322 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
323 fnsigmas[ipart][ipid]=999.;
325 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
327 for(Int_t i = 0; i != 2; ++i)
328 for(Int_t j = 0; j != 2; ++j)
329 for(Int_t iC = 0; iC < 9; iC++){
330 fMeanQ[iC][i][j] = 0.;
331 fWidthQ[iC][i][j] = 1.;
332 fMeanQv3[iC][i][j] = 0.;
333 fWidthQv3[iC][i][j] = 1.;
337 //________________________________________________________________________
338 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
339 :AliAnalysisTaskSE(name),
343 fCentralityMethod("V0A"),
344 fPPVsMultUtils(kFALSE),
346 fRequestEventPlane(kFALSE),
347 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
352 fSharedClusterCut(-1),
353 fSharedTPCmapCut(-1),
354 fSharedfraction_Pair_cut(-1),
356 skipParticlesAbove(0),
358 ffilltrigassoUNID(kFALSE),
359 ffilltrigUNIDassoID(kFALSE),
360 ffilltrigIDassoUNID(kTRUE),
361 ffilltrigIDassoID(kFALSE),
362 ffilltrigIDassoIDMCTRUTH(kFALSE),
363 fMaxNofMixingTracks(50000),
364 fPtOrderMCTruth(kTRUE),
365 fPtOrderDataReco(kTRUE),
366 fWeightPerEvent(kFALSE),
367 fTriggerSpeciesSelection(kFALSE),
368 fAssociatedSpeciesSelection(kFALSE),
369 fRandomizeReactionPlane(kFALSE),
370 fTriggerSpecies(SpPion),
371 fAssociatedSpecies(SpPion),
374 fSelectHighestPtTrig(kFALSE),
375 fcontainPIDtrig(kTRUE),
376 fcontainPIDasso(kFALSE),
378 frejectPileUp(kFALSE),
383 fselectprimaryTruth(kTRUE),
384 fonlyprimarydatareco(kFALSE),
387 ffillhistQAReco(kFALSE),
388 ffillhistQATruth(kFALSE),
389 kTrackVariablesPair(0),
422 fhistJetTrigestimate(0),
423 fTwoTrackDistancePtdip(0x0),
424 fTwoTrackDistancePtdipmix(0x0),
425 fCentralityCorrelation(0x0),
426 fHistVZEROAGainEqualizationMap(0),
427 fHistVZEROCGainEqualizationMap(0),
428 fHistVZEROChannelGainEqualizationMap(0),
429 fCentralityWeights(0),
432 fHistEQVZEROvsTPCmultiplicity(0x0),
433 fHistEQVZEROAvsTPCmultiplicity(0x0),
434 fHistEQVZEROCvsTPCmultiplicity(0x0),
435 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
436 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
437 fHistVZEROCvsVZEROAmultiplicity(0x0),
438 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
439 fHistVZEROSignal(0x0),
440 fHistEventPlaneTruth(0x0),
441 fHistPsiMinusPhi(0x0),
456 gReactionPlane(999.),
484 fControlConvResoncances(0),
499 fCorrelatonTruthPrimary(0),
500 fCorrelatonTruthPrimarymix(0),
506 fTHnCorrIDUNIDmix(0),
508 fTHnTrigcountMCTruthPrim(0),
511 fAnalysisType("AOD"),
513 ftwoTrackEfficiencyCutDataReco(kTRUE),
514 fTwoTrackCutMinRadius(0.8),
515 fTwoTrackCutMaxRadius(2.5),
516 twoTrackEfficiencyCutValue(0.02),
522 fRequestTOFPID(kTRUE),
523 fPIDType(NSigmaTPCTOF),
524 fFIllPIDQAHistos(kTRUE),
527 fdiffPIDcutvalues(kFALSE),
532 fHighPtKaonNSigmaPID(-1),
533 fHighPtKaonSigma(3.5),
534 fUseExclusiveNSigma(kFALSE),
535 fRemoveTracksT0Fill(kFALSE),
537 fTriggerSelectCharge(0),
538 fAssociatedSelectCharge(0),
539 fTriggerRestrictEta(-1),
540 fEtaOrdering(kFALSE),
541 fCutConversions(kFALSE),
542 fCutResonances(kFALSE),
543 fRejectResonanceDaughters(-1),
545 fInjectedSignals(kFALSE),
546 fRemoveWeakDecays(kFALSE),
547 fRemoveDuplicates(kFALSE),
548 fapplyTrigefficiency(kFALSE),
549 fapplyAssoefficiency(kFALSE),
550 ffillefficiency(kFALSE),
551 fmesoneffrequired(kFALSE),
552 fkaonprotoneffrequired(kFALSE),
555 // The last in the above list should not have a comma after it
559 for ( Int_t i = 0; i < 16; i++) {
563 for ( Int_t i = 0; i < 6; i++ ){
564 fTrackHistEfficiency[i] = NULL;
565 effcorection[i]=NULL;
569 for ( Int_t i = 0; i < 2; i++ ){
570 fTwoTrackDistancePt[i]=NULL;
571 fTwoTrackDistancePtmix[i]=NULL;
574 for(Int_t ipart=0;ipart<NSpecies;ipart++)
575 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
576 fnsigmas[ipart][ipid]=999.;
578 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
580 for(Int_t i = 0; i != 2; ++i)
581 for(Int_t j = 0; j != 2; ++j)
582 for(Int_t iC = 0; iC < 9; iC++){
583 fMeanQ[iC][i][j] = 0.;
584 fWidthQ[iC][i][j] = 1.;
585 fMeanQv3[iC][i][j] = 0.;
586 fWidthQv3[iC][i][j] = 1.;
590 // Define input and output slots here (never in the dummy constructor)
591 // Input slot #0 works with a TChain - it is connected to the default input container
592 // Output slot #1 writes into a TH1 container
594 DefineOutput(1, TList::Class()); // for output list
595 DefineOutput(2, TList::Class());
599 //________________________________________________________________________
600 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
602 // Destructor. Clean-up the output list, but not the histograms that are put inside
603 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
604 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
609 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
614 if (fPID) delete fPID;
615 if (fPIDCombined) delete fPIDCombined;
618 //________________________________________________________________________
620 //////////////////////////////////////////////////////////////////////////////////////////////////
622 TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
623 // returns histo named name
624 return (TH2F*) fOutputList->FindObject(name);
627 //////////////////////////////////////////////////////////////////////////////////////////////////
629 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
633 // Puts the argument in the range [-pi/2,3 pi/2].
636 if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
637 if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
642 //________________________________________________________________________
643 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
646 // Called once (on the worker node)
647 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
648 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
649 fPID = inputHandler->GetPIDResponse();
651 //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
653 //get the efficiency correction map
655 // global switch disabling the reference
656 // (to avoid "Replacing existing TH1" if several wagons are created in train)
657 Bool_t oldStatus = TH1::AddDirectoryStatus();
658 TH1::AddDirectory(kFALSE);
660 const Int_t nPsiTOF = 10;
661 const Int_t nCentrBin = 9;
664 fOutput = new TList();
665 fOutput->SetOwner(); // IMPORTANT!
667 fOutputList = new TList;
668 fOutputList->SetOwner();
669 fOutputList->SetName("PIDQAList");
673 fList->SetName("EPQAList");
675 fEventCounter = new TH1F("fEventCounter","EventCounter", 19, 0.5,19.5);
676 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
677 fEventCounter->GetXaxis()->SetBinLabel(3,"After PileUP Cut");//only for Data
678 fEventCounter->GetXaxis()->SetBinLabel(5,"Have A Vertex");
679 fEventCounter->GetXaxis()->SetBinLabel(7,"After vertex Cut");
680 fEventCounter->GetXaxis()->SetBinLabel(9,"Getting centrality");
681 fEventCounter->GetXaxis()->SetBinLabel(11,"After centrality flattening");
682 fEventCounter->GetXaxis()->SetBinLabel(13,"Within 0-100% centrality");
683 fEventCounter->GetXaxis()->SetBinLabel(15,"Event Analyzed");
684 //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
685 fOutput->Add(fEventCounter);
687 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
688 fOutput->Add(fEtaSpectrasso);
690 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
691 fOutput->Add(fphiSpectraasso);
693 if(fSampleType=="pPb" || fSampleType=="PbPb" || fPPVsMultUtils==kTRUE){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
694 fOutput->Add(fCentralityCorrelation);
697 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")
699 TString gCentName[8] = {"V0A","V0C","V0M","V0AEq","V0CEq","V0MEq","CL1","ZNA"};
700 fHistCentStats = new TH2F("fHistCentStats",
701 "Centrality statistics;;Cent percentile",
702 8,-0.5,7.5,220,-5,105);
703 for(Int_t i = 1; i <= 8; i++){
704 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
705 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
707 fOutput->Add(fHistCentStats);
710 if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
712 fhistcentrality=new TH1F("fhistcentrality","referencemultiplicity",30001,-0.5,30000.5);
713 fOutput->Add(fhistcentrality);
716 fhistcentrality=new TH1F("fhistcentrality","centrality",220,-5,105);
717 fOutput->Add(fhistcentrality);
720 if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
722 TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
723 fHistRefmult = new TH2F("fHistRefmult",
724 "Reference multiplicity",
725 4,-0.5,3.5,10000,0,20000);
726 for(Int_t i = 1; i <= 4; i++){
727 fHistRefmult->GetXaxis()->SetBinLabel(i,gmultName[i-1].Data());
728 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
730 fOutput->Add(fHistRefmult);
732 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
733 //TPC vs EQVZERO multiplicity
734 fHistEQVZEROvsTPCmultiplicity = new TH2F("fHistEQVZEROvsTPCmultiplicity","EqVZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
735 fHistEQVZEROvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO multiplicity (a.u.)");
736 fHistEQVZEROvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
737 fOutput->Add(fHistEQVZEROvsTPCmultiplicity);
740 fHistEQVZEROAvsTPCmultiplicity = new TH2F("fHistEQVZEROAvsTPCmultiplicity","EqVZERO_A vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
741 fHistEQVZEROAvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
742 fHistEQVZEROAvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
743 fOutput->Add(fHistEQVZEROAvsTPCmultiplicity);
746 fHistEQVZEROCvsTPCmultiplicity = new TH2F("fHistEQVZEROCvsTPCmultiplicity","EqVZERO_C vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
747 fHistEQVZEROCvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
748 fHistEQVZEROCvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
749 fOutput->Add(fHistEQVZEROCvsTPCmultiplicity);
751 //EQVZERO vs VZERO multiplicity
752 fHistVZEROCvsEQVZEROCmultiplicity = new TH2F("fHistVZEROCvsEQVZEROCmultiplicity","EqVZERO_C vs VZERO_C multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
753 fHistVZEROCvsEQVZEROCmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
754 fHistVZEROCvsEQVZEROCmultiplicity->GetYaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
755 fOutput->Add(fHistVZEROCvsEQVZEROCmultiplicity);
758 fHistVZEROAvsEQVZEROAmultiplicity = new TH2F("fHistVZEROAvsEQVZEROAmultiplicity","EqVZERO_A vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
759 fHistVZEROAvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
760 fHistVZEROAvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
761 fOutput->Add(fHistVZEROAvsEQVZEROAmultiplicity);
764 //VZEROC vs VZEROA multiplicity
765 fHistVZEROCvsVZEROAmultiplicity = new TH2F("fHistVZEROCvsVZEROAmultiplicity","VZERO_C vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
766 fHistVZEROCvsVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
767 fHistVZEROCvsVZEROAmultiplicity->GetYaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
768 fOutput->Add(fHistVZEROCvsVZEROAmultiplicity);
772 //EQVZEROC vs EQVZEROA multiplicity
773 fHistEQVZEROCvsEQVZEROAmultiplicity = new TH2F("fHistEQVZEROCvsEQVZEROAmultiplicity","EqVZERO_C vs EqVZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
774 fHistEQVZEROCvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
775 fHistEQVZEROCvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
776 fOutput->Add(fHistEQVZEROCvsEQVZEROAmultiplicity);
778 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
779 fOutput->Add(fHistVZEROSignal);
783 if(fRequestEventPlane){
786 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
787 fList->Add(fHistPsiMinusPhi);
789 fEventPlanePID = new TH3F("fEventPlanePID",";centrality;eventplane;PID",20,0.0,100.0,4,-0.5,3.5,4,-0.5,3.5);
790 fList->Add(fEventPlanePID);
794 if(fCutConversions || fCutResonances)
796 fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
797 fOutput->Add(fControlConvResoncances);
800 fHistoTPCdEdx = new TH2F("fHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
801 fOutputList->Add(fHistoTPCdEdx);
802 fHistoTOFbeta = new TH2F(Form("fHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
803 fOutputList->Add(fHistoTOFbeta);
805 fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
806 fOutputList->Add(fTPCTOFPion3d);
808 fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
809 fOutputList->Add(fTPCTOFKaon3d);
811 fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
812 fOutputList->Add(fTPCTOFProton3d);
816 fPionPt = new TH1F("fPionPt","p_{T} distribution",200,0.,10.);
817 fOutputList->Add(fPionPt);
818 fPionEta= new TH1F("fPionEta","#eta distribution",360,-1.8,1.8);
819 fOutputList->Add(fPionEta);
820 fPionPhi = new TH1F("fPionPhi","#phi distribution",340,0,6.8);
821 fOutputList->Add(fPionPhi);
823 fKaonPt = new TH1F("fKaonPt","p_{T} distribution",200,0.,10.);
824 fOutputList->Add(fKaonPt);
825 fKaonEta= new TH1F("fKaonEta","#eta distribution",360,-1.8,1.8);
826 fOutputList->Add(fKaonEta);
827 fKaonPhi = new TH1F("fKaonPhi","#phi distribution",340,0,6.8);
828 fOutputList->Add(fKaonPhi);
830 fProtonPt = new TH1F("fProtonPt","p_{T} distribution",200,0.,10.);
831 fOutputList->Add(fProtonPt);
832 fProtonEta= new TH1F("fProtonEta","#eta distribution",360,-1.8,1.8);
833 fOutputList->Add(fProtonEta);
834 fProtonPhi= new TH1F("fProtonPhi","#phi distribution",340,0,6.8);
835 fOutputList->Add(fProtonPhi);
838 fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
839 fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
840 fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
841 fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);
842 fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
843 fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
844 fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
845 fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
846 fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
847 fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
848 fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
849 fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
850 fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
851 fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
852 fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
853 fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
855 for(Int_t i = 0; i < 16; i++)
857 fOutput->Add(fHistQA[i]);
860 fPriHistShare = new TH1F ("fPriHistShare","Shared clusters, primaries;#shared clusters;counts",160,0,160);
861 fOutput->Add(fPriHistShare);
863 Int_t eventplaneaxis=0;
865 if (fRequestEventPlane) eventplaneaxis=1;
867 kTrackVariablesPair=6+SetChargeAxis+eventplaneaxis;
869 if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
871 if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
873 if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8+SetChargeAxis+eventplaneaxis;
876 // two particle histograms
877 Int_t anaSteps = 1; // analysis steps
878 const char* title = "d^{2}N_{ch}/d#varphid#eta";
880 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
881 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
882 TString* axisTitlePair; // axis titles for track variables
883 axisTitlePair=new TString[kTrackVariablesPair];
885 TString defaultBinningStr;
886 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"
887 "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"
888 "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
889 "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"
890 "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
891 "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
892 "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"
893 "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
895 if(fRequestEventPlane){
896 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))
900 defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
903 defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
906 if(SetChargeAxis==2){
907 defaultBinningStr += "TrigCharge: -2.0,0.0,2.0\n"; // course
908 defaultBinningStr += "AssoCharge: -2.0,0.0,2.0\n"; // course
910 // =========================================================
911 // Customization (adopted from AliUEHistograms)
912 // =========================================================
914 TObjArray* lines = defaultBinningStr.Tokenize("\n");
915 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
917 TString line(lines->At(i)->GetName());
918 TString tag = line(0, line.Index(":")+1);
919 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
920 fBinningString += line + "\n";
922 AliInfo(Form("Using custom binning for %s", tag.Data()));
925 fBinningString += fCustomBinning;
927 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
929 // =========================================================
931 // =========================================================
933 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
934 axisTitlePair[0] = " multiplicity";
936 dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
937 axisTitlePair[1] = "v_{Z} (cm)";
939 dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
940 axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
942 dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
943 axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
945 dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
946 axisTitlePair[4] = "#Delta#eta";
948 dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
949 axisTitlePair[5] = "#Delta#varphi (rad)";
953 if(fRequestEventPlane){
954 dBinsPair[dim_val] = GetBinning(fBinningString, "eventPlane", iBinPair[dim_val]);
955 axisTitlePair[dim_val] = "#varphi - #Psi_{2} (a.u.)";
959 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
960 dBinsPair[dim_val] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val]);
961 axisTitlePair[dim_val] = "TrigCharge";
963 dBinsPair[dim_val+1] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+1]);
964 axisTitlePair[dim_val+1] = "AssoCharge";
967 if(fcontainPIDtrig && !fcontainPIDasso){
968 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
969 axisTitlePair[dim_val] = "PIDTrig";
970 if(SetChargeAxis==2){
971 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
972 axisTitlePair[dim_val+1] = "TrigCharge";
974 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
975 axisTitlePair[dim_val+2] = "AssoCharge";
979 if(!fcontainPIDtrig && fcontainPIDasso){
980 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val]);
981 axisTitlePair[dim_val] = "PIDAsso";
983 if(SetChargeAxis==2){
984 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
985 axisTitlePair[dim_val+1] = "TrigCharge";
987 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
988 axisTitlePair[dim_val+2] = "AssoCharge";
992 if(fcontainPIDtrig && fcontainPIDasso){
994 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
995 axisTitlePair[dim_val] = "PIDTrig";
997 dBinsPair[dim_val+1] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val+1]);
998 axisTitlePair[dim_val+1] = "PIDAsso";
1000 if(SetChargeAxis==2){
1001 dBinsPair[dim_val+2] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+2]);
1002 axisTitlePair[dim_val+2] = "TrigCharge";
1004 dBinsPair[dim_val+3] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+3]);
1005 axisTitlePair[dim_val+3] = "AssoCharge";
1010 Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
1012 Int_t nPteffbin = -1;
1013 Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
1016 fminPtTrig=dBinsPair[2][0];
1017 fmaxPtTrig=dBinsPair[2][iBinPair[2]];
1018 fminPtAsso=dBinsPair[3][0];
1019 fmaxPtAsso=dBinsPair[3][iBinPair[3]];
1020 fmincentmult=dBinsPair[0][0];
1021 fmaxcentmult=dBinsPair[0][iBinPair[0]];
1023 //event pool manager
1024 Int_t MaxNofEvents=1000;
1025 const Int_t NofVrtxBins=10+(1+10)*2;
1026 Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
1027 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
1028 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210};
1030 if(fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
1032 const Int_t NofCentBins=10;
1033 Double_t CentralityBins[NofCentBins+1]={0.,9.,14.,19.,26.,34.,44.,58.,80.,500.,1000.};//Is This binning is fine for pp, or we don't require them....
1034 if(fRequestEventPlane){
1035 // Event plane angle (Psi) bins
1037 Double_t* psibins = NULL;
1039 psibins = GetBinning(fBinningString, "eventPlane", nPsiBins);
1041 const Int_t nPsiBins=6;
1042 Double_t psibins[nPsiBins+1]={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()};
1043 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1044 // if(psibins) delete [] psibins;
1048 const Int_t nPsiBinsd=1;
1049 Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1050 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1052 // fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
1054 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1059 const Int_t NofCentBins=15;
1060 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
1061 if(fRequestEventPlane){
1062 // Event plane angle (Psi) bins
1064 Double_t* psibins = NULL;
1066 psibins = GetBinning(fBinningString, "eventPlane", nPsiBins);
1068 const Int_t nPsiBins=6;
1069 Double_t psibins[nPsiBins+1]={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()};
1070 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1071 // if(psibins) delete [] psibins;
1075 const Int_t nPsiBinsd=1;
1076 Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1077 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1079 //fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
1081 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1086 AliError("Event Mixing required, but Pool Manager not initialized...");
1090 //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
1091 //fmaxPtComboeff=fmaxPtTrig;
1092 //THnSparses for calculation of efficiency
1094 if((fAnalysisType =="MCAOD") && ffillefficiency) {
1097 effbin[0]=iBinPair[0];
1098 effbin[1]=iBinPair[1];
1099 effbin[2]=nPteffbin;
1101 Int_t effsteps=5;//for each species type::primMCParticles(0),primRecoTracksMatched(1),allRecoTracksMatched(2),primRecoTracksMatchedPID(3),allRecoTracksMatchedPID(4)
1102 for(Int_t jj=0;jj<6;jj++)//PID type binning
1104 if(jj==5) effsteps=3;//for unidentified particles
1105 Histrename="fTrackHistEfficiency";Histrename+=jj;
1106 fTrackHistEfficiency[jj] = new AliTHn(Histrename.Data(), "Tracking efficiency", effsteps, 4, effbin);
1107 fTrackHistEfficiency[jj]->SetBinLimits(0, dBinsPair[0]);
1108 fTrackHistEfficiency[jj]->SetVarTitle(0, "Centrality");
1109 fTrackHistEfficiency[jj]->SetBinLimits(1, dBinsPair[1]);
1110 fTrackHistEfficiency[jj]->SetVarTitle(1, "zvtx");
1111 fTrackHistEfficiency[jj]->SetBinLimits(2, Pteff);
1112 fTrackHistEfficiency[jj]->SetVarTitle(2, "p_{T} (GeV/c)");
1113 fTrackHistEfficiency[jj]->SetBinLimits(3, EtaBin);
1114 fTrackHistEfficiency[jj]->SetVarTitle(3, "#eta");
1115 fOutput->Add(fTrackHistEfficiency[jj]);
1119 //AliThns for Correlation plots(data & MC)
1121 if(ffilltrigassoUNID)
1123 fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1124 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1125 fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
1126 fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
1128 fOutput->Add(fTHnCorrUNID);
1130 fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1131 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1132 fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
1133 fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1135 fOutput->Add(fTHnCorrUNIDmix);
1138 if(ffilltrigIDassoID)
1140 fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
1141 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1142 fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
1143 fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
1145 fOutput->Add(fTHnCorrID);
1147 fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1148 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1149 fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
1150 fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
1152 fOutput->Add(fTHnCorrIDmix);
1155 if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
1157 fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1158 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1159 fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
1160 fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
1162 fOutput->Add(fTHnCorrIDUNID);
1165 fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1166 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1167 fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
1168 fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1170 fOutput->Add(fTHnCorrIDUNIDmix);
1175 //ThnSparse for Correlation plots(truth MC)
1176 if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
1178 fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
1179 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1180 fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
1181 fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
1183 fOutput->Add(fCorrelatonTruthPrimary);
1186 fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
1187 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1188 fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
1189 fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
1191 fOutput->Add(fCorrelatonTruthPrimarymix);
1194 //binning for trigger no. counting
1197 if(SetChargeAxis==2) ChargeAxis=1;
1200 Int_t dims=3+ChargeAxis+eventplaneaxis;
1201 if(fcontainPIDtrig) dims=4+ChargeAxis+eventplaneaxis;
1202 fBinst= new Int_t[dims];
1203 Double_t* dBinsTrig[dims]; // bins for track variables
1204 TString* axisTitleTrig; // axis titles for track variables
1205 axisTitleTrig=new TString[dims];
1207 for(Int_t i=0; i<3;i++)
1209 fBinst[i]=iBinPair[i];
1210 dBinsTrig[i]=dBinsPair[i];
1211 axisTitleTrig[i]=axisTitlePair[i];
1213 Int_t dim_val_trig=3;
1214 if(fRequestEventPlane){
1215 fBinst[dim_val_trig]=iBinPair[6];//if fRequestEventPlane=TRUE, dim_val already becomes 7.
1216 dBinsTrig[dim_val_trig]=dBinsPair[6];
1217 axisTitleTrig[dim_val_trig]=axisTitlePair[6];
1221 if(!fcontainPIDtrig && !fcontainPIDasso && ChargeAxis==1){
1222 fBinst[dim_val_trig]=iBinPair[dim_val];
1223 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1224 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1227 if(fcontainPIDtrig && !fcontainPIDasso){
1228 fBinst[dim_val_trig]=iBinPair[dim_val];
1229 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1230 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1232 fBinst[dim_val_trig+1]=iBinPair[dim_val+1];
1233 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+1];
1234 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+1];
1238 if(!fcontainPIDtrig && fcontainPIDasso){
1240 fBinst[dim_val_trig]=iBinPair[dim_val+1];
1241 dBinsTrig[dim_val_trig]=dBinsPair[dim_val+1];
1242 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val+1];
1246 if(fcontainPIDtrig && fcontainPIDasso){
1247 fBinst[dim_val_trig]=iBinPair[dim_val];
1248 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1249 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1251 fBinst[dim_val_trig+1]=iBinPair[dim_val+2];
1252 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+2];
1253 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+2];
1257 //ThSparse for trigger counting(data & reco MC)
1258 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
1260 fTHnTrigcount = new AliTHn("fTHnTrigcount", "fTHnTrigcount", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1261 for(Int_t i=0; i<dims;i++){
1262 fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1263 fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1265 fOutput->Add(fTHnTrigcount);
1268 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
1269 //AliTHns for trigger counting(truth MC)
1270 fTHnTrigcountMCTruthPrim = new AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1271 for(Int_t i=0; i<dims;i++){
1272 fTHnTrigcountMCTruthPrim->SetBinLimits(i, dBinsTrig[i]);
1273 fTHnTrigcountMCTruthPrim->SetVarTitle(i, axisTitleTrig[i]);
1275 fOutput->Add(fTHnTrigcountMCTruthPrim);
1278 if(fAnalysisType=="MCAOD"){
1279 if(ffillhistQATruth)
1281 MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
1282 fOutputList->Add(MCtruthpt);
1284 MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
1285 fOutputList->Add(MCtrutheta);
1287 MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
1288 fOutputList->Add(MCtruthphi);
1290 MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
1291 fOutputList->Add(MCtruthpionpt);
1293 MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
1294 fOutputList->Add(MCtruthpioneta);
1296 MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
1297 fOutputList->Add(MCtruthpionphi);
1299 MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
1300 fOutputList->Add(MCtruthkaonpt);
1302 MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
1303 fOutputList->Add(MCtruthkaoneta);
1305 MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
1306 fOutputList->Add(MCtruthkaonphi);
1308 MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
1309 fOutputList->Add(MCtruthprotonpt);
1311 MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
1312 fOutputList->Add(MCtruthprotoneta);
1314 MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
1315 fOutputList->Add(MCtruthprotonphi);
1317 fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
1318 fOutputList->Add(fPioncont);
1320 fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
1321 fOutputList->Add(fKaoncont);
1323 fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
1324 fOutputList->Add(fProtoncont);
1326 fUNIDcont=new TH2F("fUNIDcont","fUNIDcont",10,-0.5,9.5,100,0.,10.);
1327 fOutputList->Add(fUNIDcont);
1330 fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1331 fEventno->GetXaxis()->SetTitle("Centrality");
1332 fEventno->GetYaxis()->SetTitle("Z_Vtx");
1333 fOutput->Add(fEventno);
1334 fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1335 fEventnobaryon->GetXaxis()->SetTitle("Centrality");
1336 fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
1337 fOutput->Add(fEventnobaryon);
1338 fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1339 fEventnomeson->GetXaxis()->SetTitle("Centrality");
1340 fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
1341 fOutput->Add(fEventnomeson);
1343 fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
1344 fOutput->Add(fhistJetTrigestimate);
1346 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);
1347 fOutput->Add(fTwoTrackDistancePtdip);
1349 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);
1350 fOutput->Add(fTwoTrackDistancePtdipmix);
1352 TString Histttrname;
1353 for(Int_t jj=0;jj<2;jj++)// PID type binning
1355 Histttrname="fTwoTrackDistancePt";Histttrname+=jj;
1356 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);
1357 fOutput->Add(fTwoTrackDistancePt[jj]);
1359 Histttrname="fTwoTrackDistancePtmix";Histttrname+=jj;
1360 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);
1361 fOutput->Add(fTwoTrackDistancePtmix[jj]);
1364 //DefineEventPool();
1366 if(fapplyTrigefficiency || fapplyAssoefficiency)
1368 const Int_t nDimt = 4;// cent zvtx pt eta
1369 Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
1370 Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
1371 Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
1373 TString Histrexname;
1374 for(Int_t jj=0;jj<6;jj++)// PID type binning
1376 Histrexname="effcorection";Histrexname+=jj;
1377 effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
1378 effcorection[jj]->Sumw2();
1379 effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
1380 effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
1381 effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
1382 effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
1383 effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
1384 effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1385 effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
1386 effcorection[jj]->GetAxis(3)->SetTitle("#eta");
1387 fOutput->Add(effcorection[jj]);
1389 // TFile *fsifile = new TFile(fefffilename,"READ");
1391 if (TString(fefffilename).BeginsWith("alien:"))
1392 TGrid::Connect("alien:");
1393 TFile *fileT=TFile::Open(fefffilename);
1395 for(Int_t jj=0;jj<6;jj++)//type binning
1397 Nameg="effmap";Nameg+=jj;
1398 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
1399 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
1401 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
1408 //*************************************************************EP plots***********************************************//
1409 if(fRequestEventPlane){
1410 // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
1412 fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
1413 fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
1414 fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
1416 fList->Add(fHResTPCv0A2);
1417 fList->Add(fHResTPCv0C2);
1418 fList->Add(fHResv0Cv0A2);
1421 fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
1422 fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
1423 fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
1425 fList->Add(fHResTPCv0A3);
1426 fList->Add(fHResTPCv0C3);
1427 fList->Add(fHResv0Cv0A3);
1429 // MC as in the dataEP resolution (but using MC tracks)
1430 if(fAnalysisType == "MCAOD" && fV2){
1431 fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
1432 fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
1433 fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
1434 fList->Add(fHResMA2);
1435 fList->Add(fHResMC2);
1436 fList->Add(fHResAC2);
1438 if(fAnalysisType == "MCAOD" && fV3){
1439 fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
1440 fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
1441 fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
1442 fList->Add(fHResMA3);
1443 fList->Add(fHResMC3);
1444 fList->Add(fHResAC3);
1448 // V0A and V0C event plane distributions
1450 fPhiRPTPC = new TH2F("fPhiRPTPCv2","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1451 fPhiRPTPCv3 = new TH2F("fPhiRPTPCv3","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1452 fList->Add(fPhiRPTPC);
1453 fList->Add(fPhiRPTPCv3);
1455 fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1456 fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1457 fList->Add(fPhiRPv0A);
1458 fList->Add(fPhiRPv0C);
1461 fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1462 fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1463 fList->Add(fPhiRPv0Av3);
1464 fList->Add(fPhiRPv0Cv3);
1466 fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth","#phi distribution of EP MCTRUTHheader;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1467 fList->Add(fHistEventPlaneTruth);
1471 //*****************************************************PIDQA histos*****************************************************//
1475 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1476 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1479 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1480 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);
1481 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1482 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1483 fOutputList->Add(fHistoNSigma);
1488 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1489 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1492 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1493 TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
1494 Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1495 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1496 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1497 fOutputList->Add(fHistoNSigma);
1502 if(fPIDType==Bayes){//use bayesianPID
1503 fPIDCombined = new AliPIDCombined();
1504 fPIDCombined->SetDefaultTPCPriors();//****************************************Need to know about it
1506 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1509 TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
1510 Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1511 fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1512 fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
1513 fOutputList->Add(fHistoBayes);
1516 TH2F *fHistoBayesTPC=new TH2F(Form("probBayes_TPC_%d",ipart),
1517 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1518 fHistoBayesTPC->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1519 fHistoBayesTPC->GetYaxis()->SetTitle(Form("Bayes prob TPC %s",kParticleSpeciesName[ipart]));
1520 fOutputList->Add(fHistoBayesTPC);
1522 TH2F *fHistoBayesTOF=new TH2F(Form("probBayes_TOF_%d",ipart),
1523 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1524 fHistoBayesTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1525 fHistoBayesTOF->GetYaxis()->SetTitle(Form("Bayes prob TOF %s",kParticleSpeciesName[ipart]));
1526 fOutputList->Add(fHistoBayesTOF);
1528 TH2F *fHistoBayesTPCTOF=new TH2F(Form("probBayes_TPCTOF_%d",ipart),
1529 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1530 fHistoBayesTPCTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1531 fHistoBayesTPCTOF->GetYaxis()->SetTitle(Form("Bayes prob TPCTOF %s",kParticleSpeciesName[ipart]));
1532 fOutputList->Add(fHistoBayesTPCTOF);
1536 //nsigma separation power plot
1537 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1540 TH2F *Pi_Ka_sep=new TH2F(Form("Pi_Ka_sep_%d",ipid),
1541 Form("Pi_Ka separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1542 Pi_Ka_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1543 Pi_Ka_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1544 fOutputList->Add(Pi_Ka_sep);
1546 TH2F *Pi_Pr_sep=new TH2F(Form("Pi_Pr_sep_%d",ipid),
1547 Form("Pi_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1548 Pi_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1549 Pi_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1550 fOutputList->Add(Pi_Pr_sep);
1552 TH2F *Ka_Pr_sep=new TH2F(Form("Ka_Pr_sep_%d",ipid),
1553 Form("Ka_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1554 Ka_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1555 Ka_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1556 fOutputList->Add(Ka_Pr_sep);
1560 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1561 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1564 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1565 TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
1566 Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1567 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1568 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1569 fOutputList->Add(fHistoNSigma);
1574 if (fAnalysisType == "MCAOD"){
1575 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1576 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1579 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1580 TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
1581 Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1582 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1583 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1584 fOutputList->Add(fHistoNSigma);
1589 for(Int_t idet=0;idet<fNDetectors;idet++){
1590 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1592 if(idet==fTOF)maxy=1.1;
1593 TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
1594 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1595 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1596 fOutputList->Add(fHistoPID);
1599 //PID signal plot, before PID cut
1600 for(Int_t idet=0;idet<fNDetectors;idet++){
1602 if(idet==fTOF)maxy=1.1;
1603 TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
1604 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1605 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1606 fOutputList->Add(fHistoPID);
1609 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
1610 PostData(2, fOutputList);
1611 if(fRequestEventPlane) PostData(3, fList);
1612 AliInfo("Finished setting up the Output");
1614 TH1::AddDirectory(oldStatus);
1619 //-------------------------------------------------------------------------------
1620 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1623 if(fAnalysisType == "AOD") {
1627 }//AOD--analysis-----
1629 else if(fAnalysisType == "MCAOD") {
1638 //-------------------------------------------------------------------------
1639 void AliTwoParticlePIDCorr::doMCAODevent()
1641 AliVEvent *event = InputEvent();
1642 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1643 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1645 AliError("Cannot get the AOD event");
1649 // count all events(physics triggered)
1650 fEventCounter->Fill(1);
1659 gReactionPlane = 999.;
1661 // get centrality object and check quality(valid for p-Pb and Pb-Pb; coming soon for pp 7 TeV)
1662 Double_t cent_v0=-1.0;
1663 Double_t effcent=1.0;
1664 Double_t refmultReco =0.0;
1666 //check the PIDResponse handler
1669 // get mag. field required for twotrack efficiency cut
1671 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1673 //check for TClonesArray(truth track MC information)
1674 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1676 AliFatal("Error: MC particles branch not found!\n");
1680 //check for AliAODMCHeader(truth event MC information)
1681 AliAODMCHeader *header=NULL;
1682 header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1684 printf("MC header branch not found!\n");
1688 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
1689 Float_t zVtxmc =header->GetVtxZ();
1690 if(TMath::Abs(zVtxmc)>fzvtxcut) return;
1692 // For productions with injected signals, figure out above which label to skip particles/tracks
1694 if (fInjectedSignals)
1696 AliGenEventHeader* eventHeader = 0;
1701 AliFatal("fInjectedSignals set but no MC header found");
1703 headers = header->GetNCocktailHeaders();
1704 eventHeader = header->GetCocktailHeader(0);
1708 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
1709 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
1710 AliError("First event header not found. Skipping this event.");
1711 //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1714 skipParticlesAbove = eventHeader->NProduced();
1715 AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
1718 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE))
1720 //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
1721 Double_t refmultTruth = GetAcceptedEventMultiplicity(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
1722 refmultReco = GetAcceptedEventMultiplicity(aod,kFALSE);
1723 if(refmultTruth<=0 || refmultReco<=0) return;
1724 cent_v0=refmultTruth;
1727 cent_v0=GetAcceptedEventMultiplicity(aod,kFALSE); //centrality value; 2nd argument has no meaning
1730 if(cent_v0<0.) return;
1731 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1733 //get the event plane in case of PbPb
1734 if(fRequestEventPlane){
1735 gReactionPlane=GetEventPlane(aod,kTRUE,cent_v0);//get the truth event plane
1736 if(gReactionPlane==999.) return;
1740 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)
1742 //TObjArray* tracksMCtruth=0;
1743 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
1744 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
1748 //There is a small difference between zvtx and zVtxmc??????
1749 //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1750 //cout<<"***********************************************zvtx="<<zvtx<<endl;
1752 //now process the truth particles(for both efficiency & correlation function)
1753 Int_t nMCTrack = fArrayMC->GetEntriesFast();
1755 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
1756 { //MC truth track loop starts
1758 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1761 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1765 //consider only charged particles
1766 if(partMC->Charge() == 0) continue;
1768 //consider only primary particles; neglect all secondary particles including from weak decays
1769 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1772 //remove injected signals(primaries above <maxLabel>)
1773 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1776 Bool_t isduplicate=kFALSE;
1777 if (fRemoveDuplicates)
1779 for (Int_t j=iMC+1; j<nMCTrack; ++j)
1780 {//2nd trutuh loop starts
1781 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1783 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1786 if (partMC->GetLabel() == partMC2->GetLabel())
1791 }//2nd truth loop ends
1793 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1795 //give only kinematic cuts at the truth level
1796 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1797 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
1799 if(!partMC) continue;//for safety
1801 //To determine multiplicity in case of PP
1803 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1804 //only physical primary(all/unidentified)
1805 if(ffillhistQATruth)
1807 MCtruthpt->Fill(partMC->Pt());
1808 MCtrutheta->Fill(partMC->Eta());
1809 MCtruthphi->Fill(partMC->Phi());
1812 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1813 Int_t particletypeTruth=-999;
1814 if (TMath::Abs(pdgtruth)==211)
1816 particletypeTruth=SpPion;
1817 if(ffillhistQATruth)
1819 MCtruthpionpt->Fill(partMC->Pt());
1820 MCtruthpioneta->Fill(partMC->Eta());
1821 MCtruthpionphi->Fill(partMC->Phi());
1824 if (TMath::Abs(pdgtruth)==321)
1826 particletypeTruth=SpKaon;
1827 if(ffillhistQATruth)
1829 MCtruthkaonpt->Fill(partMC->Pt());
1830 MCtruthkaoneta->Fill(partMC->Eta());
1831 MCtruthkaonphi->Fill(partMC->Phi());
1834 if(TMath::Abs(pdgtruth)==2212)
1836 particletypeTruth=SpProton;
1837 if(ffillhistQATruth)
1839 MCtruthprotonpt->Fill(partMC->Pt());
1840 MCtruthprotoneta->Fill(partMC->Eta());
1841 MCtruthprotonphi->Fill(partMC->Phi());
1844 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)
1846 if(fRequestEventPlane){
1847 FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
1850 // -- Fill THnSparse for efficiency and contamination calculation
1851 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
1853 Double_t primmctruth[4] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
1856 fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
1857 if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for primary truth mesons(3)
1858 if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for primary truth kaons+protons(4)
1859 if (TMath::Abs(pdgtruth)==211) fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
1860 if (TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
1861 if(TMath::Abs(pdgtruth)==2212) fTrackHistEfficiency[2]->Fill(primmctruth,0);//for protons
1864 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1865 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1867 Short_t chargeval=0;
1868 if(partMC->Charge()>0) chargeval=1;
1869 if(partMC->Charge()<0) chargeval=-1;
1870 if(chargeval==0) continue;
1871 const TBits *clustermap=0;
1872 const TBits *sharemap=0;
1873 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth,clustermap,sharemap);
1874 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1875 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1876 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1878 }//MC truth track loop ends
1880 //*********************still in event loop
1882 if (fRandomizeReactionPlane)//only for TRuth MC??
1884 Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
1885 Double_t angle = TMath::TwoPi() * centralityDigits;
1886 AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
1887 ShiftTracks(tracksMCtruth, angle);
1891 Float_t weghtval=1.0;
1893 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1895 //Fill Correlations for MC truth particles(same event)
1896 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1897 Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1900 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
1901 if (pool2 && pool2->IsReady())
1902 {//start mixing only when pool->IsReady
1903 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1904 {//proceed only when no. of trigger particles >0 in current event
1905 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
1906 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1907 { //pool event loop start
1908 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1909 if(!bgTracks6) continue;
1910 Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1912 }// pool event loop ends mixing case
1913 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1914 } //if pool->IsReady() condition ends mixing case
1916 //still in main event loop
1919 if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1923 //still in main event loop
1925 if(tracksMCtruth) delete tracksMCtruth;
1927 //now deal with reco tracks
1930 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1932 //detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
1933 if (fSampleType=="pp_2_76" || fCentralityMethod.EndsWith("_MANUAL") || (fSampleType=="pp_7" && fPPVsMultUtils==kFALSE)) cent_v0=refmultReco;
1934 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1936 if(fRequestEventPlane){
1937 gReactionPlane = GetEventPlane(aod,kFALSE,cent_v0);//get the reconstructed event plane
1938 if(gReactionPlane==999.) return;
1942 TExMap *trackMap = new TExMap();
1944 // --- track loop for mapping matrix
1947 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1948 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
1949 AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
1950 if (!track) continue;
1951 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
1952 if(tracktype!=1) continue;
1954 if(!track) continue;//for safety
1956 Int_t gid = track->GetID();
1957 trackMap->Add(gid,itrk);
1963 //TObjArray* tracksUNID=0;
1964 TObjArray* tracksUNID = new TObjArray;
1965 tracksUNID->SetOwner(kTRUE);
1967 //TObjArray* tracksID=0;
1968 TObjArray* tracksID = new TObjArray;
1969 tracksID->SetOwner(kTRUE);
1974 Double_t trackscount=0.0;
1975 // loop over reconstructed tracks
1976 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1977 { // reconstructed track loop starts
1978 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1979 if (!track) continue;
1980 //get the corresponding MC track at the truth level (doing reco matching)
1981 AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
1982 if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
1984 //remove injected signals
1985 if(fInjectedSignals)
1987 AliAODMCParticle* mother = recomatched;
1989 while (!mother->IsPhysicalPrimary())
1990 {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
1991 if (mother->GetMother() < 0)
1997 mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
2003 Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
2006 if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
2007 }//remove injected signals
2009 if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
2011 Bool_t isduplicate2=kFALSE;
2012 if (fRemoveDuplicates)
2014 for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
2016 AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
2017 if (!track2) continue;
2018 AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
2019 if(!recomatched2) continue;
2021 if (track->GetLabel() == track2->GetLabel())
2028 if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
2030 fHistQA[11]->Fill(track->GetTPCNcls());
2031 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2033 if(tracktype==0) continue;
2034 if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
2036 if(!track) continue;//for safety
2037 //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
2039 AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2041 if(fFilterBit==128){
2042 Int_t gid1 = track->GetID();
2043 //if(gid1>=0) PIDtrack = track;
2044 PIDtrack = aod->GetTrack(trackMap->GetValue(-1-gid1));
2045 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2050 //check for eta , phi holes
2051 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2052 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2054 Int_t particletypeMC=-9999;
2056 //tag all particles as unidentified
2057 particletypeMC=unidentified;
2059 Float_t effmatrix=1.;
2061 // -- Fill THnSparse for efficiency calculation
2062 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
2063 //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)
2065 //Clone & Reduce track list(TObjArray) for unidentified particles
2066 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2068 Short_t chargeval=0;
2069 if(track->Charge()>0) chargeval=1;
2070 if(track->Charge()<0) chargeval=-1;
2071 if(chargeval==0) continue;
2072 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2073 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
2074 LRCParticlePID* copy = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2075 copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
2076 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2079 //get the pdg code of the corresponding truth particle
2080 Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
2082 Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
2083 if(ffillefficiency) {
2084 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
2085 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
2086 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
2087 if(TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions
2088 if(TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
2089 if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
2091 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
2092 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
2093 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
2094 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
2095 if( TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions
2096 if( TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
2097 if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
2101 //now start the particle identification process:)
2103 Float_t dEdx = PIDtrack->GetTPCsignal();
2104 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2106 if(HasTOFPID(PIDtrack))
2108 Double_t beta = GetBeta(PIDtrack);
2109 fHistoTOFbeta->Fill(track->Pt(), beta);
2112 //do track identification(nsigma method)
2113 particletypeMC=GetParticle(PIDtrack,fFIllPIDQAHistos);//******************************problem is here
2114 switch(TMath::Abs(pdgCode)){
2116 if(fFIllPIDQAHistos){
2117 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2118 if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2119 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
2120 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
2125 if(fFIllPIDQAHistos){
2126 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2127 if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2128 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpKaon,ipid));
2129 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
2134 if(fFIllPIDQAHistos){
2135 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2136 if((ipid!=NSigmaTPC) && (!HasTOFPID(PIDtrack)))continue;//not filling TOF and combined if no TOF PID
2137 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpPion,ipid));
2138 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
2145 //2-d TPCTOF map(for each Pt interval)
2146 if(HasTOFPID(PIDtrack)){
2147 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2148 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2149 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2152 //Pt, Eta , Phi distribution of the reconstructed identified particles
2155 if (particletypeMC==SpPion)
2157 fPionPt->Fill(track->Pt());
2158 fPionEta->Fill(track->Eta());
2159 fPionPhi->Fill(track->Phi());
2161 if (particletypeMC==SpKaon)
2163 fKaonPt->Fill(track->Pt());
2164 fKaonEta->Fill(track->Eta());
2165 fKaonPhi->Fill(track->Phi());
2167 if (particletypeMC==SpProton)
2169 fProtonPt->Fill(track->Pt());
2170 fProtonEta->Fill(track->Eta());
2171 fProtonPhi->Fill(track->Phi());
2175 //for misidentification fraction calculation(do it with tuneonPID)
2176 if(particletypeMC==SpPion )
2178 if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
2179 if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
2180 if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
2181 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
2183 if(particletypeMC==SpKaon )
2185 if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
2186 if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
2187 if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
2188 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
2190 if(particletypeMC==SpProton )
2192 if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
2193 if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
2194 if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
2195 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
2197 if(particletypeMC==SpUndefined )
2199 if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
2200 if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
2201 if(TMath::Abs(pdgCode)==2212) fUNIDcont->Fill(5.,track->Pt());
2202 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fUNIDcont->Fill(7.,track->Pt());
2205 if(particletypeMC==SpUndefined) continue;
2207 if(fRequestEventPlane){
2208 FillPIDEventPlane(cent_v0,particletypeMC,track->Phi(),gReactionPlane);
2211 //fill tracking efficiency
2214 if(particletypeMC==SpPion || particletypeMC==SpKaon)
2216 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) {
2217 fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
2218 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
2221 if(particletypeMC==SpKaon || particletypeMC==SpProton)
2223 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) {
2224 fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
2225 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
2228 if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) {
2229 fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
2230 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
2232 if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
2233 fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
2234 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
2236 if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
2237 fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
2238 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
2242 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2244 Short_t chargeval=0;
2245 if(track->Charge()>0) chargeval=1;
2246 if(track->Charge()<0) chargeval=-1;
2247 if(chargeval==0) continue;
2248 if (fapplyTrigefficiency || fapplyAssoefficiency)
2249 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
2250 LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2251 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2252 tracksID->Add(copy1);
2254 }// if(tracktype==1) condition structure ands
2256 }//reco track loop ends
2258 //*************************************************************still in event loop
2263 //fill the centrality/multiplicity distribution of the selected events
2264 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2266 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??????
2268 //count selected events having centrality betn 0-100%
2269 fEventCounter->Fill(13);
2271 //***************************************event no. counting
2272 Bool_t isbaryontrig=kFALSE;
2273 Bool_t ismesontrig=kFALSE;
2274 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2276 if(tracksID && tracksID->GetEntriesFast()>0)
2278 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2279 { //trigger loop starts
2280 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2282 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2283 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2284 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2285 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2287 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2288 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2291 //same event delte-eta, delta-phi plot
2292 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2293 {//same event calculation starts
2294 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2295 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)
2298 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2299 {//same event calculation starts
2300 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2301 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2304 //still in main event loop
2306 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2307 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2308 if (pool && pool->IsReady())
2309 {//start mixing only when pool->IsReady
2310 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2311 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2312 { //pool event loop start
2313 TObjArray* bgTracks = pool->GetEvent(jMix);
2314 if(!bgTracks) continue;
2315 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2316 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2317 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2318 Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2319 }// pool event loop ends mixing case
2321 } //if pool->IsReady() condition ends mixing case
2324 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2326 }//mixing with unidentified particles
2328 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2329 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2330 if (pool1 && pool1->IsReady())
2331 {//start mixing only when pool->IsReady
2332 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2333 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2334 { //pool event loop start
2335 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2336 if(!bgTracks2) continue;
2337 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2338 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2339 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2340 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2342 }// pool event loop ends mixing case
2343 } //if pool1->IsReady() condition ends mixing case
2347 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2349 }//mixing with identified particles
2351 //no. of events analyzed
2352 fEventCounter->Fill(15);
2355 if(tracksUNID) delete tracksUNID;
2357 if(tracksID) delete tracksID;
2360 PostData(1, fOutput);
2363 //________________________________________________________________________
2364 void AliTwoParticlePIDCorr::doAODevent()
2368 AliVEvent *event = InputEvent();
2369 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
2370 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
2372 AliError("Cannot get the AOD event");
2377 fEventCounter->Fill(1);
2379 if (!fPID) return;//this should be available with each event even if we don't do PID selection
2387 gReactionPlane = 999.;
2389 Double_t cent_v0= -999.;
2390 Double_t effcent=1.0;
2392 Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
2395 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
2396 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
2399 // check event cuts and fill event histograms and return the centrality or reference multiplicity value
2400 if((cent_v0 = GetAcceptedEventMultiplicity(aod,kFALSE)) < 0){
2403 effcent=cent_v0;//required for efficiency correction case********Extremely Important
2404 //get the event plane in case of PbPb
2405 if(fRequestEventPlane){
2406 gReactionPlane = GetEventPlane(aod,kFALSE,cent_v0);
2407 if(gReactionPlane==999.) return;
2411 TExMap *trackMap = new TExMap();
2412 // --- track loop for mapping matrix
2415 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2416 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2417 AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
2418 if (!track) continue;
2419 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
2420 if(tracktype!=1) continue;
2422 if(!track) continue;//for safety
2424 Int_t gid = track->GetID();
2425 trackMap->Add(gid,itrk);
2429 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
2430 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
2432 TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
2433 tracksID->SetOwner(kTRUE); // IMPORTANT!
2438 Bool_t fTrigPtmin1=kFALSE;
2439 Bool_t fTrigPtmin2=kFALSE;
2440 Bool_t fTrigPtJet=kFALSE;
2442 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2443 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2444 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2445 if (!track) continue;
2446 fHistQA[11]->Fill(track->GetTPCNcls());
2447 Int_t particletype=-9999;//required for PID filling
2448 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2449 if(tracktype!=1) continue;
2451 if(!track) continue;//for safety
2453 AliAODTrack *PIDtrack=track;//for PID purpose, mainly important for TPC only tracks
2455 if(fFilterBit==128){
2456 Int_t gid1 = track->GetID();
2457 //if(gid1>=0) PIDtrack = track;
2458 PIDtrack = aod->GetTrack(trackMap->GetValue(-1-gid1));
2459 if(!PIDtrack) continue;//for safety; so that each of the TPC only tracks have corresponding global track along with it
2462 //check for eta , phi holes
2463 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2464 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2468 //if no applyefficiency , set the eff factor=1.0
2469 Float_t effmatrix=1.0;
2471 //tag all particles as unidentified that passed filterbit & kinematic cuts
2472 particletype=unidentified;
2474 //To count the no. of tracks having an accepted track in a certain PT(e.g. Jet Pt) range
2475 if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
2476 if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
2477 if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
2480 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
2483 //to reduce memory consumption in pool
2484 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2486 //Clone & Reduce track list(TObjArray) for unidentified particles
2487 Short_t chargeval=0;
2488 if(track->Charge()>0) chargeval=1;
2489 if(track->Charge()<0) chargeval=-1;
2490 if(chargeval==0) continue;
2491 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2492 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
2493 LRCParticlePID* copy = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2494 copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2495 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2498 //now start the particle identificaion process:)
2500 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
2502 Float_t dEdx = PIDtrack->GetTPCsignal();
2503 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2505 //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)
2506 if(HasTOFPID(PIDtrack))
2508 Double_t beta = GetBeta(PIDtrack);
2509 fHistoTOFbeta->Fill(track->Pt(), beta);
2513 //track identification(using nsigma method)
2514 particletype=GetParticle(PIDtrack,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
2516 //2-d TPCTOF map(for each Pt interval)
2517 if(HasTOFPID(PIDtrack)){
2518 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2519 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2520 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2523 //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!!!!!
2524 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)!!!!!!!!!!!
2526 if(fRequestEventPlane){
2527 FillPIDEventPlane(cent_v0,particletype,track->Phi(),gReactionPlane);
2531 //Pt, Eta , Phi distribution of the reconstructed identified particles
2534 if (particletype==SpPion)
2536 fPionPt->Fill(track->Pt());
2537 fPionEta->Fill(track->Eta());
2538 fPionPhi->Fill(track->Phi());
2540 if (particletype==SpKaon)
2542 fKaonPt->Fill(track->Pt());
2543 fKaonEta->Fill(track->Eta());
2544 fKaonPhi->Fill(track->Phi());
2546 if (particletype==SpProton)
2548 fProtonPt->Fill(track->Pt());
2549 fProtonEta->Fill(track->Eta());
2550 fProtonPhi->Fill(track->Phi());
2554 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2556 Short_t chargeval=0;
2557 if(track->Charge()>0) chargeval=1;
2558 if(track->Charge()<0) chargeval=-1;
2559 if(chargeval==0) continue;
2560 if (fapplyTrigefficiency || fapplyAssoefficiency)
2561 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
2562 LRCParticlePID* copy1 = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix,track->GetTPCClusterMapPtr(),track->GetTPCSharedMapPtr());
2563 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2564 tracksID->Add(copy1);
2566 } //track loop ends but still in event loop
2568 if(trackscount<1.0){
2569 if(tracksUNID) delete tracksUNID;
2570 if(tracksID) delete tracksID;
2574 if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
2575 if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
2576 if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
2578 Float_t weightval=1.0;
2582 //fill the centrality/multiplicity distribution of the selected events
2583 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2585 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??????
2587 //count selected events having centrality betn 0-100%
2588 fEventCounter->Fill(13);
2590 //***************************************event no. counting
2591 Bool_t isbaryontrig=kFALSE;
2592 Bool_t ismesontrig=kFALSE;
2593 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2595 if(tracksID && tracksID->GetEntriesFast()>0)
2597 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2598 { //trigger loop starts
2599 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2601 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2602 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2603 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2604 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2606 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2607 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2611 //same event delta-eta-deltaphi plot
2613 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2614 {//same event calculation starts
2615 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2616 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)
2619 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2620 {//same event calculation starts
2621 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2622 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2625 //still in main event loop
2629 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2630 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2631 if (pool && pool->IsReady())
2632 {//start mixing only when pool->IsReady
2633 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2634 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2635 { //pool event loop start
2636 TObjArray* bgTracks = pool->GetEvent(jMix);
2637 if(!bgTracks) continue;
2638 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2639 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2640 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2641 Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2642 }// pool event loop ends mixing case
2643 } //if pool->IsReady() condition ends mixing case
2646 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2648 }//mixing with unidentified particles
2651 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2652 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2653 if (pool1 && pool1->IsReady())
2654 {//start mixing only when pool->IsReady
2655 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2656 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2657 { //pool event loop start
2658 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2659 if(!bgTracks2) continue;
2660 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2661 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2662 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2663 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2665 }// pool event loop ends mixing case
2666 } //if pool1->IsReady() condition ends mixing case
2670 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2672 }//mixing with identified particles
2675 //no. of events analyzed
2676 fEventCounter->Fill(15);
2678 //still in main event loop
2681 if(tracksUNID) delete tracksUNID;
2683 if(tracksID) delete tracksID;
2686 PostData(1, fOutput);
2688 } // *************************event loop ends******************************************//_______________________________________________________________________
2689 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
2691 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
2693 TObjArray* tracksClone = new TObjArray;
2694 tracksClone->SetOwner(kTRUE);
2696 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2698 LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
2699 LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval(),particle->GetTPCPadMap(),particle->GetTPCSharedMap());
2700 copy100->SetUniqueID(particle->GetUniqueID());
2701 tracksClone->Add(copy100);
2707 //--------------------------------------------------------------------------------
2708 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)
2711 //before calling this function check that either trackstrig & tracksasso are available
2713 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
2714 TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
2715 TArrayF eta(input->GetEntriesFast());
2716 for (Int_t i=0; i<input->GetEntriesFast(); i++)
2717 eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
2720 Int_t jmax=trackstrig->GetEntriesFast();
2722 jmax=tracksasso->GetEntriesFast();
2724 // identify K, Lambda candidates and flag those particles
2725 // a TObject bit is used for this
2726 const UInt_t kResonanceDaughterFlag = 1 << 14;
2727 if (fRejectResonanceDaughters > 0)
2729 Double_t resonanceMass = -1;
2730 Double_t massDaughter1 = -1;
2731 Double_t massDaughter2 = -1;
2732 const Double_t interval = 0.02;
2733 switch (fRejectResonanceDaughters)
2735 case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
2736 case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
2737 case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
2738 default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
2741 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2742 trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2743 for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
2744 tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2746 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2748 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2749 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
2751 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2753 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
2754 if (trig->IsEqual(asso)) continue;
2756 if (trig->Charge() * asso->Charge() > 0) continue;
2758 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2760 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
2762 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2764 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
2766 trig->SetBit(kResonanceDaughterFlag);
2767 asso->SetBit(kResonanceDaughterFlag);
2769 // Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
2776 //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
2778 Float_t TriggerPtMin=fminPtTrig;
2779 Float_t TriggerPtMax=fmaxPtTrig;
2780 Int_t HighestPtTriggerIndx=-99999;
2781 TH1* triggerWeighting = 0;
2783 if(fSelectHighestPtTrig || fWeightPerEvent)//**************add this data member to the constructor
2785 if (fWeightPerEvent)
2788 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID) axis = fTHnTrigcount->GetGrid(0)->GetGrid()->GetAxis(2);
2789 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) axis = fTHnTrigcountMCTruthPrim->GetGrid(0)->GetGrid()->GetAxis(2);
2790 triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
2792 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2793 { //trigger loop starts(highest Pt trigger selection)
2794 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2796 Float_t trigpt=trig->Pt();
2797 //to avoid overflow qnd underflow
2798 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2799 Int_t particlepidtrig=trig->getparticle();
2800 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2802 Float_t trigeta=trig->Eta();
2804 // some optimization
2805 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2808 if (fOnlyOneEtaSide != 0)
2810 if (fOnlyOneEtaSide * trigeta < 0)
2813 if (fTriggerSelectCharge != 0)
2814 if (trig->Charge() * fTriggerSelectCharge < 0)
2817 if (fRejectResonanceDaughters > 0)
2818 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2820 if(fSelectHighestPtTrig){
2821 if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
2823 HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
2824 TriggerPtMin=trigpt;
2828 if (fWeightPerEvent) triggerWeighting->Fill(trigpt);
2830 }//trigger loop ends(highest Pt trigger selection)
2832 }//******************(fSelectHighestPtTrig || fWeightPerEvent) condition ends
2835 //two particle correlation filling
2836 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2837 { //trigger loop starts
2838 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2840 Float_t trigpt=trig->Pt();
2841 //to avoid overflow qnd underflow
2842 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2843 Int_t particlepidtrig=trig->getparticle();
2844 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2846 Float_t trigeta=trig->Eta();
2848 // some optimization
2849 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2852 if (fOnlyOneEtaSide != 0)
2854 if (fOnlyOneEtaSide * trigeta < 0)
2857 if (fTriggerSelectCharge != 0)
2858 if (trig->Charge() * fTriggerSelectCharge < 0)
2861 if (fRejectResonanceDaughters > 0)
2862 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2864 if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
2865 if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
2868 Float_t trigphi=trig->Phi();
2869 Float_t trackefftrig=1.0;
2870 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
2872 // Event plane (determine psi bin)
2873 Double_t gPsiMinusPhi = 0.;
2874 Double_t gPsiMinusPhiBin = -10.;
2875 if(fRequestEventPlane){
2876 gPsiMinusPhi = TMath::Abs(trigphi - ReactionPlane);
2877 //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)
2878 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2879 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
2880 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2881 gPsiMinusPhiBin = 0.0;
2883 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2884 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2885 gPsiMinusPhiBin = 0.0;
2888 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
2889 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
2890 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
2891 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
2892 gPsiMinusPhiBin = 1.0;
2894 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
2895 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
2896 gPsiMinusPhiBin = 2.0;
2899 gPsiMinusPhiBin = 3.0;
2901 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
2904 //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
2907 Int_t eventplaneAxis=0;
2908 if(fRequestEventPlane) eventplaneAxis=1;
2909 if(fcontainPIDtrig && SetChargeAxis==0) dim=4+eventplaneAxis;
2910 if(!fcontainPIDtrig && SetChargeAxis==2) dim=4+eventplaneAxis;
2911 if(fcontainPIDtrig && SetChargeAxis==2) dim=5+eventplaneAxis;
2912 trigval= new Double_t[dim];
2915 trigval[2] = trigpt;
2917 if(fRequestEventPlane){
2918 trigval[3] = gPsiMinusPhiBin;
2919 if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = particlepidtrig;
2920 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
2921 if(fcontainPIDtrig && SetChargeAxis==2) {
2922 trigval[4] = particlepidtrig;
2923 trigval[5] = trig->Charge();
2927 if(!fRequestEventPlane){
2928 if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = particlepidtrig;
2929 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
2930 if(fcontainPIDtrig && SetChargeAxis==2) {
2931 trigval[3] = particlepidtrig;
2932 trigval[4] = trig->Charge();
2938 if (fWeightPerEvent)
2940 // leads effectively to a filling of one entry per filled trigger particle pT bin
2941 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(trigval[2]);
2942 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
2943 trackefftrig *= triggerWeighting->GetBinContent(weightBin);
2947 //trigger particle counting for both same and mixed event case;;;;;step=0->same event case;;;;;step=1->mixed event case;;;;;;
2948 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
2949 if(fillup=="trigassoUNID" ) {
2950 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2951 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2954 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
2955 if(fillup=="trigassoUNID" )
2957 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2958 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2961 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
2962 if(fillup=="trigUNIDassoID")
2964 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2965 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2968 //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
2969 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
2970 if(fillup=="trigIDassoID")
2972 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2973 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2976 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
2977 if(fillup=="trigIDassoUNID" )
2979 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2980 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2983 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
2984 if(fillup=="trigIDassoID")
2986 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2987 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2991 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!!!!
2992 if(mixcase==kFALSE) fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig);
2993 if(mixcase==kTRUE && firstTime) fTHnTrigcountMCTruthPrim->Fill(trigval,1,1.0/trackefftrig);
2996 //asso loop starts within trigger loop
2997 for(Int_t j=0;j<jmax;j++)
2999 LRCParticlePID *asso=0;
3001 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
3003 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
3007 //to avoid overflow and underflow
3008 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
3010 //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
3012 if(!tracksasso && j==i) continue;
3014 // 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)
3015 // if (tracksasso && trig->IsEqual(asso)) continue;
3017 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
3020 if (asso->Pt() >= trig->Pt()) continue;
3022 Int_t particlepidasso=asso->getparticle();
3023 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
3026 if (fAssociatedSelectCharge != 0)
3027 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
3029 if (fSelectCharge > 0)
3032 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
3036 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
3042 if (trigeta < 0 && asso->Eta() < trigeta)
3044 if (trigeta > 0 && asso->Eta() > trigeta)
3048 if (fRejectResonanceDaughters > 0)
3049 if (asso->TestBit(kResonanceDaughterFlag))
3051 // Printf("Skipped j=%d", j);
3056 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
3058 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3062 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
3064 fControlConvResoncances->Fill(0.0, mass);
3066 if (mass < 0.04*0.04)
3072 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3074 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
3076 const Float_t kK0smass = 0.4976;
3078 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
3080 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
3082 fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
3084 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
3090 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3092 Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
3093 Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3095 const Float_t kLambdaMass = 1.115;
3097 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
3099 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
3101 fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
3103 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3106 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
3108 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3110 fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
3112 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3117 if (twoTrackEfficiencyCut)
3119 // the variables & cuthave been developed by the HBT group
3120 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
3121 Float_t phi1 = trig->Phi();
3122 Float_t pt1 = trig->Pt();
3123 Float_t charge1 = trig->Charge();
3124 Float_t phi2 = asso->Phi();
3125 Float_t pt2 = asso->Pt();
3126 Float_t charge2 = asso->Charge();
3128 Float_t deta= trigeta - eta[j];
3131 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
3134 // check first boundaries to see if is worth to loop and find the minimum
3135 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMinRadius, bSign);
3136 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMaxRadius, bSign);
3138 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
3140 Float_t dphistarminabs = 1e5;
3141 Float_t dphistarmin = 1e5;
3143 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
3145 for (Double_t rad=0.8; rad<2.51; rad+=0.01)
3147 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
3149 Float_t dphistarabs = TMath::Abs(dphistar);
3151 if (dphistarabs < dphistarminabs)
3153 dphistarmin = dphistar;
3154 dphistarminabs = dphistarabs;
3157 if(mixcase==kFALSE) fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3158 if(mixcase==kTRUE) fTwoTrackDistancePtmix[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3160 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
3162 // 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);
3165 if(mixcase==kFALSE) fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3166 if(mixcase==kTRUE) fTwoTrackDistancePtmix[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3172 //pair sharedfraction cut(only between the trig and asso track)
3173 if(fillup!="trigIDassoIDMCTRUTH")//******************************************NOT for TRUTH MC particles
3175 if(fSharedfraction_Pair_cut>=0){
3176 Bool_t passsharedfractionpaircut=CalculateSharedFraction(trig->GetTPCPadMap(),asso->GetTPCPadMap(),trig->GetTPCSharedMap(),asso->GetTPCSharedMap());
3177 if(!passsharedfractionpaircut) continue;
3180 Float_t weightperevent=weight;
3181 Float_t trackeffasso=1.0;
3182 if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
3183 //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
3184 Float_t deleta=trigeta-eta[j];
3185 Float_t delphi=PhiRange(trigphi-asso->Phi());
3187 Float_t delpt=trigpt-asso->Pt();
3188 //fill it with/without two track efficiency cut
3189 if(mixcase==kFALSE) fTwoTrackDistancePtdip->Fill(deleta, delphi, TMath::Abs(delpt));//for same event
3190 if(mixcase==kTRUE) fTwoTrackDistancePtdipmix->Fill(deleta, delphi, TMath::Abs(delpt));//for mixed event
3192 //here get the two particle efficiency correction factor
3193 Float_t effweight=trackefftrig*trackeffasso*weightperevent;
3194 // if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
3196 Int_t dimused=kTrackVariablesPair+eventplaneAxis;
3197 vars= new Double_t[dimused];
3206 if(fRequestEventPlane)
3208 vars[6]=gPsiMinusPhiBin;
3212 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
3213 vars[dimension]=trig->Charge();
3214 vars[dimension+1]=asso->Charge();
3216 if(fcontainPIDtrig && !fcontainPIDasso){
3217 vars[dimension]=particlepidtrig;
3218 if(SetChargeAxis==2){
3219 vars[dimension+1]=trig->Charge();
3220 vars[dimension+2]=asso->Charge();
3223 if(!fcontainPIDtrig && fcontainPIDasso){
3224 vars[dimension]=particlepidasso;
3225 if(SetChargeAxis==2){
3226 vars[dimension+1]=trig->Charge();
3227 vars[dimension+2]=asso->Charge();
3230 if(fcontainPIDtrig && fcontainPIDasso){
3231 vars[dimension]=particlepidtrig;
3232 vars[dimension+1]=particlepidasso;
3233 if(SetChargeAxis==2){
3234 vars[dimension+2]=trig->Charge();
3235 vars[dimension+3]=asso->Charge();
3239 if (fWeightPerEvent)
3241 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
3242 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3243 effweight *= triggerWeighting->GetBinContent(weightBin);
3247 //Fill Correlation ThnSparses
3248 if(fillup=="trigassoUNID")
3250 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,0,1.0/effweight);
3251 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
3253 if(fillup=="trigIDassoID")
3255 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,0,1.0/effweight);
3256 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
3258 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
3259 {//in this case effweight should be 1 always
3260 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
3261 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
3263 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
3265 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
3266 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
3272 }//trigger loop ends
3274 if (triggerWeighting)
3276 delete triggerWeighting;
3277 triggerWeighting = 0;
3281 //------------------------------------------------------------------------------------------------
3282 Bool_t AliTwoParticlePIDCorr:: CalculateSharedFraction(const TBits *triggerPadMap,const TBits *assocPadMap,const TBits *triggerShareMap,const TBits *assocShareMap)
3283 {//source code-AliFemtoShareQualityPairCut.cxx
3285 Double_t nofsharedhits=0;
3287 for(UInt_t imap=0;imap< (triggerPadMap->GetNbits() );imap++)
3289 //if they are in same pad
3290 //cout<<triggerPadMap->TestBitNumber(imap)<<" "<< assocPadMap->TestBitNumber(imap)<<endl;
3291 if (triggerPadMap->TestBitNumber(imap) &&
3292 assocPadMap->TestBitNumber(imap))
3295 //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
3296 if (triggerShareMap->TestBitNumber(imap) &&
3297 assocShareMap->TestBitNumber(imap))
3299 //cout<<triggerShareMap->TestBitNumber(imap)<<" "<<assocShareMap->TestBitNumber(imap)<<endl;
3316 //cout<< (triggerPadMap->TestBitNumber(imap) || assocPadMap->TestBitNumber(imap))<<endl;
3317 else if (triggerPadMap->TestBitNumber(imap) ||
3318 assocPadMap->TestBitNumber(imap)) {
3319 // One track has a hit, the other does not
3322 //cout<<"No hits :"<<nofhits<<endl;
3330 Double_t SharedFraction=0.0;
3331 if(nofhits>0) SharedFraction=(nofsharedhits/nofhits);
3333 //cout<<"Fraction shared hits :"<<SharedFraction<<endl;
3335 if(SharedFraction>fSharedfraction_Pair_cut) return kFALSE;
3341 //________________________________________________________________________________________________
3342 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
3344 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
3346 Float_t effvalue=1.;
3348 if(parpid==unidentified)
3350 effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
3351 effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
3352 effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
3353 effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
3354 effvalue=effcorection[5]->GetBinContent(effVars);
3356 if(parpid==SpPion || parpid==SpKaon)
3358 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3360 effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
3361 effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
3362 effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
3363 effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
3364 effvalue=effcorection[3]->GetBinContent(effVars);
3369 effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
3370 effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
3371 effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
3372 effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
3373 effvalue=effcorection[0]->GetBinContent(effVars);
3378 effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
3379 effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
3380 effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
3381 effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
3382 effvalue=effcorection[1]->GetBinContent(effVars);
3387 if(parpid==SpProton)
3389 effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
3390 effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
3391 effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
3392 effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
3393 effvalue=effcorection[2]->GetBinContent(effVars);
3396 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3398 if(parpid==SpProton || parpid==SpKaon)
3400 effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
3401 effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
3402 effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
3403 effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
3404 effvalue=effcorection[4]->GetBinContent(effVars);
3407 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
3408 if(effvalue==0.) effvalue=1.;
3413 //---------------------------------------------------------------------------------
3415 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield, Bool_t fill)
3418 if(!track) return 0;
3419 Bool_t trackOK = track->TestFilterBit(fFilterBit);
3420 if(!trackOK) return 0;
3421 if (fTrackStatus != 0 && !CheckTrack(track)) return 0;
3422 //select only primary traks(for data & reco MC tracks)
3423 if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
3424 if(track->Charge()==0) return 0;
3425 if (fill) fHistQA[12]->Fill(track->GetTPCNcls());
3428 dz = track->ZAtDCA();
3429 if (fill) fHistQA[6]->Fill(dxy);
3430 if (fill) fHistQA[7]->Fill(dz);
3431 Float_t chi2ndf = track->Chi2perNDF();
3432 if (fill) fHistQA[13]->Fill(chi2ndf);
3433 // Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
3434 Float_t nCrossedRowsTPC = track->GetTPCNCrossedRows();
3435 if (fill) fHistQA[14]->Fill(nCrossedRowsTPC);
3436 //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
3437 if (track->GetTPCNclsF()>0) {
3438 Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
3439 if (fill) fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
3442 Float_t pt=track->Pt();
3443 if(pt< fminPt || pt> fmaxPt) return 0;
3444 if(TMath::Abs(track->Eta())> fmaxeta) return 0;
3445 if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
3446 //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required???
3448 if (fdcacut && fDCAXYCut)
3455 AliAODTrack* clone =(AliAODTrack*) track->Clone();
3456 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
3461 // Printf("%f", ((AliAODTrack*)part)->DCA());
3462 // Printf("%f", pos[0]);
3463 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
3467 if (fSharedClusterCut >= 0)
3469 Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
3470 if (frac > fSharedClusterCut)
3474 // Rejects tracks with shared clusters after filling a control histogram
3475 // This overload is used for primaries
3477 // Get the shared maps
3478 const TBits sharedMap = track->GetTPCSharedMap();
3479 // Fill a control histogram
3480 fPriHistShare->Fill(sharedMap.CountBits());
3482 // Reject shared clusters
3483 if (fSharedTPCmapCut >= 0)
3485 if((sharedMap.CountBits()) >= 1) return 0;// Bad track, has too many shared clusters!
3488 if (fill) fHistQA[8]->Fill(pt);
3489 if (fill) fHistQA[9]->Fill(track->Eta());
3490 if (fill) fHistQA[10]->Fill(track->Phi());
3493 //________________________________________________________________________________
3494 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
3496 //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 )
3497 Float_t pt=track->Pt();
3499 //plot the separation power
3501 Double_t bethe[AliPID::kSPECIES]={0.};
3502 Double_t sigma_TPC[AliPID::kSPECIES]={0.};
3504 Double_t Pi_Ka_sep[NSigmaPIDType+1]={0.};
3505 Double_t Pi_Pr_sep[NSigmaPIDType+1]={0.};
3506 Double_t Ka_Pr_sep[NSigmaPIDType+1]={0.};
3509 Double_t ptpc = track->GetTPCmomentum();
3510 Int_t dEdxN = track->GetTPCsignalN();
3511 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
3512 bethe[ipart] = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
3513 //Double_t diff = dEdx - bethe;
3514 sigma_TPC[ipart] = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
3515 //nSigma[ipart] = diff / sigma;
3517 Pi_Ka_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kKaon])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kKaon])/2.0);
3518 Pi_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kProton])/2.0);
3519 Ka_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kKaon]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kKaon]+sigma_TPC[AliPID::kProton])/2.0);
3522 Double_t sigma_TOF[AliPID::kSPECIES]={0.};
3524 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3526 Double_t timei[AliPID::kSPECIES];
3527 track->GetIntegratedTimes(timei);
3528 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { sigma_TOF[ipart]= fPID->GetTOFResponse().GetExpectedSigma(track->P(), timei[ipart], AliPID::ParticleMass(ipart));}
3529 Pi_Ka_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kKaon])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kKaon])/2.0);
3530 Pi_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kProton])/2.0);
3531 Ka_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kKaon]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kKaon]+sigma_TOF[AliPID::kProton])/2.0);
3533 Pi_Ka_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Ka_sep[NSigmaTPC]*Pi_Ka_sep[NSigmaTPC]+Pi_Ka_sep[NSigmaTOF]*Pi_Ka_sep[NSigmaTOF]);
3534 Pi_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Pr_sep[NSigmaTPC]*Pi_Pr_sep[NSigmaTPC]+Pi_Pr_sep[NSigmaTOF]*Pi_Pr_sep[NSigmaTOF]);
3535 Ka_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Ka_Pr_sep[NSigmaTPC]*Ka_Pr_sep[NSigmaTPC]+Ka_Pr_sep[NSigmaTOF]*Ka_Pr_sep[NSigmaTOF]);
3539 //fill separation power histograms
3540 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3542 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3543 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3544 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3545 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3546 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3547 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3549 if(HasTOFPID(track) && pt>fPtTOFPIDmin && ipid!=0){
3550 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3551 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3552 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3553 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3554 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3555 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3560 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
3561 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
3562 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
3563 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
3565 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
3566 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
3568 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3571 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
3572 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
3573 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
3575 nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
3576 nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
3577 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
3583 // if TOF is missing and below fPtTOFPID only the TPC information is used
3584 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
3585 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
3586 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
3590 //set data member fnsigmas
3591 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
3592 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
3593 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
3595 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
3596 fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
3597 fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
3598 fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
3600 //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)
3601 fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
3602 fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
3603 fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
3606 //Fill NSigma SeparationPlot
3607 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3608 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3609 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3610 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
3611 h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
3617 //----------------------------------------------------------------------------
3618 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
3620 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3621 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
3622 //get the identity of the particle with the minimum Nsigma
3623 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3626 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3627 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
3628 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
3631 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3632 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
3633 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
3635 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3636 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3637 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3638 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3640 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3641 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3642 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3643 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3648 if(fdiffPIDcutvalues){
3649 if(track->Pt()<=4) fNSigmaPID=fPIDCutval1;
3650 if(track->Pt()>4 && track->Pt()<=6) fNSigmaPID=fPIDCutval2;
3651 if(track->Pt()>6 && track->Pt()<=8) fNSigmaPID=fPIDCutval3;
3652 if(track->Pt()>8) fNSigmaPID=fPIDCutval4;
3655 // guess the particle based on the smaller nsigma (within fNSigmaPID)
3656 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
3658 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
3659 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)
3661 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3662 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3663 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
3664 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
3669 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
3671 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3672 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3673 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
3674 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
3679 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
3681 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3682 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3683 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
3684 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
3690 // else, return undefined
3696 //------------------------------------------------------------------------------------------
3697 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
3698 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3700 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
3702 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
3704 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
3707 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
3709 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3712 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3713 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
3714 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
3717 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3718 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
3719 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
3721 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3722 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3723 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3724 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3726 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3727 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3728 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3729 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3733 // Actually the tracks in the overlapping region(in TPC-TOF nSigma plane) will be ignored
3735 if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
3736 if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
3737 if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
3742 //fill NSigma distr for double counting
3743 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3744 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
3745 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3746 if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)))continue;//not filling TOF and combined if no TOF PID
3747 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
3748 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
3755 return fHasDoubleCounting;
3758 //////////////////////////////////////////////////////////////////////////////////////////////////
3760 Bool_t* AliTwoParticlePIDCorr::GetAllCompatibleIdentitiesNSigma(AliAODTrack * trk,Bool_t FIllQAHistos){
3761 //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
3762 Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
3763 IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
3767 //////////////////////////////////////////////////////////////////////////////////////////////////
3769 UInt_t AliTwoParticlePIDCorr::CalcPIDCombined(AliAODTrack *track, Int_t detMask, Double_t* prob) const{
3771 // Bayesian PID calculation
3773 for(Int_t i=0;i<AliPID::kSPECIES;i++)
3777 fPIDCombined->SetDetectorMask(detMask);
3779 return fPIDCombined->ComputeProbabilities((AliAODTrack*)track, fPID, prob);
3782 //////////////////////////////////////////////////////////////////////////////////////////////////
3784 Int_t AliTwoParticlePIDCorr::GetIDBayes(AliAODTrack * trk, Bool_t FIllQAHistos){
3786 Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
3789 //Filling of Probability histos
3790 Double_t probTPC[AliPID::kSPECIES]={0.};
3791 Double_t probTOF[AliPID::kSPECIES]={0.};
3792 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
3794 UInt_t detUsedTPC = 0;
3795 UInt_t detUsedTOF = 0;
3796 UInt_t detUsedTPCTOF = 0;
3798 //get the TPC probability
3799 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
3800 detUsedTPC = fPIDCombined->ComputeProbabilities(trk, fPID, probTPC);
3801 if(detUsedTPC == AliPIDResponse::kDetTPC)
3803 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3805 TH2F *h=GetHistogram2D(Form("probBayes_TPC_%d",ipart));
3806 if(ipart==0) h->Fill(trk->Pt(),probTPC[AliPID::kPion]);
3807 if(ipart==1) h->Fill(trk->Pt(),probTPC[AliPID::kKaon]);
3808 if(ipart==2) h->Fill(trk->Pt(),probTPC[AliPID::kProton]);
3812 //get the TOF probability
3813 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
3814 detUsedTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTOF);
3815 if(detUsedTOF == AliPIDResponse::kDetTOF)
3817 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3818 TH2F *h=GetHistogram2D(Form("probBayes_TOF_%d",ipart));
3819 if(ipart==0) h->Fill(trk->Pt(),probTOF[AliPID::kPion]);
3820 if(ipart==1) h->Fill(trk->Pt(),probTOF[AliPID::kKaon]);
3821 if(ipart==2) h->Fill(trk->Pt(),probTOF[AliPID::kProton]);
3825 //get the TPC-TOF probability
3826 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
3827 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTPCTOF);
3828 if(detUsedTPCTOF == (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))
3830 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3831 TH2F *h=GetHistogram2D(Form("probBayes_TPCTOF_%d",ipart));
3832 if(ipart==0) h->Fill(trk->Pt(),probTPCTOF[AliPID::kPion]);
3833 if(ipart==1) h->Fill(trk->Pt(),probTPCTOF[AliPID::kKaon]);
3834 if(ipart==2) h->Fill(trk->Pt(),probTPCTOF[AliPID::kProton]);
3839 Double_t probBayes[AliPID::kSPECIES];
3842 if(HasTOFPID(trk) && trk->Pt()>fPtTOFPIDmin){//use TOF information
3843 detUsed = CalcPIDCombined(trk, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
3844 if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return SpUndefined;//check that TPC and TOF are used
3846 detUsed = CalcPIDCombined(trk,AliPIDResponse::kDetTPC, probBayes);
3847 if(detUsed != AliPIDResponse::kDetTPC)return SpUndefined;//check that TPC is used
3850 //the probability has to be normalized to one, we check it
3852 for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
3853 if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
3854 AliFatal("Bayesian probability not normalized to one");
3857 if(fdiffPIDcutvalues){
3858 if(trk->Pt()<=4) fBayesCut=fPIDCutval1;
3859 if(trk->Pt()>4 && trk->Pt()<=6) fBayesCut=fPIDCutval2;
3860 if(trk->Pt()>6 && trk->Pt()<=8) fBayesCut=fPIDCutval3;
3861 if(trk->Pt()>8) fBayesCut=fPIDCutval4;
3865 //probabilities are normalized to one, if the cut is above .5 there is no problem
3866 if(probBayes[AliPID::kPion]>fBayesCut && IDs[SpPion]==1){
3867 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpPion));
3868 h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
3871 else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[SpKaon]==1){
3872 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpKaon));
3873 h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
3876 else if(probBayes[AliPID::kProton]>fBayesCut && IDs[SpProton]==1){
3877 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpProton));
3878 h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
3887 //////////////////////////////////////////////////////////////////////////////////////////////////
3888 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
3889 //return the specie according to the minimum nsigma value
3890 //no double counting, this has to be evaluated using CheckDoubleCounting()
3892 Int_t ID=SpUndefined;
3894 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
3898 if(fPIDType==Bayes){//use bayesianPID
3901 AliFatal("PIDCombined object has to be set in the steering macro");
3904 ID = GetIDBayes(trk,FIllQAHistos);
3906 }else{ //use nsigma PID
3908 ID=FindMinNSigma(trk,FIllQAHistos);
3909 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
3911 HasDC=GetDoubleCounting(trk,FIllQAHistos);
3912 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3913 if(HasDC[ipart]==kTRUE) ID = SpUndefined;
3917 //Fill PID signal plot
3918 if(ID != SpUndefined){
3919 for(Int_t idet=0;idet<fNDetectors;idet++){
3920 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
3921 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3922 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3923 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3926 //Fill PID signal plot without cuts
3927 for(Int_t idet=0;idet<fNDetectors;idet++){
3928 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
3929 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3930 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3931 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3937 //-------------------------------------------------------------------------------------
3939 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
3942 AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
3943 if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
3944 //ULong_t status=track->GetStatus();
3945 //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
3946 //if (track->GetTPCsignal() <= 0.) return kFALSE;
3947 // 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.
3951 //___________________________________________________________
3954 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
3956 // check TOF matched track
3957 //ULong_t status=track->GetStatus();
3958 //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
3959 AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
3960 if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
3961 if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
3962 //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
3963 //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
3964 // if (probMis > 0.01) return kFALSE;
3965 if(fRemoveTracksT0Fill)
3967 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
3968 if (startTimeMask < 0)return kFALSE;
3973 //________________________________________________________________________
3974 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
3976 //it is called only when TOF PID is available
3977 //TOF beta calculation
3978 Double_t tofTime=track->GetTOFsignal();
3980 Double_t c=TMath::C()*1.E-9;// m/ns
3981 Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
3982 Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
3983 tofTime -= startTime; // subtract startTime to the signal
3984 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
3990 Double_t p = track->P();
3991 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
3993 track->GetIntegratedTimes(timei);
3994 return timei[0]/time;
3997 //------------------------------------------------------------------------------------------------------
3999 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)
4001 // calculate inv mass squared
4002 // same can be achieved, but with more computing time with
4003 /*TLorentzVector photon, p1, p2;
4004 p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
4005 p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
4009 Float_t tantheta1 = 1e10;
4011 if (eta1 < -1e-10 || eta1 > 1e-10)
4012 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
4014 Float_t tantheta2 = 1e10;
4015 if (eta2 < -1e-10 || eta2 > 1e-10)
4016 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
4018 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4019 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4021 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 ) ) );
4025 //---------------------------------------------------------------------------------
4027 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)
4029 // calculate inv mass squared approximately
4031 Float_t tantheta1 = 1e10;
4033 if (eta1 < -1e-10 || eta1 > 1e-10)
4035 Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
4036 tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4039 Float_t tantheta2 = 1e10;
4040 if (eta2 < -1e-10 || eta2 > 1e-10)
4042 Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
4043 tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
4046 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
4047 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
4050 Float_t deltaPhi = TMath::Abs(phi1 - phi2);
4051 while (deltaPhi > TMath::TwoPi())
4052 deltaPhi -= TMath::TwoPi();
4053 if (deltaPhi > TMath::Pi())
4054 deltaPhi = TMath::TwoPi() - deltaPhi;
4056 Float_t cosDeltaPhi = 0;
4057 if (deltaPhi < TMath::Pi()/3)
4058 cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
4059 else if (deltaPhi < 2*TMath::Pi()/3)
4060 cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
4062 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
4064 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
4066 // Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
4070 //--------------------------------------------------------------------------------
4071 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)
4074 // calculates dphistar
4077 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
4079 static const Double_t kPi = TMath::Pi();
4082 // if (dphistar > 2 * kPi)
4083 // dphistar -= 2 * kPi;
4084 // if (dphistar < -2 * kPi)
4085 // dphistar += 2 * kPi;
4088 dphistar = kPi * 2 - dphistar;
4089 if (dphistar < -kPi)
4090 dphistar = -kPi * 2 - dphistar;
4091 if (dphistar > kPi) // might look funny but is needed
4092 dphistar = kPi * 2 - dphistar;
4096 //_________________________________________________________________________
4098 void AliTwoParticlePIDCorr ::DefineEventPool()
4100 Int_t MaxNofEvents=1000;
4101 const Int_t NofVrtxBins=10+(1+10)*2;
4102 Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
4103 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
4104 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210
4106 //default values are for centrality
4107 Int_t NofCentBins=15;
4108 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
4110 if(fCentralityMethod.EndsWith("_MANUAL"))
4112 Int_t NofCentBins=9;
4113 CentralityBins[NofCentBins+1]={0.,9.,14.,19.,26.,34.,44.,58.,80.,500.};//Is This binning is fine for pp, or we don't require them....
4115 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
4120 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
4122 //if(!fPoolMgr) return kFALSE;
4127 //------------------------------------------------------------------------
4128 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
4130 // This method is a copy from AliUEHist::GetBinning
4131 // takes the binning from <configuration> identified by <tag>
4132 // configuration syntax example:
4133 // 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
4136 // returns bin edges which have to be deleted by the caller
4138 TString config(configuration);
4139 TObjArray* lines = config.Tokenize("\n");
4140 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
4142 TString line(lines->At(i)->GetName());
4143 if (line.BeginsWith(TString(tag) + ":"))
4145 line.Remove(0, strlen(tag) + 1);
4146 line.ReplaceAll(" ", "");
4147 TObjArray* binning = line.Tokenize(",");
4148 Double_t* bins = new Double_t[binning->GetEntriesFast()];
4149 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
4150 bins[j] = TString(binning->At(j)->GetName()).Atof();
4152 nBins = binning->GetEntriesFast() - 1;
4161 AliFatal(Form("Tag %s not found in %s", tag, configuration));
4165 //____________________________________________________________________
4167 Bool_t AliTwoParticlePIDCorr::CheckTrack(AliAODTrack * part)
4169 // check if the track status flags are set
4171 UInt_t status=((AliAODTrack*)part)->GetStatus();
4172 if ((status & fTrackStatus) == fTrackStatus)
4176 //________________________________________________________________________
4178 Bool_t AliTwoParticlePIDCorr::AcceptEventCentralityWeight(Double_t centrality)
4180 // rejects "randomly" events such that the centrality gets flat
4181 // uses fCentralityWeights histogram
4183 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
4185 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
4186 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
4188 Bool_t result = kFALSE;
4189 if (centralityDigits < weight)
4192 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
4197 //____________________________________________________________________
4198 void AliTwoParticlePIDCorr::ShiftTracks(TObjArray* tracks, Double_t angle)
4200 // shifts the phi angle of all tracks by angle
4201 // 0 <= angle <= 2pi
4203 for (Int_t i=0; i<tracks->GetEntriesFast(); ++i)
4205 LRCParticlePID *part=(LRCParticlePID*)(tracks->UncheckedAt(i));
4207 Double_t newAngle = part->Phi() + angle;
4208 if (newAngle >= TMath::TwoPi())
4209 newAngle -= TMath::TwoPi();
4211 part->SetPhi(newAngle);
4216 //________________________________________________________________________
4217 void AliTwoParticlePIDCorr::SetVZEROCalibrationFile(const char* filename,const char* lhcPeriod) {
4218 //Function to setup the VZERO gain equalization
4219 //============Get the equilization map============//
4220 TFile *calibrationFile = TFile::Open(filename);
4221 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
4222 Printf("No calibration file found!!!");
4226 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
4228 Printf("Calibration TList not found!!!");
4232 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
4233 if(!fHistVZEROAGainEqualizationMap) {
4234 Printf("VZERO-A calibration object not found!!!");
4237 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
4238 if(!fHistVZEROCGainEqualizationMap) {
4239 Printf("VZERO-C calibration object not found!!!");
4243 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
4244 if(!fHistVZEROChannelGainEqualizationMap) {
4245 Printf("VZERO channel calibration object not found!!!");
4250 //________________________________________________________________________
4251 Double_t AliTwoParticlePIDCorr::GetChannelEqualizationFactor(Int_t run,Int_t channel) {
4253 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4255 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
4256 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4257 if(gRunNumber == run)
4258 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
4264 //________________________________________________________________________
4265 Double_t AliTwoParticlePIDCorr::GetEqualizationFactor(Int_t run, const char* side) {
4267 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4269 TString gVZEROSide = side;
4270 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
4271 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4272 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
4273 if(gRunNumber == run) {
4274 if(gVZEROSide == "A")
4275 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
4276 else if(gVZEROSide == "C")
4277 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
4283 //________________________________________________________________________
4284 Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliAODEvent *event){
4285 //Function that returns the reference multiplicity from AODs (data or reco MC, Not for Truth)
4286 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
4287 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
4288 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
4290 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
4292 Printf("ERROR: AOD header not available");
4295 Int_t gRunNumber = header->GetRunNumber();
4296 Float_t bSign1=header->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
4299 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++)
4300 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
4301 AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
4302 if (!track) continue;
4303 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
4304 if(tracktype!=1) continue;
4306 if(!track) continue;//for safety
4308 gRefMultiplicityTPC += 1.0;
4312 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
4313 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
4314 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
4315 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
4318 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
4319 else if(iChannel >= 32)
4320 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
4323 //Equalization of gain
4324 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
4326 gRefMultiplicityVZEROA /= gFactorA;
4327 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
4329 gRefMultiplicityVZEROC /= gFactorC;
4330 if((gFactorA != 0)&&(gFactorC != 0))
4331 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
4334 //EQVZERO vs TPC multiplicity
4335 fHistEQVZEROvsTPCmultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
4336 fHistEQVZEROAvsTPCmultiplicity->Fill(gRefMultiplicityVZEROA,gRefMultiplicityTPC);
4337 fHistEQVZEROCvsTPCmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityTPC);
4339 //EQVZERO vs VZERO multiplicity
4340 fHistVZEROCvsEQVZEROCmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),gRefMultiplicityVZEROC);
4341 fHistVZEROAvsEQVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0A(),gRefMultiplicityVZEROA);
4343 //VZEROC vs VZEROA multiplicity
4344 fHistVZEROCvsVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),event->GetVZEROData()->GetMTotV0A());
4346 //EQVZEROC vs EQVZEROA multiplicity
4347 fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
4350 if(fCentralityMethod == "TRACKS_MANUAL") {
4351 gRefMultiplicity = gRefMultiplicityTPC;
4352 fHistRefmult->Fill(3.,gRefMultiplicityTPC);
4354 else if(fCentralityMethod == "V0M_MANUAL"){
4355 gRefMultiplicity = gRefMultiplicityVZERO;
4356 fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
4359 else if(fCentralityMethod == "V0A_MANUAL"){
4360 gRefMultiplicity = gRefMultiplicityVZEROA;
4361 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
4364 else if(fCentralityMethod == "V0C_MANUAL"){
4365 gRefMultiplicity = gRefMultiplicityVZEROC;
4366 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
4369 gRefMultiplicity = gRefMultiplicityTPC;
4371 return gRefMultiplicity;
4374 //-------------------------------------------------------------------------------------------------------
4375 Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliAODEvent *event, Bool_t truth){
4377 if(!event) return -1;
4378 // get centrality object and check quality
4379 Double_t cent_v0=-1;
4380 Double_t nooftrackstruth=0;
4381 Bool_t shift_to_TRACKS_MANUAL=kFALSE;//in case of wrong setting automatic shift to Tracks_Manual method
4383 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
4386 if(fSampleType=="pp_7" && fPPVsMultUtils)
4387 {//for pp 7 TeV case only using Alianalysisutils class
4388 if(fAnalysisUtils) cent_v0 = fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,fCentralityMethod);
4390 fHistCentStats->Fill(0.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0A"));
4391 fHistCentStats->Fill(1.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0C"));
4392 fHistCentStats->Fill(2.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0M"));
4393 fHistCentStats->Fill(3.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0AEq"));//only available for LHC10d at present (Quantile info)
4394 fHistCentStats->Fill(4.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0CEq"));//only available for LHC10d at present (Quantile info)
4395 fHistCentStats->Fill(5.,fAnalysisUtils->GetMultiplicityPercentile((AliVEvent*)event,"V0MEq"));//only available for LHC10d at present (Quantile info)
4398 if(fSampleType=="pPb" || fSampleType=="PbPb")
4400 AliCentrality *centralityObj=0;
4401 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
4402 if(!header) return -1;
4403 centralityObj = header->GetCentralityP();
4404 // if (centrality->GetQuality() != 0) return ;
4406 fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
4407 fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
4408 fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
4409 fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
4410 fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
4411 fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
4413 fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
4414 fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
4416 cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
4420 else shift_to_TRACKS_MANUAL=kTRUE;
4422 }//centralitymethod condition
4424 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
4426 if(!truth){//for data or RecoMC
4427 cent_v0 = GetReferenceMultiplicityVZEROFromAOD(event);
4428 }//for data or RecoMC
4430 if(truth && (fAnalysisType == "MCAOD")){//condition for TRUTH case
4431 //check for TClonesArray(truth track MC information)
4432 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
4434 //AliFatal("Error: MC particles branch not found!\n");
4437 //now process the truth particles(for both efficiency & correlation function)
4438 Int_t nMCTrack = fArrayMC->GetEntriesFast();
4440 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
4441 {//MC truth track loop starts
4443 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
4446 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
4450 //consider only charged particles
4451 if(partMC->Charge() == 0) continue;
4453 //consider only primary particles; neglect all secondary particles including from weak decays
4454 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
4457 //remove injected signals(primaries above <maxLabel>)
4458 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
4461 Bool_t isduplicate=kFALSE;
4462 if (fRemoveDuplicates)
4464 for (Int_t j=iMC+1; j<nMCTrack; ++j)
4465 {//2nd trutuh loop starts
4466 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
4468 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
4471 if (partMC->GetLabel() == partMC2->GetLabel())
4476 }//2nd truth loop ends
4478 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
4481 if (fCentralityMethod=="V0M_MANUAL") {
4482 if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8) continue;
4483 if (partMC->Eta() < -3.7 || partMC->Eta() > -1.7) continue;
4485 else if (fCentralityMethod=="V0A_MANUAL") {
4486 if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8) continue;}
4487 else if (fCentralityMethod=="V0C_MANUAL") {
4488 if(partMC->Eta() > -1.7 || partMC->Eta() < -3.7) continue;}
4489 else if (fCentralityMethod=="TRACKS_MANUAL") {
4490 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4491 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
4493 else{//basically returns the tracks manual case
4494 //give only kinematic cuts at the truth level
4495 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4496 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
4499 //To determine multiplicity in case of PP
4500 nooftrackstruth+= 1;;
4502 }//truth track loop ends
4503 cent_v0=nooftrackstruth;
4505 }//condition for TRUTH case
4507 }//end of MANUAL method
4509 else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC
4511 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4515 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4518 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
4519 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
4520 AliError("Event header not found. Skipping this event.");
4524 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4527 if (collGeometry) cent_v0 = collGeometry->ImpactParameter();
4529 }//end of Impact parameter method
4536 //-----------------------------------------------------------------------------------------
4537 Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliAODEvent *aod,Bool_t truth){
4538 //do the event selection(zvtx, pileup, centrality/multiplicity cut) and then return the value of the centrality of that event
4541 Float_t gRefMultiplicity = -1.;
4543 // check first event in chunk (is not needed for new reconstructions)
4544 if(fCheckFirstEventInChunk){
4545 AliAnalysisUtils ut;
4546 if(ut.IsFirstEventInChunk(aod))
4551 AliAnalysisUtils ut;
4552 ut.SetUseMVPlpSelection(kTRUE);
4553 ut.SetUseOutOfBunchPileUp(kTRUE);
4554 if(ut.IsPileUpEvent(aod))
4558 //count events after pileup selection
4559 fEventCounter->Fill(3);
4561 //vertex selection(is it fine for PP?)
4562 if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4563 trkVtx = aod->GetPrimaryVertex();
4564 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
4565 TString vtxTtl = trkVtx->GetTitle();
4566 if (!vtxTtl.Contains("VertexerTracks")) return -1;
4567 zvtx = trkVtx->GetZ();
4568 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
4569 if (!spdVtx || spdVtx->GetNContributors()<=0) return -1;
4570 TString vtxTyp = spdVtx->GetTitle();
4571 Double_t cov[6]={0};
4572 spdVtx->GetCovarianceMatrix(cov);
4573 Double_t zRes = TMath::Sqrt(cov[5]);
4574 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return -1;
4575 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return -1;
4577 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
4578 Int_t nVertex = aod->GetNumberOfVertices();
4580 trkVtx = aod->GetPrimaryVertex();
4581 Int_t nTracksPrim = trkVtx->GetNContributors();
4582 zvtx = trkVtx->GetZ();
4583 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
4584 // Reject TPC only vertex
4585 TString name(trkVtx->GetName());
4586 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return -1;
4588 // Select a quality vertex by number of tracks?
4589 if( nTracksPrim < fnTracksVertex ) {
4590 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
4593 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
4594 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
4596 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
4601 else if(fVertextype==0){//default case
4602 trkVtx = aod->GetPrimaryVertex();
4603 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
4604 zvtx = trkVtx->GetZ();
4606 trkVtx->GetCovarianceMatrix(fCov);
4607 if(fCov[5] == 0) return -1;//proper vertex resolution
4610 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
4611 return -1;//as there is no proper sample type
4614 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
4616 //count events having a proper vertex
4617 fEventCounter->Fill(5);
4619 if (TMath::Abs(zvtx) > fzvtxcut) return -1;
4621 //count events after vertex cut
4622 fEventCounter->Fill(7);
4625 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4627 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
4629 //get the centrality or multiplicity
4630 if(truth) {gRefMultiplicity = GetRefMultiOrCentrality(aod,kTRUE);}//kTRUE-->for Truth case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
4632 else {gRefMultiplicity = GetRefMultiOrCentrality(aod,kFALSE);}//kFALSE-->for data and RecoMc case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
4634 if(gRefMultiplicity<0) return -1;
4636 // take events only within the multiplicity class mentioned in the custom binning
4637 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
4639 //count events having proper centrality/ref multiplicity
4640 fEventCounter->Fill(9);
4643 // centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
4644 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
4646 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
4650 //count events after rejection due to centrality weighting
4651 fEventCounter->Fill(11);
4653 return gRefMultiplicity;
4656 //--------------------------------------------------------------------------------------------------------
4657 Float_t AliTwoParticlePIDCorr::GetEventPlane(AliAODEvent *event,Bool_t truth, Double_t v0Centr)
4659 // Get the event plane
4660 //reset Q vector info
4662 Int_t run = event->GetRunNumber();
4665 // Load the calibrations run dependent
4666 if(! fIsAfter2011) OpenInfoCalbration(run);
4672 if (v0Centr < 80){ // analysis only for 0-80% centrality classes
4674 if(v0Centr < 5) iC = 0;
4675 else if(v0Centr < 10) iC = 1;
4676 else if(v0Centr < 20) iC = 2;
4677 else if(v0Centr < 30) iC = 3;
4678 else if(v0Centr < 40) iC = 4;
4679 else if(v0Centr < 50) iC = 5;
4680 else if(v0Centr < 60) iC = 6;
4681 else if(v0Centr < 70) iC = 7;
4687 //reset Q vector info
4688 Double_t Qxa2 = 0, Qya2 = 0;
4689 Double_t Qxc2 = 0, Qyc2 = 0;
4690 Double_t Qxa3 = 0, Qya3 = 0;
4691 Double_t Qxc3 = 0, Qyc3 = 0;
4694 //MC: from reaction plane
4697 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4699 evplaneMC = header->GetReactionPlaneAngle();//[0, 360]
4700 //make it within [-pi/2,pi/2] to make it general
4701 if(evplaneMC > TMath::Pi()/2 && evplaneMC <= TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
4702 else if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
4703 fHistEventPlaneTruth->Fill(iC,evplaneMC);
4705 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4709 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4711 if (collGeometry){//get the reaction plane from MC header
4712 gReactionPlane = collGeometry->ReactionPlaneAngle();//[0,180]
4716 //taken from vnv0 code(get the TPC, V0A, V0C event plane using truth tracks)
4717 TClonesArray *mcArray = NULL;
4718 mcArray = (TClonesArray*)event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
4720 Float_t QxMCv2[3] = {0,0,0};
4721 Float_t QyMCv2[3] = {0,0,0};
4722 Float_t QxMCv3[3] = {0,0,0};
4723 Float_t QyMCv3[3] = {0,0,0};
4724 Float_t EvPlaneMCV2[3] = {0,0,0};
4725 Float_t EvPlaneMCV3[3] = {0,0,0};
4726 Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
4727 Float_t etaMax[3] = {4.88,-1.8,0.8};
4729 // analysis on MC tracks
4730 Int_t nMCtrack = mcArray->GetEntries() ;
4732 // EP computation with MC tracks
4733 for(Int_t iT=0;iT < nMCtrack;iT++){
4734 AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
4735 if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
4737 Float_t eta = mctr->Eta();
4738 for(Int_t iD=0;iD<3;iD++){
4739 if(eta > etaMin[iD] && eta < etaMax[iD]){
4740 Float_t phi = mctr->Phi();
4741 QxMCv2[iD] += TMath::Cos(2*phi);
4742 QyMCv2[iD] += TMath::Sin(2*phi);
4743 QxMCv3[iD] += TMath::Cos(3*phi);
4744 QyMCv3[iD] += TMath::Sin(3*phi);
4749 EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
4750 EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
4751 EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
4752 fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
4753 fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
4754 fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
4755 fgPsi2v0aMC = EvPlaneMCV2[0];
4756 fgPsi2v0cMC = EvPlaneMCV2[1];
4757 fgPsi2tpcMC = EvPlaneMCV2[2];
4760 EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
4761 EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
4762 EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
4763 fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
4764 fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
4765 fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
4766 fgPsi3v0aMC = EvPlaneMCV3[0];
4767 fgPsi3v0cMC = EvPlaneMCV3[1];
4768 fgPsi3tpcMC = EvPlaneMCV3[2];
4775 Int_t nAODTracks = event->GetNumberOfTracks();
4777 // TPC EP needed for resolution studies (TPC subevent)
4778 //AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
4779 //Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
4780 Double_t Qx2 = 0, Qy2 = 0;
4781 Double_t Qx3 = 0, Qy3 = 0;
4783 for(Int_t iT = 0; iT < nAODTracks; iT++) {
4785 AliAODTrack* aodTrack = event->GetTrack(iT);
4791 Bool_t trkFlag = aodTrack->TestFilterBit(1);
4793 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
4796 Double_t b[2] = {-99., -99.};
4797 Double_t bCov[3] = {-99., -99., -99.};
4800 AliAODTrack param(*aodTrack);
4801 if (!param.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov)){
4805 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
4808 Qx2 += TMath::Cos(2*aodTrack->Phi());
4809 Qy2 += TMath::Sin(2*aodTrack->Phi());
4810 Qx3 += TMath::Cos(3*aodTrack->Phi());
4811 Qy3 += TMath::Sin(3*aodTrack->Phi());
4815 Float_t evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
4816 Float_t evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
4818 fgPsi2tpc = evPlAng2;
4819 fgPsi3tpc = evPlAng3;
4821 fPhiRPTPC->Fill(iC,evPlAng2);
4822 fPhiRPTPCv3->Fill(iC,evPlAng3);
4827 AliAODVZERO* aodV0 = event->GetVZEROData();
4829 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
4830 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
4831 Float_t multv0 = aodV0->GetMultiplicity(iv0);
4834 if(fAnalysisType == "AOD"){//not for reco MC tracks, only for real data
4835 if (iv0 < 32){ // V0C
4836 Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4837 Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4838 Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4839 Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4841 Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4842 Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4843 Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4844 Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4848 if (iv0 < 32){ // V0C
4849 Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4850 Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4851 Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4852 Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4854 Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4855 Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4856 Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4857 Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4862 //grab for each centrality the proper histo with the Qx and Qy to do the recentering
4863 Double_t Qxamean2 = fMeanQ[iCcal][1][0];
4864 Double_t Qxarms2 = fWidthQ[iCcal][1][0];
4865 Double_t Qyamean2 = fMeanQ[iCcal][1][1];
4866 Double_t Qyarms2 = fWidthQ[iCcal][1][1];
4867 Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
4868 Double_t Qxarms3 = fWidthQv3[iCcal][1][0];
4869 Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
4870 Double_t Qyarms3 = fWidthQv3[iCcal][1][1];
4872 Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
4873 Double_t Qxcrms2 = fWidthQ[iCcal][0][0];
4874 Double_t Qycmean2 = fMeanQ[iCcal][0][1];
4875 Double_t Qycrms2 = fWidthQ[iCcal][0][1];
4876 Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
4877 Double_t Qxcrms3 = fWidthQv3[iCcal][0][0];
4878 Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
4879 Double_t Qycrms3 = fWidthQv3[iCcal][0][1];
4881 Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
4882 Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
4883 Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
4884 Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
4885 Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
4886 Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
4887 Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
4888 Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
4890 //to calculate 2nd order event plane with v0M
4891 Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
4892 /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
4893 Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
4894 /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
4896 //here the calculated event plane is within -Pi to +Pi(delete it , no use here , only for definition)
4897 Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
4898 Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
4899 Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
4903 Float_t evPlAngV0ACor2=999.;
4904 Float_t evPlAngV0CCor2=999.;
4905 Float_t evPlAngV0ACor3=999.;
4906 Float_t evPlAngV0CCor3=999.;
4909 if(fAnalysisType == "AOD"){
4910 evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
4911 evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
4912 evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
4913 evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
4916 evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
4917 evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
4918 evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
4919 evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
4923 AliEventplane *ep = event->GetEventplane();
4924 evPlAngV0ACor2 = ep->GetEventplane("V0A", event, 2);
4925 evPlAngV0CCor2 = ep->GetEventplane("V0C", event, 2);
4926 evPlAngV0ACor3 = ep->GetEventplane("V0A", event, 3);
4927 evPlAngV0CCor3 = ep->GetEventplane("V0C", event, 3);
4930 fgPsi2v0a = evPlAngV0ACor2;
4931 fgPsi2v0c = evPlAngV0CCor2;
4932 fgPsi3v0a = evPlAngV0ACor3;
4933 fgPsi3v0c = evPlAngV0CCor3;
4935 // Fill EP distribution histograms evPlAng
4937 fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
4938 fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
4940 fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
4941 fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
4943 // Fill histograms needed for resolution evaluation
4944 fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
4945 fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
4946 fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
4948 fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
4949 fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
4950 fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
4954 Float_t gVZEROEventPlane = -10.;
4955 Float_t gReactionPlane = -10.;
4956 Double_t qxTot = 0.0, qyTot = 0.0;
4958 AliEventplane *ep = event->GetEventplane();
4960 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
4961 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
4962 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
4963 gReactionPlane = gVZEROEventPlane;
4967 //return gReactionPlane;
4969 //make the final 2nd order event plane within 0 to Pi
4970 //using data and reco tracks only
4971 if(fgPsi2v0a!=999. && fgPsi2v0a < 0.) fgPsi2v0a += TMath::Pi();
4972 if(fgPsi2v0c!=999. && fgPsi2v0c < 0.) fgPsi2v0c += TMath::Pi();
4973 if(fgPsi2tpc!=999. && fgPsi2tpc < 0.) fgPsi2tpc += TMath::Pi();
4974 //using truth tracks only
4975 if(evplaneMC!=999. && evplaneMC < 0.) evplaneMC += TMath::Pi();
4976 if(fgPsi2v0aMC!=999. && fgPsi2v0aMC < 0.) fgPsi2v0aMC += TMath::Pi();
4977 if(fgPsi2v0cMC!=999. && fgPsi2v0cMC < 0.) fgPsi2v0cMC += TMath::Pi();
4978 if(fgPsi2tpcMC!=999. && fgPsi2tpcMC < 0.) fgPsi2tpcMC += TMath::Pi();
4979 //for the time being leave the 3rd order event planes within -pi/3 t0 +pi/3
4981 if(truth){//for truth MC
4982 if(fV2 && fEPdet=="header")gReactionPlane=evplaneMC;
4983 if(fV2 && fEPdet=="V0A")gReactionPlane=fgPsi2v0aMC;
4984 if(fV2 && fEPdet=="V0C")gReactionPlane=fgPsi2v0cMC;
4985 if(fV2 && fEPdet=="TPC")gReactionPlane=fgPsi2tpcMC;
4987 if(fV3 && fEPdet=="V0A")gReactionPlane=fgPsi3v0aMC;
4988 if(fV3 && fEPdet=="V0C")gReactionPlane=fgPsi3v0cMC;
4989 if(fV3 && fEPdet=="TPC")gReactionPlane=fgPsi3tpcMC;
4991 else{//for data and recoMC
4992 if(fV2 && fEPdet=="V0A")gReactionPlane=fgPsi2v0a;
4993 if(fV2 && fEPdet=="V0C")gReactionPlane=fgPsi2v0c;
4994 if(fV2 && fEPdet=="TPC")gReactionPlane=fgPsi2tpc;
4996 if(fV3 && fEPdet=="V0A")gReactionPlane=fgPsi3v0a;
4997 if(fV3 && fEPdet=="V0C")gReactionPlane=fgPsi3v0c;
4998 if(fV3 && fEPdet=="TPC")gReactionPlane=fgPsi3tpc;
5002 } //centrality cut condition
5004 return gReactionPlane;
5006 //------------------------------------------------------------------------------------------------------------------
5007 void AliTwoParticlePIDCorr::OpenInfoCalbration(Int_t run){
5008 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
5009 TFile *foadb = TFile::Open(oadbfilename.Data());
5012 printf("OADB file %s cannot be opened\n",oadbfilename.Data());
5016 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
5018 printf("OADB object hMultV0BefCorr is not available in the file\n");
5022 if(!(cont->GetObject(run))){
5023 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
5026 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
5028 TF1 *fpol0 = new TF1("fpol0","pol0");
5029 fMultV0->Fit(fpol0,"","",0,31);
5030 fV0Cpol = fpol0->GetParameter(0);
5031 fMultV0->Fit(fpol0,"","",32,64);
5032 fV0Apol = fpol0->GetParameter(0);
5034 for(Int_t iside=0;iside<2;iside++){
5035 for(Int_t icoord=0;icoord<2;icoord++){
5036 for(Int_t i=0;i < 9;i++){
5038 if(iside==0 && icoord==0)
5039 snprintf(namecont,100,"hQxc2_%i",i);
5040 else if(iside==1 && icoord==0)
5041 snprintf(namecont,100,"hQxa2_%i",i);
5042 else if(iside==0 && icoord==1)
5043 snprintf(namecont,100,"hQyc2_%i",i);
5044 else if(iside==1 && icoord==1)
5045 snprintf(namecont,100,"hQya2_%i",i);
5047 cont = (AliOADBContainer*) foadb->Get(namecont);
5049 printf("OADB object %s is not available in the file\n",namecont);
5053 if(!(cont->GetObject(run))){
5054 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5057 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5058 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5061 if(iside==0 && icoord==0)
5062 snprintf(namecont,100,"hQxc3_%i",i);
5063 else if(iside==1 && icoord==0)
5064 snprintf(namecont,100,"hQxa3_%i",i);
5065 else if(iside==0 && icoord==1)
5066 snprintf(namecont,100,"hQyc3_%i",i);
5067 else if(iside==1 && icoord==1)
5068 snprintf(namecont,100,"hQya3_%i",i);
5070 cont = (AliOADBContainer*) foadb->Get(namecont);
5072 printf("OADB object %s is not available in the file\n",namecont);
5076 if(!(cont->GetObject(run))){
5077 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
5080 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
5081 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
5087 //____________________________________________________________________
5088 void AliTwoParticlePIDCorr::FillPIDEventPlane(Double_t centrality,Int_t par,Float_t trigphi,Float_t fReactionPlane)
5091 // Event plane (determine psi bin)
5092 Double_t gPsiMinusPhi = 0.;
5093 Double_t gPsiMinusPhiBin = -10.;
5094 if(fRequestEventPlane){
5095 gPsiMinusPhi = TMath::Abs(trigphi - fReactionPlane);
5097 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
5098 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
5099 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
5100 gPsiMinusPhiBin = 0.0;
5102 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
5103 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
5104 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
5105 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
5106 gPsiMinusPhiBin = 1.0;
5108 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
5109 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
5110 gPsiMinusPhiBin = 2.0;
5113 gPsiMinusPhiBin = 3.0;
5115 fEventPlanePID->Fill(centrality,gPsiMinusPhiBin,(Float_t)par);
5119 //----------------------------------------------------------
5120 void AliTwoParticlePIDCorr::Terminate(Option_t *)
5122 // Draw result to screen, or perform fitting, normalizations
5123 // Called once at the end of the query
5124 fOutput = dynamic_cast<TList*> (GetOutputData(1));
5125 if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
5129 //------------------------------------------------------------------