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"),
95 fRequestEventPlane(kFALSE),
96 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
101 fSharedClusterCut(-1),
103 skipParticlesAbove(0),
105 ffilltrigassoUNID(kFALSE),
106 ffilltrigUNIDassoID(kFALSE),
107 ffilltrigIDassoUNID(kTRUE),
108 ffilltrigIDassoID(kFALSE),
109 ffilltrigIDassoIDMCTRUTH(kFALSE),
110 fMaxNofMixingTracks(50000),
111 fPtOrderMCTruth(kTRUE),
112 fPtOrderDataReco(kTRUE),
113 fWeightPerEvent(kFALSE),
114 fTriggerSpeciesSelection(kFALSE),
115 fAssociatedSpeciesSelection(kFALSE),
116 fRandomizeReactionPlane(kFALSE),
117 fTriggerSpecies(SpPion),
118 fAssociatedSpecies(SpPion),
121 fSelectHighestPtTrig(kFALSE),
122 fcontainPIDtrig(kTRUE),
123 fcontainPIDasso(kFALSE),
125 frejectPileUp(kFALSE),
130 fselectprimaryTruth(kTRUE),
131 fonlyprimarydatareco(kFALSE),
134 ffillhistQAReco(kFALSE),
135 ffillhistQATruth(kFALSE),
136 kTrackVariablesPair(0),
168 fhistJetTrigestimate(0),
169 fTwoTrackDistancePtdip(0x0),
170 fTwoTrackDistancePtdipmix(0x0),
171 fCentralityCorrelation(0x0),
172 fHistVZEROAGainEqualizationMap(0),
173 fHistVZEROCGainEqualizationMap(0),
174 fHistVZEROChannelGainEqualizationMap(0),
175 fCentralityWeights(0),
178 fHistEQVZEROvsTPCmultiplicity(0x0),
179 fHistEQVZEROAvsTPCmultiplicity(0x0),
180 fHistEQVZEROCvsTPCmultiplicity(0x0),
181 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
182 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
183 fHistVZEROCvsVZEROAmultiplicity(0x0),
184 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
185 fHistVZEROSignal(0x0),
186 fHistEventPlaneTruth(0x0),
187 fHistPsiMinusPhi(0x0),
202 gReactionPlane(999.),
230 fControlConvResoncances(0),
245 fCorrelatonTruthPrimary(0),
246 fCorrelatonTruthPrimarymix(0),
252 fTHnCorrIDUNIDmix(0),
254 fTHnTrigcountMCTruthPrim(0),
257 fAnalysisType("AOD"),
259 ftwoTrackEfficiencyCutDataReco(kTRUE),
260 twoTrackEfficiencyCutValue(0.02),
266 fRequestTOFPID(kTRUE),
267 fPIDType(NSigmaTPCTOF),
268 fFIllPIDQAHistos(kTRUE),
271 fdiffPIDcutvalues(kFALSE),
276 fHighPtKaonNSigmaPID(-1),
277 fHighPtKaonSigma(3.5),
278 fUseExclusiveNSigma(kFALSE),
279 fRemoveTracksT0Fill(kFALSE),
281 fTriggerSelectCharge(0),
282 fAssociatedSelectCharge(0),
283 fTriggerRestrictEta(-1),
284 fEtaOrdering(kFALSE),
285 fCutConversions(kFALSE),
286 fCutResonances(kFALSE),
287 fRejectResonanceDaughters(-1),
289 fInjectedSignals(kFALSE),
290 fRemoveWeakDecays(kFALSE),
291 fRemoveDuplicates(kFALSE),
292 fapplyTrigefficiency(kFALSE),
293 fapplyAssoefficiency(kFALSE),
294 ffillefficiency(kFALSE),
295 fmesoneffrequired(kFALSE),
296 fkaonprotoneffrequired(kFALSE),
301 for ( Int_t i = 0; i < 16; i++) {
305 for ( Int_t i = 0; i < 6; i++ ){
306 fTrackHistEfficiency[i] = NULL;
307 effcorection[i]=NULL;
310 for ( Int_t i = 0; i < 2; i++ ){
311 fTwoTrackDistancePt[i]=NULL;
312 fTwoTrackDistancePtmix[i]=NULL;
315 for(Int_t ipart=0;ipart<NSpecies;ipart++)
316 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
317 fnsigmas[ipart][ipid]=999.;
319 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
321 for(Int_t i = 0; i != 2; ++i)
322 for(Int_t j = 0; j != 2; ++j)
323 for(Int_t iC = 0; iC < 9; iC++){
324 fMeanQ[iC][i][j] = 0.;
325 fWidthQ[iC][i][j] = 1.;
326 fMeanQv3[iC][i][j] = 0.;
327 fWidthQv3[iC][i][j] = 1.;
331 //________________________________________________________________________
332 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
333 :AliAnalysisTaskSE(name),
337 fCentralityMethod("V0A"),
339 fRequestEventPlane(kFALSE),
340 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
345 fSharedClusterCut(-1),
347 skipParticlesAbove(0),
349 ffilltrigassoUNID(kFALSE),
350 ffilltrigUNIDassoID(kFALSE),
351 ffilltrigIDassoUNID(kTRUE),
352 ffilltrigIDassoID(kFALSE),
353 ffilltrigIDassoIDMCTRUTH(kFALSE),
354 fMaxNofMixingTracks(50000),
355 fPtOrderMCTruth(kTRUE),
356 fPtOrderDataReco(kTRUE),
357 fWeightPerEvent(kFALSE),
358 fTriggerSpeciesSelection(kFALSE),
359 fAssociatedSpeciesSelection(kFALSE),
360 fRandomizeReactionPlane(kFALSE),
361 fTriggerSpecies(SpPion),
362 fAssociatedSpecies(SpPion),
365 fSelectHighestPtTrig(kFALSE),
366 fcontainPIDtrig(kTRUE),
367 fcontainPIDasso(kFALSE),
369 frejectPileUp(kFALSE),
374 fselectprimaryTruth(kTRUE),
375 fonlyprimarydatareco(kFALSE),
378 ffillhistQAReco(kFALSE),
379 ffillhistQATruth(kFALSE),
380 kTrackVariablesPair(0),
412 fhistJetTrigestimate(0),
413 fTwoTrackDistancePtdip(0x0),
414 fTwoTrackDistancePtdipmix(0x0),
415 fCentralityCorrelation(0x0),
416 fHistVZEROAGainEqualizationMap(0),
417 fHistVZEROCGainEqualizationMap(0),
418 fHistVZEROChannelGainEqualizationMap(0),
419 fCentralityWeights(0),
422 fHistEQVZEROvsTPCmultiplicity(0x0),
423 fHistEQVZEROAvsTPCmultiplicity(0x0),
424 fHistEQVZEROCvsTPCmultiplicity(0x0),
425 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
426 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
427 fHistVZEROCvsVZEROAmultiplicity(0x0),
428 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
429 fHistVZEROSignal(0x0),
430 fHistEventPlaneTruth(0x0),
431 fHistPsiMinusPhi(0x0),
446 gReactionPlane(999.),
474 fControlConvResoncances(0),
489 fCorrelatonTruthPrimary(0),
490 fCorrelatonTruthPrimarymix(0),
496 fTHnCorrIDUNIDmix(0),
498 fTHnTrigcountMCTruthPrim(0),
501 fAnalysisType("AOD"),
503 ftwoTrackEfficiencyCutDataReco(kTRUE),
504 twoTrackEfficiencyCutValue(0.02),
510 fRequestTOFPID(kTRUE),
511 fPIDType(NSigmaTPCTOF),
512 fFIllPIDQAHistos(kTRUE),
515 fdiffPIDcutvalues(kFALSE),
520 fHighPtKaonNSigmaPID(-1),
521 fHighPtKaonSigma(3.5),
522 fUseExclusiveNSigma(kFALSE),
523 fRemoveTracksT0Fill(kFALSE),
525 fTriggerSelectCharge(0),
526 fAssociatedSelectCharge(0),
527 fTriggerRestrictEta(-1),
528 fEtaOrdering(kFALSE),
529 fCutConversions(kFALSE),
530 fCutResonances(kFALSE),
531 fRejectResonanceDaughters(-1),
533 fInjectedSignals(kFALSE),
534 fRemoveWeakDecays(kFALSE),
535 fRemoveDuplicates(kFALSE),
536 fapplyTrigefficiency(kFALSE),
537 fapplyAssoefficiency(kFALSE),
538 ffillefficiency(kFALSE),
539 fmesoneffrequired(kFALSE),
540 fkaonprotoneffrequired(kFALSE),
543 // The last in the above list should not have a comma after it
547 for ( Int_t i = 0; i < 16; i++) {
551 for ( Int_t i = 0; i < 6; i++ ){
552 fTrackHistEfficiency[i] = NULL;
553 effcorection[i]=NULL;
557 for ( Int_t i = 0; i < 2; i++ ){
558 fTwoTrackDistancePt[i]=NULL;
559 fTwoTrackDistancePtmix[i]=NULL;
562 for(Int_t ipart=0;ipart<NSpecies;ipart++)
563 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
564 fnsigmas[ipart][ipid]=999.;
566 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
568 for(Int_t i = 0; i != 2; ++i)
569 for(Int_t j = 0; j != 2; ++j)
570 for(Int_t iC = 0; iC < 9; iC++){
571 fMeanQ[iC][i][j] = 0.;
572 fWidthQ[iC][i][j] = 1.;
573 fMeanQv3[iC][i][j] = 0.;
574 fWidthQv3[iC][i][j] = 1.;
578 // Define input and output slots here (never in the dummy constructor)
579 // Input slot #0 works with a TChain - it is connected to the default input container
580 // Output slot #1 writes into a TH1 container
582 DefineOutput(1, TList::Class()); // for output list
583 DefineOutput(2, TList::Class());
587 //________________________________________________________________________
588 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
590 // Destructor. Clean-up the output list, but not the histograms that are put inside
591 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
592 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
597 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
602 if (fPID) delete fPID;
603 if (fPIDCombined) delete fPIDCombined;
606 //________________________________________________________________________
608 //////////////////////////////////////////////////////////////////////////////////////////////////
610 TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
611 // returns histo named name
612 return (TH2F*) fOutputList->FindObject(name);
615 //////////////////////////////////////////////////////////////////////////////////////////////////
617 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
621 // Puts the argument in the range [-pi/2,3 pi/2].
624 if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
625 if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
630 //________________________________________________________________________
631 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
634 // Called once (on the worker node)
635 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
636 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
637 fPID = inputHandler->GetPIDResponse();
639 //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
641 //get the efficiency correction map
643 // global switch disabling the reference
644 // (to avoid "Replacing existing TH1" if several wagons are created in train)
645 Bool_t oldStatus = TH1::AddDirectoryStatus();
646 TH1::AddDirectory(kFALSE);
648 const Int_t nPsiTOF = 10;
649 const Int_t nCentrBin = 9;
652 fOutput = new TList();
653 fOutput->SetOwner(); // IMPORTANT!
655 fOutputList = new TList;
656 fOutputList->SetOwner();
657 fOutputList->SetName("PIDQAList");
661 fList->SetName("EPQAList");
663 fEventCounter = new TH1F("fEventCounter","EventCounter", 19, 0.5,19.5);
664 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
665 fEventCounter->GetXaxis()->SetBinLabel(3,"After PileUP Cut");//only for Data
666 fEventCounter->GetXaxis()->SetBinLabel(5,"Have A Vertex");
667 fEventCounter->GetXaxis()->SetBinLabel(7,"After vertex Cut");
668 fEventCounter->GetXaxis()->SetBinLabel(9,"Getting centrality");
669 fEventCounter->GetXaxis()->SetBinLabel(11,"After centrality flattening");
670 fEventCounter->GetXaxis()->SetBinLabel(13,"Within 0-100% centrality");
671 fEventCounter->GetXaxis()->SetBinLabel(15,"Event Analyzed");
672 //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
673 fOutput->Add(fEventCounter);
675 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
676 fOutput->Add(fEtaSpectrasso);
678 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
679 fOutput->Add(fphiSpectraasso);
681 if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
682 fOutput->Add(fCentralityCorrelation);
685 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")
687 TString gCentName[8] = {"V0A","V0C","V0M","V0AEq","V0CEq","V0MEq","CL1","ZNA"};
688 fHistCentStats = new TH2F("fHistCentStats",
689 "Centrality statistics;;Cent percentile",
690 8,-0.5,7.5,220,-5,105);
691 for(Int_t i = 1; i <= 8; i++){
692 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
693 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
695 fOutput->Add(fHistCentStats);
698 if(fCentralityMethod.EndsWith("_MANUAL"))
700 fhistcentrality=new TH1F("fhistcentrality","referencemultiplicity",30001,-0.5,30000.5);
701 fOutput->Add(fhistcentrality);
704 fhistcentrality=new TH1F("fhistcentrality","centrality",220,-5,105);
705 fOutput->Add(fhistcentrality);
708 if(fCentralityMethod.EndsWith("_MANUAL"))
710 TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
711 fHistRefmult = new TH2F("fHistRefmult",
712 "Reference multiplicity",
713 4,-0.5,3.5,10000,0,20000);
714 for(Int_t i = 1; i <= 4; i++){
715 fHistRefmult->GetXaxis()->SetBinLabel(i,gmultName[i-1].Data());
716 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
718 fOutput->Add(fHistRefmult);
720 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
721 //TPC vs EQVZERO multiplicity
722 fHistEQVZEROvsTPCmultiplicity = new TH2F("fHistEQVZEROvsTPCmultiplicity","EqVZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
723 fHistEQVZEROvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO multiplicity (a.u.)");
724 fHistEQVZEROvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
725 fOutput->Add(fHistEQVZEROvsTPCmultiplicity);
728 fHistEQVZEROAvsTPCmultiplicity = new TH2F("fHistEQVZEROAvsTPCmultiplicity","EqVZERO_A vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
729 fHistEQVZEROAvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
730 fHistEQVZEROAvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
731 fOutput->Add(fHistEQVZEROAvsTPCmultiplicity);
734 fHistEQVZEROCvsTPCmultiplicity = new TH2F("fHistEQVZEROCvsTPCmultiplicity","EqVZERO_C vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
735 fHistEQVZEROCvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
736 fHistEQVZEROCvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
737 fOutput->Add(fHistEQVZEROCvsTPCmultiplicity);
739 //EQVZERO vs VZERO multiplicity
740 fHistVZEROCvsEQVZEROCmultiplicity = new TH2F("fHistVZEROCvsEQVZEROCmultiplicity","EqVZERO_C vs VZERO_C multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
741 fHistVZEROCvsEQVZEROCmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
742 fHistVZEROCvsEQVZEROCmultiplicity->GetYaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
743 fOutput->Add(fHistVZEROCvsEQVZEROCmultiplicity);
746 fHistVZEROAvsEQVZEROAmultiplicity = new TH2F("fHistVZEROAvsEQVZEROAmultiplicity","EqVZERO_A vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
747 fHistVZEROAvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
748 fHistVZEROAvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
749 fOutput->Add(fHistVZEROAvsEQVZEROAmultiplicity);
752 //VZEROC vs VZEROA multiplicity
753 fHistVZEROCvsVZEROAmultiplicity = new TH2F("fHistVZEROCvsVZEROAmultiplicity","VZERO_C vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
754 fHistVZEROCvsVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
755 fHistVZEROCvsVZEROAmultiplicity->GetYaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
756 fOutput->Add(fHistVZEROCvsVZEROAmultiplicity);
760 //EQVZEROC vs EQVZEROA multiplicity
761 fHistEQVZEROCvsEQVZEROAmultiplicity = new TH2F("fHistEQVZEROCvsEQVZEROAmultiplicity","EqVZERO_C vs EqVZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
762 fHistEQVZEROCvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
763 fHistEQVZEROCvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
764 fOutput->Add(fHistEQVZEROCvsEQVZEROAmultiplicity);
766 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
767 fOutput->Add(fHistVZEROSignal);
771 if(fRequestEventPlane){
774 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
775 fList->Add(fHistPsiMinusPhi);
777 fEventPlanePID = new TH3F("fEventPlanePID",";centrality;eventplane;PID",20,0.0,100.0,4,-0.5,3.5,4,-0.5,3.5);
778 fList->Add(fEventPlanePID);
782 if(fCutConversions || fCutResonances)
784 fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
785 fOutput->Add(fControlConvResoncances);
788 fHistoTPCdEdx = new TH2F("fHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
789 fOutputList->Add(fHistoTPCdEdx);
790 fHistoTOFbeta = new TH2F(Form("fHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
791 fOutputList->Add(fHistoTOFbeta);
793 fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
794 fOutputList->Add(fTPCTOFPion3d);
796 fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
797 fOutputList->Add(fTPCTOFKaon3d);
799 fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
800 fOutputList->Add(fTPCTOFProton3d);
804 fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);
805 fOutputList->Add(fPionPt);
806 fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);
807 fOutputList->Add(fPionEta);
808 fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);
809 fOutputList->Add(fPionPhi);
811 fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);
812 fOutputList->Add(fKaonPt);
813 fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);
814 fOutputList->Add(fKaonEta);
815 fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);
816 fOutputList->Add(fKaonPhi);
818 fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);
819 fOutputList->Add(fProtonPt);
820 fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);
821 fOutputList->Add(fProtonEta);
822 fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);
823 fOutputList->Add(fProtonPhi);
826 fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
827 fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
828 fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
829 fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);
830 fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
831 fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
832 fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
833 fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
834 fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
835 fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
836 fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
837 fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
838 fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
839 fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
840 fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
841 fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
843 for(Int_t i = 0; i < 16; i++)
845 fOutput->Add(fHistQA[i]);
848 Int_t eventplaneaxis=0;
850 if (fRequestEventPlane) eventplaneaxis=1;
852 kTrackVariablesPair=6+SetChargeAxis+eventplaneaxis;
854 if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
856 if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
858 if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8+SetChargeAxis+eventplaneaxis;
861 // two particle histograms
862 Int_t anaSteps = 1; // analysis steps
863 const char* title = "d^{2}N_{ch}/d#varphid#eta";
865 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
866 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
867 TString* axisTitlePair; // axis titles for track variables
868 axisTitlePair=new TString[kTrackVariablesPair];
870 TString defaultBinningStr;
871 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"
872 "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"
873 "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
874 "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"
875 "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
876 "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
877 "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"
878 "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
880 if(fRequestEventPlane){
881 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))
885 defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
888 defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
891 if(SetChargeAxis==2){
892 defaultBinningStr += "TrigCharge: -2.0,0.0,2.0\n"; // course
893 defaultBinningStr += "AssoCharge: -2.0,0.0,2.0\n"; // course
895 // =========================================================
896 // Customization (adopted from AliUEHistograms)
897 // =========================================================
899 TObjArray* lines = defaultBinningStr.Tokenize("\n");
900 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
902 TString line(lines->At(i)->GetName());
903 TString tag = line(0, line.Index(":")+1);
904 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
905 fBinningString += line + "\n";
907 AliInfo(Form("Using custom binning for %s", tag.Data()));
910 fBinningString += fCustomBinning;
912 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
914 // =========================================================
916 // =========================================================
918 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
919 axisTitlePair[0] = " multiplicity";
921 dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
922 axisTitlePair[1] = "v_{Z} (cm)";
924 dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
925 axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
927 dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
928 axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
930 dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
931 axisTitlePair[4] = "#Delta#eta";
933 dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
934 axisTitlePair[5] = "#Delta#varphi (rad)";
938 if(fRequestEventPlane){
939 dBinsPair[dim_val] = GetBinning(fBinningString, "eventPlane", iBinPair[dim_val]);
940 axisTitlePair[dim_val] = "#varphi - #Psi_{2} (a.u.)";
944 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
945 dBinsPair[dim_val] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val]);
946 axisTitlePair[dim_val] = "TrigCharge";
948 dBinsPair[dim_val+1] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+1]);
949 axisTitlePair[dim_val+1] = "AssoCharge";
952 if(fcontainPIDtrig && !fcontainPIDasso){
953 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
954 axisTitlePair[dim_val] = "PIDTrig";
955 if(SetChargeAxis==2){
956 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
957 axisTitlePair[dim_val+1] = "TrigCharge";
959 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
960 axisTitlePair[dim_val+2] = "AssoCharge";
964 if(!fcontainPIDtrig && fcontainPIDasso){
965 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val]);
966 axisTitlePair[dim_val] = "PIDAsso";
968 if(SetChargeAxis==2){
969 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
970 axisTitlePair[dim_val+1] = "TrigCharge";
972 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
973 axisTitlePair[dim_val+2] = "AssoCharge";
977 if(fcontainPIDtrig && fcontainPIDasso){
979 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
980 axisTitlePair[dim_val] = "PIDTrig";
982 dBinsPair[dim_val+1] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val+1]);
983 axisTitlePair[dim_val+1] = "PIDAsso";
985 if(SetChargeAxis==2){
986 dBinsPair[dim_val+2] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+2]);
987 axisTitlePair[dim_val+2] = "TrigCharge";
989 dBinsPair[dim_val+3] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+3]);
990 axisTitlePair[dim_val+3] = "AssoCharge";
995 Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
997 Int_t nPteffbin = -1;
998 Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
1001 fminPtTrig=dBinsPair[2][0];
1002 fmaxPtTrig=dBinsPair[2][iBinPair[2]];
1003 fminPtAsso=dBinsPair[3][0];
1004 fmaxPtAsso=dBinsPair[3][iBinPair[3]];
1005 fmincentmult=dBinsPair[0][0];
1006 fmaxcentmult=dBinsPair[0][iBinPair[0]];
1008 //event pool manager
1009 Int_t MaxNofEvents=1000;
1010 const Int_t NofVrtxBins=10+(1+10)*2;
1011 Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
1012 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
1013 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210};
1015 if(fCentralityMethod.EndsWith("_MANUAL"))
1017 const Int_t NofCentBins=10;
1018 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....
1019 if(fRequestEventPlane){
1020 // Event plane angle (Psi) bins
1022 Double_t* psibins = NULL;
1024 psibins = GetBinning(fBinningString, "eventPlane", nPsiBins);
1026 const Int_t nPsiBins=6;
1027 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()};
1028 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1029 // if(psibins) delete [] psibins;
1033 const Int_t nPsiBinsd=1;
1034 Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1035 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1037 // fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
1039 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1044 const Int_t NofCentBins=15;
1045 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
1046 if(fRequestEventPlane){
1047 // Event plane angle (Psi) bins
1049 Double_t* psibins = NULL;
1051 psibins = GetBinning(fBinningString, "eventPlane", nPsiBins);
1053 const Int_t nPsiBins=6;
1054 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()};
1055 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1056 // if(psibins) delete [] psibins;
1060 const Int_t nPsiBinsd=1;
1061 Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1062 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1064 //fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
1066 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1071 AliError("Event Mixing required, but Pool Manager not initialized...");
1075 //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
1076 //fmaxPtComboeff=fmaxPtTrig;
1077 //THnSparses for calculation of efficiency
1079 if((fAnalysisType =="MCAOD") && ffillefficiency) {
1082 effbin[0]=iBinPair[0];
1083 effbin[1]=iBinPair[1];
1084 effbin[2]=nPteffbin;
1086 Int_t effsteps=5;//for each species type::primMCParticles(0),primRecoTracksMatched(1),allRecoTracksMatched(2),primRecoTracksMatchedPID(3),allRecoTracksMatchedPID(4)
1087 for(Int_t jj=0;jj<6;jj++)//PID type binning
1089 if(jj==5) effsteps=3;//for unidentified particles
1090 Histrename="fTrackHistEfficiency";Histrename+=jj;
1091 fTrackHistEfficiency[jj] = new AliTHn(Histrename.Data(), "Tracking efficiency", effsteps, 4, effbin);
1092 fTrackHistEfficiency[jj]->SetBinLimits(0, dBinsPair[0]);
1093 fTrackHistEfficiency[jj]->SetVarTitle(0, "Centrality");
1094 fTrackHistEfficiency[jj]->SetBinLimits(1, dBinsPair[1]);
1095 fTrackHistEfficiency[jj]->SetVarTitle(1, "zvtx");
1096 fTrackHistEfficiency[jj]->SetBinLimits(2, Pteff);
1097 fTrackHistEfficiency[jj]->SetVarTitle(2, "p_{T} (GeV/c)");
1098 fTrackHistEfficiency[jj]->SetBinLimits(3, EtaBin);
1099 fTrackHistEfficiency[jj]->SetVarTitle(3, "#eta");
1100 fOutput->Add(fTrackHistEfficiency[jj]);
1104 //AliThns for Correlation plots(data & MC)
1106 if(ffilltrigassoUNID)
1108 fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1109 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1110 fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
1111 fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
1113 fOutput->Add(fTHnCorrUNID);
1115 fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1116 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1117 fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
1118 fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1120 fOutput->Add(fTHnCorrUNIDmix);
1123 if(ffilltrigIDassoID)
1125 fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
1126 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1127 fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
1128 fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
1130 fOutput->Add(fTHnCorrID);
1132 fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1133 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1134 fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
1135 fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
1137 fOutput->Add(fTHnCorrIDmix);
1140 if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
1142 fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1143 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1144 fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
1145 fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
1147 fOutput->Add(fTHnCorrIDUNID);
1150 fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1151 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1152 fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
1153 fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1155 fOutput->Add(fTHnCorrIDUNIDmix);
1160 //ThnSparse for Correlation plots(truth MC)
1161 if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
1163 fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
1164 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1165 fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
1166 fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
1168 fOutput->Add(fCorrelatonTruthPrimary);
1171 fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
1172 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1173 fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
1174 fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
1176 fOutput->Add(fCorrelatonTruthPrimarymix);
1179 //binning for trigger no. counting
1182 if(SetChargeAxis==2) ChargeAxis=1;
1185 Int_t dims=3+ChargeAxis+eventplaneaxis;
1186 if(fcontainPIDtrig) dims=4+ChargeAxis+eventplaneaxis;
1187 fBinst= new Int_t[dims];
1188 Double_t* dBinsTrig[dims]; // bins for track variables
1189 TString* axisTitleTrig; // axis titles for track variables
1190 axisTitleTrig=new TString[dims];
1192 for(Int_t i=0; i<3;i++)
1194 fBinst[i]=iBinPair[i];
1195 dBinsTrig[i]=dBinsPair[i];
1196 axisTitleTrig[i]=axisTitlePair[i];
1198 Int_t dim_val_trig=3;
1199 if(fRequestEventPlane){
1200 fBinst[dim_val_trig]=iBinPair[6];//if fRequestEventPlane=TRUE, dim_val already becomes 7.
1201 dBinsTrig[dim_val_trig]=dBinsPair[6];
1202 axisTitleTrig[dim_val_trig]=axisTitlePair[6];
1206 if(!fcontainPIDtrig && !fcontainPIDasso && ChargeAxis==1){
1207 fBinst[dim_val_trig]=iBinPair[dim_val];
1208 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1209 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1212 if(fcontainPIDtrig && !fcontainPIDasso){
1213 fBinst[dim_val_trig]=iBinPair[dim_val];
1214 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1215 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1217 fBinst[dim_val_trig+1]=iBinPair[dim_val+1];
1218 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+1];
1219 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+1];
1223 if(!fcontainPIDtrig && fcontainPIDasso){
1225 fBinst[dim_val_trig]=iBinPair[dim_val+1];
1226 dBinsTrig[dim_val_trig]=dBinsPair[dim_val+1];
1227 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val+1];
1231 if(fcontainPIDtrig && fcontainPIDasso){
1232 fBinst[dim_val_trig]=iBinPair[dim_val];
1233 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1234 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1236 fBinst[dim_val_trig+1]=iBinPair[dim_val+2];
1237 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+2];
1238 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+2];
1242 //ThSparse for trigger counting(data & reco MC)
1243 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
1245 fTHnTrigcount = new AliTHn("fTHnTrigcount", "fTHnTrigcount", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1246 for(Int_t i=0; i<dims;i++){
1247 fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1248 fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1250 fOutput->Add(fTHnTrigcount);
1253 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
1254 //AliTHns for trigger counting(truth MC)
1255 fTHnTrigcountMCTruthPrim = new AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1256 for(Int_t i=0; i<dims;i++){
1257 fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1258 fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1260 fOutput->Add(fTHnTrigcountMCTruthPrim);
1263 if(fAnalysisType=="MCAOD"){
1264 if(ffillhistQATruth)
1266 MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
1267 fOutputList->Add(MCtruthpt);
1269 MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
1270 fOutputList->Add(MCtrutheta);
1272 MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
1273 fOutputList->Add(MCtruthphi);
1275 MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
1276 fOutputList->Add(MCtruthpionpt);
1278 MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
1279 fOutputList->Add(MCtruthpioneta);
1281 MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
1282 fOutputList->Add(MCtruthpionphi);
1284 MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
1285 fOutputList->Add(MCtruthkaonpt);
1287 MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
1288 fOutputList->Add(MCtruthkaoneta);
1290 MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
1291 fOutputList->Add(MCtruthkaonphi);
1293 MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
1294 fOutputList->Add(MCtruthprotonpt);
1296 MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
1297 fOutputList->Add(MCtruthprotoneta);
1299 MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
1300 fOutputList->Add(MCtruthprotonphi);
1302 fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
1303 fOutputList->Add(fPioncont);
1305 fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
1306 fOutputList->Add(fKaoncont);
1308 fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
1309 fOutputList->Add(fProtoncont);
1311 fUNIDcont=new TH2F("fUNIDcont","fUNIDcont",10,-0.5,9.5,100,0.,10.);
1312 fOutputList->Add(fUNIDcont);
1315 fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1316 fEventno->GetXaxis()->SetTitle("Centrality");
1317 fEventno->GetYaxis()->SetTitle("Z_Vtx");
1318 fOutput->Add(fEventno);
1319 fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1320 fEventnobaryon->GetXaxis()->SetTitle("Centrality");
1321 fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
1322 fOutput->Add(fEventnobaryon);
1323 fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1324 fEventnomeson->GetXaxis()->SetTitle("Centrality");
1325 fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
1326 fOutput->Add(fEventnomeson);
1328 fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
1329 fOutput->Add(fhistJetTrigestimate);
1331 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);
1332 fOutput->Add(fTwoTrackDistancePtdip);
1334 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);
1335 fOutput->Add(fTwoTrackDistancePtdipmix);
1337 TString Histttrname;
1338 for(Int_t jj=0;jj<2;jj++)// PID type binning
1340 Histttrname="fTwoTrackDistancePt";Histttrname+=jj;
1341 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);
1342 fOutput->Add(fTwoTrackDistancePt[jj]);
1344 Histttrname="fTwoTrackDistancePtmix";Histttrname+=jj;
1345 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);
1346 fOutput->Add(fTwoTrackDistancePtmix[jj]);
1349 //DefineEventPool();
1351 if(fapplyTrigefficiency || fapplyAssoefficiency)
1353 const Int_t nDimt = 4;// cent zvtx pt eta
1354 Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
1355 Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
1356 Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
1358 TString Histrexname;
1359 for(Int_t jj=0;jj<6;jj++)// PID type binning
1361 Histrexname="effcorection";Histrexname+=jj;
1362 effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
1363 effcorection[jj]->Sumw2();
1364 effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
1365 effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
1366 effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
1367 effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
1368 effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
1369 effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1370 effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
1371 effcorection[jj]->GetAxis(3)->SetTitle("#eta");
1372 fOutput->Add(effcorection[jj]);
1374 // TFile *fsifile = new TFile(fefffilename,"READ");
1376 if (TString(fefffilename).BeginsWith("alien:"))
1377 TGrid::Connect("alien:");
1378 TFile *fileT=TFile::Open(fefffilename);
1380 for(Int_t jj=0;jj<6;jj++)//type binning
1382 Nameg="effmap";Nameg+=jj;
1383 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
1384 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
1386 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
1393 //*************************************************************EP plots***********************************************//
1394 if(fRequestEventPlane){
1395 // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
1397 fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
1398 fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
1399 fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
1401 fList->Add(fHResTPCv0A2);
1402 fList->Add(fHResTPCv0C2);
1403 fList->Add(fHResv0Cv0A2);
1406 fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
1407 fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
1408 fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
1410 fList->Add(fHResTPCv0A3);
1411 fList->Add(fHResTPCv0C3);
1412 fList->Add(fHResv0Cv0A3);
1414 // MC as in the dataEP resolution (but using MC tracks)
1415 if(fAnalysisType == "MCAOD" && fV2){
1416 fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
1417 fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
1418 fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
1419 fList->Add(fHResMA2);
1420 fList->Add(fHResMC2);
1421 fList->Add(fHResAC2);
1423 if(fAnalysisType == "MCAOD" && fV3){
1424 fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
1425 fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
1426 fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
1427 fList->Add(fHResMA3);
1428 fList->Add(fHResMC3);
1429 fList->Add(fHResAC3);
1433 // V0A and V0C event plane distributions
1435 fPhiRPTPC = new TH2F("fPhiRPTPCv2","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1436 fPhiRPTPCv3 = new TH2F("fPhiRPTPCv3","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1437 fList->Add(fPhiRPTPC);
1438 fList->Add(fPhiRPTPCv3);
1440 fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1441 fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1442 fList->Add(fPhiRPv0A);
1443 fList->Add(fPhiRPv0C);
1446 fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1447 fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1448 fList->Add(fPhiRPv0Av3);
1449 fList->Add(fPhiRPv0Cv3);
1451 fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth","#phi distribution of EP MCTRUTHheader;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1452 fList->Add(fHistEventPlaneTruth);
1456 //*****************************************************PIDQA histos*****************************************************//
1460 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1461 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1464 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1465 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);
1466 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1467 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1468 fOutputList->Add(fHistoNSigma);
1473 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1474 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1477 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1478 TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
1479 Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1480 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1481 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1482 fOutputList->Add(fHistoNSigma);
1487 if(fPIDType==Bayes){//use bayesianPID
1488 fPIDCombined = new AliPIDCombined();
1489 fPIDCombined->SetDefaultTPCPriors();//****************************************Need to know about it
1491 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1494 TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
1495 Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1496 fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1497 fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
1498 fOutputList->Add(fHistoBayes);
1501 TH2F *fHistoBayesTPC=new TH2F(Form("probBayes_TPC_%d",ipart),
1502 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1503 fHistoBayesTPC->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1504 fHistoBayesTPC->GetYaxis()->SetTitle(Form("Bayes prob TPC %s",kParticleSpeciesName[ipart]));
1505 fOutputList->Add(fHistoBayesTPC);
1507 TH2F *fHistoBayesTOF=new TH2F(Form("probBayes_TOF_%d",ipart),
1508 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1509 fHistoBayesTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1510 fHistoBayesTOF->GetYaxis()->SetTitle(Form("Bayes prob TOF %s",kParticleSpeciesName[ipart]));
1511 fOutputList->Add(fHistoBayesTOF);
1513 TH2F *fHistoBayesTPCTOF=new TH2F(Form("probBayes_TPCTOF_%d",ipart),
1514 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1515 fHistoBayesTPCTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1516 fHistoBayesTPCTOF->GetYaxis()->SetTitle(Form("Bayes prob TPCTOF %s",kParticleSpeciesName[ipart]));
1517 fOutputList->Add(fHistoBayesTPCTOF);
1521 //nsigma separation power plot
1522 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1525 TH2F *Pi_Ka_sep=new TH2F(Form("Pi_Ka_sep_%d",ipid),
1526 Form("Pi_Ka separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1527 Pi_Ka_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1528 Pi_Ka_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1529 fOutputList->Add(Pi_Ka_sep);
1531 TH2F *Pi_Pr_sep=new TH2F(Form("Pi_Pr_sep_%d",ipid),
1532 Form("Pi_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1533 Pi_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1534 Pi_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1535 fOutputList->Add(Pi_Pr_sep);
1537 TH2F *Ka_Pr_sep=new TH2F(Form("Ka_Pr_sep_%d",ipid),
1538 Form("Ka_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1539 Ka_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1540 Ka_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1541 fOutputList->Add(Ka_Pr_sep);
1545 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1546 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1549 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1550 TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
1551 Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1552 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1553 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1554 fOutputList->Add(fHistoNSigma);
1559 if (fAnalysisType == "MCAOD"){
1560 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1561 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1564 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1565 TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
1566 Form("n#sigma for MC %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 for(Int_t idet=0;idet<fNDetectors;idet++){
1575 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1577 if(idet==fTOF)maxy=1.1;
1578 TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
1579 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1580 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1581 fOutputList->Add(fHistoPID);
1584 //PID signal plot, before PID cut
1585 for(Int_t idet=0;idet<fNDetectors;idet++){
1587 if(idet==fTOF)maxy=1.1;
1588 TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
1589 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1590 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1591 fOutputList->Add(fHistoPID);
1594 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
1595 PostData(2, fOutputList);
1596 if(fRequestEventPlane) PostData(3, fList);
1597 AliInfo("Finished setting up the Output");
1599 TH1::AddDirectory(oldStatus);
1604 //-------------------------------------------------------------------------------
1605 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1608 if(fAnalysisType == "AOD") {
1612 }//AOD--analysis-----
1614 else if(fAnalysisType == "MCAOD") {
1623 //-------------------------------------------------------------------------
1624 void AliTwoParticlePIDCorr::doMCAODevent()
1626 AliVEvent *event = InputEvent();
1627 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1628 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1630 AliError("Cannot get the AOD event");
1634 // count all events(physics triggered)
1635 fEventCounter->Fill(1);
1644 gReactionPlane = 999.;
1646 // get centrality object and check quality(valid for p-Pb and Pb-Pb; coming soon for pp 7 TeV)
1647 Double_t cent_v0=-1.0;
1648 Double_t effcent=1.0;
1649 Double_t refmultReco =0.0;
1651 //check the PIDResponse handler
1654 // get mag. field required for twotrack efficiency cut
1656 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1658 //check for TClonesArray(truth track MC information)
1659 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1661 AliFatal("Error: MC particles branch not found!\n");
1665 //check for AliAODMCHeader(truth event MC information)
1666 AliAODMCHeader *header=NULL;
1667 header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1669 printf("MC header branch not found!\n");
1673 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
1674 Float_t zVtxmc =header->GetVtxZ();
1675 if(TMath::Abs(zVtxmc)>fzvtxcut) return;
1677 // For productions with injected signals, figure out above which label to skip particles/tracks
1679 if (fInjectedSignals)
1681 AliGenEventHeader* eventHeader = 0;
1686 AliFatal("fInjectedSignals set but no MC header found");
1688 headers = header->GetNCocktailHeaders();
1689 eventHeader = header->GetCocktailHeader(0);
1693 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
1694 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
1695 AliError("First event header not found. Skipping this event.");
1696 //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1699 skipParticlesAbove = eventHeader->NProduced();
1700 AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
1703 if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL"))
1705 //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
1706 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
1707 refmultReco = GetAcceptedEventMultiplicity(aod,kFALSE);
1708 if(refmultTruth<=0 || refmultReco<=0) return;
1709 cent_v0=refmultTruth;
1712 cent_v0=GetAcceptedEventMultiplicity(aod,kTRUE); //centrality value; 2nd argument has no meaning
1713 if(cent_v0<0.) return;
1716 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1718 //get the event plane in case of PbPb
1719 if(fRequestEventPlane){
1720 gReactionPlane=GetEventPlane(aod,kTRUE,cent_v0);//get the truth event plane
1721 if(gReactionPlane==999.) return;
1725 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)
1727 //TObjArray* tracksMCtruth=0;
1728 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
1729 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
1733 //There is a small difference between zvtx and zVtxmc??????
1734 //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1735 //cout<<"***********************************************zvtx="<<zvtx<<endl;
1737 //now process the truth particles(for both efficiency & correlation function)
1738 Int_t nMCTrack = fArrayMC->GetEntriesFast();
1740 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
1741 { //MC truth track loop starts
1743 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1746 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1750 //consider only charged particles
1751 if(partMC->Charge() == 0) continue;
1753 //consider only primary particles; neglect all secondary particles including from weak decays
1754 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1757 //remove injected signals(primaries above <maxLabel>)
1758 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1761 Bool_t isduplicate=kFALSE;
1762 if (fRemoveDuplicates)
1764 for (Int_t j=iMC+1; j<nMCTrack; ++j)
1765 {//2nd trutuh loop starts
1766 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1768 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1771 if (partMC->GetLabel() == partMC2->GetLabel())
1776 }//2nd truth loop ends
1778 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1780 //give only kinematic cuts at the truth level
1781 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1782 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
1784 if(!partMC) continue;//for safety
1786 //To determine multiplicity in case of PP
1788 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1789 //only physical primary(all/unidentified)
1790 if(ffillhistQATruth)
1792 MCtruthpt->Fill(partMC->Pt());
1793 MCtrutheta->Fill(partMC->Eta());
1794 MCtruthphi->Fill(partMC->Phi());
1797 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1798 Int_t particletypeTruth=-999;
1799 if (TMath::Abs(pdgtruth)==211)
1801 particletypeTruth=SpPion;
1802 if(ffillhistQATruth)
1804 MCtruthpionpt->Fill(partMC->Pt());
1805 MCtruthpioneta->Fill(partMC->Eta());
1806 MCtruthpionphi->Fill(partMC->Phi());
1809 if (TMath::Abs(pdgtruth)==321)
1811 particletypeTruth=SpKaon;
1812 if(ffillhistQATruth)
1814 MCtruthkaonpt->Fill(partMC->Pt());
1815 MCtruthkaoneta->Fill(partMC->Eta());
1816 MCtruthkaonphi->Fill(partMC->Phi());
1819 if(TMath::Abs(pdgtruth)==2212)
1821 particletypeTruth=SpProton;
1822 if(ffillhistQATruth)
1824 MCtruthprotonpt->Fill(partMC->Pt());
1825 MCtruthprotoneta->Fill(partMC->Eta());
1826 MCtruthprotonphi->Fill(partMC->Phi());
1829 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)
1831 if(fRequestEventPlane){
1832 FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
1835 // -- Fill THnSparse for efficiency and contamination calculation
1836 if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) 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
1838 Double_t primmctruth[4] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
1841 fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
1842 if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for primary truth mesons(3)
1843 if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for primary truth kaons+protons(4)
1844 if (TMath::Abs(pdgtruth)==211) fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
1845 if (TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
1846 if(TMath::Abs(pdgtruth)==2212) fTrackHistEfficiency[2]->Fill(primmctruth,0);//for protons
1849 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1850 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1852 Short_t chargeval=0;
1853 if(partMC->Charge()>0) chargeval=1;
1854 if(partMC->Charge()<0) chargeval=-1;
1855 if(chargeval==0) continue;
1856 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);
1857 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1858 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1859 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1861 }//MC truth track loop ends
1863 //*********************still in event loop
1865 if (fSampleType=="PbPb"){
1866 if (fRandomizeReactionPlane)//only for TRuth MC??
1868 Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
1869 Double_t angle = TMath::TwoPi() * centralityDigits;
1870 AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
1871 ShiftTracks(tracksMCtruth, angle);
1875 Float_t weghtval=1.0;
1877 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1879 //Fill Correlations for MC truth particles(same event)
1880 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1881 Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1884 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
1885 if (pool2 && pool2->IsReady())
1886 {//start mixing only when pool->IsReady
1887 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1888 {//proceed only when no. of trigger particles >0 in current event
1889 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
1890 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1891 { //pool event loop start
1892 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1893 if(!bgTracks6) continue;
1894 Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1896 }// pool event loop ends mixing case
1897 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1898 } //if pool->IsReady() condition ends mixing case
1900 //still in main event loop
1903 if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1907 //still in main event loop
1909 if(tracksMCtruth) delete tracksMCtruth;
1911 //now deal with reco tracks
1913 //detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
1914 if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) cent_v0=refmultReco;
1915 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1917 if(fRequestEventPlane){
1918 gReactionPlane = GetEventPlane(aod,kFALSE,cent_v0);//get the reconstructed event plane
1919 if(gReactionPlane==999.) return;
1922 //TObjArray* tracksUNID=0;
1923 TObjArray* tracksUNID = new TObjArray;
1924 tracksUNID->SetOwner(kTRUE);
1926 //TObjArray* tracksID=0;
1927 TObjArray* tracksID = new TObjArray;
1928 tracksID->SetOwner(kTRUE);
1931 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1934 Double_t trackscount=0.0;
1935 // loop over reconstructed tracks
1936 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1937 { // reconstructed track loop starts
1938 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1939 if (!track) continue;
1940 //get the corresponding MC track at the truth level (doing reco matching)
1941 AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
1942 if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
1944 //remove injected signals
1945 if(fInjectedSignals)
1947 AliAODMCParticle* mother = recomatched;
1949 while (!mother->IsPhysicalPrimary())
1950 {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
1951 if (mother->GetMother() < 0)
1957 mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1963 Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
1966 if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
1967 }//remove injected signals
1969 if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
1971 Bool_t isduplicate2=kFALSE;
1972 if (fRemoveDuplicates)
1974 for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
1976 AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
1977 if (!track2) continue;
1978 AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
1979 if(!recomatched2) continue;
1981 if (track->GetLabel() == track2->GetLabel())
1988 if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
1990 fHistQA[11]->Fill(track->GetTPCNcls());
1991 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
1993 if(tracktype==0) continue;
1994 if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
1996 if(!track) continue;//for safety
1997 //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
2000 //check for eta , phi holes
2001 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2002 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2004 Int_t particletypeMC=-9999;
2006 //tag all particles as unidentified
2007 particletypeMC=unidentified;
2009 Float_t effmatrix=1.;
2011 // -- Fill THnSparse for efficiency calculation
2012 if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) 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
2013 //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)
2015 //Clone & Reduce track list(TObjArray) for unidentified particles
2016 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2018 Short_t chargeval=0;
2019 if(track->Charge()>0) chargeval=1;
2020 if(track->Charge()<0) chargeval=-1;
2021 if(chargeval==0) continue;
2022 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2023 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
2024 LRCParticlePID* copy = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2025 copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
2026 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2029 //get the pdg code of the corresponding truth particle
2030 Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
2032 Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
2033 if(ffillefficiency) {
2034 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
2035 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
2036 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
2037 if(TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions
2038 if(TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
2039 if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
2041 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
2042 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
2043 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
2044 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
2045 if( TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions
2046 if( TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
2047 if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
2051 //now start the particle identification process:)
2053 Float_t dEdx = track->GetTPCsignal();
2054 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2056 if(HasTOFPID(track))
2058 Double_t beta = GetBeta(track);
2059 fHistoTOFbeta->Fill(track->Pt(), beta);
2062 //do track identification(nsigma method)
2063 particletypeMC=GetParticle(track,fFIllPIDQAHistos);//******************************problem is here
2065 switch(TMath::Abs(pdgCode)){
2067 if(fFIllPIDQAHistos){
2068 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2069 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
2070 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
2071 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
2076 if(fFIllPIDQAHistos){
2077 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2078 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
2079 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpKaon,ipid));
2080 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
2085 if(fFIllPIDQAHistos){
2086 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2087 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
2088 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpPion,ipid));
2089 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
2096 //2-d TPCTOF map(for each Pt interval)
2097 if(HasTOFPID(track)){
2098 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2099 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2100 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2103 //Pt, Eta , Phi distribution of the reconstructed identified particles
2106 if (particletypeMC==SpPion)
2108 fPionPt->Fill(track->Pt());
2109 fPionEta->Fill(track->Eta());
2110 fPionPhi->Fill(track->Phi());
2112 if (particletypeMC==SpKaon)
2114 fKaonPt->Fill(track->Pt());
2115 fKaonEta->Fill(track->Eta());
2116 fKaonPhi->Fill(track->Phi());
2118 if (particletypeMC==SpProton)
2120 fProtonPt->Fill(track->Pt());
2121 fProtonEta->Fill(track->Eta());
2122 fProtonPhi->Fill(track->Phi());
2126 //for misidentification fraction calculation(do it with tuneonPID)
2127 if(particletypeMC==SpPion )
2129 if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
2130 if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
2131 if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
2132 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
2134 if(particletypeMC==SpKaon )
2136 if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
2137 if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
2138 if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
2139 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
2141 if(particletypeMC==SpProton )
2143 if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
2144 if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
2145 if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
2146 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
2148 if(particletypeMC==SpUndefined )
2150 if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
2151 if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
2152 if(TMath::Abs(pdgCode)==2212) fUNIDcont->Fill(5.,track->Pt());
2153 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fUNIDcont->Fill(7.,track->Pt());
2156 if(particletypeMC==SpUndefined) continue;
2158 if(fRequestEventPlane){
2159 FillPIDEventPlane(cent_v0,particletypeMC,track->Phi(),gReactionPlane);
2162 //fill tracking efficiency
2165 if(particletypeMC==SpPion || particletypeMC==SpKaon)
2167 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) {
2168 fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
2169 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
2172 if(particletypeMC==SpKaon || particletypeMC==SpProton)
2174 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) {
2175 fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
2176 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
2179 if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) {
2180 fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
2181 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
2183 if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
2184 fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
2185 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
2187 if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
2188 fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
2189 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
2193 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2195 Short_t chargeval=0;
2196 if(track->Charge()>0) chargeval=1;
2197 if(track->Charge()<0) chargeval=-1;
2198 if(chargeval==0) continue;
2199 if (fapplyTrigefficiency || fapplyAssoefficiency)
2200 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
2201 LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2202 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2203 tracksID->Add(copy1);
2205 }// if(tracktype==1) condition structure ands
2207 }//reco track loop ends
2209 //*************************************************************still in event loop
2214 //fill the centrality/multiplicity distribution of the selected events
2215 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2217 if (fSampleType=="pPb" || fSampleType=="PbPb") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2219 //count selected events having centrality betn 0-100%
2220 fEventCounter->Fill(13);
2222 //***************************************event no. counting
2223 Bool_t isbaryontrig=kFALSE;
2224 Bool_t ismesontrig=kFALSE;
2225 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2227 if(tracksID && tracksID->GetEntriesFast()>0)
2229 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2230 { //trigger loop starts
2231 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2233 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2234 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2235 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2236 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2238 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2239 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2242 //same event delte-eta, delta-phi plot
2243 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2244 {//same event calculation starts
2245 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2246 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)
2249 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2250 {//same event calculation starts
2251 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)
2252 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2255 //still in main event loop
2257 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2258 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2259 if (pool && pool->IsReady())
2260 {//start mixing only when pool->IsReady
2261 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2262 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2263 { //pool event loop start
2264 TObjArray* bgTracks = pool->GetEvent(jMix);
2265 if(!bgTracks) continue;
2266 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2267 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2268 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2269 Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2270 }// pool event loop ends mixing case
2272 } //if pool->IsReady() condition ends mixing case
2275 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2277 }//mixing with unidentified particles
2279 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2280 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2281 if (pool1 && pool1->IsReady())
2282 {//start mixing only when pool->IsReady
2283 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2284 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2285 { //pool event loop start
2286 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2287 if(!bgTracks2) continue;
2288 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2289 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2290 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2291 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2293 }// pool event loop ends mixing case
2294 } //if pool1->IsReady() condition ends mixing case
2298 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2300 }//mixing with identified particles
2302 //no. of events analyzed
2303 fEventCounter->Fill(15);
2306 if(tracksUNID) delete tracksUNID;
2308 if(tracksID) delete tracksID;
2311 PostData(1, fOutput);
2314 //________________________________________________________________________
2315 void AliTwoParticlePIDCorr::doAODevent()
2319 AliVEvent *event = InputEvent();
2320 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
2321 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
2323 AliError("Cannot get the AOD event");
2328 fEventCounter->Fill(1);
2330 if (!fPID) return;//this should be available with each event even if we don't do PID selection
2338 gReactionPlane = 999.;
2340 Double_t cent_v0= -999.;
2341 Double_t effcent=1.0;
2343 Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
2346 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
2347 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
2350 // check event cuts and fill event histograms and return the centrality or reference multiplicity value
2351 if((cent_v0 = GetAcceptedEventMultiplicity(aod,kFALSE)) < 0){
2355 //get the event plane in case of PbPb
2356 if(fRequestEventPlane){
2357 gReactionPlane = GetEventPlane(aod,kFALSE,cent_v0);
2358 if(gReactionPlane==999.) return;
2361 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
2362 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
2364 TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
2365 tracksID->SetOwner(kTRUE); // IMPORTANT!
2370 Bool_t fTrigPtmin1=kFALSE;
2371 Bool_t fTrigPtmin2=kFALSE;
2372 Bool_t fTrigPtJet=kFALSE;
2374 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2375 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2376 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2377 if (!track) continue;
2378 fHistQA[11]->Fill(track->GetTPCNcls());
2379 Int_t particletype=-9999;//required for PID filling
2380 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2381 if(tracktype!=1) continue;
2383 if(!track) continue;//for safety
2385 //check for eta , phi holes
2386 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2387 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2391 //if no applyefficiency , set the eff factor=1.0
2392 Float_t effmatrix=1.0;
2394 //tag all particles as unidentified that passed filterbit & kinematic cuts
2395 particletype=unidentified;
2397 //To count the no. of tracks having an accepted track in a certain PT(e.g. Jet Pt) range
2398 if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
2399 if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
2400 if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
2403 if (fSampleType=="pp") 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
2406 //to reduce memory consumption in pool
2407 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2409 //Clone & Reduce track list(TObjArray) for unidentified particles
2410 Short_t chargeval=0;
2411 if(track->Charge()>0) chargeval=1;
2412 if(track->Charge()<0) chargeval=-1;
2413 if(chargeval==0) continue;
2414 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2415 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
2416 LRCParticlePID* copy = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2417 copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2418 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2421 //now start the particle identificaion process:)
2423 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
2425 Float_t dEdx = track->GetTPCsignal();
2426 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2428 //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)
2429 if(HasTOFPID(track))
2431 Double_t beta = GetBeta(track);
2432 fHistoTOFbeta->Fill(track->Pt(), beta);
2436 //track identification(using nsigma method)
2437 particletype=GetParticle(track,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
2440 //2-d TPCTOF map(for each Pt interval)
2441 if(HasTOFPID(track)){
2442 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2443 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2444 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2447 //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!!!!!
2448 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)!!!!!!!!!!!
2450 if(fRequestEventPlane){
2451 FillPIDEventPlane(cent_v0,particletype,track->Phi(),gReactionPlane);
2455 //Pt, Eta , Phi distribution of the reconstructed identified particles
2458 if (particletype==SpPion)
2460 fPionPt->Fill(track->Pt());
2461 fPionEta->Fill(track->Eta());
2462 fPionPhi->Fill(track->Phi());
2464 if (particletype==SpKaon)
2466 fKaonPt->Fill(track->Pt());
2467 fKaonEta->Fill(track->Eta());
2468 fKaonPhi->Fill(track->Phi());
2470 if (particletype==SpProton)
2472 fProtonPt->Fill(track->Pt());
2473 fProtonEta->Fill(track->Eta());
2474 fProtonPhi->Fill(track->Phi());
2478 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2480 Short_t chargeval=0;
2481 if(track->Charge()>0) chargeval=1;
2482 if(track->Charge()<0) chargeval=-1;
2483 if(chargeval==0) continue;
2484 if (fapplyTrigefficiency || fapplyAssoefficiency)
2485 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
2486 LRCParticlePID* copy1 = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2487 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2488 tracksID->Add(copy1);
2490 } //track loop ends but still in event loop
2492 if(trackscount<1.0){
2493 if(tracksUNID) delete tracksUNID;
2494 if(tracksID) delete tracksID;
2498 if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
2499 if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
2500 if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
2502 Float_t weightval=1.0;
2506 //fill the centrality/multiplicity distribution of the selected events
2507 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2509 if(fSampleType=="pPb" || fSampleType=="PbPb") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2511 //count selected events having centrality betn 0-100%
2512 fEventCounter->Fill(13);
2514 //***************************************event no. counting
2515 Bool_t isbaryontrig=kFALSE;
2516 Bool_t ismesontrig=kFALSE;
2517 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2519 if(tracksID && tracksID->GetEntriesFast()>0)
2521 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2522 { //trigger loop starts
2523 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2525 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2526 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2527 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2528 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2530 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2531 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2535 //same event delta-eta-deltaphi plot
2537 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2538 {//same event calculation starts
2539 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2540 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)
2543 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2544 {//same event calculation starts
2545 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)
2546 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2549 //still in main event loop
2553 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2554 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2555 if (pool && pool->IsReady())
2556 {//start mixing only when pool->IsReady
2557 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2558 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2559 { //pool event loop start
2560 TObjArray* bgTracks = pool->GetEvent(jMix);
2561 if(!bgTracks) continue;
2562 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2563 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2564 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2565 Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2566 }// pool event loop ends mixing case
2567 } //if pool->IsReady() condition ends mixing case
2570 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2572 }//mixing with unidentified particles
2575 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2576 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2577 if (pool1 && pool1->IsReady())
2578 {//start mixing only when pool->IsReady
2579 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2580 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2581 { //pool event loop start
2582 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2583 if(!bgTracks2) continue;
2584 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2585 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2586 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2587 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2589 }// pool event loop ends mixing case
2590 } //if pool1->IsReady() condition ends mixing case
2594 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2596 }//mixing with identified particles
2599 //no. of events analyzed
2600 fEventCounter->Fill(15);
2602 //still in main event loop
2605 if(tracksUNID) delete tracksUNID;
2607 if(tracksID) delete tracksID;
2610 PostData(1, fOutput);
2612 } // *************************event loop ends******************************************//_______________________________________________________________________
2613 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
2615 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
2617 TObjArray* tracksClone = new TObjArray;
2618 tracksClone->SetOwner(kTRUE);
2620 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2622 LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
2623 LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval());
2624 copy100->SetUniqueID(particle->GetUniqueID());
2625 tracksClone->Add(copy100);
2631 //--------------------------------------------------------------------------------
2632 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)
2635 //before calling this function check that either trackstrig & tracksasso are available
2637 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
2638 TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
2639 TArrayF eta(input->GetEntriesFast());
2640 for (Int_t i=0; i<input->GetEntriesFast(); i++)
2641 eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
2644 Int_t jmax=trackstrig->GetEntriesFast();
2646 jmax=tracksasso->GetEntriesFast();
2648 // identify K, Lambda candidates and flag those particles
2649 // a TObject bit is used for this
2650 const UInt_t kResonanceDaughterFlag = 1 << 14;
2651 if (fRejectResonanceDaughters > 0)
2653 Double_t resonanceMass = -1;
2654 Double_t massDaughter1 = -1;
2655 Double_t massDaughter2 = -1;
2656 const Double_t interval = 0.02;
2657 switch (fRejectResonanceDaughters)
2659 case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
2660 case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
2661 case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
2662 default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
2665 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2666 trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2667 for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
2668 tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2670 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2672 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2673 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
2675 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2677 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
2678 if (trig->IsEqual(asso)) continue;
2680 if (trig->Charge() * asso->Charge() > 0) continue;
2682 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2684 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
2686 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2688 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
2690 trig->SetBit(kResonanceDaughterFlag);
2691 asso->SetBit(kResonanceDaughterFlag);
2693 // Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
2700 //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
2702 Float_t TriggerPtMin=fminPtTrig;
2703 Float_t TriggerPtMax=fmaxPtTrig;
2704 Int_t HighestPtTriggerIndx=-99999;
2705 TH1* triggerWeighting = 0;
2707 if(fSelectHighestPtTrig || fWeightPerEvent)//**************add this data member to the constructor
2709 if (fWeightPerEvent)
2712 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID) axis = fTHnTrigcount->GetGrid(0)->GetGrid()->GetAxis(2);
2713 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) axis = fTHnTrigcountMCTruthPrim->GetGrid(0)->GetGrid()->GetAxis(2);
2714 triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
2716 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2717 { //trigger loop starts(highest Pt trigger selection)
2718 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2720 Float_t trigpt=trig->Pt();
2721 //to avoid overflow qnd underflow
2722 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2723 Int_t particlepidtrig=trig->getparticle();
2724 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2726 Float_t trigeta=trig->Eta();
2728 // some optimization
2729 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2732 if (fOnlyOneEtaSide != 0)
2734 if (fOnlyOneEtaSide * trigeta < 0)
2737 if (fTriggerSelectCharge != 0)
2738 if (trig->Charge() * fTriggerSelectCharge < 0)
2741 if (fRejectResonanceDaughters > 0)
2742 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2744 if(fSelectHighestPtTrig){
2745 if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
2747 HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
2748 TriggerPtMin=trigpt;
2752 if (fWeightPerEvent) triggerWeighting->Fill(trigpt);
2754 }//trigger loop ends(highest Pt trigger selection)
2756 }//******************(fSelectHighestPtTrig || fWeightPerEvent) condition ends
2759 //two particle correlation filling
2760 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2761 { //trigger loop starts
2762 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2764 Float_t trigpt=trig->Pt();
2765 //to avoid overflow qnd underflow
2766 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2767 Int_t particlepidtrig=trig->getparticle();
2768 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2770 Float_t trigeta=trig->Eta();
2772 // some optimization
2773 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2776 if (fOnlyOneEtaSide != 0)
2778 if (fOnlyOneEtaSide * trigeta < 0)
2781 if (fTriggerSelectCharge != 0)
2782 if (trig->Charge() * fTriggerSelectCharge < 0)
2785 if (fRejectResonanceDaughters > 0)
2786 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2788 if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
2789 if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
2792 Float_t trigphi=trig->Phi();
2793 Float_t trackefftrig=1.0;
2794 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
2796 // Event plane (determine psi bin)
2797 Double_t gPsiMinusPhi = 0.;
2798 Double_t gPsiMinusPhiBin = -10.;
2799 if(fRequestEventPlane){
2800 gPsiMinusPhi = TMath::Abs(trigphi - ReactionPlane);
2801 //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)
2802 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2803 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
2804 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2805 gPsiMinusPhiBin = 0.0;
2807 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2808 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2809 gPsiMinusPhiBin = 0.0;
2812 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
2813 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
2814 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
2815 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
2816 gPsiMinusPhiBin = 1.0;
2818 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
2819 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
2820 gPsiMinusPhiBin = 2.0;
2823 gPsiMinusPhiBin = 3.0;
2825 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
2828 //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
2831 Int_t eventplaneAxis=0;
2832 if(fRequestEventPlane) eventplaneAxis=1;
2833 if(fcontainPIDtrig && SetChargeAxis==0) dim=4+eventplaneAxis;
2834 if(!fcontainPIDtrig && SetChargeAxis==2) dim=4+eventplaneAxis;
2835 if(fcontainPIDtrig && SetChargeAxis==2) dim=5+eventplaneAxis;
2836 trigval= new Double_t[dim];
2839 trigval[2] = trigpt;
2841 if(fRequestEventPlane){
2842 trigval[3] = gPsiMinusPhiBin;
2843 if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = particlepidtrig;
2844 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
2845 if(fcontainPIDtrig && SetChargeAxis==2) {
2846 trigval[4] = particlepidtrig;
2847 trigval[5] = trig->Charge();
2851 if(!fRequestEventPlane){
2852 if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = particlepidtrig;
2853 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
2854 if(fcontainPIDtrig && SetChargeAxis==2) {
2855 trigval[3] = particlepidtrig;
2856 trigval[4] = trig->Charge();
2862 if (fWeightPerEvent)
2864 // leads effectively to a filling of one entry per filled trigger particle pT bin
2865 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(trigval[2]);
2866 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
2867 trackefftrig *= triggerWeighting->GetBinContent(weightBin);
2871 //trigger particle counting for both same and mixed event case;;;;;step=0->same event case;;;;;step=1->mixed event case;;;;;;
2872 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
2873 if(fillup=="trigassoUNID" ) {
2874 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2875 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2878 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
2879 if(fillup=="trigassoUNID" )
2881 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2882 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2885 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
2886 if(fillup=="trigUNIDassoID")
2888 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2889 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2892 //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
2893 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
2894 if(fillup=="trigIDassoID")
2896 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2897 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2900 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
2901 if(fillup=="trigIDassoUNID" )
2903 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2904 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2907 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
2908 if(fillup=="trigIDassoID")
2910 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2911 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2915 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!!!!
2916 if(mixcase==kFALSE) fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig);
2917 if(mixcase==kTRUE && firstTime) fTHnTrigcountMCTruthPrim->Fill(trigval,1,1.0/trackefftrig);
2920 //asso loop starts within trigger loop
2921 for(Int_t j=0;j<jmax;j++)
2923 LRCParticlePID *asso=0;
2925 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
2927 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2931 //to avoid overflow and underflow
2932 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
2934 //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
2936 if(!tracksasso && j==i) continue;
2938 // 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 Pi 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)
2939 // if (tracksasso && trig->IsEqual(asso)) continue;
2941 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
2944 if (asso->Pt() >= trig->Pt()) continue;
2946 Int_t particlepidasso=asso->getparticle();
2947 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
2950 if (fAssociatedSelectCharge != 0)
2951 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
2953 if (fSelectCharge > 0)
2956 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
2960 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
2966 if (trigeta < 0 && asso->Eta() < trigeta)
2968 if (trigeta > 0 && asso->Eta() > trigeta)
2972 if (fRejectResonanceDaughters > 0)
2973 if (asso->TestBit(kResonanceDaughterFlag))
2975 // Printf("Skipped j=%d", j);
2980 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
2982 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2986 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2988 fControlConvResoncances->Fill(0.0, mass);
2990 if (mass < 0.04*0.04)
2996 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2998 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
3000 const Float_t kK0smass = 0.4976;
3002 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
3004 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
3006 fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
3008 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
3014 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3016 Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
3017 Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3019 const Float_t kLambdaMass = 1.115;
3021 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
3023 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
3025 fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
3027 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3030 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
3032 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3034 fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
3036 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3041 if (twoTrackEfficiencyCut)
3043 // the variables & cuthave been developed by the HBT group
3044 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
3045 Float_t phi1 = trig->Phi();
3046 Float_t pt1 = trig->Pt();
3047 Float_t charge1 = trig->Charge();
3048 Float_t phi2 = asso->Phi();
3049 Float_t pt2 = asso->Pt();
3050 Float_t charge2 = asso->Charge();
3052 Float_t deta= trigeta - eta[j];
3055 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
3058 // check first boundaries to see if is worth to loop and find the minimum
3059 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
3060 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
3062 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
3064 Float_t dphistarminabs = 1e5;
3065 Float_t dphistarmin = 1e5;
3067 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
3069 for (Double_t rad=0.8; rad<2.51; rad+=0.01)
3071 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
3073 Float_t dphistarabs = TMath::Abs(dphistar);
3075 if (dphistarabs < dphistarminabs)
3077 dphistarmin = dphistar;
3078 dphistarminabs = dphistarabs;
3081 if(mixcase==kFALSE) fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3082 if(mixcase==kTRUE) fTwoTrackDistancePtmix[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3084 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
3086 // 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);
3089 if(mixcase==kFALSE) fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3090 if(mixcase==kTRUE) fTwoTrackDistancePtmix[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3096 Float_t weightperevent=weight;
3097 Float_t trackeffasso=1.0;
3098 if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
3099 //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
3100 Float_t deleta=trigeta-eta[j];
3101 Float_t delphi=PhiRange(trigphi-asso->Phi());
3103 Float_t delpt=trigpt-asso->Pt();
3104 //fill it with/without two track efficiency cut
3105 if(mixcase==kFALSE) fTwoTrackDistancePtdip->Fill(deleta, delphi, TMath::Abs(delpt));//for same event
3106 if(mixcase==kTRUE) fTwoTrackDistancePtdipmix->Fill(deleta, delphi, TMath::Abs(delpt));//for mixed event
3108 //here get the two particle efficiency correction factor
3109 Float_t effweight=trackefftrig*trackeffasso*weightperevent;
3110 // if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
3112 Int_t dimused=kTrackVariablesPair+eventplaneAxis;
3113 vars= new Double_t[dimused];
3122 if(fRequestEventPlane)
3124 vars[6]=gPsiMinusPhiBin;
3128 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
3129 vars[dimension]=trig->Charge();
3130 vars[dimension+1]=asso->Charge();
3132 if(fcontainPIDtrig && !fcontainPIDasso){
3133 vars[dimension]=particlepidtrig;
3134 if(SetChargeAxis==2){
3135 vars[dimension+1]=trig->Charge();
3136 vars[dimension+2]=asso->Charge();
3139 if(!fcontainPIDtrig && fcontainPIDasso){
3140 vars[dimension]=particlepidasso;
3141 if(SetChargeAxis==2){
3142 vars[dimension+1]=trig->Charge();
3143 vars[dimension+2]=asso->Charge();
3146 if(fcontainPIDtrig && fcontainPIDasso){
3147 vars[dimension]=particlepidtrig;
3148 vars[dimension+1]=particlepidasso;
3149 if(SetChargeAxis==2){
3150 vars[dimension+2]=trig->Charge();
3151 vars[dimension+3]=asso->Charge();
3155 if (fWeightPerEvent)
3157 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
3158 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3159 effweight *= triggerWeighting->GetBinContent(weightBin);
3163 //Fill Correlation ThnSparses
3164 if(fillup=="trigassoUNID")
3166 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,0,1.0/effweight);
3167 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
3169 if(fillup=="trigIDassoID")
3171 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,0,1.0/effweight);
3172 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
3174 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
3175 {//in this case effweight should be 1 always
3176 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
3177 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
3179 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
3181 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
3182 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
3188 }//trigger loop ends
3190 if (triggerWeighting)
3192 delete triggerWeighting;
3193 triggerWeighting = 0;
3197 //--------------------------------------------------------------------------------
3198 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
3200 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
3202 Float_t effvalue=1.;
3204 if(parpid==unidentified)
3206 effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
3207 effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
3208 effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
3209 effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
3210 effvalue=effcorection[5]->GetBinContent(effVars);
3212 if(parpid==SpPion || parpid==SpKaon)
3214 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3216 effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
3217 effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
3218 effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
3219 effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
3220 effvalue=effcorection[3]->GetBinContent(effVars);
3225 effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
3226 effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
3227 effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
3228 effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
3229 effvalue=effcorection[0]->GetBinContent(effVars);
3234 effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
3235 effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
3236 effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
3237 effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
3238 effvalue=effcorection[1]->GetBinContent(effVars);
3243 if(parpid==SpProton)
3245 effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
3246 effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
3247 effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
3248 effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
3249 effvalue=effcorection[2]->GetBinContent(effVars);
3252 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3254 if(parpid==SpProton || parpid==SpKaon)
3256 effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
3257 effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
3258 effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
3259 effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
3260 effvalue=effcorection[4]->GetBinContent(effVars);
3263 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
3264 if(effvalue==0.) effvalue=1.;
3269 //---------------------------------------------------------------------------------
3271 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield, Bool_t fill)
3274 if(!track) return 0;
3275 Bool_t trackOK = track->TestFilterBit(fFilterBit);
3276 if(!trackOK) return 0;
3277 if (fTrackStatus != 0 && !CheckTrack(track)) return 0;
3278 //select only primary traks(for data & reco MC tracks)
3279 if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
3280 if(track->Charge()==0) return 0;
3281 if (fill) fHistQA[12]->Fill(track->GetTPCNcls());
3284 dz = track->ZAtDCA();
3285 if (fill) fHistQA[6]->Fill(dxy);
3286 if (fill) fHistQA[7]->Fill(dz);
3287 Float_t chi2ndf = track->Chi2perNDF();
3288 if (fill) fHistQA[13]->Fill(chi2ndf);
3289 // Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
3290 Float_t nCrossedRowsTPC = track->GetTPCNCrossedRows();
3291 if (fill) fHistQA[14]->Fill(nCrossedRowsTPC);
3292 //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
3293 if (track->GetTPCNclsF()>0) {
3294 Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
3295 if (fill) fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
3298 Float_t pt=track->Pt();
3299 if(pt< fminPt || pt> fmaxPt) return 0;
3300 if(TMath::Abs(track->Eta())> fmaxeta) return 0;
3301 if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
3302 //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required???
3304 if (fdcacut && fDCAXYCut)
3311 AliAODTrack* clone =(AliAODTrack*) track->Clone();
3312 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
3317 // Printf("%f", ((AliAODTrack*)part)->DCA());
3318 // Printf("%f", pos[0]);
3319 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
3323 if (fSharedClusterCut >= 0)
3325 Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
3326 if (frac > fSharedClusterCut)
3329 if (fill) fHistQA[8]->Fill(pt);
3330 if (fill) fHistQA[9]->Fill(track->Eta());
3331 if (fill) fHistQA[10]->Fill(track->Phi());
3334 //________________________________________________________________________________
3335 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
3337 //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 )
3338 Float_t pt=track->Pt();
3340 //plot the separation power
3342 Double_t bethe[AliPID::kSPECIES]={0.};
3343 Double_t sigma_TPC[AliPID::kSPECIES]={0.};
3345 Double_t Pi_Ka_sep[NSigmaPIDType+1]={0.};
3346 Double_t Pi_Pr_sep[NSigmaPIDType+1]={0.};
3347 Double_t Ka_Pr_sep[NSigmaPIDType+1]={0.};
3350 Double_t ptpc = track->GetTPCmomentum();
3351 Int_t dEdxN = track->GetTPCsignalN();
3352 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
3353 bethe[ipart] = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
3354 //Double_t diff = dEdx - bethe;
3355 sigma_TPC[ipart] = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
3356 //nSigma[ipart] = diff / sigma;
3358 Pi_Ka_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kKaon])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kKaon])/2.0);
3359 Pi_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kProton])/2.0);
3360 Ka_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kKaon]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kKaon]+sigma_TPC[AliPID::kProton])/2.0);
3363 Double_t sigma_TOF[AliPID::kSPECIES]={0.};
3365 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3367 Double_t timei[AliPID::kSPECIES];
3368 track->GetIntegratedTimes(timei);
3369 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { sigma_TOF[ipart]= fPID->GetTOFResponse().GetExpectedSigma(track->P(), timei[ipart], AliPID::ParticleMass(ipart));}
3370 Pi_Ka_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kKaon])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kKaon])/2.0);
3371 Pi_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kProton])/2.0);
3372 Ka_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kKaon]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kKaon]+sigma_TOF[AliPID::kProton])/2.0);
3374 Pi_Ka_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Ka_sep[NSigmaTPC]*Pi_Ka_sep[NSigmaTPC]+Pi_Ka_sep[NSigmaTOF]*Pi_Ka_sep[NSigmaTOF]);
3375 Pi_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Pr_sep[NSigmaTPC]*Pi_Pr_sep[NSigmaTPC]+Pi_Pr_sep[NSigmaTOF]*Pi_Pr_sep[NSigmaTOF]);
3376 Ka_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Ka_Pr_sep[NSigmaTPC]*Ka_Pr_sep[NSigmaTPC]+Ka_Pr_sep[NSigmaTOF]*Ka_Pr_sep[NSigmaTOF]);
3380 //fill separation power histograms
3381 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3383 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3384 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3385 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3386 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3387 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3388 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3390 if(HasTOFPID(track) && pt>fPtTOFPIDmin && ipid!=0){
3391 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3392 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3393 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3394 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3395 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3396 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3401 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
3402 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
3403 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
3404 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
3406 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
3407 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
3409 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3412 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
3413 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
3414 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
3416 nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
3417 nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
3418 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
3424 // if TOF is missing and below fPtTOFPID only the TPC information is used
3425 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
3426 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
3427 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
3431 //set data member fnsigmas
3432 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
3433 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
3434 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
3436 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
3437 fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
3438 fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
3439 fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
3441 //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)
3442 fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
3443 fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
3444 fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
3447 //Fill NSigma SeparationPlot
3448 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3449 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3450 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3451 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
3452 h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
3458 //----------------------------------------------------------------------------
3459 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
3461 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3462 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
3463 //get the identity of the particle with the minimum Nsigma
3464 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3467 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3468 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
3469 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
3472 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3473 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
3474 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
3476 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3477 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3478 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3479 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3481 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3482 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3483 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3484 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3489 if(fdiffPIDcutvalues){
3490 if(track->Pt()<=4) fNSigmaPID=fPIDCutval1;
3491 if(track->Pt()>4 && track->Pt()<=6) fNSigmaPID=fPIDCutval2;
3492 if(track->Pt()>6 && track->Pt()<=8) fNSigmaPID=fPIDCutval3;
3493 if(track->Pt()>8) fNSigmaPID=fPIDCutval4;
3496 // guess the particle based on the smaller nsigma (within fNSigmaPID)
3497 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
3499 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
3500 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)
3502 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3503 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3504 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
3505 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
3510 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
3512 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3513 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3514 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
3515 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
3520 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
3522 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3523 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3524 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
3525 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
3531 // else, return undefined
3537 //------------------------------------------------------------------------------------------
3538 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
3539 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3541 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
3543 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
3545 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
3548 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
3550 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3553 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3554 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
3555 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
3558 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3559 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
3560 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
3562 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3563 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3564 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3565 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3567 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3568 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3569 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3570 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3574 // Actually the tracks in the overlapping region(in TPC-TOF nSigma plane) will be ignored
3576 if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
3577 if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
3578 if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
3583 //fill NSigma distr for double counting
3584 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3585 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
3586 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3587 if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)))continue;//not filling TOF and combined if no TOF PID
3588 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
3589 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
3596 return fHasDoubleCounting;
3599 //////////////////////////////////////////////////////////////////////////////////////////////////
3601 Bool_t* AliTwoParticlePIDCorr::GetAllCompatibleIdentitiesNSigma(AliAODTrack * trk,Bool_t FIllQAHistos){
3602 //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
3603 Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
3604 IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
3608 //////////////////////////////////////////////////////////////////////////////////////////////////
3610 UInt_t AliTwoParticlePIDCorr::CalcPIDCombined(AliAODTrack *track, Int_t detMask, Double_t* prob) const{
3612 // Bayesian PID calculation
3614 for(Int_t i=0;i<AliPID::kSPECIES;i++)
3618 fPIDCombined->SetDetectorMask(detMask);
3620 return fPIDCombined->ComputeProbabilities((AliAODTrack*)track, fPID, prob);
3623 //////////////////////////////////////////////////////////////////////////////////////////////////
3625 Int_t AliTwoParticlePIDCorr::GetIDBayes(AliAODTrack * trk, Bool_t FIllQAHistos){
3627 Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
3630 //Filling of Probability histos
3631 Double_t probTPC[AliPID::kSPECIES]={0.};
3632 Double_t probTOF[AliPID::kSPECIES]={0.};
3633 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
3635 UInt_t detUsedTPC = 0;
3636 UInt_t detUsedTOF = 0;
3637 UInt_t detUsedTPCTOF = 0;
3639 //get the TPC probability
3640 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
3641 detUsedTPC = fPIDCombined->ComputeProbabilities(trk, fPID, probTPC);
3642 if(detUsedTPC == AliPIDResponse::kDetTPC)
3644 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3646 TH2F *h=GetHistogram2D(Form("probBayes_TPC_%d",ipart));
3647 if(ipart==0) h->Fill(trk->Pt(),probTPC[AliPID::kPion]);
3648 if(ipart==1) h->Fill(trk->Pt(),probTPC[AliPID::kKaon]);
3649 if(ipart==2) h->Fill(trk->Pt(),probTPC[AliPID::kProton]);
3653 //get the TOF probability
3654 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
3655 detUsedTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTOF);
3656 if(detUsedTOF == AliPIDResponse::kDetTOF)
3658 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3659 TH2F *h=GetHistogram2D(Form("probBayes_TOF_%d",ipart));
3660 if(ipart==0) h->Fill(trk->Pt(),probTOF[AliPID::kPion]);
3661 if(ipart==1) h->Fill(trk->Pt(),probTOF[AliPID::kKaon]);
3662 if(ipart==2) h->Fill(trk->Pt(),probTOF[AliPID::kProton]);
3666 //get the TPC-TOF probability
3667 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
3668 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTPCTOF);
3669 if(detUsedTPCTOF == (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))
3671 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3672 TH2F *h=GetHistogram2D(Form("probBayes_TPCTOF_%d",ipart));
3673 if(ipart==0) h->Fill(trk->Pt(),probTPCTOF[AliPID::kPion]);
3674 if(ipart==1) h->Fill(trk->Pt(),probTPCTOF[AliPID::kKaon]);
3675 if(ipart==2) h->Fill(trk->Pt(),probTPCTOF[AliPID::kProton]);
3680 Double_t probBayes[AliPID::kSPECIES];
3683 if(HasTOFPID(trk) && trk->Pt()>fPtTOFPIDmin){//use TOF information
3684 detUsed = CalcPIDCombined(trk, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
3685 if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return SpUndefined;//check that TPC and TOF are used
3687 detUsed = CalcPIDCombined(trk,AliPIDResponse::kDetTPC, probBayes);
3688 if(detUsed != AliPIDResponse::kDetTPC)return SpUndefined;//check that TPC is used
3691 //the probability has to be normalized to one, we check it
3693 for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
3694 if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
3695 AliFatal("Bayesian probability not normalized to one");
3698 if(fdiffPIDcutvalues){
3699 if(trk->Pt()<=4) fBayesCut=fPIDCutval1;
3700 if(trk->Pt()>4 && trk->Pt()<=6) fBayesCut=fPIDCutval2;
3701 if(trk->Pt()>6 && trk->Pt()<=8) fBayesCut=fPIDCutval3;
3702 if(trk->Pt()>8) fBayesCut=fPIDCutval4;
3706 //probabilities are normalized to one, if the cut is above .5 there is no problem
3707 if(probBayes[AliPID::kPion]>fBayesCut && IDs[SpPion]==1){
3708 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpPion));
3709 h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
3712 else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[SpKaon]==1){
3713 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpKaon));
3714 h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
3717 else if(probBayes[AliPID::kProton]>fBayesCut && IDs[SpProton]==1){
3718 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpProton));
3719 h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
3728 //////////////////////////////////////////////////////////////////////////////////////////////////
3729 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
3730 //return the specie according to the minimum nsigma value
3731 //no double counting, this has to be evaluated using CheckDoubleCounting()
3733 Int_t ID=SpUndefined;
3735 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
3739 if(fPIDType==Bayes){//use bayesianPID
3742 AliFatal("PIDCombined object has to be set in the steering macro");
3745 ID = GetIDBayes(trk,FIllQAHistos);
3747 }else{ //use nsigma PID
3749 ID=FindMinNSigma(trk,FIllQAHistos);
3750 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
3752 HasDC=GetDoubleCounting(trk,FIllQAHistos);
3753 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3754 if(HasDC[ipart]==kTRUE) ID = SpUndefined;
3758 //Fill PID signal plot
3759 if(ID != SpUndefined){
3760 for(Int_t idet=0;idet<fNDetectors;idet++){
3761 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
3762 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3763 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3764 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3767 //Fill PID signal plot without cuts
3768 for(Int_t idet=0;idet<fNDetectors;idet++){
3769 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
3770 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3771 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3772 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3778 //-------------------------------------------------------------------------------------
3780 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
3783 AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
3784 if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
3785 //ULong_t status=track->GetStatus();
3786 //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
3787 //if (track->GetTPCsignal() <= 0.) return kFALSE;
3788 // 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.
3792 //___________________________________________________________
3795 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
3797 // check TOF matched track
3798 //ULong_t status=track->GetStatus();
3799 //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
3800 AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
3801 if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
3802 if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
3803 //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
3804 //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
3805 // if (probMis > 0.01) return kFALSE;
3806 if(fRemoveTracksT0Fill)
3808 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
3809 if (startTimeMask < 0)return kFALSE;
3814 //________________________________________________________________________
3815 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
3817 //it is called only when TOF PID is available
3818 //TOF beta calculation
3819 Double_t tofTime=track->GetTOFsignal();
3821 Double_t c=TMath::C()*1.E-9;// m/ns
3822 Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
3823 Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
3824 tofTime -= startTime; // subtract startTime to the signal
3825 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
3831 Double_t p = track->P();
3832 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
3834 track->GetIntegratedTimes(timei);
3835 return timei[0]/time;
3838 //------------------------------------------------------------------------------------------------------
3840 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)
3842 // calculate inv mass squared
3843 // same can be achieved, but with more computing time with
3844 /*TLorentzVector photon, p1, p2;
3845 p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
3846 p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
3850 Float_t tantheta1 = 1e10;
3852 if (eta1 < -1e-10 || eta1 > 1e-10)
3853 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
3855 Float_t tantheta2 = 1e10;
3856 if (eta2 < -1e-10 || eta2 > 1e-10)
3857 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
3859 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3860 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3862 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 ) ) );
3866 //---------------------------------------------------------------------------------
3868 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)
3870 // calculate inv mass squared approximately
3872 Float_t tantheta1 = 1e10;
3874 if (eta1 < -1e-10 || eta1 > 1e-10)
3876 Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
3877 tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3880 Float_t tantheta2 = 1e10;
3881 if (eta2 < -1e-10 || eta2 > 1e-10)
3883 Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
3884 tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3887 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3888 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3891 Float_t deltaPhi = TMath::Abs(phi1 - phi2);
3892 while (deltaPhi > TMath::TwoPi())
3893 deltaPhi -= TMath::TwoPi();
3894 if (deltaPhi > TMath::Pi())
3895 deltaPhi = TMath::TwoPi() - deltaPhi;
3897 Float_t cosDeltaPhi = 0;
3898 if (deltaPhi < TMath::Pi()/3)
3899 cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
3900 else if (deltaPhi < 2*TMath::Pi()/3)
3901 cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
3903 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
3905 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
3907 // Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
3911 //--------------------------------------------------------------------------------
3912 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)
3915 // calculates dphistar
3918 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
3920 static const Double_t kPi = TMath::Pi();
3923 // if (dphistar > 2 * kPi)
3924 // dphistar -= 2 * kPi;
3925 // if (dphistar < -2 * kPi)
3926 // dphistar += 2 * kPi;
3929 dphistar = kPi * 2 - dphistar;
3930 if (dphistar < -kPi)
3931 dphistar = -kPi * 2 - dphistar;
3932 if (dphistar > kPi) // might look funny but is needed
3933 dphistar = kPi * 2 - dphistar;
3937 //_________________________________________________________________________
3939 void AliTwoParticlePIDCorr ::DefineEventPool()
3941 Int_t MaxNofEvents=1000;
3942 const Int_t NofVrtxBins=10+(1+10)*2;
3943 Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
3944 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
3945 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210
3947 //default values are for centrality
3948 Int_t NofCentBins=15;
3949 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
3951 if(fCentralityMethod.EndsWith("_MANUAL"))
3953 Int_t NofCentBins=9;
3954 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....
3956 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
3961 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
3963 //if(!fPoolMgr) return kFALSE;
3968 //------------------------------------------------------------------------
3969 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
3971 // This method is a copy from AliUEHist::GetBinning
3972 // takes the binning from <configuration> identified by <tag>
3973 // configuration syntax example:
3974 // 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
3977 // returns bin edges which have to be deleted by the caller
3979 TString config(configuration);
3980 TObjArray* lines = config.Tokenize("\n");
3981 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
3983 TString line(lines->At(i)->GetName());
3984 if (line.BeginsWith(TString(tag) + ":"))
3986 line.Remove(0, strlen(tag) + 1);
3987 line.ReplaceAll(" ", "");
3988 TObjArray* binning = line.Tokenize(",");
3989 Double_t* bins = new Double_t[binning->GetEntriesFast()];
3990 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
3991 bins[j] = TString(binning->At(j)->GetName()).Atof();
3993 nBins = binning->GetEntriesFast() - 1;
4002 AliFatal(Form("Tag %s not found in %s", tag, configuration));
4006 //____________________________________________________________________
4008 Bool_t AliTwoParticlePIDCorr::CheckTrack(AliAODTrack * part)
4010 // check if the track status flags are set
4012 UInt_t status=((AliAODTrack*)part)->GetStatus();
4013 if ((status & fTrackStatus) == fTrackStatus)
4017 //________________________________________________________________________
4019 Bool_t AliTwoParticlePIDCorr::AcceptEventCentralityWeight(Double_t centrality)
4021 // rejects "randomly" events such that the centrality gets flat
4022 // uses fCentralityWeights histogram
4024 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
4026 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
4027 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
4029 Bool_t result = kFALSE;
4030 if (centralityDigits < weight)
4033 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
4038 //____________________________________________________________________
4039 void AliTwoParticlePIDCorr::ShiftTracks(TObjArray* tracks, Double_t angle)
4041 // shifts the phi angle of all tracks by angle
4042 // 0 <= angle <= 2pi
4044 for (Int_t i=0; i<tracks->GetEntriesFast(); ++i)
4046 LRCParticlePID *part=(LRCParticlePID*)(tracks->UncheckedAt(i));
4048 Double_t newAngle = part->Phi() + angle;
4049 if (newAngle >= TMath::TwoPi())
4050 newAngle -= TMath::TwoPi();
4052 part->SetPhi(newAngle);
4057 //________________________________________________________________________
4058 void AliTwoParticlePIDCorr::SetVZEROCalibrationFile(const char* filename,const char* lhcPeriod) {
4059 //Function to setup the VZERO gain equalization
4060 //============Get the equilization map============//
4061 TFile *calibrationFile = TFile::Open(filename);
4062 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
4063 Printf("No calibration file found!!!");
4067 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
4069 Printf("Calibration TList not found!!!");
4073 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
4074 if(!fHistVZEROAGainEqualizationMap) {
4075 Printf("VZERO-A calibration object not found!!!");
4078 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
4079 if(!fHistVZEROCGainEqualizationMap) {
4080 Printf("VZERO-C calibration object not found!!!");
4084 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
4085 if(!fHistVZEROChannelGainEqualizationMap) {
4086 Printf("VZERO channel calibration object not found!!!");
4091 //________________________________________________________________________
4092 Double_t AliTwoParticlePIDCorr::GetChannelEqualizationFactor(Int_t run,Int_t channel) {
4094 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4096 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
4097 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4098 if(gRunNumber == run)
4099 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
4105 //________________________________________________________________________
4106 Double_t AliTwoParticlePIDCorr::GetEqualizationFactor(Int_t run, const char* side) {
4108 if(!fHistVZEROAGainEqualizationMap) return 1.0;
4110 TString gVZEROSide = side;
4111 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
4112 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4113 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
4114 if(gRunNumber == run) {
4115 if(gVZEROSide == "A")
4116 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
4117 else if(gVZEROSide == "C")
4118 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
4124 //________________________________________________________________________
4125 Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliAODEvent *event){
4126 //Function that returns the reference multiplicity from AODs (data or reco MC, Not for Truth)
4127 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
4128 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
4129 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
4131 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
4133 Printf("ERROR: AOD header not available");
4136 Int_t gRunNumber = header->GetRunNumber();
4137 Float_t bSign1=header->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
4140 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++)
4141 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
4142 AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
4143 if (!track) continue;
4144 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
4145 if(tracktype!=1) continue;
4147 if(!track) continue;//for safety
4149 gRefMultiplicityTPC += 1.0;
4153 if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
4154 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
4155 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
4156 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
4159 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
4160 else if(iChannel >= 32)
4161 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
4164 //Equalization of gain
4165 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
4167 gRefMultiplicityVZEROA /= gFactorA;
4168 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
4170 gRefMultiplicityVZEROC /= gFactorC;
4171 if((gFactorA != 0)&&(gFactorC != 0))
4172 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
4175 //EQVZERO vs TPC multiplicity
4176 fHistEQVZEROvsTPCmultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
4177 fHistEQVZEROAvsTPCmultiplicity->Fill(gRefMultiplicityVZEROA,gRefMultiplicityTPC);
4178 fHistEQVZEROCvsTPCmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityTPC);
4180 //EQVZERO vs VZERO multiplicity
4181 fHistVZEROCvsEQVZEROCmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),gRefMultiplicityVZEROC);
4182 fHistVZEROAvsEQVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0A(),gRefMultiplicityVZEROA);
4184 //VZEROC vs VZEROA multiplicity
4185 fHistVZEROCvsVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),event->GetVZEROData()->GetMTotV0A());
4187 //EQVZEROC vs EQVZEROA multiplicity
4188 fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
4191 if(fCentralityMethod == "TRACKS_MANUAL") {
4192 gRefMultiplicity = gRefMultiplicityTPC;
4193 fHistRefmult->Fill(3.,gRefMultiplicityTPC);
4195 else if(fCentralityMethod == "V0M_MANUAL"){
4196 gRefMultiplicity = gRefMultiplicityVZERO;
4197 fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
4200 else if(fCentralityMethod == "V0A_MANUAL"){
4201 gRefMultiplicity = gRefMultiplicityVZEROA;
4202 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
4205 else if(fCentralityMethod == "V0C_MANUAL"){
4206 gRefMultiplicity = gRefMultiplicityVZEROC;
4207 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
4211 return gRefMultiplicity;
4214 //-------------------------------------------------------------------------------------------------------
4215 Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliAODEvent *event, Bool_t truth){
4217 if(!event) return -1;
4218 // get centrality object and check quality
4219 Double_t cent_v0=-1;
4220 Double_t nooftrackstruth=0;
4222 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
4224 AliCentrality *centralityObj=0;
4225 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
4226 if(!header) return -1;
4227 centralityObj = header->GetCentralityP();
4228 // if (centrality->GetQuality() != 0) return ;
4232 fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
4233 fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
4234 fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
4235 if(fSampleType=="pp") fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
4236 if(fSampleType=="pp") fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
4237 if(fSampleType=="pp") fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
4239 if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
4240 if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
4242 cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
4245 }//centralitymethod condition
4247 else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL")//data or RecoMc and also for TRUTH
4249 if(!truth){//for data or RecoMC
4250 cent_v0 = GetReferenceMultiplicityVZEROFromAOD(event);
4251 }//for data or RecoMC
4253 if(truth){//condition for TRUTH case
4254 //check for TClonesArray(truth track MC information)
4255 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
4257 //AliFatal("Error: MC particles branch not found!\n");
4260 //now process the truth particles(for both efficiency & correlation function)
4261 Int_t nMCTrack = fArrayMC->GetEntriesFast();
4263 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
4264 {//MC truth track loop starts
4266 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
4269 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
4273 //consider only charged particles
4274 if(partMC->Charge() == 0) continue;
4276 //consider only primary particles; neglect all secondary particles including from weak decays
4277 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
4280 //remove injected signals(primaries above <maxLabel>)
4281 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
4284 Bool_t isduplicate=kFALSE;
4285 if (fRemoveDuplicates)
4287 for (Int_t j=iMC+1; j<nMCTrack; ++j)
4288 {//2nd trutuh loop starts
4289 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
4291 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
4294 if (partMC->GetLabel() == partMC2->GetLabel())
4299 }//2nd truth loop ends
4301 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
4304 if (fCentralityMethod=="V0M_MANUAL") {
4305 if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8) continue;
4306 if (partMC->Eta() < -3.7 || partMC->Eta() > -1.7) continue;
4308 else if (fCentralityMethod=="V0A_MANUAL") {
4309 if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8) continue;}
4310 else if (fCentralityMethod=="V0C_MANUAL") {
4311 if(partMC->Eta() > -1.7 || partMC->Eta() < -3.7) continue;}
4312 else if (fCentralityMethod=="TRACKS_MANUAL") {
4313 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4314 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
4316 else{//basically returns the tracks manual case
4317 //give only kinematic cuts at the truth level
4318 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4319 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
4322 //To determine multiplicity in case of PP
4323 nooftrackstruth+= 1;;
4325 }//truth track loop ends
4326 cent_v0=nooftrackstruth;
4328 }//condition for TRUTH case
4330 }//end of MANUAL method
4332 else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC
4334 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4338 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4341 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
4342 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
4343 AliError("Event header not found. Skipping this event.");
4347 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4350 if (collGeometry) cent_v0 = collGeometry->ImpactParameter();
4352 }//end of Impact parameter method
4359 //-----------------------------------------------------------------------------------------
4360 Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliAODEvent *aod,Bool_t truth){
4361 //do the event selection(zvtx, pileup, centrality/multiplicity cut) and then return the value of the centrality of that event
4364 Float_t gRefMultiplicity = -1.;
4366 // check first event in chunk (is not needed for new reconstructions)
4367 if(fCheckFirstEventInChunk){
4368 AliAnalysisUtils ut;
4369 if(ut.IsFirstEventInChunk(aod))
4374 AliAnalysisUtils ut;
4375 ut.SetUseMVPlpSelection(kTRUE);
4376 ut.SetUseOutOfBunchPileUp(kTRUE);
4377 if(ut.IsPileUpEvent(aod))
4381 //count events after pileup selection
4382 fEventCounter->Fill(3);
4384 //vertex selection(is it fine for PP?)
4385 if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4386 trkVtx = aod->GetPrimaryVertex();
4387 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
4388 TString vtxTtl = trkVtx->GetTitle();
4389 if (!vtxTtl.Contains("VertexerTracks")) return -1;
4390 zvtx = trkVtx->GetZ();
4391 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
4392 if (!spdVtx || spdVtx->GetNContributors()<=0) return -1;
4393 TString vtxTyp = spdVtx->GetTitle();
4394 Double_t cov[6]={0};
4395 spdVtx->GetCovarianceMatrix(cov);
4396 Double_t zRes = TMath::Sqrt(cov[5]);
4397 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return -1;
4398 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return -1;
4400 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
4401 Int_t nVertex = aod->GetNumberOfVertices();
4403 trkVtx = aod->GetPrimaryVertex();
4404 Int_t nTracksPrim = trkVtx->GetNContributors();
4405 zvtx = trkVtx->GetZ();
4406 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
4407 // Reject TPC only vertex
4408 TString name(trkVtx->GetName());
4409 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return -1;
4411 // Select a quality vertex by number of tracks?
4412 if( nTracksPrim < fnTracksVertex ) {
4413 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
4416 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
4417 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
4419 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
4424 else if(fVertextype==0){//default case
4425 trkVtx = aod->GetPrimaryVertex();
4426 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
4427 zvtx = trkVtx->GetZ();
4429 trkVtx->GetCovarianceMatrix(fCov);
4430 if(fCov[5] == 0) return -1;//proper vertex resolution
4433 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
4434 return -1;//as there is no proper sample type
4437 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
4439 //count events having a proper vertex
4440 fEventCounter->Fill(5);
4442 if (TMath::Abs(zvtx) > fzvtxcut) return -1;
4444 //count events after vertex cut
4445 fEventCounter->Fill(7);
4448 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4450 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
4452 //get the centrality or multiplicity
4453 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)
4455 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)
4457 if(gRefMultiplicity<0) return -1;
4459 // take events only within the multiplicity class mentioned in the custom binning
4460 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
4462 //count events having proper centrality/ref multiplicity
4463 fEventCounter->Fill(9);
4466 // centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
4467 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
4469 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
4473 //count events after rejection due to centrality weighting
4474 fEventCounter->Fill(11);
4476 return gRefMultiplicity;
4479 //--------------------------------------------------------------------------------------------------------
4480 Float_t AliTwoParticlePIDCorr::GetEventPlane(AliAODEvent *event,Bool_t truth, Double_t v0Centr)
4482 // Get the event plane
4483 //reset Q vector info
4485 Int_t run = event->GetRunNumber();
4488 // Load the calibrations run dependent
4489 if(! fIsAfter2011) OpenInfoCalbration(run);
4495 if (v0Centr < 80){ // analysis only for 0-80% centrality classes
4497 if(v0Centr < 5) iC = 0;
4498 else if(v0Centr < 10) iC = 1;
4499 else if(v0Centr < 20) iC = 2;
4500 else if(v0Centr < 30) iC = 3;
4501 else if(v0Centr < 40) iC = 4;
4502 else if(v0Centr < 50) iC = 5;
4503 else if(v0Centr < 60) iC = 6;
4504 else if(v0Centr < 70) iC = 7;
4510 //reset Q vector info
4511 Double_t Qxa2 = 0, Qya2 = 0;
4512 Double_t Qxc2 = 0, Qyc2 = 0;
4513 Double_t Qxa3 = 0, Qya3 = 0;
4514 Double_t Qxc3 = 0, Qyc3 = 0;
4517 //MC: from reaction plane
4520 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4522 evplaneMC = header->GetReactionPlaneAngle();//[0, 360]
4523 //make it within [-pi/2,pi/2] to make it general
4524 if(evplaneMC > TMath::Pi()/2 && evplaneMC <= TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
4525 else if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
4526 fHistEventPlaneTruth->Fill(iC,evplaneMC);
4528 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4532 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4534 if (collGeometry){//get the reaction plane from MC header
4535 gReactionPlane = collGeometry->ReactionPlaneAngle();//[0,180]
4539 //taken from vnv0 code(get the TPC, V0A, V0C event plane using truth tracks)
4540 TClonesArray *mcArray = NULL;
4541 mcArray = (TClonesArray*)event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
4543 Float_t QxMCv2[3] = {0,0,0};
4544 Float_t QyMCv2[3] = {0,0,0};
4545 Float_t QxMCv3[3] = {0,0,0};
4546 Float_t QyMCv3[3] = {0,0,0};
4547 Float_t EvPlaneMCV2[3] = {0,0,0};
4548 Float_t EvPlaneMCV3[3] = {0,0,0};
4549 Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
4550 Float_t etaMax[3] = {4.88,-1.8,0.8};
4552 // analysis on MC tracks
4553 Int_t nMCtrack = mcArray->GetEntries() ;
4555 // EP computation with MC tracks
4556 for(Int_t iT=0;iT < nMCtrack;iT++){
4557 AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
4558 if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
4560 Float_t eta = mctr->Eta();
4561 for(Int_t iD=0;iD<3;iD++){
4562 if(eta > etaMin[iD] && eta < etaMax[iD]){
4563 Float_t phi = mctr->Phi();
4564 QxMCv2[iD] += TMath::Cos(2*phi);
4565 QyMCv2[iD] += TMath::Sin(2*phi);
4566 QxMCv3[iD] += TMath::Cos(3*phi);
4567 QyMCv3[iD] += TMath::Sin(3*phi);
4572 EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
4573 EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
4574 EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
4575 fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
4576 fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
4577 fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
4578 fgPsi2v0aMC = EvPlaneMCV2[0];
4579 fgPsi2v0cMC = EvPlaneMCV2[1];
4580 fgPsi2tpcMC = EvPlaneMCV2[2];
4583 EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
4584 EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
4585 EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
4586 fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
4587 fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
4588 fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
4589 fgPsi3v0aMC = EvPlaneMCV3[0];
4590 fgPsi3v0cMC = EvPlaneMCV3[1];
4591 fgPsi3tpcMC = EvPlaneMCV3[2];
4598 Int_t nAODTracks = event->GetNumberOfTracks();
4600 // TPC EP needed for resolution studies (TPC subevent)
4601 //AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
4602 //Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
4603 Double_t Qx2 = 0, Qy2 = 0;
4604 Double_t Qx3 = 0, Qy3 = 0;
4606 for(Int_t iT = 0; iT < nAODTracks; iT++) {
4608 AliAODTrack* aodTrack = event->GetTrack(iT);
4614 Bool_t trkFlag = aodTrack->TestFilterBit(1);
4616 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
4619 Double_t b[2] = {-99., -99.};
4620 Double_t bCov[3] = {-99., -99., -99.};
4623 AliAODTrack param(*aodTrack);
4624 if (!param.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov)){
4628 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
4631 Qx2 += TMath::Cos(2*aodTrack->Phi());
4632 Qy2 += TMath::Sin(2*aodTrack->Phi());
4633 Qx3 += TMath::Cos(3*aodTrack->Phi());
4634 Qy3 += TMath::Sin(3*aodTrack->Phi());
4638 Float_t evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
4639 Float_t evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
4641 fgPsi2tpc = evPlAng2;
4642 fgPsi3tpc = evPlAng3;
4644 fPhiRPTPC->Fill(iC,evPlAng2);
4645 fPhiRPTPCv3->Fill(iC,evPlAng3);
4650 AliAODVZERO* aodV0 = event->GetVZEROData();
4652 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
4653 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
4654 Float_t multv0 = aodV0->GetMultiplicity(iv0);
4657 if(fAnalysisType == "AOD"){//not for reco MC tracks, only for real data
4658 if (iv0 < 32){ // V0C
4659 Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4660 Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4661 Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4662 Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4664 Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4665 Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4666 Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4667 Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4671 if (iv0 < 32){ // V0C
4672 Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4673 Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4674 Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4675 Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4677 Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4678 Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4679 Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4680 Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4685 //grab for each centrality the proper histo with the Qx and Qy to do the recentering
4686 Double_t Qxamean2 = fMeanQ[iCcal][1][0];
4687 Double_t Qxarms2 = fWidthQ[iCcal][1][0];
4688 Double_t Qyamean2 = fMeanQ[iCcal][1][1];
4689 Double_t Qyarms2 = fWidthQ[iCcal][1][1];
4690 Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
4691 Double_t Qxarms3 = fWidthQv3[iCcal][1][0];
4692 Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
4693 Double_t Qyarms3 = fWidthQv3[iCcal][1][1];
4695 Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
4696 Double_t Qxcrms2 = fWidthQ[iCcal][0][0];
4697 Double_t Qycmean2 = fMeanQ[iCcal][0][1];
4698 Double_t Qycrms2 = fWidthQ[iCcal][0][1];
4699 Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
4700 Double_t Qxcrms3 = fWidthQv3[iCcal][0][0];
4701 Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
4702 Double_t Qycrms3 = fWidthQv3[iCcal][0][1];
4704 Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
4705 Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
4706 Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
4707 Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
4708 Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
4709 Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
4710 Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
4711 Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
4713 //to calculate 2nd order event plane with v0M
4714 Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
4715 /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
4716 Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
4717 /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
4719 //here the calculated event plane is within -Pi to +Pi(delete it , no use here , only for definition)
4720 Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
4721 Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
4722 Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
4726 Float_t evPlAngV0ACor2=999.;
4727 Float_t evPlAngV0CCor2=999.;
4728 Float_t evPlAngV0ACor3=999.;
4729 Float_t evPlAngV0CCor3=999.;
4732 if(fAnalysisType == "AOD"){
4733 evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
4734 evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
4735 evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
4736 evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
4739 evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
4740 evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
4741 evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
4742 evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
4746 AliEventplane *ep = event->GetEventplane();
4747 evPlAngV0ACor2 = ep->GetEventplane("V0A", event, 2);
4748 evPlAngV0CCor2 = ep->GetEventplane("V0C", event, 2);
4749 evPlAngV0ACor3 = ep->GetEventplane("V0A", event, 3);
4750 evPlAngV0CCor3 = ep->GetEventplane("V0C", event, 3);
4753 fgPsi2v0a = evPlAngV0ACor2;
4754 fgPsi2v0c = evPlAngV0CCor2;
4755 fgPsi3v0a = evPlAngV0ACor3;
4756 fgPsi3v0c = evPlAngV0CCor3;
4758 // Fill EP distribution histograms evPlAng
4760 fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
4761 fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
4763 fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
4764 fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
4766 // Fill histograms needed for resolution evaluation
4767 fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
4768 fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
4769 fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
4771 fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
4772 fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
4773 fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
4777 Float_t gVZEROEventPlane = -10.;
4778 Float_t gReactionPlane = -10.;
4779 Double_t qxTot = 0.0, qyTot = 0.0;
4781 AliEventplane *ep = event->GetEventplane();
4783 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
4784 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
4785 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
4786 gReactionPlane = gVZEROEventPlane;
4790 //return gReactionPlane;
4792 //make the final 2nd order event plane within 0 to Pi
4793 //using data and reco tracks only
4794 if(fgPsi2v0a!=999. && fgPsi2v0a < 0.) fgPsi2v0a += TMath::Pi();
4795 if(fgPsi2v0c!=999. && fgPsi2v0c < 0.) fgPsi2v0c += TMath::Pi();
4796 if(fgPsi2tpc!=999. && fgPsi2tpc < 0.) fgPsi2tpc += TMath::Pi();
4797 //using truth tracks only
4798 if(evplaneMC!=999. && evplaneMC < 0.) evplaneMC += TMath::Pi();
4799 if(fgPsi2v0aMC!=999. && fgPsi2v0aMC < 0.) fgPsi2v0aMC += TMath::Pi();
4800 if(fgPsi2v0cMC!=999. && fgPsi2v0cMC < 0.) fgPsi2v0cMC += TMath::Pi();
4801 if(fgPsi2tpcMC!=999. && fgPsi2tpcMC < 0.) fgPsi2tpcMC += TMath::Pi();
4802 //for the time being leave the 3rd order event planes within -pi/3 t0 +pi/3
4804 if(truth){//for truth MC
4805 if(fV2 && fEPdet=="header")gReactionPlane=evplaneMC;
4806 if(fV2 && fEPdet=="V0A")gReactionPlane=fgPsi2v0aMC;
4807 if(fV2 && fEPdet=="V0C")gReactionPlane=fgPsi2v0cMC;
4808 if(fV2 && fEPdet=="TPC")gReactionPlane=fgPsi2tpcMC;
4810 if(fV3 && fEPdet=="V0A")gReactionPlane=fgPsi3v0aMC;
4811 if(fV3 && fEPdet=="V0C")gReactionPlane=fgPsi3v0cMC;
4812 if(fV3 && fEPdet=="TPC")gReactionPlane=fgPsi3tpcMC;
4814 else{//for data and recoMC
4815 if(fV2 && fEPdet=="V0A")gReactionPlane=fgPsi2v0a;
4816 if(fV2 && fEPdet=="V0C")gReactionPlane=fgPsi2v0c;
4817 if(fV2 && fEPdet=="TPC")gReactionPlane=fgPsi2tpc;
4819 if(fV3 && fEPdet=="V0A")gReactionPlane=fgPsi3v0a;
4820 if(fV3 && fEPdet=="V0C")gReactionPlane=fgPsi3v0c;
4821 if(fV3 && fEPdet=="TPC")gReactionPlane=fgPsi3tpc;
4825 } //centrality cut condition
4827 return gReactionPlane;
4829 //------------------------------------------------------------------------------------------------------------------
4830 void AliTwoParticlePIDCorr::OpenInfoCalbration(Int_t run){
4831 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
4832 TFile *foadb = TFile::Open(oadbfilename.Data());
4835 printf("OADB file %s cannot be opened\n",oadbfilename.Data());
4839 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
4841 printf("OADB object hMultV0BefCorr is not available in the file\n");
4845 if(!(cont->GetObject(run))){
4846 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
4849 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
4851 TF1 *fpol0 = new TF1("fpol0","pol0");
4852 fMultV0->Fit(fpol0,"","",0,31);
4853 fV0Cpol = fpol0->GetParameter(0);
4854 fMultV0->Fit(fpol0,"","",32,64);
4855 fV0Apol = fpol0->GetParameter(0);
4857 for(Int_t iside=0;iside<2;iside++){
4858 for(Int_t icoord=0;icoord<2;icoord++){
4859 for(Int_t i=0;i < 9;i++){
4861 if(iside==0 && icoord==0)
4862 snprintf(namecont,100,"hQxc2_%i",i);
4863 else if(iside==1 && icoord==0)
4864 snprintf(namecont,100,"hQxa2_%i",i);
4865 else if(iside==0 && icoord==1)
4866 snprintf(namecont,100,"hQyc2_%i",i);
4867 else if(iside==1 && icoord==1)
4868 snprintf(namecont,100,"hQya2_%i",i);
4870 cont = (AliOADBContainer*) foadb->Get(namecont);
4872 printf("OADB object %s is not available in the file\n",namecont);
4876 if(!(cont->GetObject(run))){
4877 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
4880 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
4881 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
4884 if(iside==0 && icoord==0)
4885 snprintf(namecont,100,"hQxc3_%i",i);
4886 else if(iside==1 && icoord==0)
4887 snprintf(namecont,100,"hQxa3_%i",i);
4888 else if(iside==0 && icoord==1)
4889 snprintf(namecont,100,"hQyc3_%i",i);
4890 else if(iside==1 && icoord==1)
4891 snprintf(namecont,100,"hQya3_%i",i);
4893 cont = (AliOADBContainer*) foadb->Get(namecont);
4895 printf("OADB object %s is not available in the file\n",namecont);
4899 if(!(cont->GetObject(run))){
4900 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
4903 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
4904 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
4910 //____________________________________________________________________
4911 void AliTwoParticlePIDCorr::FillPIDEventPlane(Double_t centrality,Int_t par,Float_t trigphi,Float_t fReactionPlane)
4914 // Event plane (determine psi bin)
4915 Double_t gPsiMinusPhi = 0.;
4916 Double_t gPsiMinusPhiBin = -10.;
4917 if(fRequestEventPlane){
4918 gPsiMinusPhi = TMath::Abs(trigphi - fReactionPlane);
4920 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
4921 (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
4922 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
4923 gPsiMinusPhiBin = 0.0;
4925 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
4926 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
4927 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
4928 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
4929 gPsiMinusPhiBin = 1.0;
4931 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
4932 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
4933 gPsiMinusPhiBin = 2.0;
4936 gPsiMinusPhiBin = 3.0;
4938 fEventPlanePID->Fill(centrality,gPsiMinusPhiBin,(Float_t)par);
4942 //----------------------------------------------------------
4943 void AliTwoParticlePIDCorr::Terminate(Option_t *)
4945 // Draw result to screen, or perform fitting, normalizations
4946 // Called once at the end of the query
4947 fOutput = dynamic_cast<TList*> (GetOutputData(1));
4948 if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
4952 //------------------------------------------------------------------