1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 #include "AliTwoParticlePIDCorr.h"
17 #include "AliVParticle.h"
29 #include "AliCentrality.h"
30 #include "Riostream.h"
33 #include "AliCFContainer.h"
35 #include "THnSparse.h"
39 #include "AliESDpid.h"
40 #include "AliAODpidUtil.h"
41 #include <AliPIDResponse.h>
43 #include <AliAnalysisManager.h>
44 #include <AliInputEventHandler.h>
45 #include "AliAODInputHandler.h"
47 #include "AliAnalysisTaskSE.h"
48 #include "AliAnalysisManager.h"
49 #include "AliCentrality.h"
51 #include "AliVEvent.h"
52 #include "AliAODEvent.h"
53 #include "AliAODTrack.h"
54 #include "AliVTrack.h"
56 #include "THnSparse.h"
58 #include "AliAODMCHeader.h"
59 #include "AliAODMCParticle.h"
60 #include "AliMCEventHandler.h"
61 #include "AliMCEvent.h"
62 #include "AliMCParticle.h"
63 #include "TParticle.h"
64 #include <TDatabasePDG.h>
65 #include <TParticlePDG.h>
67 #include "AliGenCocktailEventHeader.h"
68 #include "AliGenEventHeader.h"
70 #include "AliEventPoolManager.h"
71 //#include "AliAnalysisUtils.h"
72 using namespace AliPIDNameSpace;
75 ClassImp(AliTwoParticlePIDCorr)
76 ClassImp(LRCParticlePID)
78 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
79 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
80 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
82 //________________________________________________________________________
83 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
86 fCentralityMethod("V0A"),
88 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
92 fSharedClusterCut(-1),
95 ffilltrigassoUNID(kFALSE),
96 ffilltrigUNIDassoID(kFALSE),
97 ffilltrigIDassoUNID(kTRUE),
98 ffilltrigIDassoID(kFALSE),
99 ffilltrigIDassoIDMCTRUTH(kFALSE),
100 fMaxNofMixingTracks(50000),
101 fPtOrderMCTruth(kTRUE),
102 fPtOrderDataReco(kTRUE),
103 fTriggerSpeciesSelection(kFALSE),
104 fAssociatedSpeciesSelection(kFALSE),
105 fTriggerSpecies(SpPion),
106 fAssociatedSpecies(SpPion),
109 fSelectHighestPtTrig(kFALSE),
110 fcontainPIDtrig(kTRUE),
111 fcontainPIDasso(kFALSE),
112 //frejectPileUp(kFALSE),
117 fminprotonsigmacut(-6.0),
118 fmaxprotonsigmacut(-3.0),
119 fminpionsigmacut(0.0),
120 fmaxpionsigmacut(4.0),
121 fselectprimaryTruth(kTRUE),
122 fonlyprimarydatareco(kFALSE),
125 ffillhistQAReco(kFALSE),
126 ffillhistQATruth(kFALSE),
127 kTrackVariablesPair(0),
156 fhistJetTrigestimate(0),
157 fCentralityCorrelation(0x0),
172 fCorrelatonTruthPrimary(0),
173 fCorrelatonTruthPrimarymix(0),
179 fTHnCorrIDUNIDmix(0),
181 fTHnTrigcountMCTruthPrim(0),
184 fAnalysisType("AOD"),
186 twoTrackEfficiencyCutValue(0.02),
187 //fControlConvResoncances(0),
192 fRequestTOFPID(kTRUE),
193 fPIDType(NSigmaTPCTOF),
194 fFIllPIDQAHistos(kTRUE),
196 fUseExclusiveNSigma(kFALSE),
197 fRemoveTracksT0Fill(kFALSE),
199 fTriggerSelectCharge(0),
200 fAssociatedSelectCharge(0),
201 fTriggerRestrictEta(-1),
202 fEtaOrdering(kFALSE),
203 fCutConversions(kFALSE),
204 fCutResonances(kFALSE),
205 fRejectResonanceDaughters(-1),
207 fInjectedSignals(kFALSE),
208 fRemoveWeakDecays(kFALSE),
209 fRemoveDuplicates(kFALSE),
210 fapplyTrigefficiency(kFALSE),
211 fapplyAssoefficiency(kFALSE),
212 ffillefficiency(kFALSE),
213 fmesoneffrequired(kFALSE),
214 fkaonprotoneffrequired(kFALSE),
215 //fAnalysisUtils(0x0),
219 for ( Int_t i = 0; i < 16; i++) {
223 for ( Int_t i = 0; i < 6; i++ ){
224 fTHnrecomatchedallPid[i] = NULL;
225 fTHngenprimPidTruth[i] = NULL;
226 effcorection[i]=NULL;
231 for(Int_t ipart=0;ipart<NSpecies;ipart++)
232 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
233 fnsigmas[ipart][ipid]=999.;
235 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
238 //________________________________________________________________________
239 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
240 :AliAnalysisTaskSE(name),
242 fCentralityMethod("V0A"),
244 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
248 fSharedClusterCut(-1),
251 ffilltrigassoUNID(kFALSE),
252 ffilltrigUNIDassoID(kFALSE),
253 ffilltrigIDassoUNID(kTRUE),
254 ffilltrigIDassoID(kFALSE),
255 ffilltrigIDassoIDMCTRUTH(kFALSE),
256 fMaxNofMixingTracks(50000),
257 fPtOrderMCTruth(kTRUE),
258 fPtOrderDataReco(kTRUE),
259 fTriggerSpeciesSelection(kFALSE),
260 fAssociatedSpeciesSelection(kFALSE),
261 fTriggerSpecies(SpPion),
262 fAssociatedSpecies(SpPion),
265 fSelectHighestPtTrig(kFALSE),
266 fcontainPIDtrig(kTRUE),
267 fcontainPIDasso(kFALSE),
268 // frejectPileUp(kFALSE),
273 fminprotonsigmacut(-6.0),
274 fmaxprotonsigmacut(-3.0),
275 fminpionsigmacut(0.0),
276 fmaxpionsigmacut(4.0),
277 fselectprimaryTruth(kTRUE),
278 fonlyprimarydatareco(kFALSE),
281 ffillhistQAReco(kFALSE),
282 ffillhistQATruth(kFALSE),
283 kTrackVariablesPair(0),
312 fhistJetTrigestimate(0),
313 fCentralityCorrelation(0x0),
328 fCorrelatonTruthPrimary(0),
329 fCorrelatonTruthPrimarymix(0),
335 fTHnCorrIDUNIDmix(0),
337 fTHnTrigcountMCTruthPrim(0),
340 fAnalysisType("AOD"),
342 twoTrackEfficiencyCutValue(0.02),
343 //fControlConvResoncances(0),
348 fRequestTOFPID(kTRUE),
349 fPIDType(NSigmaTPCTOF),
350 fFIllPIDQAHistos(kTRUE),
352 fUseExclusiveNSigma(kFALSE),
353 fRemoveTracksT0Fill(kFALSE),
355 fTriggerSelectCharge(0),
356 fAssociatedSelectCharge(0),
357 fTriggerRestrictEta(-1),
358 fEtaOrdering(kFALSE),
359 fCutConversions(kFALSE),
360 fCutResonances(kFALSE),
361 fRejectResonanceDaughters(-1),
363 fInjectedSignals(kFALSE),
364 fRemoveWeakDecays(kFALSE),
365 fRemoveDuplicates(kFALSE),
366 fapplyTrigefficiency(kFALSE),
367 fapplyAssoefficiency(kFALSE),
368 ffillefficiency(kFALSE),
369 fmesoneffrequired(kFALSE),
370 fkaonprotoneffrequired(kFALSE),
371 //fAnalysisUtils(0x0),
375 for ( Int_t i = 0; i < 16; i++) {
379 for ( Int_t i = 0; i < 6; i++ ){
380 fTHnrecomatchedallPid[i] = NULL;
381 fTHngenprimPidTruth[i] = NULL;
382 effcorection[i]=NULL;
387 for(Int_t ipart=0;ipart<NSpecies;ipart++)
388 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
389 fnsigmas[ipart][ipid]=999.;
391 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
393 // The last in the above list should not have a comma after it
396 // Define input and output slots here (never in the dummy constructor)
397 // Input slot #0 works with a TChain - it is connected to the default input container
398 // Output slot #1 writes into a TH1 container
400 DefineOutput(1, TList::Class()); // for output list
401 DefineOutput(2, TList::Class());
405 //________________________________________________________________________
406 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
408 // Destructor. Clean-up the output list, but not the histograms that are put inside
409 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
410 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
415 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
420 if (fPID) delete fPID;
423 //________________________________________________________________________
425 //////////////////////////////////////////////////////////////////////////////////////////////////
427 TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
428 // returns histo named name
429 return (TH2F*) fOutputList->FindObject(name);
432 //////////////////////////////////////////////////////////////////////////////////////////////////
434 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
438 // Puts the argument in the range [-pi/2,3 pi/2].
441 if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
442 if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
447 //________________________________________________________________________
448 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
451 // Called once (on the worker node)
452 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
453 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
454 fPID = inputHandler->GetPIDResponse();
456 //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
458 //get the efficiency correction map
460 // global switch disabling the reference
461 // (to avoid "Replacing existing TH1" if several wagons are created in train)
462 Bool_t oldStatus = TH1::AddDirectoryStatus();
463 TH1::AddDirectory(kFALSE);
465 fOutput = new TList();
466 fOutput->SetOwner(); // IMPORTANT!
468 fOutputList = new TList;
469 fOutputList->SetOwner();
470 fOutputList->SetName("PIDQAList");
473 Int_t centmultbins=10;
474 Double_t centmultmin=0.0;
475 Double_t centmultmax=100.0;
476 if(fSampleType=="pPb" || fSampleType=="PbPb")
482 if(fSampleType=="pp")
489 fhistcentrality=new TH1F("fhistcentrality","fhistcentrality",centmultbins*4,centmultmin,centmultmax);
490 fOutput->Add(fhistcentrality);
492 fEventCounter = new TH1F("fEventCounter","EventCounter", 10, 0.5,10.5);
493 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
494 fEventCounter->GetXaxis()->SetBinLabel(2,"Have a vertex");
495 fEventCounter->GetXaxis()->SetBinLabel(5,"After vertex Cut");
496 fEventCounter->GetXaxis()->SetBinLabel(6,"Within 0-100% centrality");
497 fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
498 //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
499 fOutput->Add(fEventCounter);
501 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
502 fOutput->Add(fEtaSpectrasso);
504 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
505 fOutput->Add(fphiSpectraasso);
507 if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
508 fOutput->Add(fCentralityCorrelation);
511 fHistoTPCdEdx = new TH2F("hHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
512 fOutputList->Add(fHistoTPCdEdx);
513 fHistoTOFbeta = new TH2F(Form("hHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
514 fOutputList->Add(fHistoTOFbeta);
516 fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
517 fOutputList->Add(fTPCTOFPion3d);
519 fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
520 fOutputList->Add(fTPCTOFKaon3d);
522 fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
523 fOutputList->Add(fTPCTOFProton3d);
527 fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);
528 fOutputList->Add(fPionPt);
529 fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);
530 fOutputList->Add(fPionEta);
531 fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);
532 fOutputList->Add(fPionPhi);
534 fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);
535 fOutputList->Add(fKaonPt);
536 fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);
537 fOutputList->Add(fKaonEta);
538 fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);
539 fOutputList->Add(fKaonPhi);
541 fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);
542 fOutputList->Add(fProtonPt);
543 fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);
544 fOutputList->Add(fProtonEta);
545 fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);
546 fOutputList->Add(fProtonPhi);
549 fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
550 fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
551 fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
552 fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);
553 fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
554 fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
555 fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
556 fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
557 fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
558 fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
559 fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
560 fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
561 fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
562 fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
563 fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
564 fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
566 for(Int_t i = 0; i < 16; i++)
568 fOutput->Add(fHistQA[i]);
571 kTrackVariablesPair=6 ;
573 if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7;
575 if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7;
577 if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8;
580 // two particle histograms
581 Int_t anaSteps = 1; // analysis steps
582 const char* title = "d^{2}N_{ch}/d#varphid#eta";
584 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
585 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
586 TString* axisTitlePair; // axis titles for track variables
587 axisTitlePair=new TString[kTrackVariablesPair];
589 TString defaultBinningStr;
590 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"
591 "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"
592 "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
593 "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"
594 "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
595 "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
596 "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"
597 "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
600 defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
603 defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
605 // =========================================================
606 // Customization (adopted from AliUEHistograms)
607 // =========================================================
609 TObjArray* lines = defaultBinningStr.Tokenize("\n");
610 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
612 TString line(lines->At(i)->GetName());
613 TString tag = line(0, line.Index(":")+1);
614 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
615 fBinningString += line + "\n";
617 AliInfo(Form("Using custom binning for %s", tag.Data()));
620 fBinningString += fCustomBinning;
622 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
624 // =========================================================
626 // =========================================================
628 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
629 axisTitlePair[0] = " multiplicity";
631 dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
632 axisTitlePair[1] = "v_{Z} (cm)";
634 dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
635 axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
637 dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
638 axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
640 dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
641 axisTitlePair[4] = "#Delta#eta";
643 dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
644 axisTitlePair[5] = "#Delta#varphi (rad)";
646 if(fcontainPIDtrig && !fcontainPIDasso){
647 dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
648 axisTitlePair[6] = "PIDTrig";
651 if(!fcontainPIDtrig && fcontainPIDasso){
652 dBinsPair[6] = GetBinning(fBinningString, "PIDAsso", iBinPair[6]);
653 axisTitlePair[6] = "PIDAsso";
656 if(fcontainPIDtrig && fcontainPIDasso){
658 dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
659 axisTitlePair[6] = "PIDTrig";
661 dBinsPair[7] = GetBinning(fBinningString, "PIDAsso", iBinPair[7]);
662 axisTitlePair[7] = "PIDAsso";
666 Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
668 Int_t nPteffbin = -1;
669 Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
672 fminPtTrig=dBinsPair[2][0];
673 fmaxPtTrig=dBinsPair[2][iBinPair[2]];
674 fminPtAsso=dBinsPair[3][0];
675 fmaxPtAsso=dBinsPair[3][iBinPair[3]];
677 //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
678 //fmaxPtComboeff=fmaxPtTrig;
679 //THnSparses for calculation of efficiency
681 if((fAnalysisType =="MCAOD") && ffillefficiency) {
682 const Int_t nDim = 4;// cent zvtx pt eta
683 Int_t fBinsCh[nDim] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//**********************change it
684 Double_t fMinCh[nDim] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
685 Double_t fMaxCh[nDim] = { dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
688 for(Int_t jj=0;jj<6;jj++)//PID type binning
690 Histrename="fTHnrecomatchedallPid";Histrename+=jj;
691 fTHnrecomatchedallPid[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);
692 fTHnrecomatchedallPid[jj]->Sumw2();
693 fTHnrecomatchedallPid[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
694 fTHnrecomatchedallPid[jj]->GetAxis(0)->SetTitle("Centrality");
695 fTHnrecomatchedallPid[jj]->GetAxis(1)->Set(iBinPair[1],dBinsPair[1]);
696 fTHnrecomatchedallPid[jj]->GetAxis(1)->SetTitle("zvtx");
697 fTHnrecomatchedallPid[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
698 fTHnrecomatchedallPid[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
699 fTHnrecomatchedallPid[jj]->GetAxis(3)->Set(nEtaBin,EtaBin);
700 fTHnrecomatchedallPid[jj]->GetAxis(3)->SetTitle("#eta");
701 fOutput->Add(fTHnrecomatchedallPid[jj]);
703 Histrename="fTHngenprimPidTruth";Histrename+=jj;
704 fTHngenprimPidTruth[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);
705 fTHngenprimPidTruth[jj]->Sumw2();
706 fTHngenprimPidTruth[jj]->GetAxis(0)->Set(iBinPair[0],dBinsPair[0]);
707 fTHngenprimPidTruth[jj]->GetAxis(0)->SetTitle("Centrality");
708 fTHngenprimPidTruth[jj]->GetAxis(1)->Set(iBinPair[1], dBinsPair[1]);
709 fTHngenprimPidTruth[jj]->GetAxis(1)->SetTitle("zvtx");
710 fTHngenprimPidTruth[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
711 fTHngenprimPidTruth[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
712 fTHngenprimPidTruth[jj]->GetAxis(3)->Set(nEtaBin, EtaBin);
713 fTHngenprimPidTruth[jj]->GetAxis(3)->SetTitle("#eta");
714 fOutput->Add(fTHngenprimPidTruth[jj]);
718 //AliThns for Correlation plots(data & MC)
720 if(ffilltrigassoUNID)
722 fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
723 for (Int_t j=0; j<kTrackVariablesPair; j++) {
724 fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
725 fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
727 fOutput->Add(fTHnCorrUNID);
729 fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
730 for (Int_t j=0; j<kTrackVariablesPair; j++) {
731 fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
732 fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
734 fOutput->Add(fTHnCorrUNIDmix);
737 if(ffilltrigIDassoID)
739 fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
740 for (Int_t j=0; j<kTrackVariablesPair; j++) {
741 fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
742 fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
744 fOutput->Add(fTHnCorrID);
746 fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
747 for (Int_t j=0; j<kTrackVariablesPair; j++) {
748 fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
749 fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
751 fOutput->Add(fTHnCorrIDmix);
754 if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
756 fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
757 for (Int_t j=0; j<kTrackVariablesPair; j++) {
758 fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
759 fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
761 fOutput->Add(fTHnCorrIDUNID);
764 fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
765 for (Int_t j=0; j<kTrackVariablesPair; j++) {
766 fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
767 fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
769 fOutput->Add(fTHnCorrIDUNIDmix);
774 //ThnSparse for Correlation plots(truth MC)
775 if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
777 fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
778 for (Int_t j=0; j<kTrackVariablesPair; j++) {
779 fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
780 fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
782 fOutput->Add(fCorrelatonTruthPrimary);
785 fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
786 for (Int_t j=0; j<kTrackVariablesPair; j++) {
787 fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
788 fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
790 fOutput->Add(fCorrelatonTruthPrimarymix);
793 //binning for trigger no. counting
797 if(fcontainPIDtrig) dims=4;
798 fBinst= new Int_t[dims];
799 for(Int_t i=0; i<3;i++)
801 fBinst[i]=iBinPair[i];
804 fBinst[3]=iBinPair[6];
807 //ThSparse for trigger counting(data & reco MC)
808 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
810 fTHnTrigcount = new AliTHn("fTHnTrigcount", "fTHnTrigcount", anaSteps, dims, fBinst);
811 for(Int_t i=0; i<3;i++){
812 fTHnTrigcount->SetBinLimits(i, dBinsPair[i]);
813 fTHnTrigcount->SetVarTitle(i, axisTitlePair[i]);
817 fTHnTrigcount->SetBinLimits(3, dBinsPair[6]);
818 fTHnTrigcount->SetVarTitle(3, axisTitlePair[6]);
820 fOutput->Add(fTHnTrigcount);
823 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
824 //AliTHns for trigger counting(truth MC)
825 fTHnTrigcountMCTruthPrim = new AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", anaSteps, dims, fBinst);
826 for(Int_t i=0; i<3;i++){
827 fTHnTrigcountMCTruthPrim->SetBinLimits(i, dBinsPair[i]);
828 fTHnTrigcountMCTruthPrim->SetVarTitle(i, axisTitlePair[i]);
832 fTHnTrigcountMCTruthPrim->SetBinLimits(3, dBinsPair[6]);
833 fTHnTrigcountMCTruthPrim->SetVarTitle(3, axisTitlePair[6]);
835 fOutput->Add(fTHnTrigcountMCTruthPrim);
838 if(fAnalysisType=="MCAOD"){
841 MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
842 fOutputList->Add(MCtruthpt);
844 MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
845 fOutputList->Add(MCtrutheta);
847 MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
848 fOutputList->Add(MCtruthphi);
850 MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
851 fOutputList->Add(MCtruthpionpt);
853 MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
854 fOutputList->Add(MCtruthpioneta);
856 MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
857 fOutputList->Add(MCtruthpionphi);
859 MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
860 fOutputList->Add(MCtruthkaonpt);
862 MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
863 fOutputList->Add(MCtruthkaoneta);
865 MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
866 fOutputList->Add(MCtruthkaonphi);
868 MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
869 fOutputList->Add(MCtruthprotonpt);
871 MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
872 fOutputList->Add(MCtruthprotoneta);
874 MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
875 fOutputList->Add(MCtruthprotonphi);
877 fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
878 fOutputList->Add(fPioncont);
880 fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
881 fOutputList->Add(fKaoncont);
883 fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
884 fOutputList->Add(fProtoncont);
887 fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
888 fEventno->GetXaxis()->SetTitle("Centrality");
889 fEventno->GetYaxis()->SetTitle("Z_Vtx");
890 fOutput->Add(fEventno);
891 fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
892 fEventnobaryon->GetXaxis()->SetTitle("Centrality");
893 fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
894 fOutput->Add(fEventnobaryon);
895 fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
896 fEventnomeson->GetXaxis()->SetTitle("Centrality");
897 fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
898 fOutput->Add(fEventnomeson);
900 fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
901 fOutput->Add(fhistJetTrigestimate);
907 if(fapplyTrigefficiency || fapplyAssoefficiency)
909 const Int_t nDimt = 4;// cent zvtx pt eta
910 Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
911 Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
912 Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
915 for(Int_t jj=0;jj<6;jj++)// PID type binning
917 Histrexname="effcorection";Histrexname+=jj;
918 effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
919 effcorection[jj]->Sumw2();
920 effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
921 effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
922 effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
923 effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
924 effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
925 effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
926 effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
927 effcorection[jj]->GetAxis(3)->SetTitle("#eta");
928 fOutput->Add(effcorection[jj]);
930 // TFile *fsifile = new TFile(fefffilename,"READ");
932 if (TString(fefffilename).BeginsWith("alien:"))
933 TGrid::Connect("alien:");
934 TFile *fileT=TFile::Open(fefffilename);
936 for(Int_t jj=0;jj<6;jj++)//type binning
938 Nameg="effmap";Nameg+=jj;
939 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
940 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
942 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
949 //fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
950 // fOutput->Add(fControlConvResoncances);
956 //*****************************************************PIDQA histos*****************************************************//
960 for(Int_t ipart=0;ipart<NSpecies;ipart++){
961 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
964 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
965 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);
966 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
967 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
968 fOutputList->Add(fHistoNSigma);
973 for(Int_t ipart=0;ipart<NSpecies;ipart++){
974 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
977 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
978 TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
979 Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
980 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
981 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
982 fOutputList->Add(fHistoNSigma);
987 if(fUseExclusiveNSigma) {
988 for(Int_t ipart=0;ipart<NSpecies;ipart++){
989 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
992 if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
993 TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
994 Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
995 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
996 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
997 fOutputList->Add(fHistoNSigma);
1002 if (fAnalysisType == "MCAOD"){
1003 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1004 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1007 if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1008 TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
1009 Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1010 fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1011 fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1012 fOutputList->Add(fHistoNSigma);
1017 for(Int_t idet=0;idet<fNDetectors;idet++){
1018 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1020 if(idet==fTOF)maxy=1.1;
1021 TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
1022 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1023 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1024 fOutputList->Add(fHistoPID);
1027 //PID signal plot, before PID cut
1028 for(Int_t idet=0;idet<fNDetectors;idet++){
1030 if(idet==fTOF)maxy=1.1;
1031 TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
1032 fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1033 fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1034 fOutputList->Add(fHistoPID);
1037 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
1038 PostData(2, fOutputList);
1040 AliInfo("Finished setting up the Output");
1042 TH1::AddDirectory(oldStatus);
1047 //-------------------------------------------------------------------------------
1048 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1051 if(fAnalysisType == "AOD") {
1055 }//AOD--analysis-----
1057 else if(fAnalysisType == "MCAOD") {
1066 //-------------------------------------------------------------------------
1067 void AliTwoParticlePIDCorr::doMCAODevent()
1069 AliVEvent *event = InputEvent();
1070 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1071 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1073 AliError("Cannot get the AOD event");
1077 // count all events(physics triggered)
1078 fEventCounter->Fill(1);
1080 // get centrality object and check quality(valid for p-Pb and Pb-Pb)
1081 Double_t cent_v0=0.0;
1083 if(fSampleType=="pPb" || fSampleType=="PbPb")
1085 AliCentrality *centrality=0;
1087 centrality = aod->GetHeader()->GetCentralityP();
1088 // if (centrality->GetQuality() != 0) return ;
1092 cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
1100 //check the PIDResponse handler
1103 // get mag. field required for twotrack efficiency cut
1105 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1107 //check for TClonesArray(truth track MC information)
1108 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1110 AliFatal("Error: MC particles branch not found!\n");
1114 //check for AliAODMCHeader(truth event MC information)
1115 AliAODMCHeader *header=NULL;
1116 header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1118 printf("MC header branch not found!\n");
1122 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
1123 Float_t zVtxmc =header->GetVtxZ();
1124 if(TMath::Abs(zVtxmc)>fzvtxcut) return;
1126 // For productions with injected signals, figure out above which label to skip particles/tracks
1127 Int_t skipParticlesAbove = 0;
1129 if (fInjectedSignals)
1131 AliGenEventHeader* eventHeader = 0;
1136 AliFatal("fInjectedSignals set but no MC header found");
1138 headers = header->GetNCocktailHeaders();
1139 eventHeader = header->GetCocktailHeader(0);
1143 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
1144 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
1145 AliError("First event header not found. Skipping this event.");
1146 //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1149 skipParticlesAbove = eventHeader->NProduced();
1150 AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
1153 //vertex selection(is it fine for PP?)
1154 if ( fVertextype==1){
1155 trkVtx = aod->GetPrimaryVertex();
1156 if (!trkVtx || trkVtx->GetNContributors()<=0) return;
1157 TString vtxTtl = trkVtx->GetTitle();
1158 if (!vtxTtl.Contains("VertexerTracks")) return;
1159 zvtx = trkVtx->GetZ();
1160 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
1161 if (!spdVtx || spdVtx->GetNContributors()<=0) return;
1162 TString vtxTyp = spdVtx->GetTitle();
1163 Double_t cov[6]={0};
1164 spdVtx->GetCovarianceMatrix(cov);
1165 Double_t zRes = TMath::Sqrt(cov[5]);
1166 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1167 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1169 else if(fVertextype==2) {//for pp and pb-pb case
1170 Int_t nVertex = aod->GetNumberOfVertices();
1172 trkVtx = aod->GetPrimaryVertex();
1173 Int_t nTracksPrim = trkVtx->GetNContributors();
1174 zvtx = trkVtx->GetZ();
1175 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1176 // Reject TPC only vertex
1177 TString name(trkVtx->GetName());
1178 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
1180 // Select a quality vertex by number of tracks?
1181 if( nTracksPrim < fnTracksVertex ) {
1182 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1185 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1186 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1188 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1193 else if(fVertextype==0){//default case
1194 trkVtx = aod->GetPrimaryVertex();
1195 if (!trkVtx || trkVtx->GetNContributors()<=0) return;//proper number of contributors
1196 zvtx = trkVtx->GetZ();
1198 trkVtx->GetCovarianceMatrix(fCov);
1199 if(fCov[5] == 0) return;//proper vertex resolution
1202 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
1203 return;//as there is no proper sample type
1207 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
1209 //count events having a proper vertex
1210 fEventCounter->Fill(2);
1212 if (TMath::Abs(zvtx) > fzvtxcut) return;
1214 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut for trkVtx only
1216 //now we have events passed physics trigger, centrality,eventzvtx cut
1218 //count events after vertex cut
1219 fEventCounter->Fill(5);
1221 if(!aod) return; //for safety
1223 if (fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
1225 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)
1227 Double_t cent_v0_truth=0.0;
1229 //TObjArray* tracksMCtruth=0;
1230 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
1231 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
1235 //There is a small difference between zvtx and zVtxmc??????
1236 //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1237 //cout<<"***********************************************zvtx="<<zvtx<<endl;
1239 //now process the truth particles(for both efficiency & correlation function)
1240 Int_t nMCTrack = fArrayMC->GetEntriesFast();
1242 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
1243 { //MC truth track loop starts
1245 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1248 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1252 //consider only charged particles
1253 if(partMC->Charge() == 0) continue;
1255 //consider only primary particles; neglect all secondary particles including from weak decays
1256 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1259 //remove injected signals(primaries above <maxLabel>)
1260 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1263 Bool_t isduplicate=kFALSE;
1264 if (fRemoveDuplicates)
1266 for (Int_t j=iMC+1; j<nMCTrack; ++j)
1267 {//2nd trutuh loop starts
1268 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1270 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1273 if (partMC->GetLabel() == partMC2->GetLabel())
1278 }//2nd truth loop ends
1280 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1282 //give only kinematic cuts at the truth level
1283 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1284 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
1286 if(!partMC) continue;//for safety
1288 //To determine multiplicity in case of PP
1290 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1291 //only physical primary(all/unidentified)
1292 if(ffillhistQATruth)
1294 MCtruthpt->Fill(partMC->Pt());
1295 MCtrutheta->Fill(partMC->Eta());
1296 MCtruthphi->Fill(partMC->Phi());
1299 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1300 Int_t particletypeTruth=-999;
1301 if (TMath::Abs(pdgtruth)==211)
1303 particletypeTruth=SpPion;
1304 if(ffillhistQATruth)
1306 MCtruthpionpt->Fill(partMC->Pt());
1307 MCtruthpioneta->Fill(partMC->Eta());
1308 MCtruthpionphi->Fill(partMC->Phi());
1311 if (TMath::Abs(pdgtruth)==321)
1313 particletypeTruth=SpKaon;
1314 if(ffillhistQATruth)
1316 MCtruthkaonpt->Fill(partMC->Pt());
1317 MCtruthkaoneta->Fill(partMC->Eta());
1318 MCtruthkaonphi->Fill(partMC->Phi());
1321 if(TMath::Abs(pdgtruth)==2212)
1323 particletypeTruth=SpProton;
1324 if(ffillhistQATruth)
1326 MCtruthprotonpt->Fill(partMC->Pt());
1327 MCtruthprotoneta->Fill(partMC->Eta());
1328 MCtruthprotonphi->Fill(partMC->Phi());
1331 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)
1333 // -- Fill THnSparse for efficiency and contamination calculation
1334 if (fSampleType=="pp") cent_v0=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
1336 Double_t primmctruth[4] = {cent_v0, zVtxmc,partMC->Pt(), partMC->Eta()};
1339 fTHngenprimPidTruth[5]->Fill(primmctruth);//for all primary truth particles(4)
1340 if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[3]->Fill(primmctruth);//for primary truth mesons(3)
1341 if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[4]->Fill(primmctruth);//for primary truth kaons+protons(4)
1342 if (TMath::Abs(pdgtruth)==211) fTHngenprimPidTruth[0]->Fill(primmctruth);//for pions
1343 if (TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[1]->Fill(primmctruth);//for kaons
1344 if(TMath::Abs(pdgtruth)==2212) fTHngenprimPidTruth[2]->Fill(primmctruth);//for protons
1347 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1348 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1350 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,partMC->Charge(),partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);
1351 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1352 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1353 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1355 }//MC truth track loop ends
1357 //*********************still in event loop
1358 Float_t weghtval=1.0;
1360 if (fSampleType=="pp") cent_v0_truth=nooftrackstruth;
1361 else cent_v0_truth=cent_v0;//the notation cent_v0 is reserved for reco/data case
1363 //now cent_v0_truth should be used for all correlation function calculation
1365 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1367 //Fill Correlations for MC truth particles(same event)
1368 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1369 Fillcorrelation(tracksMCtruth,0,cent_v0_truth,zVtxmc,weghtval,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1372 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0_truth, zVtxmc+200);
1373 if (pool2 && pool2->IsReady())
1374 {//start mixing only when pool->IsReady
1375 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1376 {//proceed only when no. of trigger particles >0 in current event
1377 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
1378 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1379 { //pool event loop start
1380 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1381 if(!bgTracks6) continue;
1382 Fillcorrelation(tracksMCtruth,bgTracks6,cent_v0_truth,zVtxmc,nmix,bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1384 }// pool event loop ends mixing case
1385 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1386 } //if pool->IsReady() condition ends mixing case
1388 //still in main event loop
1391 if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1395 //still in main event loop
1397 if(tracksMCtruth) delete tracksMCtruth;
1399 //now deal with reco tracks
1400 //TObjArray* tracksUNID=0;
1401 TObjArray* tracksUNID = new TObjArray;
1402 tracksUNID->SetOwner(kTRUE);
1404 //TObjArray* tracksID=0;
1405 TObjArray* tracksID = new TObjArray;
1406 tracksID->SetOwner(kTRUE);
1409 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1411 Double_t trackscount=0.0;
1413 // loop over reconstructed tracks
1414 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1415 { // reconstructed track loop starts
1416 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1417 if (!track) continue;
1418 //get the corresponding MC track at the truth level (doing reco matching)
1419 AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
1420 if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
1422 //remove injected signals
1423 if(fInjectedSignals)
1425 AliAODMCParticle* mother = recomatched;
1427 while (!mother->IsPhysicalPrimary())
1428 {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
1429 if (mother->GetMother() < 0)
1435 mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1441 Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
1444 if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
1445 }//remove injected signals
1447 if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
1449 Bool_t isduplicate2=kFALSE;
1450 if (fRemoveDuplicates)
1452 for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
1454 AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
1455 if (!track2) continue;
1456 AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
1457 if(!recomatched2) continue;
1459 if (track->GetLabel() == track2->GetLabel())
1466 if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
1468 fHistQA[11]->Fill(track->GetTPCNcls());
1469 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
1471 if(tracktype==0) continue;
1472 if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
1474 if(!track) continue;//for safety
1475 //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
1478 //check for eta , phi holes
1479 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1480 fphiSpectraasso->Fill(track->Phi(),track->Pt());
1482 Int_t particletypeMC=-9999;
1484 //tag all particles as unidentified
1485 particletypeMC=unidentified;
1487 Float_t effmatrix=1.;
1489 // -- Fill THnSparse for efficiency calculation
1490 if (fSampleType=="pp") cent_v0=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important
1491 //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)
1493 //Clone & Reduce track list(TObjArray) for unidentified particles
1494 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1496 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
1497 effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);
1498 LRCParticlePID* copy = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1499 copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
1500 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
1502 Double_t allrecomatchedpid[4] = {cent_v0, zVtxmc,recomatched->Pt(), recomatched->Eta()};
1503 if(ffillefficiency) fTHnrecomatchedallPid[5]->Fill(allrecomatchedpid);//for all
1505 //now start the particle identification process:)
1507 //get the pdg code of the corresponding truth particle
1508 Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
1510 switch(TMath::Abs(pdgCode)){
1512 if(fFIllPIDQAHistos){
1513 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1514 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)) && track->Pt()<fPtTOFPIDmin)continue;//not filling TOF and combined if no TOF PID
1515 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
1516 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
1521 if(fFIllPIDQAHistos){
1522 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1523 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)) && track->Pt()<fPtTOFPIDmin)continue;//not filling TOF and combined if no TOF PID
1524 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpKaon,ipid));
1525 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
1530 if(fFIllPIDQAHistos){
1531 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1532 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)) && track->Pt()<fPtTOFPIDmin)continue;//not filling TOF and combined if no TOF PID
1533 TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpPion,ipid));
1534 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
1541 Float_t dEdx = track->GetTPCsignal();
1542 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1544 if(HasTOFPID(track))
1546 Double_t beta = GetBeta(track);
1547 fHistoTOFbeta->Fill(track->Pt(), beta);
1550 //do track identification(nsigma method)
1551 particletypeMC=GetParticle(track,fFIllPIDQAHistos);//******************************problem is here
1553 //2-d TPCTOF map(for each Pt interval)
1554 if(HasTOFPID(track)){
1555 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
1556 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
1557 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
1559 if(particletypeMC==SpUndefined) continue;
1561 //Pt, Eta , Phi distribution of the reconstructed identified particles
1564 if (particletypeMC==SpPion)
1566 fPionPt->Fill(track->Pt());
1567 fPionEta->Fill(track->Eta());
1568 fPionPhi->Fill(track->Phi());
1570 if (particletypeMC==SpKaon)
1572 fKaonPt->Fill(track->Pt());
1573 fKaonEta->Fill(track->Eta());
1574 fKaonPhi->Fill(track->Phi());
1576 if (particletypeMC==SpProton)
1578 fProtonPt->Fill(track->Pt());
1579 fProtonEta->Fill(track->Eta());
1580 fProtonPhi->Fill(track->Phi());
1584 //for misidentification fraction calculation(do it with tuneonPID)
1585 if(particletypeMC==SpPion )
1587 if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
1588 if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
1589 if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
1590 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
1592 if(particletypeMC==SpKaon )
1594 if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
1595 if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
1596 if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
1597 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
1599 if(particletypeMC==SpProton )
1601 if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
1602 if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
1603 if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
1604 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
1607 //fill tracking efficiency
1610 if(particletypeMC==SpPion || particletypeMC==SpKaon)
1612 if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[3]->Fill(allrecomatchedpid);//for mesons
1614 if(particletypeMC==SpKaon || particletypeMC==SpProton)
1616 if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[4]->Fill(allrecomatchedpid);//for kaons+protons
1618 if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) fTHnrecomatchedallPid[0]->Fill(allrecomatchedpid);//for pions
1619 if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[1]->Fill(allrecomatchedpid);//for kaons
1620 if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[2]->Fill(allrecomatchedpid);//for protons
1623 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1625 if (fapplyTrigefficiency || fapplyAssoefficiency)
1626 effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
1627 LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1628 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1629 tracksID->Add(copy1);
1631 }// if(tracktype==1) condition structure ands
1633 }//reco track loop ends
1635 //*************************************************************still in event loop
1637 //same event delta-eta-deltaphi plot
1638 if(fSampleType=="pp") cent_v0=trackscount;//multiplicity
1642 //fill the centrality/multiplicity distribution of the selected events
1643 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
1645 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??????
1647 //count selected events having centrality betn 0-100%
1648 fEventCounter->Fill(6);
1650 //***************************************event no. counting
1651 Bool_t isbaryontrig=kFALSE;
1652 Bool_t ismesontrig=kFALSE;
1653 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
1655 if(tracksID && tracksID->GetEntriesFast()>0)
1657 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
1658 { //trigger loop starts
1659 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
1661 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
1662 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
1663 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
1664 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
1666 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
1667 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
1670 //same event delte-eta, delta-phi plot
1671 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
1672 {//same event calculation starts
1673 if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,weghtval,bSign,fPtOrderDataReco,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
1674 if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,weghtval,bSign,fPtOrderDataReco,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
1677 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
1678 {//same event calculation starts
1679 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,weghtval,bSign,fPtOrderDataReco,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
1680 if(ffilltrigIDassoID) Fillcorrelation(tracksID,0,cent_v0,zvtx,weghtval,bSign,fPtOrderDataReco,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
1683 //still in main event loop
1685 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
1686 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
1687 if (pool && pool->IsReady())
1688 {//start mixing only when pool->IsReady
1689 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
1690 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
1691 { //pool event loop start
1692 TObjArray* bgTracks = pool->GetEvent(jMix);
1693 if(!bgTracks) continue;
1694 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
1695 Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,nmix1,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
1696 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
1697 Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,nmix1,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
1698 }// pool event loop ends mixing case
1700 } //if pool->IsReady() condition ends mixing case
1703 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
1705 }//mixing with unidentified particles
1707 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
1708 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
1709 if (pool1 && pool1->IsReady())
1710 {//start mixing only when pool->IsReady
1711 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
1712 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
1713 { //pool event loop start
1714 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
1715 if(!bgTracks2) continue;
1716 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
1717 Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
1718 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
1719 Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,nmix2,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
1721 }// pool event loop ends mixing case
1722 } //if pool1->IsReady() condition ends mixing case
1726 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
1728 }//mixing with identified particles
1730 //no. of events analyzed
1731 fEventCounter->Fill(7);
1734 if(tracksUNID) delete tracksUNID;
1736 if(tracksID) delete tracksID;
1739 PostData(1, fOutput);
1742 //________________________________________________________________________
1743 void AliTwoParticlePIDCorr::doAODevent()
1747 AliVEvent *event = InputEvent();
1748 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1749 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1751 AliError("Cannot get the AOD event");
1756 fEventCounter->Fill(1);
1758 // get centrality object and check quality
1761 if(fSampleType=="pPb" || fSampleType=="PbPb")
1763 AliCentrality *centrality=0;
1765 centrality = aod->GetHeader()->GetCentralityP();
1766 // if (centrality->GetQuality() != 0) return ;
1770 cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
1778 Float_t bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
1779 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
1781 // Pileup selection ************************************************
1782 // if(frejectPileUp) //applicable only for TPC only tracks,not for hybrid tracks?.
1784 //if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(aod)) return;
1787 if (!fPID) return;//this should be available with each event even if we don't do PID selection
1789 //vertex selection(is it fine for PP?)
1790 if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
1791 trkVtx = aod->GetPrimaryVertex();
1792 if (!trkVtx || trkVtx->GetNContributors()<=0) return;
1793 TString vtxTtl = trkVtx->GetTitle();
1794 if (!vtxTtl.Contains("VertexerTracks")) return;
1795 zvtx = trkVtx->GetZ();
1796 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
1797 if (!spdVtx || spdVtx->GetNContributors()<=0) return;
1798 TString vtxTyp = spdVtx->GetTitle();
1799 Double_t cov[6]={0};
1800 spdVtx->GetCovarianceMatrix(cov);
1801 Double_t zRes = TMath::Sqrt(cov[5]);
1802 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1803 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1805 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
1806 Int_t nVertex = aod->GetNumberOfVertices();
1808 trkVtx = aod->GetPrimaryVertex();
1809 Int_t nTracksPrim = trkVtx->GetNContributors();
1810 zvtx = trkVtx->GetZ();
1811 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1812 // Reject TPC only vertex
1813 TString name(trkVtx->GetName());
1814 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
1816 // Select a quality vertex by number of tracks?
1817 if( nTracksPrim < fnTracksVertex ) {
1818 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1821 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1822 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1824 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1829 else if(fVertextype==0){//default case
1830 trkVtx = aod->GetPrimaryVertex();
1831 if (!trkVtx || trkVtx->GetNContributors()<=0) return;//proper number of contributors
1832 zvtx = trkVtx->GetZ();
1834 trkVtx->GetCovarianceMatrix(fCov);
1835 if(fCov[5] == 0) return;//proper vertex resolution
1838 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
1839 return;//as there is no proper sample type
1842 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
1844 //count events having a proper vertex
1845 fEventCounter->Fill(2);
1847 if (TMath::Abs(zvtx) > fzvtxcut) return;
1849 //count events after vertex cut
1850 fEventCounter->Fill(5);
1853 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
1855 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
1858 if(!aod) return; //for safety
1860 if(fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
1862 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
1863 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
1865 TObjArray* tracksID= new TObjArray;//only pions, kaons,protonsi.e. after doing the PID selection
1866 tracksID->SetOwner(kTRUE); // IMPORTANT!
1868 Double_t trackscount=0.0;//in case of pp this will give the multiplicity after the track loop(only for unidentified particles that pass the filterbit and kinematic cuts)
1872 Bool_t fTrigPtmin1=kFALSE;
1873 Bool_t fTrigPtmin2=kFALSE;
1874 Bool_t fTrigPtJet=kFALSE;
1876 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1877 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
1878 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1879 if (!track) continue;
1880 fHistQA[11]->Fill(track->GetTPCNcls());
1881 Int_t particletype=-9999;//required for PID filling
1882 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
1883 if(tracktype!=1) continue;
1885 if(!track) continue;//for safety
1887 //check for eta , phi holes
1888 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1889 fphiSpectraasso->Fill(track->Phi(),track->Pt());
1893 //if no applyefficiency , set the eff factor=1.0
1894 Float_t effmatrix=1.0;
1896 //tag all particles as unidentified that passed filterbit & kinematic cuts
1897 particletype=unidentified;
1900 if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
1901 if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
1902 if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
1905 if (fSampleType=="pp") cent_v0=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
1908 //to reduce memory consumption in pool
1909 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
1911 //Clone & Reduce track list(TObjArray) for unidentified particles
1912 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
1913 effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);
1914 LRCParticlePID* copy = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1915 copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1916 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
1919 //now start the particle identificaion process:)
1921 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
1923 Float_t dEdx = track->GetTPCsignal();
1924 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1926 //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)
1927 if(HasTOFPID(track))
1929 Double_t beta = GetBeta(track);
1930 fHistoTOFbeta->Fill(track->Pt(), beta);
1934 //track identification(using nsigma method)
1935 particletype=GetParticle(track,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
1938 //2-d TPCTOF map(for each Pt interval)
1939 if(HasTOFPID(track)){
1940 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
1941 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
1942 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
1945 //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!!!!!
1946 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)!!!!!!!!!!!
1948 //Pt, Eta , Phi distribution of the reconstructed identified particles
1951 if (particletype==SpPion)
1953 fPionPt->Fill(track->Pt());
1954 fPionEta->Fill(track->Eta());
1955 fPionPhi->Fill(track->Phi());
1957 if (particletype==SpKaon)
1959 fKaonPt->Fill(track->Pt());
1960 fKaonEta->Fill(track->Eta());
1961 fKaonPhi->Fill(track->Phi());
1963 if (particletype==SpProton)
1965 fProtonPt->Fill(track->Pt());
1966 fProtonEta->Fill(track->Eta());
1967 fProtonPhi->Fill(track->Phi());
1971 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
1973 if (fapplyTrigefficiency || fapplyAssoefficiency)
1974 effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles; Bool_t mesoneffrequired=kFALSE
1975 LRCParticlePID* copy1 = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1976 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1977 tracksID->Add(copy1);
1979 } //track loop ends but still in event loop
1981 if(trackscount<1.0){
1982 if(tracksUNID) delete tracksUNID;
1983 if(tracksID) delete tracksID;
1987 if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
1988 if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
1989 if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
1991 Float_t weightval=1.0;
1994 // cout<<fminPtAsso<<"***"<<fmaxPtAsso<<"*********************"<<fminPtTrig<<"***"<<fmaxPtTrig<<"*****************"<<kTrackVariablesPair<<endl;
1995 if(fSampleType=="pp") cent_v0=trackscount;//multiplicity
1997 //fill the centrality/multiplicity distribution of the selected events
1998 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2000 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??????
2002 //count selected events having centrality betn 0-100%
2003 fEventCounter->Fill(6);
2005 //***************************************event no. counting
2006 Bool_t isbaryontrig=kFALSE;
2007 Bool_t ismesontrig=kFALSE;
2008 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2010 if(tracksID && tracksID->GetEntriesFast()>0)
2012 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2013 { //trigger loop starts
2014 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2016 if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2017 Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2018 if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2019 if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2021 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2022 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2026 //same event delta-eta-deltaphi plot
2028 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2029 {//same event calculation starts
2030 if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,weightval,bSign,fPtOrderDataReco,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2031 if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,weightval,bSign,fPtOrderDataReco,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
2034 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2035 {//same event calculation starts
2036 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,weightval,bSign,fPtOrderDataReco,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2037 if(ffilltrigIDassoID) Fillcorrelation(tracksID,0,cent_v0,zvtx,weightval,bSign,fPtOrderDataReco,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2040 //still in main event loop
2044 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2045 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
2046 if (pool && pool->IsReady())
2047 {//start mixing only when pool->IsReady
2048 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2049 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2050 { //pool event loop start
2051 TObjArray* bgTracks = pool->GetEvent(jMix);
2052 if(!bgTracks) continue;
2053 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2054 Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,nmix1,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
2055 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2056 Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,nmix1,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2057 }// pool event loop ends mixing case
2059 } //if pool->IsReady() condition ends mixing case
2062 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2064 }//mixing with unidentified particles
2067 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2068 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
2069 if (pool1 && pool1->IsReady())
2070 {//start mixing only when pool->IsReady
2071 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2072 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2073 { //pool event loop start
2074 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2075 if(!bgTracks2) continue;
2076 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2077 Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2078 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2079 Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,nmix2,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
2081 }// pool event loop ends mixing case
2082 } //if pool1->IsReady() condition ends mixing case
2086 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2088 }//mixing with identified particles
2091 //no. of events analyzed
2092 fEventCounter->Fill(7);
2094 //still in main event loop
2097 if(tracksUNID) delete tracksUNID;
2099 if(tracksID) delete tracksID;
2102 PostData(1, fOutput);
2104 } // *************************event loop ends******************************************//_______________________________________________________________________
2105 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
2107 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
2109 TObjArray* tracksClone = new TObjArray;
2110 tracksClone->SetOwner(kTRUE);
2112 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2114 LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
2115 LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval());
2116 copy100->SetUniqueID(particle->GetUniqueID());
2117 tracksClone->Add(copy100);
2123 //--------------------------------------------------------------------------------
2124 void AliTwoParticlePIDCorr::Fillcorrelation(TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t weight,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)
2127 //before calling this function check that either trackstrig & tracksasso are available
2129 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
2130 TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
2131 TArrayF eta(input->GetEntriesFast());
2132 for (Int_t i=0; i<input->GetEntriesFast(); i++)
2133 eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
2136 Int_t jmax=trackstrig->GetEntriesFast();
2138 jmax=tracksasso->GetEntriesFast();
2140 // identify K, Lambda candidates and flag those particles
2141 // a TObject bit is used for this
2142 const UInt_t kResonanceDaughterFlag = 1 << 14;
2143 if (fRejectResonanceDaughters > 0)
2145 Double_t resonanceMass = -1;
2146 Double_t massDaughter1 = -1;
2147 Double_t massDaughter2 = -1;
2148 const Double_t interval = 0.02;
2149 switch (fRejectResonanceDaughters)
2151 case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
2152 case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
2153 case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
2154 default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
2157 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2158 trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2159 for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
2160 tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2162 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2164 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2165 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
2167 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2169 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
2170 if (trig->IsEqual(asso)) continue;
2172 if (trig->Charge() * asso->Charge() > 0) continue;
2174 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2176 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
2178 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2180 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
2182 trig->SetBit(kResonanceDaughterFlag);
2183 asso->SetBit(kResonanceDaughterFlag);
2185 // Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
2192 //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
2194 Float_t TriggerPtMin=fminPtTrig;
2195 Float_t TriggerPtMax=fmaxPtTrig;
2196 Int_t HighestPtTriggerIndx=-99999;
2198 if(fSelectHighestPtTrig)//**************add this data member to the constructor
2200 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2201 { //trigger loop starts(highest Pt trigger selection)
2202 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2204 Float_t trigpt=trig->Pt();
2205 //to avoid overflow qnd underflow
2206 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2207 Int_t particlepidtrig=trig->getparticle();
2208 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2210 Float_t trigeta=trig->Eta();
2212 // some optimization
2213 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2216 if (fOnlyOneEtaSide != 0)
2218 if (fOnlyOneEtaSide * trigeta < 0)
2221 if (fTriggerSelectCharge != 0)
2222 if (trig->Charge() * fTriggerSelectCharge < 0)
2225 if (fRejectResonanceDaughters > 0)
2226 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2228 if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
2230 //TriggerPt=trigpt;//*********think about it
2232 HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
2233 TriggerPtMin=trigpt;
2236 }//trigger loop ends(highest Pt trigger selection)
2238 }//******************SelectHighestPtTrig condition ends
2241 //two particle correlation filling
2242 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2243 { //trigger loop starts
2244 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2246 Float_t trigpt=trig->Pt();
2247 //to avoid overflow qnd underflow
2248 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2249 Int_t particlepidtrig=trig->getparticle();
2250 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2252 Float_t trigeta=trig->Eta();
2254 // some optimization
2255 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2258 if (fOnlyOneEtaSide != 0)
2260 if (fOnlyOneEtaSide * trigeta < 0)
2263 if (fTriggerSelectCharge != 0)
2264 if (trig->Charge() * fTriggerSelectCharge < 0)
2267 if (fRejectResonanceDaughters > 0)
2268 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2270 if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
2271 if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
2274 Float_t trigphi=trig->Phi();
2275 Float_t trackefftrig=1.0;
2276 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
2277 //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
2280 if(fcontainPIDtrig) dim=4;
2281 trigval= new Double_t[dim];
2284 trigval[2] = trigpt;
2285 if(fcontainPIDtrig) trigval[3] = particlepidtrig;
2287 //filling only for same event case(mixcase=kFALSE)
2289 //trigger particle counting only when mixcase=kFALSE i.e for same event correlation function calculation
2292 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
2293 if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2295 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
2296 if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2298 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
2299 if(fillup=="trigUNIDassoID") fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2301 //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
2302 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
2303 if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2305 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
2306 if(fillup=="trigIDassoUNID" ) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2308 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
2309 if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2312 if(fillup=="trigIDassoIDMCTRUTH") fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig); //In truth MC case "Unidentified" means any particle other than pion,kaon or proton and no efficiency correction(default value 1.0)************************be careful!!!!
2315 //asso loop starts within trigger loop
2316 for(Int_t j=0;j<jmax;j++)
2318 LRCParticlePID *asso=0;
2320 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
2322 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2326 //to avoid overflow qnd underflow
2327 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
2329 //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
2331 if(!tracksasso && j==i) continue;
2333 // 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)
2334 // if (tracksasso && trig->IsEqual(asso)) continue;
2336 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
2339 if (asso->Pt() >= trig->Pt()) continue;
2341 Int_t particlepidasso=asso->getparticle();
2342 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
2345 if (fAssociatedSelectCharge != 0)
2346 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
2348 if (fSelectCharge > 0)
2351 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
2355 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
2361 if (trigeta < 0 && asso->Eta() < trigeta)
2363 if (trigeta > 0 && asso->Eta() > trigeta)
2367 if (fRejectResonanceDaughters > 0)
2368 if (asso->TestBit(kResonanceDaughterFlag))
2370 // Printf("Skipped j=%d", j);
2375 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
2377 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2381 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2383 //fControlConvResoncances->Fill(0.0, mass);
2385 if (mass < 0.04*0.04)
2391 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2393 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
2395 const Float_t kK0smass = 0.4976;
2397 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
2399 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
2401 //fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
2403 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
2409 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2411 Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
2412 Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2414 const Float_t kLambdaMass = 1.115;
2416 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
2418 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
2420 //fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
2422 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2425 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
2427 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2429 //fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
2431 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2436 if (twoTrackEfficiencyCut)
2438 // the variables & cuthave been developed by the HBT group
2439 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
2440 Float_t phi1 = trig->Phi();
2441 Float_t pt1 = trig->Pt();
2442 Float_t charge1 = trig->Charge();
2443 Float_t phi2 = asso->Phi();
2444 Float_t pt2 = asso->Pt();
2445 Float_t charge2 = asso->Charge();
2447 Float_t deta= trigeta - eta[j];
2450 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
2453 // check first boundaries to see if is worth to loop and find the minimum
2454 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
2455 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
2457 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
2459 Float_t dphistarminabs = 1e5;
2460 Float_t dphistarmin = 1e5;
2462 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
2464 for (Double_t rad=0.8; rad<2.51; rad+=0.01)
2466 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
2468 Float_t dphistarabs = TMath::Abs(dphistar);
2470 if (dphistarabs < dphistarminabs)
2472 dphistarmin = dphistar;
2473 dphistarminabs = dphistarabs;
2477 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
2479 // 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);
2482 //fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
2488 Float_t weightperevent=weight;
2489 Float_t trackeffasso=1.0;
2490 if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
2491 //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
2492 Float_t deleta=trigeta-eta[j];
2493 Float_t delphi=PhiRange(trigphi-asso->Phi());
2495 //here get the two particle efficiency correction factor
2496 Float_t effweight=trackefftrig*trackeffasso*weightperevent;
2497 // if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
2499 vars= new Double_t[kTrackVariablesPair];
2506 if(fcontainPIDtrig && !fcontainPIDasso) vars[6]=particlepidtrig;
2507 if(!fcontainPIDtrig && fcontainPIDasso) vars[6]=particlepidasso;
2508 if(fcontainPIDtrig && fcontainPIDasso){
2509 vars[6]=particlepidtrig;
2510 vars[7]=particlepidasso;
2513 //Fill Correlation ThnSparses
2514 if(fillup=="trigassoUNID")
2516 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,0,1.0/effweight);
2517 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
2519 if(fillup=="trigIDassoID")
2521 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,0,1.0/effweight);
2522 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
2524 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
2525 {//in this case effweight should be 1 always
2526 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
2527 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
2529 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
2531 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
2532 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
2538 }//trigger loop ends
2542 //--------------------------------------------------------------------------------
2543 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
2545 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
2547 Float_t effvalue=1.;
2549 if(parpid==unidentified)
2551 effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
2552 effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
2553 effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
2554 effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
2555 effvalue=effcorection[5]->GetBinContent(effVars);
2557 if(parpid==SpPion || parpid==SpKaon)
2559 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2561 effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
2562 effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
2563 effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
2564 effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
2565 effvalue=effcorection[3]->GetBinContent(effVars);
2570 effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
2571 effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
2572 effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
2573 effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
2574 effvalue=effcorection[0]->GetBinContent(effVars);
2579 effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
2580 effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
2581 effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
2582 effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
2583 effvalue=effcorection[1]->GetBinContent(effVars);
2588 if(parpid==SpProton)
2590 effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
2591 effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
2592 effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
2593 effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
2594 effvalue=effcorection[2]->GetBinContent(effVars);
2597 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2599 if(parpid==SpProton || parpid==SpKaon)
2601 effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
2602 effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
2603 effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
2604 effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
2605 effvalue=effcorection[4]->GetBinContent(effVars);
2608 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
2609 if(effvalue==0.) effvalue=1.;
2614 //-----------------------------------------------------------------------
2616 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield)
2619 if(!track) return 0;
2620 Bool_t trackOK = track->TestFilterBit(fFilterBit);
2621 if(!trackOK) return 0;
2622 //select only primary traks(for data & reco MC tracks)
2623 if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
2624 if(track->Charge()==0) return 0;
2625 fHistQA[12]->Fill(track->GetTPCNcls());
2628 dz = track->ZAtDCA();
2629 fHistQA[6]->Fill(dxy);
2630 fHistQA[7]->Fill(dz);
2631 Float_t chi2ndf = track->Chi2perNDF();
2632 fHistQA[13]->Fill(chi2ndf);
2633 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
2634 fHistQA[14]->Fill(nCrossedRowsTPC);
2635 //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
2636 if (track->GetTPCNclsF()>0) {
2637 Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
2638 fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
2641 Float_t pt=track->Pt();
2642 if(pt< fminPt || pt> fmaxPt) return 0;
2643 if(TMath::Abs(track->Eta())> fmaxeta) return 0;
2644 if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
2645 //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required
2647 if (fdcacut && fDCAXYCut)
2654 AliAODTrack* clone =(AliAODTrack*) track->Clone();
2655 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
2660 // Printf("%f", ((AliAODTrack*)part)->DCA());
2661 // Printf("%f", pos[0]);
2662 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
2666 if (fSharedClusterCut >= 0)
2668 Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
2669 if (frac > fSharedClusterCut)
2672 fHistQA[8]->Fill(pt);
2673 fHistQA[9]->Fill(track->Eta());
2674 fHistQA[10]->Fill(track->Phi());
2677 //________________________________________________________________________________
2678 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
2680 //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 )
2681 Float_t pt=track->Pt();
2683 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
2684 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
2685 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
2686 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
2688 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
2689 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
2691 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
2694 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
2695 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
2696 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
2698 nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
2699 nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
2700 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
2706 // if TOF is missing and below fPtTOFPID only the TPC information is used
2707 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
2708 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
2709 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
2713 //set data member fnsigmas
2714 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
2715 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
2716 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
2718 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
2719 fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
2720 fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
2721 fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
2723 //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)
2724 fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
2725 fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
2726 fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
2729 //Fill NSigma SeparationPlot
2730 for(Int_t ipart=0;ipart<NSpecies;ipart++){
2731 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2732 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)) && track->Pt()<fPtTOFPIDmin)continue;//not filling TOF and combined if no TOF PID
2733 TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
2734 h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
2740 //----------------------------------------------------------------------------
2741 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
2743 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
2744 if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && track->Pt()<=fPtTOFPIDmax && (!HasTOFPID(track)) )return SpUndefined;//so any track having Pt>0.6 && Pt<=4.0 Gev withot having proper TOF response will be defined as SpUndefined
2745 //get the identity of the particle with the minimum Nsigma
2746 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
2749 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
2750 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
2751 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
2754 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
2755 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
2756 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
2758 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
2759 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
2760 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
2761 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
2765 if(track->Pt()<=fPtTOFPIDmax) {
2766 // guess the particle based on the smaller nsigma (within fNSigmaPID)
2767 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
2768 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
2770 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2771 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)) && (track->Pt()<fPtTOFPIDmin))continue;//not filling TOF and combined if no TOF PID
2772 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
2773 h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
2778 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
2780 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2781 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)) && (track->Pt()<fPtTOFPIDmin))continue;//not filling TOF and combined if no TOF PID
2782 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
2783 h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
2788 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
2790 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2791 if((ipid!=NSigmaTPC) && (!HasTOFPID(track)) && (track->Pt()<fPtTOFPIDmin))continue;//not filling TOF and combined if no TOF PID
2792 TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
2793 h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
2799 // else, return undefined
2802 else {//asymmetric nsigma cut around pion's mean position for tracks having Pt>4.0 Gev
2803 if(fminprotonsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxprotonsigmacut) return SpProton;
2804 if(fminpionsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxpionsigmacut) return SpPion;
2805 // else, return undefined(here we don't detect kaons)
2811 //------------------------------------------------------------------------------------------
2812 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
2813 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
2815 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
2817 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
2819 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
2822 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
2824 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
2827 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
2828 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
2829 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
2832 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
2833 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
2834 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
2836 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
2837 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
2838 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
2839 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
2843 //there is chance of overlapping only for particles having pt below 4.0 GEv
2845 if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
2846 if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
2847 if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
2852 //fill NSigma distr for double counting
2853 for(Int_t ipart=0;ipart<NSpecies;ipart++){
2854 if(fHasDoubleCounting[ipart]){
2855 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2856 if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)) && (trk->Pt()<fPtTOFPIDmin))continue;//not filling TOF and combined if no TOF PID
2857 TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
2858 h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
2865 return fHasDoubleCounting;
2868 //////////////////////////////////////////////////////////////////////////////////////////////////
2869 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
2870 //return the specie according to the minimum nsigma value
2871 //no double counting, this has to be evaluated using CheckDoubleCounting()
2873 Int_t ID=SpUndefined;
2875 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
2877 ID=FindMinNSigma(trk,FIllQAHistos);
2879 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
2881 HasDC=GetDoubleCounting(trk,FIllQAHistos);
2882 for(Int_t ipart=0;ipart<NSpecies;ipart++){
2883 if(HasDC[ipart]==kTRUE) ID = SpUndefined;
2887 //Fill PID signal plot
2888 if(ID != SpUndefined){
2889 for(Int_t idet=0;idet<fNDetectors;idet++){
2890 TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
2891 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
2892 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
2893 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
2896 //Fill PID signal plot without cuts
2897 for(Int_t idet=0;idet<fNDetectors;idet++){
2898 TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
2899 if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
2900 if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
2901 if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
2907 //-------------------------------------------------------------------------------------
2909 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
2912 AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
2913 if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
2914 //ULong_t status=track->GetStatus();
2915 //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
2916 //if (track->GetTPCsignal() <= 0.) return kFALSE;
2917 // 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.
2921 //___________________________________________________________
2924 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
2926 // check TOF matched track
2927 //ULong_t status=track->GetStatus();
2928 //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
2929 AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
2930 if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
2931 if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
2932 //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
2933 //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
2934 // if (probMis > 0.01) return kFALSE;
2935 if(fRemoveTracksT0Fill)
2937 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
2938 if (startTimeMask < 0)return kFALSE;
2943 //________________________________________________________________________
2944 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
2946 //it is called only when TOF PID is available
2947 //TOF beta calculation
2948 Double_t tofTime=track->GetTOFsignal();
2950 Double_t c=TMath::C()*1.E-9;// m/ns
2951 Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
2952 Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
2953 tofTime -= startTime; // subtract startTime to the signal
2954 Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
2960 Double_t p = track->P();
2961 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
2963 track->GetIntegratedTimes(timei);
2964 return timei[0]/time;
2967 //------------------------------------------------------------------------------------------------------
2969 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)
2971 // calculate inv mass squared
2972 // same can be achieved, but with more computing time with
2973 /*TLorentzVector photon, p1, p2;
2974 p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
2975 p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
2979 Float_t tantheta1 = 1e10;
2981 if (eta1 < -1e-10 || eta1 > 1e-10)
2982 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
2984 Float_t tantheta2 = 1e10;
2985 if (eta2 < -1e-10 || eta2 > 1e-10)
2986 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
2988 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
2989 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
2991 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 ) ) );
2995 //---------------------------------------------------------------------------------
2997 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)
2999 // calculate inv mass squared approximately
3001 Float_t tantheta1 = 1e10;
3003 if (eta1 < -1e-10 || eta1 > 1e-10)
3005 Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
3006 tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3009 Float_t tantheta2 = 1e10;
3010 if (eta2 < -1e-10 || eta2 > 1e-10)
3012 Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
3013 tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3016 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3017 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3020 Float_t deltaPhi = TMath::Abs(phi1 - phi2);
3021 while (deltaPhi > TMath::TwoPi())
3022 deltaPhi -= TMath::TwoPi();
3023 if (deltaPhi > TMath::Pi())
3024 deltaPhi = TMath::TwoPi() - deltaPhi;
3026 Float_t cosDeltaPhi = 0;
3027 if (deltaPhi < TMath::Pi()/3)
3028 cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
3029 else if (deltaPhi < 2*TMath::Pi()/3)
3030 cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
3032 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
3034 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
3036 // Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
3040 //--------------------------------------------------------------------------------
3041 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)
3044 // calculates dphistar
3047 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
3049 static const Double_t kPi = TMath::Pi();
3052 // if (dphistar > 2 * kPi)
3053 // dphistar -= 2 * kPi;
3054 // if (dphistar < -2 * kPi)
3055 // dphistar += 2 * kPi;
3058 dphistar = kPi * 2 - dphistar;
3059 if (dphistar < -kPi)
3060 dphistar = -kPi * 2 - dphistar;
3061 if (dphistar > kPi) // might look funny but is needed
3062 dphistar = kPi * 2 - dphistar;
3066 //_________________________________________________________________________
3067 void AliTwoParticlePIDCorr ::DefineEventPool()
3069 Int_t MaxNofEvents=1000;
3070 const Int_t NofVrtxBins=10+(1+10)*2;
3071 Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
3072 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
3073 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210
3075 if (fSampleType=="pp"){
3076 const Int_t NofCentBins=4;
3077 Double_t CentralityBins[NofCentBins+1]={0.,20., 40.,60.,200.1};//Is This binning is fine for pp, or we don't require them....
3078 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
3080 if (fSampleType=="pPb" || fSampleType=="PbPb")
3082 const Int_t NofCentBins=15;
3083 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
3084 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
3086 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
3088 //if(!fPoolMgr) return kFALSE;
3092 //------------------------------------------------------------------------
3093 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
3095 // This method is a copy from AliUEHist::GetBinning
3096 // takes the binning from <configuration> identified by <tag>
3097 // configuration syntax example:
3098 // 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
3101 // returns bin edges which have to be deleted by the caller
3103 TString config(configuration);
3104 TObjArray* lines = config.Tokenize("\n");
3105 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
3107 TString line(lines->At(i)->GetName());
3108 if (line.BeginsWith(TString(tag) + ":"))
3110 line.Remove(0, strlen(tag) + 1);
3111 line.ReplaceAll(" ", "");
3112 TObjArray* binning = line.Tokenize(",");
3113 Double_t* bins = new Double_t[binning->GetEntriesFast()];
3114 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
3115 bins[j] = TString(binning->At(j)->GetName()).Atof();
3117 nBins = binning->GetEntriesFast() - 1;
3126 AliFatal(Form("Tag %s not found in %s", tag, configuration));
3129 //________________________________________________________________________
3130 void AliTwoParticlePIDCorr::Terminate(Option_t *)
3132 // Draw result to screen, or perform fitting, normalizations
3133 // Called once at the end of the query
3134 fOutput = dynamic_cast<TList*> (GetOutputData(1));
3135 if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
3139 //------------------------------------------------------------------