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 **************************************************************************/
15 //Source code::dphicorrelations, TaskBFpsi, AliHelperPID
17 #include "AliTwoParticlePIDCorr.h"
18 #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"
73 #include "AliEventPoolManager.h"
74 #include "AliAnalysisUtils.h"
75 using namespace AliPIDNameSpace;
78 ClassImp(AliTwoParticlePIDCorr)
79 ClassImp(LRCParticlePID)
81 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
82 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
83 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
85 //________________________________________________________________________
86 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
90 fCentralityMethod("V0A"),
92 fRequestEventPlane(kFALSE),
93 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
98 fSharedClusterCut(-1),
100 skipParticlesAbove(0),
102 ffilltrigassoUNID(kFALSE),
103 ffilltrigUNIDassoID(kFALSE),
104 ffilltrigIDassoUNID(kTRUE),
105 ffilltrigIDassoID(kFALSE),
106 ffilltrigIDassoIDMCTRUTH(kFALSE),
107 fMaxNofMixingTracks(50000),
108 fPtOrderMCTruth(kTRUE),
109 fPtOrderDataReco(kTRUE),
110 fWeightPerEvent(kFALSE),
111 fTriggerSpeciesSelection(kFALSE),
112 fAssociatedSpeciesSelection(kFALSE),
113 fRandomizeReactionPlane(kFALSE),
114 fTriggerSpecies(SpPion),
115 fAssociatedSpecies(SpPion),
118 fSelectHighestPtTrig(kFALSE),
119 fcontainPIDtrig(kTRUE),
120 fcontainPIDasso(kFALSE),
122 frejectPileUp(kFALSE),
127 fselectprimaryTruth(kTRUE),
128 fonlyprimarydatareco(kFALSE),
131 ffillhistQAReco(kFALSE),
132 ffillhistQATruth(kFALSE),
133 kTrackVariablesPair(0),
165 fhistJetTrigestimate(0),
166 fCentralityCorrelation(0x0),
167 fHistVZEROAGainEqualizationMap(0),
168 fHistVZEROCGainEqualizationMap(0),
169 fHistVZEROChannelGainEqualizationMap(0),
170 fCentralityWeights(0),
173 fHistEQVZEROvsTPCmultiplicity(0x0),
174 fHistEQVZEROAvsTPCmultiplicity(0x0),
175 fHistEQVZEROCvsTPCmultiplicity(0x0),
176 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
177 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
178 fHistVZEROCvsVZEROAmultiplicity(0x0),
179 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
180 fHistVZEROSignal(0x0),
181 fHistEventPlaneReco(0x0),
182 fHistEventPlaneTruth(0x0),
183 fHistPsiMinusPhi(0x0),
184 fControlConvResoncances(0),
199 fCorrelatonTruthPrimary(0),
200 fCorrelatonTruthPrimarymix(0),
206 fTHnCorrIDUNIDmix(0),
208 fTHnTrigcountMCTruthPrim(0),
211 fAnalysisType("AOD"),
213 ftwoTrackEfficiencyCutDataReco(kTRUE),
214 twoTrackEfficiencyCutValue(0.02),
220 fRequestTOFPID(kTRUE),
221 fPIDType(NSigmaTPCTOF),
222 fFIllPIDQAHistos(kTRUE),
225 fdiffPIDcutvalues(kFALSE),
230 fHighPtKaonNSigmaPID(-1),
231 fHighPtKaonSigma(3.5),
232 fUseExclusiveNSigma(kFALSE),
233 fRemoveTracksT0Fill(kFALSE),
235 fTriggerSelectCharge(0),
236 fAssociatedSelectCharge(0),
237 fTriggerRestrictEta(-1),
238 fEtaOrdering(kFALSE),
239 fCutConversions(kFALSE),
240 fCutResonances(kFALSE),
241 fRejectResonanceDaughters(-1),
243 fInjectedSignals(kFALSE),
244 fRemoveWeakDecays(kFALSE),
245 fRemoveDuplicates(kFALSE),
246 fapplyTrigefficiency(kFALSE),
247 fapplyAssoefficiency(kFALSE),
248 ffillefficiency(kFALSE),
249 fmesoneffrequired(kFALSE),
250 fkaonprotoneffrequired(kFALSE),
255 for ( Int_t i = 0; i < 16; i++) {
259 for ( Int_t i = 0; i < 6; i++ ){
260 fTrackHistEfficiency[i] = NULL;
261 effcorection[i]=NULL;
266 for(Int_t ipart=0;ipart<NSpecies;ipart++)
267 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
268 fnsigmas[ipart][ipid]=999.;
270 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
273 //________________________________________________________________________
274 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
275 :AliAnalysisTaskSE(name),
278 fCentralityMethod("V0A"),
280 fRequestEventPlane(kFALSE),
281 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
286 fSharedClusterCut(-1),
288 skipParticlesAbove(0),
290 ffilltrigassoUNID(kFALSE),
291 ffilltrigUNIDassoID(kFALSE),
292 ffilltrigIDassoUNID(kTRUE),
293 ffilltrigIDassoID(kFALSE),
294 ffilltrigIDassoIDMCTRUTH(kFALSE),
295 fMaxNofMixingTracks(50000),
296 fPtOrderMCTruth(kTRUE),
297 fPtOrderDataReco(kTRUE),
298 fWeightPerEvent(kFALSE),
299 fTriggerSpeciesSelection(kFALSE),
300 fAssociatedSpeciesSelection(kFALSE),
301 fRandomizeReactionPlane(kFALSE),
302 fTriggerSpecies(SpPion),
303 fAssociatedSpecies(SpPion),
306 fSelectHighestPtTrig(kFALSE),
307 fcontainPIDtrig(kTRUE),
308 fcontainPIDasso(kFALSE),
310 frejectPileUp(kFALSE),
315 fselectprimaryTruth(kTRUE),
316 fonlyprimarydatareco(kFALSE),
319 ffillhistQAReco(kFALSE),
320 ffillhistQATruth(kFALSE),
321 kTrackVariablesPair(0),
353 fhistJetTrigestimate(0),
354 fCentralityCorrelation(0x0),
355 fHistVZEROAGainEqualizationMap(0),
356 fHistVZEROCGainEqualizationMap(0),
357 fHistVZEROChannelGainEqualizationMap(0),
358 fCentralityWeights(0),
361 fHistEQVZEROvsTPCmultiplicity(0x0),
362 fHistEQVZEROAvsTPCmultiplicity(0x0),
363 fHistEQVZEROCvsTPCmultiplicity(0x0),
364 fHistVZEROCvsEQVZEROCmultiplicity(0x0),
365 fHistVZEROAvsEQVZEROAmultiplicity(0x0),
366 fHistVZEROCvsVZEROAmultiplicity(0x0),
367 fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
368 fHistVZEROSignal(0x0),
369 fHistEventPlaneReco(0x0),
370 fHistEventPlaneTruth(0x0),
371 fHistPsiMinusPhi(0x0),
372 fControlConvResoncances(0),
387 fCorrelatonTruthPrimary(0),
388 fCorrelatonTruthPrimarymix(0),
394 fTHnCorrIDUNIDmix(0),
396 fTHnTrigcountMCTruthPrim(0),
399 fAnalysisType("AOD"),
401 ftwoTrackEfficiencyCutDataReco(kTRUE),
402 twoTrackEfficiencyCutValue(0.02),
408 fRequestTOFPID(kTRUE),
409 fPIDType(NSigmaTPCTOF),
410 fFIllPIDQAHistos(kTRUE),
413 fdiffPIDcutvalues(kFALSE),
418 fHighPtKaonNSigmaPID(-1),
419 fHighPtKaonSigma(3.5),
420 fUseExclusiveNSigma(kFALSE),
421 fRemoveTracksT0Fill(kFALSE),
423 fTriggerSelectCharge(0),
424 fAssociatedSelectCharge(0),
425 fTriggerRestrictEta(-1),
426 fEtaOrdering(kFALSE),
427 fCutConversions(kFALSE),
428 fCutResonances(kFALSE),
429 fRejectResonanceDaughters(-1),
431 fInjectedSignals(kFALSE),
432 fRemoveWeakDecays(kFALSE),
433 fRemoveDuplicates(kFALSE),
434 fapplyTrigefficiency(kFALSE),
435 fapplyAssoefficiency(kFALSE),
436 ffillefficiency(kFALSE),
437 fmesoneffrequired(kFALSE),
438 fkaonprotoneffrequired(kFALSE),
443 for ( Int_t i = 0; i < 16; i++) {
447 for ( Int_t i = 0; i < 6; i++ ){
448 fTrackHistEfficiency[i] = NULL;
449 effcorection[i]=NULL;
453 for(Int_t ipart=0;ipart<NSpecies;ipart++)
454 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
455 fnsigmas[ipart][ipid]=999.;
457 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
459 // The last in the above list should not have a comma after it
462 // Define input and output slots here (never in the dummy constructor)
463 // Input slot #0 works with a TChain - it is connected to the default input container
464 // Output slot #1 writes into a TH1 container
466 DefineOutput(1, TList::Class()); // for output list
467 DefineOutput(2, TList::Class());
471 //________________________________________________________________________
472 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
474 // Destructor. Clean-up the output list, but not the histograms that are put inside
475 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
476 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
481 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
486 if (fPID) delete fPID;
487 if (fPIDCombined) delete fPIDCombined;
490 //________________________________________________________________________
492 //////////////////////////////////////////////////////////////////////////////////////////////////
494 TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
495 // returns histo named name
496 return (TH2F*) fOutputList->FindObject(name);
499 //////////////////////////////////////////////////////////////////////////////////////////////////
501 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
505 // Puts the argument in the range [-pi/2,3 pi/2].
508 if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
509 if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
514 //________________________________________________________________________
515 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
518 // Called once (on the worker node)
519 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
520 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
521 fPID = inputHandler->GetPIDResponse();
523 //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
525 //get the efficiency correction map
527 // global switch disabling the reference
528 // (to avoid "Replacing existing TH1" if several wagons are created in train)
529 Bool_t oldStatus = TH1::AddDirectoryStatus();
530 TH1::AddDirectory(kFALSE);
532 fOutput = new TList();
533 fOutput->SetOwner(); // IMPORTANT!
535 fOutputList = new TList;
536 fOutputList->SetOwner();
537 fOutputList->SetName("PIDQAList");
539 fEventCounter = new TH1F("fEventCounter","EventCounter", 19, 0.5,19.5);
540 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
541 fEventCounter->GetXaxis()->SetBinLabel(3,"After PileUP Cut");//only for Data
542 fEventCounter->GetXaxis()->SetBinLabel(5,"Have A Vertex");
543 fEventCounter->GetXaxis()->SetBinLabel(7,"After vertex Cut");
544 fEventCounter->GetXaxis()->SetBinLabel(9,"After centrality flattening");
545 fEventCounter->GetXaxis()->SetBinLabel(11,"Within 0-100% centrality");
546 fEventCounter->GetXaxis()->SetBinLabel(13,"Event Analyzed");
547 //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
548 fOutput->Add(fEventCounter);
550 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
551 fOutput->Add(fEtaSpectrasso);
553 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
554 fOutput->Add(fphiSpectraasso);
556 if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
557 fOutput->Add(fCentralityCorrelation);
560 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")
562 TString gCentName[8] = {"V0A","V0C","V0M","V0AEq","V0CEq","V0MEq","CL1","ZNA"};
563 fHistCentStats = new TH2F("fHistCentStats",
564 "Centrality statistics;;Cent percentile",
565 8,-0.5,7.5,220,-5,105);
566 for(Int_t i = 1; i <= 8; i++){
567 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
568 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
570 fOutput->Add(fHistCentStats);
573 if(fCentralityMethod.EndsWith("_MANUAL"))
575 fhistcentrality=new TH1F("fhistcentrality","referencemultiplicity",30001,-0.5,30000.5);
576 fOutput->Add(fhistcentrality);
579 fhistcentrality=new TH1F("fhistcentrality","centrality",220,-5,105);
580 fOutput->Add(fhistcentrality);
583 if(fCentralityMethod.EndsWith("_MANUAL"))
585 TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
586 fHistRefmult = new TH2F("fHistRefmult",
587 "Reference multiplicity",
588 4,-0.5,3.5,10000,0,20000);
589 for(Int_t i = 1; i <= 4; i++){
590 fHistRefmult->GetXaxis()->SetBinLabel(i,gmultName[i-1].Data());
591 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
593 fOutput->Add(fHistRefmult);
596 //TPC vs EQVZERO multiplicity
597 fHistEQVZEROvsTPCmultiplicity = new TH2F("fHistEQVZEROvsTPCmultiplicity","EqVZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
598 fHistEQVZEROvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO multiplicity (a.u.)");
599 fHistEQVZEROvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
600 fOutput->Add(fHistEQVZEROvsTPCmultiplicity);
603 fHistEQVZEROAvsTPCmultiplicity = new TH2F("fHistEQVZEROAvsTPCmultiplicity","EqVZERO_A vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
604 fHistEQVZEROAvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
605 fHistEQVZEROAvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
606 fOutput->Add(fHistEQVZEROAvsTPCmultiplicity);
609 fHistEQVZEROCvsTPCmultiplicity = new TH2F("fHistEQVZEROCvsTPCmultiplicity","EqVZERO_C vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
610 fHistEQVZEROCvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
611 fHistEQVZEROCvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
612 fOutput->Add(fHistEQVZEROCvsTPCmultiplicity);
614 //EQVZERO vs VZERO multiplicity
615 fHistVZEROCvsEQVZEROCmultiplicity = new TH2F("fHistVZEROCvsEQVZEROCmultiplicity","EqVZERO_C vs VZERO_C multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
616 fHistVZEROCvsEQVZEROCmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
617 fHistVZEROCvsEQVZEROCmultiplicity->GetYaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
618 fOutput->Add(fHistVZEROCvsEQVZEROCmultiplicity);
621 fHistVZEROAvsEQVZEROAmultiplicity = new TH2F("fHistVZEROAvsEQVZEROAmultiplicity","EqVZERO_A vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
622 fHistVZEROAvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
623 fHistVZEROAvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
624 fOutput->Add(fHistVZEROAvsEQVZEROAmultiplicity);
627 //VZEROC vs VZEROA multiplicity
628 fHistVZEROCvsVZEROAmultiplicity = new TH2F("fHistVZEROCvsVZEROAmultiplicity","VZERO_C vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
629 fHistVZEROCvsVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
630 fHistVZEROCvsVZEROAmultiplicity->GetYaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
631 fOutput->Add(fHistVZEROCvsVZEROAmultiplicity);
635 //EQVZEROC vs EQVZEROA multiplicity
636 fHistEQVZEROCvsEQVZEROAmultiplicity = new TH2F("fHistEQVZEROCvsEQVZEROAmultiplicity","EqVZERO_C vs EqVZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
637 fHistEQVZEROCvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
638 fHistEQVZEROCvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
639 fOutput->Add(fHistEQVZEROCvsEQVZEROAmultiplicity);
641 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
642 fOutput->Add(fHistVZEROSignal);
645 if(fSampleType=="PbPb" && fRequestEventPlane){
647 fHistEventPlaneReco = new TH2F("fHistEventPlaneReco",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
648 fOutput->Add(fHistEventPlaneReco);
651 fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
652 fOutput->Add(fHistEventPlaneTruth);
654 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
655 fOutput->Add(fHistPsiMinusPhi);
659 if(fCutConversions || fCutResonances)
661 fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
662 fOutput->Add(fControlConvResoncances);
665 fHistoTPCdEdx = new TH2F("fHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
666 fOutputList->Add(fHistoTPCdEdx);
667 fHistoTOFbeta = new TH2F(Form("fHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
668 fOutputList->Add(fHistoTOFbeta);
670 fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
671 fOutputList->Add(fTPCTOFPion3d);
673 fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
674 fOutputList->Add(fTPCTOFKaon3d);
676 fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
677 fOutputList->Add(fTPCTOFProton3d);
681 fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);
682 fOutputList->Add(fPionPt);
683 fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);
684 fOutputList->Add(fPionEta);
685 fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);
686 fOutputList->Add(fPionPhi);
688 fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);
689 fOutputList->Add(fKaonPt);
690 fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);
691 fOutputList->Add(fKaonEta);
692 fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);
693 fOutputList->Add(fKaonPhi);
695 fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);
696 fOutputList->Add(fProtonPt);
697 fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);
698 fOutputList->Add(fProtonEta);
699 fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);
700 fOutputList->Add(fProtonPhi);
703 fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
704 fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
705 fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
706 fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);
707 fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
708 fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
709 fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
710 fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
711 fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
712 fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
713 fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
714 fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
715 fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
716 fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
717 fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
718 fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
720 for(Int_t i = 0; i < 16; i++)
722 fOutput->Add(fHistQA[i]);
725 Int_t eventplaneaxis=0;
727 if (fRequestEventPlane) eventplaneaxis=1;
729 kTrackVariablesPair=6+SetChargeAxis+eventplaneaxis;
731 if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
733 if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
735 if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8+SetChargeAxis+eventplaneaxis;
738 // two particle histograms
739 Int_t anaSteps = 1; // analysis steps
740 const char* title = "d^{2}N_{ch}/d#varphid#eta";
742 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
743 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
744 TString* axisTitlePair; // axis titles for track variables
745 axisTitlePair=new TString[kTrackVariablesPair];
747 TString defaultBinningStr;
748 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"
749 "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"
750 "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
751 "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"
752 "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
753 "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
754 "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"
755 "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
757 if(fRequestEventPlane){
758 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))
762 defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
765 defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
768 if(SetChargeAxis==2){
769 defaultBinningStr += "TrigCharge: -2.0,0.0,2.0\n"; // course
770 defaultBinningStr += "AssoCharge: -2.0,0.0,2.0\n"; // course
772 // =========================================================
773 // Customization (adopted from AliUEHistograms)
774 // =========================================================
776 TObjArray* lines = defaultBinningStr.Tokenize("\n");
777 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
779 TString line(lines->At(i)->GetName());
780 TString tag = line(0, line.Index(":")+1);
781 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
782 fBinningString += line + "\n";
784 AliInfo(Form("Using custom binning for %s", tag.Data()));
787 fBinningString += fCustomBinning;
789 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
791 // =========================================================
793 // =========================================================
795 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
796 axisTitlePair[0] = " multiplicity";
798 dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
799 axisTitlePair[1] = "v_{Z} (cm)";
801 dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
802 axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
804 dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
805 axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
807 dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
808 axisTitlePair[4] = "#Delta#eta";
810 dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
811 axisTitlePair[5] = "#Delta#varphi (rad)";
815 if(fRequestEventPlane){
816 dBinsPair[dim_val] = GetBinning(fBinningString, "eventPlane", iBinPair[dim_val]);
817 axisTitlePair[dim_val] = "#varphi - #Psi_{2} (a.u.)";
821 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
822 dBinsPair[dim_val] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val]);
823 axisTitlePair[dim_val] = "TrigCharge";
825 dBinsPair[dim_val+1] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+1]);
826 axisTitlePair[dim_val+1] = "AssoCharge";
829 if(fcontainPIDtrig && !fcontainPIDasso){
830 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
831 axisTitlePair[dim_val] = "PIDTrig";
832 if(SetChargeAxis==2){
833 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
834 axisTitlePair[dim_val+1] = "TrigCharge";
836 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
837 axisTitlePair[dim_val+2] = "AssoCharge";
841 if(!fcontainPIDtrig && fcontainPIDasso){
842 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val]);
843 axisTitlePair[dim_val] = "PIDAsso";
845 if(SetChargeAxis==2){
846 dBinsPair[dim_val+1] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
847 axisTitlePair[dim_val+1] = "TrigCharge";
849 dBinsPair[dim_val+2] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
850 axisTitlePair[dim_val+2] = "AssoCharge";
854 if(fcontainPIDtrig && fcontainPIDasso){
856 dBinsPair[dim_val] = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
857 axisTitlePair[dim_val] = "PIDTrig";
859 dBinsPair[dim_val+1] = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val+1]);
860 axisTitlePair[dim_val+1] = "PIDAsso";
862 if(SetChargeAxis==2){
863 dBinsPair[dim_val+2] = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+2]);
864 axisTitlePair[dim_val+2] = "TrigCharge";
866 dBinsPair[dim_val+3] = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+3]);
867 axisTitlePair[dim_val+3] = "AssoCharge";
872 Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
874 Int_t nPteffbin = -1;
875 Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
878 fminPtTrig=dBinsPair[2][0];
879 fmaxPtTrig=dBinsPair[2][iBinPair[2]];
880 fminPtAsso=dBinsPair[3][0];
881 fmaxPtAsso=dBinsPair[3][iBinPair[3]];
882 fmincentmult=dBinsPair[0][0];
883 fmaxcentmult=dBinsPair[0][iBinPair[0]];
885 //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
886 //fmaxPtComboeff=fmaxPtTrig;
887 //THnSparses for calculation of efficiency
889 if((fAnalysisType =="MCAOD") && ffillefficiency) {
892 effbin[0]=iBinPair[0];
893 effbin[1]=iBinPair[1];
896 Int_t effsteps=5;//for each species type::primMCParticles(0),primRecoTracksMatched(1),allRecoTracksMatched(2),primRecoTracksMatchedPID(3),allRecoTracksMatchedPID(4)
897 for(Int_t jj=0;jj<6;jj++)//PID type binning
899 if(jj==5) effsteps=3;//for unidentified particles
900 Histrename="fTrackHistEfficiency";Histrename+=jj;
901 fTrackHistEfficiency[jj] = new AliTHn(Histrename.Data(), "Tracking efficiency", effsteps, 4, effbin);
902 fTrackHistEfficiency[jj]->SetBinLimits(0, dBinsPair[0]);
903 fTrackHistEfficiency[jj]->SetVarTitle(0, "Centrality");
904 fTrackHistEfficiency[jj]->SetBinLimits(1, dBinsPair[1]);
905 fTrackHistEfficiency[jj]->SetVarTitle(1, "zvtx");
906 fTrackHistEfficiency[jj]->SetBinLimits(2, Pteff);
907 fTrackHistEfficiency[jj]->SetVarTitle(2, "p_{T} (GeV/c)");
908 fTrackHistEfficiency[jj]->SetBinLimits(3, EtaBin);
909 fTrackHistEfficiency[jj]->SetVarTitle(3, "#eta");
910 fOutput->Add(fTrackHistEfficiency[jj]);
914 //AliThns for Correlation plots(data & MC)
916 if(ffilltrigassoUNID)
918 fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
919 for (Int_t j=0; j<kTrackVariablesPair; j++) {
920 fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
921 fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
923 fOutput->Add(fTHnCorrUNID);
925 fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
926 for (Int_t j=0; j<kTrackVariablesPair; j++) {
927 fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
928 fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
930 fOutput->Add(fTHnCorrUNIDmix);
933 if(ffilltrigIDassoID)
935 fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
936 for (Int_t j=0; j<kTrackVariablesPair; j++) {
937 fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
938 fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
940 fOutput->Add(fTHnCorrID);
942 fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
943 for (Int_t j=0; j<kTrackVariablesPair; j++) {
944 fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
945 fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
947 fOutput->Add(fTHnCorrIDmix);
950 if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
952 fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
953 for (Int_t j=0; j<kTrackVariablesPair; j++) {
954 fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
955 fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
957 fOutput->Add(fTHnCorrIDUNID);
960 fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
961 for (Int_t j=0; j<kTrackVariablesPair; j++) {
962 fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
963 fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
965 fOutput->Add(fTHnCorrIDUNIDmix);
970 //ThnSparse for Correlation plots(truth MC)
971 if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
973 fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
974 for (Int_t j=0; j<kTrackVariablesPair; j++) {
975 fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
976 fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
978 fOutput->Add(fCorrelatonTruthPrimary);
981 fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
982 for (Int_t j=0; j<kTrackVariablesPair; j++) {
983 fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
984 fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
986 fOutput->Add(fCorrelatonTruthPrimarymix);
989 //binning for trigger no. counting
992 if(SetChargeAxis==2) ChargeAxis=1;
995 Int_t dims=3+ChargeAxis+eventplaneaxis;
996 if(fcontainPIDtrig) dims=4+ChargeAxis+eventplaneaxis;
997 fBinst= new Int_t[dims];
998 Double_t* dBinsTrig[dims]; // bins for track variables
999 TString* axisTitleTrig; // axis titles for track variables
1000 axisTitleTrig=new TString[dims];
1002 for(Int_t i=0; i<3;i++)
1004 fBinst[i]=iBinPair[i];
1005 dBinsTrig[i]=dBinsPair[i];
1006 axisTitleTrig[i]=axisTitlePair[i];
1008 Int_t dim_val_trig=3;
1009 if(fRequestEventPlane){
1010 fBinst[dim_val_trig]=iBinPair[6];//if fRequestEventPlane=TRUE, dim_val already becomes 7.
1011 dBinsTrig[dim_val_trig]=dBinsPair[6];
1012 axisTitleTrig[dim_val_trig]=axisTitlePair[6];
1016 if(!fcontainPIDtrig && !fcontainPIDasso && ChargeAxis==1){
1017 fBinst[dim_val_trig]=iBinPair[dim_val];
1018 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1019 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1022 if(fcontainPIDtrig && !fcontainPIDasso){
1023 fBinst[dim_val_trig]=iBinPair[dim_val];
1024 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1025 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1027 fBinst[dim_val_trig+1]=iBinPair[dim_val+1];
1028 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+1];
1029 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+1];
1033 if(!fcontainPIDtrig && fcontainPIDasso){
1035 fBinst[dim_val_trig]=iBinPair[dim_val+1];
1036 dBinsTrig[dim_val_trig]=dBinsPair[dim_val+1];
1037 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val+1];
1041 if(fcontainPIDtrig && fcontainPIDasso){
1042 fBinst[dim_val_trig]=iBinPair[dim_val];
1043 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1044 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1046 fBinst[dim_val_trig+1]=iBinPair[dim_val+2];
1047 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+2];
1048 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+2];
1052 //ThSparse for trigger counting(data & reco MC)
1053 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
1055 fTHnTrigcount = new AliTHn("fTHnTrigcount", "fTHnTrigcount", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1056 for(Int_t i=0; i<dims;i++){
1057 fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1058 fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1060 fOutput->Add(fTHnTrigcount);
1063 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
1064 //AliTHns for trigger counting(truth MC)
1065 fTHnTrigcountMCTruthPrim = new AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1066 for(Int_t i=0; i<dims;i++){
1067 fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1068 fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1070 fOutput->Add(fTHnTrigcountMCTruthPrim);
1073 if(fAnalysisType=="MCAOD"){
1074 if(ffillhistQATruth)
1076 MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
1077 fOutputList->Add(MCtruthpt);
1079 MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
1080 fOutputList->Add(MCtrutheta);
1082 MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
1083 fOutputList->Add(MCtruthphi);
1085 MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
1086 fOutputList->Add(MCtruthpionpt);
1088 MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
1089 fOutputList->Add(MCtruthpioneta);
1091 MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
1092 fOutputList->Add(MCtruthpionphi);
1094 MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
1095 fOutputList->Add(MCtruthkaonpt);
1097 MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
1098 fOutputList->Add(MCtruthkaoneta);
1100 MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
1101 fOutputList->Add(MCtruthkaonphi);
1103 MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
1104 fOutputList->Add(MCtruthprotonpt);
1106 MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
1107 fOutputList->Add(MCtruthprotoneta);
1109 MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
1110 fOutputList->Add(MCtruthprotonphi);
1112 fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
1113 fOutputList->Add(fPioncont);
1115 fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
1116 fOutputList->Add(fKaoncont);
1118 fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
1119 fOutputList->Add(fProtoncont);
1121 fUNIDcont=new TH2F("fUNIDcont","fUNIDcont",10,-0.5,9.5,100,0.,10.);
1122 fOutputList->Add(fUNIDcont);
1125 fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1126 fEventno->GetXaxis()->SetTitle("Centrality");
1127 fEventno->GetYaxis()->SetTitle("Z_Vtx");
1128 fOutput->Add(fEventno);
1129 fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1130 fEventnobaryon->GetXaxis()->SetTitle("Centrality");
1131 fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
1132 fOutput->Add(fEventnobaryon);
1133 fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1134 fEventnomeson->GetXaxis()->SetTitle("Centrality");
1135 fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
1136 fOutput->Add(fEventnomeson);
1138 fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
1139 fOutput->Add(fhistJetTrigestimate);
1145 if(fapplyTrigefficiency || fapplyAssoefficiency)
1147 const Int_t nDimt = 4;// cent zvtx pt eta
1148 Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
1149 Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
1150 Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
1152 TString Histrexname;
1153 for(Int_t jj=0;jj<6;jj++)// PID type binning
1155 Histrexname="effcorection";Histrexname+=jj;
1156 effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
1157 effcorection[jj]->Sumw2();
1158 effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
1159 effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
1160 effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
1161 effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
1162 effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
1163 effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1164 effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
1165 effcorection[jj]->GetAxis(3)->SetTitle("#eta");
1166 fOutput->Add(effcorection[jj]);
1168 // TFile *fsifile = new TFile(fefffilename,"READ");
1170 if (TString(fefffilename).BeginsWith("alien:"))
1171 TGrid::Connect("alien:");
1172 TFile *fileT=TFile::Open(fefffilename);
1174 for(Int_t jj=0;jj<6;jj++)//type binning
1176 Nameg="effmap";Nameg+=jj;
1177 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
1178 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
1180 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
1187 //*****************************************************PIDQA histos*****************************************************//
1191 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1192 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1195 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1196 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);
1197 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1198 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1199 fOutputList->Add(fHistoNSigma);
1204 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1205 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1208 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1209 TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
1210 Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1211 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1212 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1213 fOutputList->Add(fHistoNSigma);
1218 if(fPIDType==Bayes){//use bayesianPID
1219 fPIDCombined = new AliPIDCombined();
1220 fPIDCombined->SetDefaultTPCPriors();//****************************************Need to know about it
1222 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1225 TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
1226 Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1227 fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1228 fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
1229 fOutputList->Add(fHistoBayes);
1232 TH2F *fHistoBayesTPC=new TH2F(Form("probBayes_TPC_%d",ipart),
1233 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1234 fHistoBayesTPC->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1235 fHistoBayesTPC->GetYaxis()->SetTitle(Form("Bayes prob TPC %s",kParticleSpeciesName[ipart]));
1236 fOutputList->Add(fHistoBayesTPC);
1238 TH2F *fHistoBayesTOF=new TH2F(Form("probBayes_TOF_%d",ipart),
1239 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1240 fHistoBayesTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1241 fHistoBayesTOF->GetYaxis()->SetTitle(Form("Bayes prob TOF %s",kParticleSpeciesName[ipart]));
1242 fOutputList->Add(fHistoBayesTOF);
1244 TH2F *fHistoBayesTPCTOF=new TH2F(Form("probBayes_TPCTOF_%d",ipart),
1245 Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1246 fHistoBayesTPCTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1247 fHistoBayesTPCTOF->GetYaxis()->SetTitle(Form("Bayes prob TPCTOF %s",kParticleSpeciesName[ipart]));
1248 fOutputList->Add(fHistoBayesTPCTOF);
1252 //nsigma separation power plot
1253 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1256 TH2F *Pi_Ka_sep=new TH2F(Form("Pi_Ka_sep_%d",ipid),
1257 Form("Pi_Ka separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1258 Pi_Ka_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1259 Pi_Ka_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1260 fOutputList->Add(Pi_Ka_sep);
1262 TH2F *Pi_Pr_sep=new TH2F(Form("Pi_Pr_sep_%d",ipid),
1263 Form("Pi_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1264 Pi_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1265 Pi_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1266 fOutputList->Add(Pi_Pr_sep);
1268 TH2F *Ka_Pr_sep=new TH2F(Form("Ka_Pr_sep_%d",ipid),
1269 Form("Ka_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1270 Ka_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1271 Ka_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1272 fOutputList->Add(Ka_Pr_sep);
1276 if(fUseExclusiveNSigma) {
1277 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1278 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1281 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1282 TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
1283 Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1284 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1285 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1286 fOutputList->Add(fHistoNSigma);
1291 if (fAnalysisType == "MCAOD"){
1292 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1293 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1296 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1297 TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
1298 Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1299 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1300 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1301 fOutputList->Add(fHistoNSigma);
1306 for(Int_t idet=0;idet<fNDetectors;idet++){
1307 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1309 if(idet==fTOF)maxy=1.1;
1310 TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
1311 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1312 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1313 fOutputList->Add(fHistoPID);
1316 //PID signal plot, before PID cut
1317 for(Int_t idet=0;idet<fNDetectors;idet++){
1319 if(idet==fTOF)maxy=1.1;
1320 TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
1321 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1322 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1323 fOutputList->Add(fHistoPID);
1326 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
1327 PostData(2, fOutputList);
1329 AliInfo("Finished setting up the Output");
1331 TH1::AddDirectory(oldStatus);
1336 //-------------------------------------------------------------------------------
1337 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1340 if(fAnalysisType == "AOD") {
1344 }//AOD--analysis-----
1346 else if(fAnalysisType == "MCAOD") {
1355 //-------------------------------------------------------------------------
1356 void AliTwoParticlePIDCorr::doMCAODevent()
1358 AliVEvent *event = InputEvent();
1359 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1360 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1362 AliError("Cannot get the AOD event");
1366 // count all events(physics triggered)
1367 fEventCounter->Fill(1);
1369 // get centrality object and check quality(valid for p-Pb and Pb-Pb; coming soon for pp 7 TeV)
1370 Double_t cent_v0=-1.0;
1371 Double_t effcent=1.0;
1372 Double_t refmultReco =0.0;
1373 Double_t gReactionPlane = -1.;
1375 //check the PIDResponse handler
1378 // get mag. field required for twotrack efficiency cut
1380 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1382 //check for TClonesArray(truth track MC information)
1383 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1385 AliFatal("Error: MC particles branch not found!\n");
1389 //check for AliAODMCHeader(truth event MC information)
1390 AliAODMCHeader *header=NULL;
1391 header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1393 printf("MC header branch not found!\n");
1397 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
1398 Float_t zVtxmc =header->GetVtxZ();
1399 if(TMath::Abs(zVtxmc)>fzvtxcut) return;
1401 // For productions with injected signals, figure out above which label to skip particles/tracks
1403 if (fInjectedSignals)
1405 AliGenEventHeader* eventHeader = 0;
1410 AliFatal("fInjectedSignals set but no MC header found");
1412 headers = header->GetNCocktailHeaders();
1413 eventHeader = header->GetCocktailHeader(0);
1417 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
1418 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
1419 AliError("First event header not found. Skipping this event.");
1420 //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1423 skipParticlesAbove = eventHeader->NProduced();
1424 AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
1427 if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL"))
1429 //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
1430 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
1431 refmultReco = GetAcceptedEventMultiplicity(aod,kFALSE);
1432 if(refmultTruth<=0 || refmultReco<=0) return;
1433 cent_v0=refmultTruth;
1436 cent_v0=GetAcceptedEventMultiplicity(aod,kTRUE); //centrality value; 2nd argument has no meaning
1437 if(cent_v0<0.) return;
1440 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1442 //get the event plane in case of PbPb
1443 if(fSampleType=="PbPb"){
1444 if(fRequestEventPlane){
1445 gReactionPlane = GetEventPlane(aod,kTRUE);//get the truth event plane
1446 fHistEventPlaneTruth->Fill(gReactionPlane,cent_v0);
1447 if(gReactionPlane < 0) return;
1451 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)
1453 //TObjArray* tracksMCtruth=0;
1454 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
1455 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
1459 //There is a small difference between zvtx and zVtxmc??????
1460 //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1461 //cout<<"***********************************************zvtx="<<zvtx<<endl;
1463 //now process the truth particles(for both efficiency & correlation function)
1464 Int_t nMCTrack = fArrayMC->GetEntriesFast();
1466 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
1467 { //MC truth track loop starts
1469 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1472 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1476 //consider only charged particles
1477 if(partMC->Charge() == 0) continue;
1479 //consider only primary particles; neglect all secondary particles including from weak decays
1480 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1483 //remove injected signals(primaries above <maxLabel>)
1484 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1487 Bool_t isduplicate=kFALSE;
1488 if (fRemoveDuplicates)
1490 for (Int_t j=iMC+1; j<nMCTrack; ++j)
1491 {//2nd trutuh loop starts
1492 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1494 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1497 if (partMC->GetLabel() == partMC2->GetLabel())
1502 }//2nd truth loop ends
1504 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1506 //give only kinematic cuts at the truth level
1507 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1508 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
1510 if(!partMC) continue;//for safety
1512 //To determine multiplicity in case of PP
1514 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1515 //only physical primary(all/unidentified)
1516 if(ffillhistQATruth)
1518 MCtruthpt->Fill(partMC->Pt());
1519 MCtrutheta->Fill(partMC->Eta());
1520 MCtruthphi->Fill(partMC->Phi());
1523 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1524 Int_t particletypeTruth=-999;
1525 if (TMath::Abs(pdgtruth)==211)
1527 particletypeTruth=SpPion;
1528 if(ffillhistQATruth)
1530 MCtruthpionpt->Fill(partMC->Pt());
1531 MCtruthpioneta->Fill(partMC->Eta());
1532 MCtruthpionphi->Fill(partMC->Phi());
1535 if (TMath::Abs(pdgtruth)==321)
1537 particletypeTruth=SpKaon;
1538 if(ffillhistQATruth)
1540 MCtruthkaonpt->Fill(partMC->Pt());
1541 MCtruthkaoneta->Fill(partMC->Eta());
1542 MCtruthkaonphi->Fill(partMC->Phi());
1545 if(TMath::Abs(pdgtruth)==2212)
1547 particletypeTruth=SpProton;
1548 if(ffillhistQATruth)
1550 MCtruthprotonpt->Fill(partMC->Pt());
1551 MCtruthprotoneta->Fill(partMC->Eta());
1552 MCtruthprotonphi->Fill(partMC->Phi());
1555 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)
1557 // -- Fill THnSparse for efficiency and contamination calculation
1558 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
1560 Double_t primmctruth[4] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
1563 fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
1564 if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for primary truth mesons(3)
1565 if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for primary truth kaons+protons(4)
1566 if (TMath::Abs(pdgtruth)==211) fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
1567 if (TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
1568 if(TMath::Abs(pdgtruth)==2212) fTrackHistEfficiency[2]->Fill(primmctruth,0);//for protons
1571 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1572 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1574 Short_t chargeval=0;
1575 if(partMC->Charge()>0) chargeval=1;
1576 if(partMC->Charge()<0) chargeval=-1;
1577 if(chargeval==0) continue;
1578 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);
1579 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1580 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1581 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1583 }//MC truth track loop ends
1585 //*********************still in event loop
1587 if (fSampleType=="PbPb"){
1588 if (fRandomizeReactionPlane)//only for TRuth MC??
1590 Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
1591 Double_t angle = TMath::TwoPi() * centralityDigits;
1592 AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
1593 ShiftTracks(tracksMCtruth, angle);
1597 Float_t weghtval=1.0;
1599 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1601 //Fill Correlations for MC truth particles(same event)
1602 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1603 Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1606 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200);
1607 if (pool2 && pool2->IsReady())
1608 {//start mixing only when pool->IsReady
1609 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1610 {//proceed only when no. of trigger particles >0 in current event
1611 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
1612 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1613 { //pool event loop start
1614 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1615 if(!bgTracks6) continue;
1616 Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1618 }// pool event loop ends mixing case
1619 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1620 } //if pool->IsReady() condition ends mixing case
1622 //still in main event loop
1625 if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1629 //still in main event loop
1631 if(tracksMCtruth) delete tracksMCtruth;
1633 //now deal with reco tracks
1635 //detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
1636 if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) cent_v0=refmultReco;
1637 effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1639 if(fSampleType=="PbPb"){
1640 if(fRequestEventPlane){
1641 gReactionPlane = GetEventPlane(aod,kFALSE);//get the reconstructed event plane
1642 fHistEventPlaneReco->Fill(gReactionPlane,cent_v0);
1643 if(gReactionPlane < 0) return;
1646 //TObjArray* tracksUNID=0;
1647 TObjArray* tracksUNID = new TObjArray;
1648 tracksUNID->SetOwner(kTRUE);
1650 //TObjArray* tracksID=0;
1651 TObjArray* tracksID = new TObjArray;
1652 tracksID->SetOwner(kTRUE);
1655 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1658 Double_t trackscount=0.0;
1659 // loop over reconstructed tracks
1660 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1661 { // reconstructed track loop starts
1662 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1663 if (!track) continue;
1664 //get the corresponding MC track at the truth level (doing reco matching)
1665 AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
1666 if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
1668 //remove injected signals
1669 if(fInjectedSignals)
1671 AliAODMCParticle* mother = recomatched;
1673 while (!mother->IsPhysicalPrimary())
1674 {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
1675 if (mother->GetMother() < 0)
1681 mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1687 Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
1690 if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
1691 }//remove injected signals
1693 if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
1695 Bool_t isduplicate2=kFALSE;
1696 if (fRemoveDuplicates)
1698 for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
1700 AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
1701 if (!track2) continue;
1702 AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
1703 if(!recomatched2) continue;
1705 if (track->GetLabel() == track2->GetLabel())
1712 if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
1714 fHistQA[11]->Fill(track->GetTPCNcls());
1715 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
1717 if(tracktype==0) continue;
1718 if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
1720 if(!track) continue;//for safety
1721 //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
1724 //check for eta , phi holes
1725 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1726 fphiSpectraasso->Fill(track->Phi(),track->Pt());
1728 Int_t particletypeMC=-9999;
1730 //tag all particles as unidentified
1731 particletypeMC=unidentified;
1733 Float_t effmatrix=1.;
1735 // -- Fill THnSparse for efficiency calculation
1736 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
1737 //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)
1739 //Clone & Reduce track list(TObjArray) for unidentified particles
1740 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1742 Short_t chargeval=0;
1743 if(track->Charge()>0) chargeval=1;
1744 if(track->Charge()<0) chargeval=-1;
1745 if(chargeval==0) continue;
1746 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
1747 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
1748 LRCParticlePID* copy = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
1749 copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
1750 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
1753 //get the pdg code of the corresponding truth particle
1754 Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
1756 Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
1757 if(ffillefficiency) {
1758 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
1759 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
1760 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
1761 if(TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions
1762 if(TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
1763 if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
1765 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
1766 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
1767 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
1768 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
1769 if( TMath::Abs(pdgCode)==211) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions
1770 if( TMath::Abs(pdgCode)==321) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
1771 if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
1775 //now start the particle identification process:)
1777 Float_t dEdx = track->GetTPCsignal();
1778 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1780 if(HasTOFPID(track))
1782 Double_t beta = GetBeta(track);
1783 fHistoTOFbeta->Fill(track->Pt(), beta);
1786 //do track identification(nsigma method)
1787 particletypeMC=GetParticle(track,fFIllPIDQAHistos);//******************************problem is here
1789 switch(TMath::Abs(pdgCode)){
1791 if(fFIllPIDQAHistos){
1792 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1793 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
1794 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
1795 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
1800 if(fFIllPIDQAHistos){
1801 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1802 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
1803 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpKaon,ipid));
1804 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
1809 if(fFIllPIDQAHistos){
1810 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1811 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
1812 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpPion,ipid));
1813 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
1820 //2-d TPCTOF map(for each Pt interval)
1821 if(HasTOFPID(track)){
1822 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
1823 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
1824 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
1827 //Pt, Eta , Phi distribution of the reconstructed identified particles
1830 if (particletypeMC==SpPion)
1832 fPionPt->Fill(track->Pt());
1833 fPionEta->Fill(track->Eta());
1834 fPionPhi->Fill(track->Phi());
1836 if (particletypeMC==SpKaon)
1838 fKaonPt->Fill(track->Pt());
1839 fKaonEta->Fill(track->Eta());
1840 fKaonPhi->Fill(track->Phi());
1842 if (particletypeMC==SpProton)
1844 fProtonPt->Fill(track->Pt());
1845 fProtonEta->Fill(track->Eta());
1846 fProtonPhi->Fill(track->Phi());
1850 //for misidentification fraction calculation(do it with tuneonPID)
1851 if(particletypeMC==SpPion )
1853 if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
1854 if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
1855 if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
1856 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
1858 if(particletypeMC==SpKaon )
1860 if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
1861 if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
1862 if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
1863 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
1865 if(particletypeMC==SpProton )
1867 if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
1868 if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
1869 if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
1870 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
1872 if(particletypeMC==SpUndefined )
1874 if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
1875 if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
1876 if(TMath::Abs(pdgCode)==2212) fUNIDcont->Fill(5.,track->Pt());
1877 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fUNIDcont->Fill(7.,track->Pt());
1880 if(particletypeMC==SpUndefined) continue;
1883 //fill tracking efficiency
1886 if(particletypeMC==SpPion || particletypeMC==SpKaon)
1888 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) {
1889 fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
1890 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
1893 if(particletypeMC==SpKaon || particletypeMC==SpProton)
1895 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) {
1896 fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
1897 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
1900 if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) {
1901 fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
1902 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
1904 if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
1905 fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
1906 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
1908 if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
1909 fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
1910 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
1914 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1916 Short_t chargeval=0;
1917 if(track->Charge()>0) chargeval=1;
1918 if(track->Charge()<0) chargeval=-1;
1919 if(chargeval==0) continue;
1920 if (fapplyTrigefficiency || fapplyAssoefficiency)
1921 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
1922 LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
1923 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1924 tracksID->Add(copy1);
1926 }// if(tracktype==1) condition structure ands
1928 }//reco track loop ends
1930 //*************************************************************still in event loop
1935 //fill the centrality/multiplicity distribution of the selected events
1936 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
1938 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??????
1940 //count selected events having centrality betn 0-100%
1941 fEventCounter->Fill(11);
1943 //***************************************event no. counting
1944 Bool_t isbaryontrig=kFALSE;
1945 Bool_t ismesontrig=kFALSE;
1946 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
1948 if(tracksID && tracksID->GetEntriesFast()>0)
1950 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
1951 { //trigger loop starts
1952 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
1954 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
1955 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
1956 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
1957 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
1959 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
1960 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
1963 //same event delte-eta, delta-phi plot
1964 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
1965 {//same event calculation starts
1966 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
1967 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)
1970 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
1971 {//same event calculation starts
1972 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)
1973 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
1976 //still in main event loop
1978 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
1979 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
1980 if (pool && pool->IsReady())
1981 {//start mixing only when pool->IsReady
1982 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
1983 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
1984 { //pool event loop start
1985 TObjArray* bgTracks = pool->GetEvent(jMix);
1986 if(!bgTracks) continue;
1987 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
1988 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
1989 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
1990 Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
1991 }// pool event loop ends mixing case
1993 } //if pool->IsReady() condition ends mixing case
1996 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
1998 }//mixing with unidentified particles
2000 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2001 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
2002 if (pool1 && pool1->IsReady())
2003 {//start mixing only when pool->IsReady
2004 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2005 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2006 { //pool event loop start
2007 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2008 if(!bgTracks2) continue;
2009 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2010 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2011 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2012 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2014 }// pool event loop ends mixing case
2015 } //if pool1->IsReady() condition ends mixing case
2019 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2021 }//mixing with identified particles
2023 //no. of events analyzed
2024 fEventCounter->Fill(13);
2027 if(tracksUNID) delete tracksUNID;
2029 if(tracksID) delete tracksID;
2032 PostData(1, fOutput);
2035 //________________________________________________________________________
2036 void AliTwoParticlePIDCorr::doAODevent()
2040 AliVEvent *event = InputEvent();
2041 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
2042 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
2044 AliError("Cannot get the AOD event");
2049 fEventCounter->Fill(1);
2051 if (!fPID) return;//this should be available with each event even if we don't do PID selection
2053 Double_t cent_v0= -999.;
2054 Double_t effcent=1.0;
2055 Double_t gReactionPlane = -1.;
2057 Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
2060 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
2061 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
2064 // check event cuts and fill event histograms and return the centrality or reference multiplicity value
2065 if((cent_v0 = GetAcceptedEventMultiplicity(aod,kFALSE)) < 0){
2069 //get the event plane in case of PbPb
2070 if(fSampleType=="PbPb"){
2071 if(fRequestEventPlane){
2072 gReactionPlane = GetEventPlane(aod,kFALSE);
2073 fHistEventPlaneReco->Fill(gReactionPlane,cent_v0);
2074 if(gReactionPlane < 0) return;
2078 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
2079 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
2081 TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
2082 tracksID->SetOwner(kTRUE); // IMPORTANT!
2087 Bool_t fTrigPtmin1=kFALSE;
2088 Bool_t fTrigPtmin2=kFALSE;
2089 Bool_t fTrigPtJet=kFALSE;
2091 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2092 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2093 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2094 if (!track) continue;
2095 fHistQA[11]->Fill(track->GetTPCNcls());
2096 Int_t particletype=-9999;//required for PID filling
2097 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2098 if(tracktype!=1) continue;
2100 if(!track) continue;//for safety
2102 //check for eta , phi holes
2103 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2104 fphiSpectraasso->Fill(track->Phi(),track->Pt());
2108 //if no applyefficiency , set the eff factor=1.0
2109 Float_t effmatrix=1.0;
2111 //tag all particles as unidentified that passed filterbit & kinematic cuts
2112 particletype=unidentified;
2114 //To count the no. of tracks having an accepted track in a certain PT(e.g. Jet Pt) range
2115 if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
2116 if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
2117 if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
2120 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
2123 //to reduce memory consumption in pool
2124 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2126 //Clone & Reduce track list(TObjArray) for unidentified particles
2127 Short_t chargeval=0;
2128 if(track->Charge()>0) chargeval=1;
2129 if(track->Charge()<0) chargeval=-1;
2130 if(chargeval==0) continue;
2131 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2132 effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
2133 LRCParticlePID* copy = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2134 copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2135 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2138 //now start the particle identificaion process:)
2140 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
2142 Float_t dEdx = track->GetTPCsignal();
2143 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2145 //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)
2146 if(HasTOFPID(track))
2148 Double_t beta = GetBeta(track);
2149 fHistoTOFbeta->Fill(track->Pt(), beta);
2153 //track identification(using nsigma method)
2154 particletype=GetParticle(track,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
2157 //2-d TPCTOF map(for each Pt interval)
2158 if(HasTOFPID(track)){
2159 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2160 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2161 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2164 //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!!!!!
2165 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)!!!!!!!!!!!
2167 //Pt, Eta , Phi distribution of the reconstructed identified particles
2170 if (particletype==SpPion)
2172 fPionPt->Fill(track->Pt());
2173 fPionEta->Fill(track->Eta());
2174 fPionPhi->Fill(track->Phi());
2176 if (particletype==SpKaon)
2178 fKaonPt->Fill(track->Pt());
2179 fKaonEta->Fill(track->Eta());
2180 fKaonPhi->Fill(track->Phi());
2182 if (particletype==SpProton)
2184 fProtonPt->Fill(track->Pt());
2185 fProtonEta->Fill(track->Eta());
2186 fProtonPhi->Fill(track->Phi());
2190 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2192 Short_t chargeval=0;
2193 if(track->Charge()>0) chargeval=1;
2194 if(track->Charge()<0) chargeval=-1;
2195 if(chargeval==0) continue;
2196 if (fapplyTrigefficiency || fapplyAssoefficiency)
2197 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
2198 LRCParticlePID* copy1 = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2199 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2200 tracksID->Add(copy1);
2202 } //track loop ends but still in event loop
2204 if(trackscount<1.0){
2205 if(tracksUNID) delete tracksUNID;
2206 if(tracksID) delete tracksID;
2210 if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
2211 if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
2212 if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
2214 Float_t weightval=1.0;
2218 //fill the centrality/multiplicity distribution of the selected events
2219 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2221 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??????
2223 //count selected events having centrality betn 0-100%
2224 fEventCounter->Fill(11);
2226 //***************************************event no. counting
2227 Bool_t isbaryontrig=kFALSE;
2228 Bool_t ismesontrig=kFALSE;
2229 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2231 if(tracksID && tracksID->GetEntriesFast()>0)
2233 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2234 { //trigger loop starts
2235 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2237 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2238 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2239 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2240 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2242 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2243 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2247 //same event delta-eta-deltaphi plot
2249 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2250 {//same event calculation starts
2251 if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2252 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)
2255 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2256 {//same event calculation starts
2257 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)
2258 if(ffilltrigIDassoID) Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2261 //still in main event loop
2265 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2266 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
2267 if (pool && pool->IsReady())
2268 {//start mixing only when pool->IsReady
2269 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2270 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2271 { //pool event loop start
2272 TObjArray* bgTracks = pool->GetEvent(jMix);
2273 if(!bgTracks) continue;
2274 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2275 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2276 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2277 Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2278 }// pool event loop ends mixing case
2279 } //if pool->IsReady() condition ends mixing case
2282 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2284 }//mixing with unidentified particles
2287 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2288 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
2289 if (pool1 && pool1->IsReady())
2290 {//start mixing only when pool->IsReady
2291 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2292 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2293 { //pool event loop start
2294 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2295 if(!bgTracks2) continue;
2296 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2297 Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2298 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2299 Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2301 }// pool event loop ends mixing case
2302 } //if pool1->IsReady() condition ends mixing case
2306 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2308 }//mixing with identified particles
2311 //no. of events analyzed
2312 fEventCounter->Fill(13);
2314 //still in main event loop
2317 if(tracksUNID) delete tracksUNID;
2319 if(tracksID) delete tracksID;
2322 PostData(1, fOutput);
2324 } // *************************event loop ends******************************************//_______________________________________________________________________
2325 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
2327 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
2329 TObjArray* tracksClone = new TObjArray;
2330 tracksClone->SetOwner(kTRUE);
2332 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2334 LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
2335 LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval());
2336 copy100->SetUniqueID(particle->GetUniqueID());
2337 tracksClone->Add(copy100);
2343 //--------------------------------------------------------------------------------
2344 void AliTwoParticlePIDCorr::Fillcorrelation(Double_t gReactionPlane,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)
2347 //before calling this function check that either trackstrig & tracksasso are available
2349 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
2350 TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
2351 TArrayF eta(input->GetEntriesFast());
2352 for (Int_t i=0; i<input->GetEntriesFast(); i++)
2353 eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
2356 Int_t jmax=trackstrig->GetEntriesFast();
2358 jmax=tracksasso->GetEntriesFast();
2360 // identify K, Lambda candidates and flag those particles
2361 // a TObject bit is used for this
2362 const UInt_t kResonanceDaughterFlag = 1 << 14;
2363 if (fRejectResonanceDaughters > 0)
2365 Double_t resonanceMass = -1;
2366 Double_t massDaughter1 = -1;
2367 Double_t massDaughter2 = -1;
2368 const Double_t interval = 0.02;
2369 switch (fRejectResonanceDaughters)
2371 case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
2372 case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
2373 case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
2374 default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
2377 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2378 trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2379 for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
2380 tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2382 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2384 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2385 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
2387 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2389 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
2390 if (trig->IsEqual(asso)) continue;
2392 if (trig->Charge() * asso->Charge() > 0) continue;
2394 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2396 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
2398 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2400 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
2402 trig->SetBit(kResonanceDaughterFlag);
2403 asso->SetBit(kResonanceDaughterFlag);
2405 // Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
2412 //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
2414 Float_t TriggerPtMin=fminPtTrig;
2415 Float_t TriggerPtMax=fmaxPtTrig;
2416 Int_t HighestPtTriggerIndx=-99999;
2417 TH1* triggerWeighting = 0;
2419 if(fSelectHighestPtTrig || fWeightPerEvent)//**************add this data member to the constructor
2421 if (fWeightPerEvent)
2424 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID) axis = fTHnTrigcount->GetGrid(0)->GetGrid()->GetAxis(2);
2425 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) axis = fTHnTrigcountMCTruthPrim->GetGrid(0)->GetGrid()->GetAxis(2);
2426 triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
2428 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2429 { //trigger loop starts(highest Pt trigger selection)
2430 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2432 Float_t trigpt=trig->Pt();
2433 //to avoid overflow qnd underflow
2434 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2435 Int_t particlepidtrig=trig->getparticle();
2436 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2438 Float_t trigeta=trig->Eta();
2440 // some optimization
2441 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2444 if (fOnlyOneEtaSide != 0)
2446 if (fOnlyOneEtaSide * trigeta < 0)
2449 if (fTriggerSelectCharge != 0)
2450 if (trig->Charge() * fTriggerSelectCharge < 0)
2453 if (fRejectResonanceDaughters > 0)
2454 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2456 if(fSelectHighestPtTrig){
2457 if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
2459 HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
2460 TriggerPtMin=trigpt;
2464 if (fWeightPerEvent) triggerWeighting->Fill(trigpt);
2466 }//trigger loop ends(highest Pt trigger selection)
2468 }//******************(fSelectHighestPtTrig || fWeightPerEvent) condition ends
2471 //two particle correlation filling
2472 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2473 { //trigger loop starts
2474 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2476 Float_t trigpt=trig->Pt();
2477 //to avoid overflow qnd underflow
2478 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2479 Int_t particlepidtrig=trig->getparticle();
2480 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2482 Float_t trigeta=trig->Eta();
2484 // some optimization
2485 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2488 if (fOnlyOneEtaSide != 0)
2490 if (fOnlyOneEtaSide * trigeta < 0)
2493 if (fTriggerSelectCharge != 0)
2494 if (trig->Charge() * fTriggerSelectCharge < 0)
2497 if (fRejectResonanceDaughters > 0)
2498 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2500 if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
2501 if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
2504 Float_t trigphi=trig->Phi();
2505 Float_t trackefftrig=1.0;
2506 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
2508 // Event plane (determine psi bin)
2509 Double_t gPsiMinusPhi = 0.;
2510 Double_t gPsiMinusPhiBin = -10.;
2511 if(fRequestEventPlane){
2512 gPsiMinusPhi = TMath::Abs(trigphi - gReactionPlane);
2514 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2515 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2516 gPsiMinusPhiBin = 0.0;
2518 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
2519 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
2520 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
2521 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
2522 gPsiMinusPhiBin = 1.0;
2524 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
2525 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
2526 gPsiMinusPhiBin = 2.0;
2529 gPsiMinusPhiBin = 3.0;
2531 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
2534 //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
2537 Int_t eventplaneAxis=0;
2538 if(fRequestEventPlane) eventplaneAxis=1;
2539 if(fcontainPIDtrig && SetChargeAxis==0) dim=4+eventplaneAxis;
2540 if(!fcontainPIDtrig && SetChargeAxis==2) dim=4+eventplaneAxis;
2541 if(fcontainPIDtrig && SetChargeAxis==2) dim=5+eventplaneAxis;
2542 trigval= new Double_t[dim];
2545 trigval[2] = trigpt;
2547 if(fRequestEventPlane){
2548 trigval[3] = gPsiMinusPhiBin;
2549 if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = particlepidtrig;
2550 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
2551 if(fcontainPIDtrig && SetChargeAxis==2) {
2552 trigval[4] = particlepidtrig;
2553 trigval[5] = trig->Charge();
2557 if(!fRequestEventPlane){
2558 if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = particlepidtrig;
2559 if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
2560 if(fcontainPIDtrig && SetChargeAxis==2) {
2561 trigval[3] = particlepidtrig;
2562 trigval[4] = trig->Charge();
2568 if (fWeightPerEvent)
2570 // leads effectively to a filling of one entry per filled trigger particle pT bin
2571 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(trigval[2]);
2572 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
2573 trackefftrig *= triggerWeighting->GetBinContent(weightBin);
2577 //trigger particle counting for both same and mixed event case;;;;;step=0->same event case;;;;;step=1->mixed event case;;;;;;
2578 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
2579 if(fillup=="trigassoUNID" ) {
2580 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2581 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2584 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
2585 if(fillup=="trigassoUNID" )
2587 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2588 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2591 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
2592 if(fillup=="trigUNIDassoID")
2594 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2595 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2598 //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
2599 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
2600 if(fillup=="trigIDassoID")
2602 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2603 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2606 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
2607 if(fillup=="trigIDassoUNID" )
2609 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2610 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2613 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
2614 if(fillup=="trigIDassoID")
2616 if(mixcase==kFALSE) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2617 if(mixcase==kTRUE && firstTime) fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2621 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!!!!
2622 if(mixcase==kFALSE) fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig);
2623 if(mixcase==kTRUE && firstTime) fTHnTrigcountMCTruthPrim->Fill(trigval,1,1.0/trackefftrig);
2626 //asso loop starts within trigger loop
2627 for(Int_t j=0;j<jmax;j++)
2629 LRCParticlePID *asso=0;
2631 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
2633 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2637 //to avoid overflow and underflow
2638 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
2640 //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
2642 if(!tracksasso && j==i) continue;
2644 // 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)
2645 // if (tracksasso && trig->IsEqual(asso)) continue;
2647 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
2650 if (asso->Pt() >= trig->Pt()) continue;
2652 Int_t particlepidasso=asso->getparticle();
2653 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
2656 if (fAssociatedSelectCharge != 0)
2657 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
2659 if (fSelectCharge > 0)
2662 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
2666 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
2672 if (trigeta < 0 && asso->Eta() < trigeta)
2674 if (trigeta > 0 && asso->Eta() > trigeta)
2678 if (fRejectResonanceDaughters > 0)
2679 if (asso->TestBit(kResonanceDaughterFlag))
2681 // Printf("Skipped j=%d", j);
2686 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
2688 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2692 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2694 fControlConvResoncances->Fill(0.0, mass);
2696 if (mass < 0.04*0.04)
2702 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2704 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
2706 const Float_t kK0smass = 0.4976;
2708 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
2710 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
2712 fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
2714 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
2720 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2722 Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
2723 Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2725 const Float_t kLambdaMass = 1.115;
2727 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
2729 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
2731 fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
2733 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2736 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
2738 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2740 fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
2742 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2747 if (twoTrackEfficiencyCut)
2749 // the variables & cuthave been developed by the HBT group
2750 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
2751 Float_t phi1 = trig->Phi();
2752 Float_t pt1 = trig->Pt();
2753 Float_t charge1 = trig->Charge();
2754 Float_t phi2 = asso->Phi();
2755 Float_t pt2 = asso->Pt();
2756 Float_t charge2 = asso->Charge();
2758 Float_t deta= trigeta - eta[j];
2761 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
2764 // check first boundaries to see if is worth to loop and find the minimum
2765 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
2766 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
2768 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
2770 Float_t dphistarminabs = 1e5;
2771 Float_t dphistarmin = 1e5;
2773 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
2775 for (Double_t rad=0.8; rad<2.51; rad+=0.01)
2777 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
2779 Float_t dphistarabs = TMath::Abs(dphistar);
2781 if (dphistarabs < dphistarminabs)
2783 dphistarmin = dphistar;
2784 dphistarminabs = dphistarabs;
2788 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
2790 // 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);
2793 //fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
2799 Float_t weightperevent=weight;
2800 Float_t trackeffasso=1.0;
2801 if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
2802 //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
2803 Float_t deleta=trigeta-eta[j];
2804 Float_t delphi=PhiRange(trigphi-asso->Phi());
2806 //here get the two particle efficiency correction factor
2807 Float_t effweight=trackefftrig*trackeffasso*weightperevent;
2808 // if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
2810 Int_t dimused=kTrackVariablesPair+eventplaneAxis;
2811 vars= new Double_t[dimused];
2820 if(fRequestEventPlane)
2822 vars[6]=gPsiMinusPhiBin;
2826 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
2827 vars[dimension]=trig->Charge();
2828 vars[dimension+1]=asso->Charge();
2830 if(fcontainPIDtrig && !fcontainPIDasso){
2831 vars[dimension]=particlepidtrig;
2832 if(SetChargeAxis==2){
2833 vars[dimension+1]=trig->Charge();
2834 vars[dimension+2]=asso->Charge();
2837 if(!fcontainPIDtrig && fcontainPIDasso){
2838 vars[dimension]=particlepidasso;
2839 if(SetChargeAxis==2){
2840 vars[dimension+1]=trig->Charge();
2841 vars[dimension+2]=asso->Charge();
2844 if(fcontainPIDtrig && fcontainPIDasso){
2845 vars[dimension]=particlepidtrig;
2846 vars[dimension+1]=particlepidasso;
2847 if(SetChargeAxis==2){
2848 vars[dimension+2]=trig->Charge();
2849 vars[dimension+3]=asso->Charge();
2853 if (fWeightPerEvent)
2855 Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
2856 // Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
2857 effweight *= triggerWeighting->GetBinContent(weightBin);
2861 //Fill Correlation ThnSparses
2862 if(fillup=="trigassoUNID")
2864 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,0,1.0/effweight);
2865 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
2867 if(fillup=="trigIDassoID")
2869 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,0,1.0/effweight);
2870 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
2872 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
2873 {//in this case effweight should be 1 always
2874 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
2875 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
2877 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
2879 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
2880 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
2886 }//trigger loop ends
2888 if (triggerWeighting)
2890 delete triggerWeighting;
2891 triggerWeighting = 0;
2895 //--------------------------------------------------------------------------------
2896 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
2898 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
2900 Float_t effvalue=1.;
2902 if(parpid==unidentified)
2904 effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
2905 effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
2906 effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
2907 effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
2908 effvalue=effcorection[5]->GetBinContent(effVars);
2910 if(parpid==SpPion || parpid==SpKaon)
2912 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2914 effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
2915 effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
2916 effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
2917 effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
2918 effvalue=effcorection[3]->GetBinContent(effVars);
2923 effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
2924 effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
2925 effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
2926 effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
2927 effvalue=effcorection[0]->GetBinContent(effVars);
2932 effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
2933 effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
2934 effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
2935 effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
2936 effvalue=effcorection[1]->GetBinContent(effVars);
2941 if(parpid==SpProton)
2943 effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
2944 effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
2945 effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
2946 effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
2947 effvalue=effcorection[2]->GetBinContent(effVars);
2950 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2952 if(parpid==SpProton || parpid==SpKaon)
2954 effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
2955 effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
2956 effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
2957 effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
2958 effvalue=effcorection[4]->GetBinContent(effVars);
2961 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
2962 if(effvalue==0.) effvalue=1.;
2967 //---------------------------------------------------------------------------------
2969 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield, Bool_t fill)
2972 if(!track) return 0;
2973 Bool_t trackOK = track->TestFilterBit(fFilterBit);
2974 if(!trackOK) return 0;
2975 if (fTrackStatus != 0 && !CheckTrack(track)) return 0;
2976 //select only primary traks(for data & reco MC tracks)
2977 if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
2978 if(track->Charge()==0) return 0;
2979 if (fill) fHistQA[12]->Fill(track->GetTPCNcls());
2982 dz = track->ZAtDCA();
2983 if (fill) fHistQA[6]->Fill(dxy);
2984 if (fill) fHistQA[7]->Fill(dz);
2985 Float_t chi2ndf = track->Chi2perNDF();
2986 if (fill) fHistQA[13]->Fill(chi2ndf);
2987 // Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
2988 Float_t nCrossedRowsTPC = track->GetTPCNCrossedRows();
2989 if (fill) fHistQA[14]->Fill(nCrossedRowsTPC);
2990 //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
2991 if (track->GetTPCNclsF()>0) {
2992 Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
2993 if (fill) fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
2996 Float_t pt=track->Pt();
2997 if(pt< fminPt || pt> fmaxPt) return 0;
2998 if(TMath::Abs(track->Eta())> fmaxeta) return 0;
2999 if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
3000 //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required???
3002 if (fdcacut && fDCAXYCut)
3009 AliAODTrack* clone =(AliAODTrack*) track->Clone();
3010 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
3015 // Printf("%f", ((AliAODTrack*)part)->DCA());
3016 // Printf("%f", pos[0]);
3017 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
3021 if (fSharedClusterCut >= 0)
3023 Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
3024 if (frac > fSharedClusterCut)
3027 if (fill) fHistQA[8]->Fill(pt);
3028 if (fill) fHistQA[9]->Fill(track->Eta());
3029 if (fill) fHistQA[10]->Fill(track->Phi());
3032 //________________________________________________________________________________
3033 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
3035 //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 )
3036 Float_t pt=track->Pt();
3038 //plot the separation power
3040 Double_t bethe[AliPID::kSPECIES]={0.};
3041 Double_t sigma_TPC[AliPID::kSPECIES]={0.};
3043 Double_t Pi_Ka_sep[NSigmaPIDType+1]={0.};
3044 Double_t Pi_Pr_sep[NSigmaPIDType+1]={0.};
3045 Double_t Ka_Pr_sep[NSigmaPIDType+1]={0.};
3048 Double_t ptpc = track->GetTPCmomentum();
3049 Int_t dEdxN = track->GetTPCsignalN();
3050 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
3051 bethe[ipart] = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
3052 //Double_t diff = dEdx - bethe;
3053 sigma_TPC[ipart] = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
3054 //nSigma[ipart] = diff / sigma;
3056 Pi_Ka_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kKaon])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kKaon])/2.0);
3057 Pi_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kProton])/2.0);
3058 Ka_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kKaon]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kKaon]+sigma_TPC[AliPID::kProton])/2.0);
3061 Double_t sigma_TOF[AliPID::kSPECIES]={0.};
3063 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3065 Double_t timei[AliPID::kSPECIES];
3066 track->GetIntegratedTimes(timei);
3067 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { sigma_TOF[ipart]= fPID->GetTOFResponse().GetExpectedSigma(track->P(), timei[ipart], AliPID::ParticleMass(ipart));}
3068 Pi_Ka_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kKaon])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kKaon])/2.0);
3069 Pi_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kProton])/2.0);
3070 Ka_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kKaon]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kKaon]+sigma_TOF[AliPID::kProton])/2.0);
3072 Pi_Ka_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Ka_sep[NSigmaTPC]*Pi_Ka_sep[NSigmaTPC]+Pi_Ka_sep[NSigmaTOF]*Pi_Ka_sep[NSigmaTOF]);
3073 Pi_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Pr_sep[NSigmaTPC]*Pi_Pr_sep[NSigmaTPC]+Pi_Pr_sep[NSigmaTOF]*Pi_Pr_sep[NSigmaTOF]);
3074 Ka_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Ka_Pr_sep[NSigmaTPC]*Ka_Pr_sep[NSigmaTPC]+Ka_Pr_sep[NSigmaTOF]*Ka_Pr_sep[NSigmaTOF]);
3078 //fill separation power histograms
3079 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3081 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3082 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3083 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3084 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3085 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3086 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3088 if(HasTOFPID(track) && pt>fPtTOFPIDmin && ipid!=0){
3089 TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3090 h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3091 TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3092 h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3093 TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3094 h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3099 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
3100 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
3101 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
3102 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
3104 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
3105 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
3107 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3110 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
3111 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
3112 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
3114 nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
3115 nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
3116 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
3122 // if TOF is missing and below fPtTOFPID only the TPC information is used
3123 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
3124 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
3125 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
3129 //set data member fnsigmas
3130 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
3131 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
3132 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
3134 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
3135 fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
3136 fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
3137 fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
3139 //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)
3140 fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
3141 fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
3142 fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
3145 //Fill NSigma SeparationPlot
3146 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3147 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3148 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3149 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
3150 h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
3156 //----------------------------------------------------------------------------
3157 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
3159 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3160 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
3161 //get the identity of the particle with the minimum Nsigma
3162 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3165 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3166 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
3167 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
3170 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3171 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
3172 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
3174 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3175 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3176 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3177 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3179 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3180 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3181 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3182 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3187 if(fdiffPIDcutvalues){
3188 if(track->Pt()<=4) fNSigmaPID=fPIDCutval1;
3189 if(track->Pt()>4 && track->Pt()<=6) fNSigmaPID=fPIDCutval2;
3190 if(track->Pt()>6 && track->Pt()<=8) fNSigmaPID=fPIDCutval3;
3191 if(track->Pt()>8) fNSigmaPID=fPIDCutval4;
3194 // guess the particle based on the smaller nsigma (within fNSigmaPID)
3195 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
3197 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
3198 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)
3200 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3201 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3202 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
3203 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
3208 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
3210 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3211 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3212 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
3213 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
3218 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
3220 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3221 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3222 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
3223 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
3229 // else, return undefined
3235 //------------------------------------------------------------------------------------------
3236 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
3237 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3239 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
3241 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
3243 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
3246 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
3248 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3251 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3252 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
3253 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
3256 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3257 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
3258 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
3260 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3261 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3262 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3263 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3265 case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3266 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3267 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
3268 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
3272 // Actually the tracks in the overlapping region(in TPC-TOF nSigma plane) will be ignored
3274 if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
3275 if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
3276 if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
3281 //fill NSigma distr for double counting
3282 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3283 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
3284 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3285 if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)))continue;//not filling TOF and combined if no TOF PID
3286 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
3287 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
3294 return fHasDoubleCounting;
3297 //////////////////////////////////////////////////////////////////////////////////////////////////
3299 Bool_t* AliTwoParticlePIDCorr::GetAllCompatibleIdentitiesNSigma(AliAODTrack * trk,Bool_t FIllQAHistos){
3300 //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
3301 Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
3302 IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
3306 //////////////////////////////////////////////////////////////////////////////////////////////////
3308 UInt_t AliTwoParticlePIDCorr::CalcPIDCombined(AliAODTrack *track, Int_t detMask, Double_t* prob) const{
3310 // Bayesian PID calculation
3312 for(Int_t i=0;i<AliPID::kSPECIES;i++)
3316 fPIDCombined->SetDetectorMask(detMask);
3318 return fPIDCombined->ComputeProbabilities((AliAODTrack*)track, fPID, prob);
3321 //////////////////////////////////////////////////////////////////////////////////////////////////
3323 Int_t AliTwoParticlePIDCorr::GetIDBayes(AliAODTrack * trk, Bool_t FIllQAHistos){
3325 Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
3328 //Filling of Probability histos
3329 Double_t probTPC[AliPID::kSPECIES]={0.};
3330 Double_t probTOF[AliPID::kSPECIES]={0.};
3331 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
3333 UInt_t detUsedTPC = 0;
3334 UInt_t detUsedTOF = 0;
3335 UInt_t detUsedTPCTOF = 0;
3337 //get the TPC probability
3338 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
3339 detUsedTPC = fPIDCombined->ComputeProbabilities(trk, fPID, probTPC);
3340 if(detUsedTPC == AliPIDResponse::kDetTPC)
3342 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3344 TH2F *h=GetHistogram2D(Form("probBayes_TPC_%d",ipart));
3345 if(ipart==0) h->Fill(trk->Pt(),probTPC[AliPID::kPion]);
3346 if(ipart==1) h->Fill(trk->Pt(),probTPC[AliPID::kKaon]);
3347 if(ipart==2) h->Fill(trk->Pt(),probTPC[AliPID::kProton]);
3351 //get the TOF probability
3352 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
3353 detUsedTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTOF);
3354 if(detUsedTOF == AliPIDResponse::kDetTOF)
3356 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3357 TH2F *h=GetHistogram2D(Form("probBayes_TOF_%d",ipart));
3358 if(ipart==0) h->Fill(trk->Pt(),probTOF[AliPID::kPion]);
3359 if(ipart==1) h->Fill(trk->Pt(),probTOF[AliPID::kKaon]);
3360 if(ipart==2) h->Fill(trk->Pt(),probTOF[AliPID::kProton]);
3364 //get the TPC-TOF probability
3365 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
3366 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTPCTOF);
3367 if(detUsedTPCTOF == (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))
3369 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3370 TH2F *h=GetHistogram2D(Form("probBayes_TPCTOF_%d",ipart));
3371 if(ipart==0) h->Fill(trk->Pt(),probTPCTOF[AliPID::kPion]);
3372 if(ipart==1) h->Fill(trk->Pt(),probTPCTOF[AliPID::kKaon]);
3373 if(ipart==2) h->Fill(trk->Pt(),probTPCTOF[AliPID::kProton]);
3378 Double_t probBayes[AliPID::kSPECIES];
3381 if(HasTOFPID(trk) && trk->Pt()>fPtTOFPIDmin){//use TOF information
3382 detUsed = CalcPIDCombined(trk, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
3383 if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return SpUndefined;//check that TPC and TOF are used
3385 detUsed = CalcPIDCombined(trk,AliPIDResponse::kDetTPC, probBayes);
3386 if(detUsed != AliPIDResponse::kDetTPC)return SpUndefined;//check that TPC is used
3389 //the probability has to be normalized to one, we check it
3391 for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
3392 if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
3393 AliFatal("Bayesian probability not normalized to one");
3396 if(fdiffPIDcutvalues){
3397 if(trk->Pt()<=4) fBayesCut=fPIDCutval1;
3398 if(trk->Pt()>4 && trk->Pt()<=6) fBayesCut=fPIDCutval2;
3399 if(trk->Pt()>6 && trk->Pt()<=8) fBayesCut=fPIDCutval3;
3400 if(trk->Pt()>8) fBayesCut=fPIDCutval4;
3404 //probabilities are normalized to one, if the cut is above .5 there is no problem
3405 if(probBayes[AliPID::kPion]>fBayesCut && IDs[SpPion]==1){
3406 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpPion));
3407 h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
3410 else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[SpKaon]==1){
3411 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpKaon));
3412 h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
3415 else if(probBayes[AliPID::kProton]>fBayesCut && IDs[SpProton]==1){
3416 TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpProton));
3417 h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
3426 //////////////////////////////////////////////////////////////////////////////////////////////////
3427 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
3428 //return the specie according to the minimum nsigma value
3429 //no double counting, this has to be evaluated using CheckDoubleCounting()
3431 Int_t ID=SpUndefined;
3433 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
3437 if(fPIDType==Bayes){//use bayesianPID
3440 AliFatal("PIDCombined object has to be set in the steering macro");
3443 ID = GetIDBayes(trk,FIllQAHistos);
3445 }else{ //use nsigma PID
3447 ID=FindMinNSigma(trk,FIllQAHistos);
3448 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
3450 HasDC=GetDoubleCounting(trk,FIllQAHistos);
3451 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3452 if(HasDC[ipart]==kTRUE) ID = SpUndefined;
3456 //Fill PID signal plot
3457 if(ID != SpUndefined){
3458 for(Int_t idet=0;idet<fNDetectors;idet++){
3459 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
3460 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3461 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3462 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3465 //Fill PID signal plot without cuts
3466 for(Int_t idet=0;idet<fNDetectors;idet++){
3467 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
3468 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3469 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3470 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3476 //-------------------------------------------------------------------------------------
3478 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
3481 AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
3482 if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
3483 //ULong_t status=track->GetStatus();
3484 //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
3485 //if (track->GetTPCsignal() <= 0.) return kFALSE;
3486 // 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.
3490 //___________________________________________________________
3493 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
3495 // check TOF matched track
3496 //ULong_t status=track->GetStatus();
3497 //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
3498 AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
3499 if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
3500 if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
3501 //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
3502 //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
3503 // if (probMis > 0.01) return kFALSE;
3504 if(fRemoveTracksT0Fill)
3506 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
3507 if (startTimeMask < 0)return kFALSE;
3512 //________________________________________________________________________
3513 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
3515 //it is called only when TOF PID is available
3516 //TOF beta calculation
3517 Double_t tofTime=track->GetTOFsignal();
3519 Double_t c=TMath::C()*1.E-9;// m/ns
3520 Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
3521 Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
3522 tofTime -= startTime; // subtract startTime to the signal
3523 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
3529 Double_t p = track->P();
3530 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
3532 track->GetIntegratedTimes(timei);
3533 return timei[0]/time;
3536 //------------------------------------------------------------------------------------------------------
3538 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)
3540 // calculate inv mass squared
3541 // same can be achieved, but with more computing time with
3542 /*TLorentzVector photon, p1, p2;
3543 p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
3544 p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
3548 Float_t tantheta1 = 1e10;
3550 if (eta1 < -1e-10 || eta1 > 1e-10)
3551 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
3553 Float_t tantheta2 = 1e10;
3554 if (eta2 < -1e-10 || eta2 > 1e-10)
3555 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
3557 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3558 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3560 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 ) ) );
3564 //---------------------------------------------------------------------------------
3566 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)
3568 // calculate inv mass squared approximately
3570 Float_t tantheta1 = 1e10;
3572 if (eta1 < -1e-10 || eta1 > 1e-10)
3574 Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
3575 tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3578 Float_t tantheta2 = 1e10;
3579 if (eta2 < -1e-10 || eta2 > 1e-10)
3581 Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
3582 tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3585 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3586 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3589 Float_t deltaPhi = TMath::Abs(phi1 - phi2);
3590 while (deltaPhi > TMath::TwoPi())
3591 deltaPhi -= TMath::TwoPi();
3592 if (deltaPhi > TMath::Pi())
3593 deltaPhi = TMath::TwoPi() - deltaPhi;
3595 Float_t cosDeltaPhi = 0;
3596 if (deltaPhi < TMath::Pi()/3)
3597 cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
3598 else if (deltaPhi < 2*TMath::Pi()/3)
3599 cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
3601 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
3603 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
3605 // Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
3609 //--------------------------------------------------------------------------------
3610 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)
3613 // calculates dphistar
3616 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
3618 static const Double_t kPi = TMath::Pi();
3621 // if (dphistar > 2 * kPi)
3622 // dphistar -= 2 * kPi;
3623 // if (dphistar < -2 * kPi)
3624 // dphistar += 2 * kPi;
3627 dphistar = kPi * 2 - dphistar;
3628 if (dphistar < -kPi)
3629 dphistar = -kPi * 2 - dphistar;
3630 if (dphistar > kPi) // might look funny but is needed
3631 dphistar = kPi * 2 - dphistar;
3635 //_________________________________________________________________________
3636 void AliTwoParticlePIDCorr ::DefineEventPool()
3638 Int_t MaxNofEvents=1000;
3639 const Int_t NofVrtxBins=10+(1+10)*2;
3640 Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
3641 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
3642 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210
3644 if (fSampleType=="pp"){
3645 const Int_t NofCentBins=4;
3646 Double_t CentralityBins[NofCentBins+1]={0.,20., 40.,60.,200.1};//Is This binning is fine for pp, or we don't require them....
3647 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
3649 if (fSampleType=="pPb" || fSampleType=="PbPb")
3651 const Int_t NofCentBins=15;
3652 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
3653 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
3655 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
3657 //if(!fPoolMgr) return kFALSE;
3661 //------------------------------------------------------------------------
3662 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
3664 // This method is a copy from AliUEHist::GetBinning
3665 // takes the binning from <configuration> identified by <tag>
3666 // configuration syntax example:
3667 // 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
3670 // returns bin edges which have to be deleted by the caller
3672 TString config(configuration);
3673 TObjArray* lines = config.Tokenize("\n");
3674 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
3676 TString line(lines->At(i)->GetName());
3677 if (line.BeginsWith(TString(tag) + ":"))
3679 line.Remove(0, strlen(tag) + 1);
3680 line.ReplaceAll(" ", "");
3681 TObjArray* binning = line.Tokenize(",");
3682 Double_t* bins = new Double_t[binning->GetEntriesFast()];
3683 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
3684 bins[j] = TString(binning->At(j)->GetName()).Atof();
3686 nBins = binning->GetEntriesFast() - 1;
3695 AliFatal(Form("Tag %s not found in %s", tag, configuration));
3699 //____________________________________________________________________
3701 Bool_t AliTwoParticlePIDCorr::CheckTrack(AliAODTrack * part)
3703 // check if the track status flags are set
3705 UInt_t status=((AliAODTrack*)part)->GetStatus();
3706 if ((status & fTrackStatus) == fTrackStatus)
3710 //________________________________________________________________________
3712 Bool_t AliTwoParticlePIDCorr::AcceptEventCentralityWeight(Double_t centrality)
3714 // rejects "randomly" events such that the centrality gets flat
3715 // uses fCentralityWeights histogram
3717 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
3719 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
3720 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
3722 Bool_t result = kFALSE;
3723 if (centralityDigits < weight)
3726 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
3731 //____________________________________________________________________
3732 void AliTwoParticlePIDCorr::ShiftTracks(TObjArray* tracks, Double_t angle)
3734 // shifts the phi angle of all tracks by angle
3735 // 0 <= angle <= 2pi
3737 for (Int_t i=0; i<tracks->GetEntriesFast(); ++i)
3739 LRCParticlePID *part=(LRCParticlePID*)(tracks->UncheckedAt(i));
3741 Double_t newAngle = part->Phi() + angle;
3742 if (newAngle >= TMath::TwoPi())
3743 newAngle -= TMath::TwoPi();
3745 part->SetPhi(newAngle);
3750 //________________________________________________________________________
3751 void AliTwoParticlePIDCorr::SetVZEROCalibrationFile(const char* filename,const char* lhcPeriod) {
3752 //Function to setup the VZERO gain equalization
3753 //============Get the equilization map============//
3754 TFile *calibrationFile = TFile::Open(filename);
3755 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
3756 Printf("No calibration file found!!!");
3760 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
3762 Printf("Calibration TList not found!!!");
3766 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
3767 if(!fHistVZEROAGainEqualizationMap) {
3768 Printf("VZERO-A calibration object not found!!!");
3771 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
3772 if(!fHistVZEROCGainEqualizationMap) {
3773 Printf("VZERO-C calibration object not found!!!");
3777 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
3778 if(!fHistVZEROChannelGainEqualizationMap) {
3779 Printf("VZERO channel calibration object not found!!!");
3784 //________________________________________________________________________
3785 Double_t AliTwoParticlePIDCorr::GetChannelEqualizationFactor(Int_t run,Int_t channel) {
3787 if(!fHistVZEROAGainEqualizationMap) return 1.0;
3789 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
3790 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
3791 if(gRunNumber == run)
3792 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
3798 //________________________________________________________________________
3799 Double_t AliTwoParticlePIDCorr::GetEqualizationFactor(Int_t run, const char* side) {
3801 if(!fHistVZEROAGainEqualizationMap) return 1.0;
3803 TString gVZEROSide = side;
3804 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
3805 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
3806 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
3807 if(gRunNumber == run) {
3808 if(gVZEROSide == "A")
3809 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
3810 else if(gVZEROSide == "C")
3811 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
3817 //________________________________________________________________________
3818 Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliAODEvent *event){
3819 //Function that returns the reference multiplicity from AODs (data or reco MC)
3820 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
3821 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
3822 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
3824 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
3826 Printf("ERROR: AOD header not available");
3829 Int_t gRunNumber = header->GetRunNumber();
3830 Float_t bSign1=header->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
3833 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++)
3834 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
3835 AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
3836 if (!track) continue;
3837 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
3838 if(tracktype!=1) continue;
3840 if(!track) continue;//for safety
3842 gRefMultiplicityTPC += 1.0;
3846 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
3847 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
3848 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
3851 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
3852 else if(iChannel >= 32)
3853 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
3856 //Equalization of gain
3857 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
3859 gRefMultiplicityVZEROA /= gFactorA;
3860 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
3862 gRefMultiplicityVZEROC /= gFactorC;
3863 if((gFactorA != 0)&&(gFactorC != 0))
3864 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
3867 //EQVZERO vs TPC multiplicity
3868 fHistEQVZEROvsTPCmultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
3869 fHistEQVZEROAvsTPCmultiplicity->Fill(gRefMultiplicityVZEROA,gRefMultiplicityTPC);
3870 fHistEQVZEROCvsTPCmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityTPC);
3872 //EQVZERO vs VZERO multiplicity
3873 fHistVZEROCvsEQVZEROCmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),gRefMultiplicityVZEROC);
3874 fHistVZEROAvsEQVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0A(),gRefMultiplicityVZEROA);
3876 //VZEROC vs VZEROA multiplicity
3877 fHistVZEROCvsVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),event->GetVZEROData()->GetMTotV0A());
3879 //EQVZEROC vs EQVZEROA multiplicity
3880 fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
3883 if(fCentralityMethod == "TRACKS_MANUAL")
3884 gRefMultiplicity = gRefMultiplicityTPC;
3885 else if(fCentralityMethod == "V0M_MANUAL")
3886 gRefMultiplicity = gRefMultiplicityVZERO;
3887 else if(fCentralityMethod == "V0A_MANUAL")
3888 gRefMultiplicity = gRefMultiplicityVZEROA;
3889 else if(fCentralityMethod == "V0C_MANUAL")
3890 gRefMultiplicity = gRefMultiplicityVZEROC;
3893 fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
3894 fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
3895 fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
3896 fHistRefmult->Fill(3.,gRefMultiplicityTPC);
3899 return gRefMultiplicity;
3902 //-------------------------------------------------------------------------------------------------------
3903 Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliAODEvent *event, Bool_t truth){
3905 if(!event) return -1;
3906 // get centrality object and check quality
3907 Double_t cent_v0=-1;
3908 Double_t nooftrackstruth=0;
3910 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
3912 AliCentrality *centralityObj=0;
3913 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
3914 if(!header) return -1;
3915 centralityObj = header->GetCentralityP();
3916 // if (centrality->GetQuality() != 0) return ;
3920 fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
3921 fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
3922 fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
3923 if(fSampleType=="pp") fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
3924 if(fSampleType=="pp") fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
3925 if(fSampleType=="pp") fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
3927 if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
3928 if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
3930 cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
3933 }//centralitymethod condition
3935 else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL")//data or RecoMc and also for TRUTH
3937 if(!truth){//for data or RecoMC
3938 cent_v0 = GetReferenceMultiplicityVZEROFromAOD(event);
3939 }//for data or RecoMC
3941 if(truth){//condition for TRUTH case
3942 //check for TClonesArray(truth track MC information)
3943 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
3945 //AliFatal("Error: MC particles branch not found!\n");
3948 //now process the truth particles(for both efficiency & correlation function)
3949 Int_t nMCTrack = fArrayMC->GetEntriesFast();
3951 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
3952 {//MC truth track loop starts
3954 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
3957 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
3961 //consider only charged particles
3962 if(partMC->Charge() == 0) continue;
3964 //consider only primary particles; neglect all secondary particles including from weak decays
3965 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
3968 //remove injected signals(primaries above <maxLabel>)
3969 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
3972 Bool_t isduplicate=kFALSE;
3973 if (fRemoveDuplicates)
3975 for (Int_t j=iMC+1; j<nMCTrack; ++j)
3976 {//2nd trutuh loop starts
3977 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
3979 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
3982 if (partMC->GetLabel() == partMC2->GetLabel())
3987 }//2nd truth loop ends
3989 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
3992 if (fCentralityMethod=="V0M_MANUAL") {
3993 if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8) continue;
3994 if (partMC->Eta() < -3.7 || partMC->Eta() > -1.7) continue;
3996 else if (fCentralityMethod=="V0A_MANUAL") {
3997 if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8) continue;}
3998 else if (fCentralityMethod=="V0C_MANUAL") {
3999 if(partMC->Eta() > -1.7 || partMC->Eta() < -3.7) continue;}
4000 else if (fCentralityMethod=="TRACKS_MANUAL") {
4001 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4002 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
4004 else{//basically returns the tracks manual case
4005 //give only kinematic cuts at the truth level
4006 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4007 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
4010 //To determine multiplicity in case of PP
4011 nooftrackstruth+= 1;;
4013 }//truth track loop ends
4014 cent_v0=nooftrackstruth;
4016 }//condition for TRUTH case
4018 }//end of MANUAL method
4020 else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC
4022 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4026 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4029 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
4030 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
4031 AliError("Event header not found. Skipping this event.");
4035 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4038 if (collGeometry) cent_v0 = collGeometry->ImpactParameter();
4040 }//end of Impact parameter method
4047 //-----------------------------------------------------------------------------------------
4048 Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliAODEvent *aod,Bool_t truth){
4052 Float_t gRefMultiplicity = -1.;
4054 // check first event in chunk (is not needed for new reconstructions)
4055 if(fCheckFirstEventInChunk){
4056 AliAnalysisUtils ut;
4057 if(ut.IsFirstEventInChunk(aod))
4062 AliAnalysisUtils ut;
4063 ut.SetUseMVPlpSelection(kTRUE);
4064 ut.SetUseOutOfBunchPileUp(kTRUE);
4065 if(ut.IsPileUpEvent(aod))
4069 //count events after pileup selection
4070 fEventCounter->Fill(3);
4072 //vertex selection(is it fine for PP?)
4073 if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4074 trkVtx = aod->GetPrimaryVertex();
4075 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
4076 TString vtxTtl = trkVtx->GetTitle();
4077 if (!vtxTtl.Contains("VertexerTracks")) return -1;
4078 zvtx = trkVtx->GetZ();
4079 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
4080 if (!spdVtx || spdVtx->GetNContributors()<=0) return -1;
4081 TString vtxTyp = spdVtx->GetTitle();
4082 Double_t cov[6]={0};
4083 spdVtx->GetCovarianceMatrix(cov);
4084 Double_t zRes = TMath::Sqrt(cov[5]);
4085 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return -1;
4086 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return -1;
4088 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
4089 Int_t nVertex = aod->GetNumberOfVertices();
4091 trkVtx = aod->GetPrimaryVertex();
4092 Int_t nTracksPrim = trkVtx->GetNContributors();
4093 zvtx = trkVtx->GetZ();
4094 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
4095 // Reject TPC only vertex
4096 TString name(trkVtx->GetName());
4097 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return -1;
4099 // Select a quality vertex by number of tracks?
4100 if( nTracksPrim < fnTracksVertex ) {
4101 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
4104 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
4105 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
4107 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
4112 else if(fVertextype==0){//default case
4113 trkVtx = aod->GetPrimaryVertex();
4114 if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
4115 zvtx = trkVtx->GetZ();
4117 trkVtx->GetCovarianceMatrix(fCov);
4118 if(fCov[5] == 0) return -1;//proper vertex resolution
4121 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
4122 return -1;//as there is no proper sample type
4125 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
4127 //count events having a proper vertex
4128 fEventCounter->Fill(5);
4130 if (TMath::Abs(zvtx) > fzvtxcut) return -1;
4132 //count events after vertex cut
4133 fEventCounter->Fill(7);
4136 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4138 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
4140 //get the centrality or multiplicity
4141 if(truth) {gRefMultiplicity = GetRefMultiOrCentrality(aod,kTRUE);}//kTRUE-->for Truth case(only meaningful in case of ref multiplicity)
4143 else {gRefMultiplicity = GetRefMultiOrCentrality(aod,kFALSE);}//kFALSE-->for data and RecoMc case(only meaningful in case of ref multiplicity)
4145 if(gRefMultiplicity<0) return -1;
4147 // take events only within the multiplicity class mentioned in the custom binning
4148 if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
4150 //count events having proper centrality/ref multiplicity
4151 fEventCounter->Fill(9);
4154 // centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
4155 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
4157 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
4161 //count events after rejection due to centrality weighting
4162 fEventCounter->Fill(11);
4164 return gRefMultiplicity;
4167 //--------------------------------------------------------------------------------------------------------
4168 Double_t AliTwoParticlePIDCorr::GetEventPlane(AliAODEvent *event,Bool_t truth){
4169 // Get the event plane
4172 Float_t gVZEROEventPlane = -10.;
4173 Float_t gReactionPlane = -10.;
4174 Double_t qxTot = 0.0, qyTot = 0.0;
4176 //MC: from reaction plane
4179 AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4182 AliGenEventHeader* eventHeader = header->GetCocktailHeader(0); // get first MC header from either ESD/AOD (including cocktail header if available)
4186 AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4188 if (collGeometry) gReactionPlane = collGeometry->ReactionPlaneAngle();
4194 AliEventplane *ep = event->GetEventplane();
4196 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
4197 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
4198 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
4199 gReactionPlane = gVZEROEventPlane;
4202 return gReactionPlane;
4210 //____________________________________________________________________
4211 void AliTwoParticlePIDCorr::Terminate(Option_t *)
4213 // Draw result to screen, or perform fitting, normalizations
4214 // Called once at the end of the query
4215 fOutput = dynamic_cast<TList*> (GetOutputData(1));
4216 if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
4220 //------------------------------------------------------------------
4223 // get centrality object and check quality
4227 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C")//for PbPb, pPb, pp7TeV(still to be introduced)
4229 AliCentrality *centralityObj=0;
4231 AliAODHeader *header = (AliAODHeader*) aod->GetHeader();
4233 centralityObj = header->GetCentralityP();
4234 // if (centrality->GetQuality() != 0) return ;
4238 if(fSampleType=="pPb" || fSampleType=="PbPb" || fSampleType=="pp") fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0M"));//only available for LHC10d at present (Quantile info)
4239 if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0A"));
4240 if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0C"));
4241 if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
4242 if(fSampleType=="pPb" || fSampleType=="PbPb") fHistCentStats->Fill(8.,centralityObj->GetCentralityPercentile("ZNA"));
4244 cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
4251 }//centralitymethod condition
4253 else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL")
4255 cent_v0 = GetReferenceMultiplicityVZEROFromAOD(aod);
4256 fHistrefMultiplicity->Fill(cent_v0);
4261 if(fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
4266 Float_t bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
4267 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
4269 // Pileup selection ************************************************
4270 if(frejectPileUp) //applicable only for TPC only tracks,not for hybrid tracks?.
4272 if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(aod)) return;
4273 //count events after PileUP cut
4274 fEventCounter->Fill(3);
4278 //vertex selection(is it fine for PP?)
4279 if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4280 trkVtx = aod->GetPrimaryVertex();
4281 if (!trkVtx || trkVtx->GetNContributors()<=0) return;
4282 TString vtxTtl = trkVtx->GetTitle();
4283 if (!vtxTtl.Contains("VertexerTracks")) return;
4284 zvtx = trkVtx->GetZ();
4285 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
4286 if (!spdVtx || spdVtx->GetNContributors()<=0) return;
4287 TString vtxTyp = spdVtx->GetTitle();
4288 Double_t cov[6]={0};
4289 spdVtx->GetCovarianceMatrix(cov);
4290 Double_t zRes = TMath::Sqrt(cov[5]);
4291 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
4292 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
4294 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
4295 Int_t nVertex = aod->GetNumberOfVertices();
4297 trkVtx = aod->GetPrimaryVertex();
4298 Int_t nTracksPrim = trkVtx->GetNContributors();
4299 zvtx = trkVtx->GetZ();
4300 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
4301 // Reject TPC only vertex
4302 TString name(trkVtx->GetName());
4303 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
4305 // Select a quality vertex by number of tracks?
4306 if( nTracksPrim < fnTracksVertex ) {
4307 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
4310 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
4311 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
4313 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
4318 else if(fVertextype==0){//default case
4319 trkVtx = aod->GetPrimaryVertex();
4320 if (!trkVtx || trkVtx->GetNContributors()<=0) return;//proper number of contributors
4321 zvtx = trkVtx->GetZ();
4323 trkVtx->GetCovarianceMatrix(fCov);
4324 if(fCov[5] == 0) return;//proper vertex resolution
4327 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
4328 return;//as there is no proper sample type
4331 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
4333 //count events having a proper vertex
4334 fEventCounter->Fill(5);
4336 if (TMath::Abs(zvtx) > fzvtxcut) return;
4338 //count events after vertex cut
4339 fEventCounter->Fill(7);
4342 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4344 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
4347 if(!aod) return; //for safety
4349 Double_t frefMult=0;
4351 //reference multiplicity for pp 7 TeV
4352 if ((fMultiplicityEstimator == "TRACKS_MANUAL") || (fMultiplicityEstimator == "V0M_MANUAL")|| (fMultiplicityEstimator == "V0A_MANUAL")||(fMultiplicityEstimator == "V0C_MANUAL")) {cent_v0=GetReferenceMultiplicityVZEROFromAOD(aod,bSign1);}
4353 else {frefMult=GetReferenceMultiplicityVZEROFromAOD(aod,bSign1);}
4357 // centrality weighting (optional for 2011 if central and semicentral triggers are used)
4358 if (fCentralityWeights && !AcceptEventCentralityWeight(cent_v0))
4360 AliInfo(Form("Rejecting event because of centrality weighting: %f", cent_v0));
4364 //count events after rejection due to centrality weighting
4365 fEventCounter->Fill(9);