]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/Correlations/DPhi/TriggerPID/AliTwoParticlePIDCorr.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
TriggerPID: Updates from Debojit
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / TriggerPID / AliTwoParticlePIDCorr.cxx
CommitLineData
c683985a 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16#include "AliTwoParticlePIDCorr.h"
17#include "AliVParticle.h"
18#include "TFormula.h"
19#include "TAxis.h"
20#include "TChain.h"
21#include "TTree.h"
22#include "TH1F.h"
23#include "TH2F.h"
24#include "TH3F.h"
25#include "TList.h"
26#include "TFile.h"
4a2a9605 27#include "TGrid.h"
c683985a 28
29#include "AliCentrality.h"
30#include "Riostream.h"
31
83cdcf92 32#include "AliTHn.h"
33#include "AliCFContainer.h"
34#include "THn.h"
35#include "THnSparse.h"
36
c683985a 37#include <TSpline.h>
38#include <AliPID.h>
39#include "AliESDpid.h"
40#include "AliAODpidUtil.h"
41#include <AliPIDResponse.h>
42
43#include <AliAnalysisManager.h>
44#include <AliInputEventHandler.h>
45#include "AliAODInputHandler.h"
46
47#include "AliAnalysisTaskSE.h"
48#include "AliAnalysisManager.h"
49#include "AliCentrality.h"
50
51#include "AliVEvent.h"
52#include "AliAODEvent.h"
53#include "AliAODTrack.h"
54#include "AliVTrack.h"
55
56#include "THnSparse.h"
57
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>
66
67#include "AliGenCocktailEventHeader.h"
68#include "AliGenEventHeader.h"
69
70#include "AliEventPoolManager.h"
71//#include "AliAnalysisUtils.h"
72using namespace AliPIDNameSpace;
73using namespace std;
74
75ClassImp(AliTwoParticlePIDCorr)
76ClassImp(LRCParticlePID)
83cdcf92 77
78const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
79const char * kDetectorName[]={"ITS","TPC","TOF"} ;
80const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
81
c683985a 82//________________________________________________________________________
83AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
84:AliAnalysisTaskSE(),
85 fOutput(0),
86 fCentralityMethod("V0A"),
87 fSampleType("pPb"),
ef349d3c 88 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
c683985a 89 trkVtx(0),
90 zvtx(0),
91 fFilterBit(768),
ef349d3c 92 fSharedClusterCut(-1),
93 fVertextype(1),
c683985a 94 fzvtxcut(10.0),
95 ffilltrigassoUNID(kFALSE),
96 ffilltrigUNIDassoID(kFALSE),
97 ffilltrigIDassoUNID(kTRUE),
98 ffilltrigIDassoID(kFALSE),
99 ffilltrigIDassoIDMCTRUTH(kFALSE),
ef349d3c 100 fMaxNofMixingTracks(50000),
0684fc6a 101 fPtOrderMCTruth(kTRUE),
102 fPtOrderDataReco(kTRUE),
c683985a 103 fTriggerSpeciesSelection(kFALSE),
104 fAssociatedSpeciesSelection(kFALSE),
105 fTriggerSpecies(SpPion),
106 fAssociatedSpecies(SpPion),
107 fCustomBinning(""),
108 fBinningString(""),
ef349d3c 109 fSelectHighestPtTrig(kFALSE),
c683985a 110 fcontainPIDtrig(kTRUE),
111 fcontainPIDasso(kFALSE),
112//frejectPileUp(kFALSE),
113 fminPt(0.2),
114 fmaxPt(10.0),
115 fmineta(-0.8),
116 fmaxeta(0.8),
117 fminprotonsigmacut(-6.0),
118 fmaxprotonsigmacut(-3.0),
119 fminpionsigmacut(0.0),
120 fmaxpionsigmacut(4.0),
121 fselectprimaryTruth(kTRUE),
122 fonlyprimarydatareco(kFALSE),
123 fdcacut(kFALSE),
124 fdcacutvalue(3.0),
125 ffillhistQAReco(kFALSE),
126 ffillhistQATruth(kFALSE),
127 kTrackVariablesPair(0),
128 fminPtTrig(0),
129 fmaxPtTrig(0),
130 fminPtComboeff(2.0),
131 fmaxPtComboeff(4.0),
132 fminPtAsso(0),
133 fmaxPtAsso(0),
134 fhistcentrality(0),
135 fEventCounter(0),
136 fEtaSpectrasso(0),
137 fphiSpectraasso(0),
138 MCtruthpt(0),
139 MCtrutheta(0),
140 MCtruthphi(0),
141 MCtruthpionpt(0),
142 MCtruthpioneta(0),
143 MCtruthpionphi(0),
144 MCtruthkaonpt(0),
145 MCtruthkaoneta(0),
146 MCtruthkaonphi(0),
147 MCtruthprotonpt(0),
148 MCtruthprotoneta(0),
149 MCtruthprotonphi(0),
150 fPioncont(0),
151 fKaoncont(0),
152 fProtoncont(0),
83cdcf92 153 fEventno(0),
154 fEventnobaryon(0),
155 fEventnomeson(0),
0684fc6a 156 fhistJetTrigestimate(0),
c683985a 157 fCentralityCorrelation(0x0),
158 fHistoTPCdEdx(0x0),
159 fHistoTOFbeta(0x0),
160 fTPCTOFPion3d(0),
161 fTPCTOFKaon3d(0),
162 fTPCTOFProton3d(0),
163 fPionPt(0),
164 fPionEta(0),
165 fPionPhi(0),
166 fKaonPt(0),
167 fKaonEta(0),
168 fKaonPhi(0),
169 fProtonPt(0),
170 fProtonEta(0),
171 fProtonPhi(0),
172 fCorrelatonTruthPrimary(0),
173 fCorrelatonTruthPrimarymix(0),
174 fTHnCorrUNID(0),
175 fTHnCorrUNIDmix(0),
176 fTHnCorrID(0),
177 fTHnCorrIDmix(0),
178 fTHnCorrIDUNID(0),
179 fTHnCorrIDUNIDmix(0),
180 fTHnTrigcount(0),
181 fTHnTrigcountMCTruthPrim(0),
182 fPoolMgr(0x0),
183 fArrayMC(0),
4a2a9605 184 fAnalysisType("AOD"),
c683985a 185 fefffilename(""),
186 twoTrackEfficiencyCutValue(0.02),
187//fControlConvResoncances(0),
188 fPID(NULL),
189 eventno(0),
190 fPtTOFPIDmin(0.6),
191 fPtTOFPIDmax(4.0),
192 fRequestTOFPID(kTRUE),
193 fPIDType(NSigmaTPCTOF),
83cdcf92 194 fFIllPIDQAHistos(kTRUE),
c683985a 195 fNSigmaPID(3),
196 fUseExclusiveNSigma(kFALSE),
197 fRemoveTracksT0Fill(kFALSE),
198fSelectCharge(0),
199fTriggerSelectCharge(0),
200fAssociatedSelectCharge(0),
201fTriggerRestrictEta(-1),
202fEtaOrdering(kFALSE),
203fCutConversions(kFALSE),
204fCutResonances(kFALSE),
205fRejectResonanceDaughters(-1),
206 fOnlyOneEtaSide(0),
207fInjectedSignals(kFALSE),
208 fRemoveWeakDecays(kFALSE),
209fRemoveDuplicates(kFALSE),
210 fapplyTrigefficiency(kFALSE),
211 fapplyAssoefficiency(kFALSE),
212 ffillefficiency(kFALSE),
213 fmesoneffrequired(kFALSE),
214 fkaonprotoneffrequired(kFALSE),
215//fAnalysisUtils(0x0),
216 fDCAXYCut(0)
217
218{
219 for ( Int_t i = 0; i < 16; i++) {
220 fHistQA[i] = NULL;
221 }
222
223 for ( Int_t i = 0; i < 6; i++ ){
224 fTHnrecomatchedallPid[i] = NULL;
225 fTHngenprimPidTruth[i] = NULL;
226 effcorection[i]=NULL;
227 //effmap[i]=NULL;
228
229 }
230
231 for(Int_t ipart=0;ipart<NSpecies;ipart++)
232 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
233 fnsigmas[ipart][ipid]=999.;
234
235 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
236
237 }
238//________________________________________________________________________
239AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
240 :AliAnalysisTaskSE(name),
241 fOutput(0),
242 fCentralityMethod("V0A"),
243 fSampleType("pPb"),
244 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
245 trkVtx(0),
246 zvtx(0),
247 fFilterBit(768),
ef349d3c 248 fSharedClusterCut(-1),
249 fVertextype(1),
c683985a 250 fzvtxcut(10.0),
251 ffilltrigassoUNID(kFALSE),
252 ffilltrigUNIDassoID(kFALSE),
253 ffilltrigIDassoUNID(kTRUE),
254 ffilltrigIDassoID(kFALSE),
255 ffilltrigIDassoIDMCTRUTH(kFALSE),
ef349d3c 256 fMaxNofMixingTracks(50000),
0684fc6a 257 fPtOrderMCTruth(kTRUE),
258 fPtOrderDataReco(kTRUE),
c683985a 259 fTriggerSpeciesSelection(kFALSE),
260 fAssociatedSpeciesSelection(kFALSE),
261 fTriggerSpecies(SpPion),
262 fAssociatedSpecies(SpPion),
263 fCustomBinning(""),
264 fBinningString(""),
ef349d3c 265 fSelectHighestPtTrig(kFALSE),
c683985a 266 fcontainPIDtrig(kTRUE),
267 fcontainPIDasso(kFALSE),
268 // frejectPileUp(kFALSE),
269 fminPt(0.2),
270 fmaxPt(10.0),
271 fmineta(-0.8),
272 fmaxeta(0.8),
273 fminprotonsigmacut(-6.0),
274 fmaxprotonsigmacut(-3.0),
275 fminpionsigmacut(0.0),
276 fmaxpionsigmacut(4.0),
277 fselectprimaryTruth(kTRUE),
278 fonlyprimarydatareco(kFALSE),
279 fdcacut(kFALSE),
280 fdcacutvalue(3.0),
281 ffillhistQAReco(kFALSE),
282 ffillhistQATruth(kFALSE),
283 kTrackVariablesPair(0),
284 fminPtTrig(0),
285 fmaxPtTrig(0),
286 fminPtComboeff(2.0),
287 fmaxPtComboeff(4.0),
288 fminPtAsso(0),
289 fmaxPtAsso(0),
290 fhistcentrality(0),
291 fEventCounter(0),
292 fEtaSpectrasso(0),
293 fphiSpectraasso(0),
294 MCtruthpt(0),
295 MCtrutheta(0),
296 MCtruthphi(0),
297 MCtruthpionpt(0),
298 MCtruthpioneta(0),
299 MCtruthpionphi(0),
300 MCtruthkaonpt(0),
301 MCtruthkaoneta(0),
302 MCtruthkaonphi(0),
303 MCtruthprotonpt(0),
304 MCtruthprotoneta(0),
305 MCtruthprotonphi(0),
306 fPioncont(0),
307 fKaoncont(0),
308 fProtoncont(0),
83cdcf92 309 fEventno(0),
310 fEventnobaryon(0),
311 fEventnomeson(0),
0684fc6a 312 fhistJetTrigestimate(0),
c683985a 313 fCentralityCorrelation(0x0),
314 fHistoTPCdEdx(0x0),
315 fHistoTOFbeta(0x0),
316 fTPCTOFPion3d(0),
317 fTPCTOFKaon3d(0),
318 fTPCTOFProton3d(0),
319 fPionPt(0),
320 fPionEta(0),
321 fPionPhi(0),
322 fKaonPt(0),
323 fKaonEta(0),
324 fKaonPhi(0),
325 fProtonPt(0),
326 fProtonEta(0),
327 fProtonPhi(0),
328 fCorrelatonTruthPrimary(0),
329 fCorrelatonTruthPrimarymix(0),
330 fTHnCorrUNID(0),
331 fTHnCorrUNIDmix(0),
332 fTHnCorrID(0),
333 fTHnCorrIDmix(0),
334 fTHnCorrIDUNID(0),
335 fTHnCorrIDUNIDmix(0),
336 fTHnTrigcount(0),
337 fTHnTrigcountMCTruthPrim(0),
338 fPoolMgr(0x0),
339 fArrayMC(0),
4a2a9605 340 fAnalysisType("AOD"),
c683985a 341 fefffilename(""),
342 twoTrackEfficiencyCutValue(0.02),
343//fControlConvResoncances(0),
344 fPID(NULL),
345 eventno(0),
346 fPtTOFPIDmin(0.6),
347 fPtTOFPIDmax(4.0),
348 fRequestTOFPID(kTRUE),
349 fPIDType(NSigmaTPCTOF),
83cdcf92 350 fFIllPIDQAHistos(kTRUE),
c683985a 351 fNSigmaPID(3),
352 fUseExclusiveNSigma(kFALSE),
353 fRemoveTracksT0Fill(kFALSE),
354fSelectCharge(0),
355fTriggerSelectCharge(0),
356fAssociatedSelectCharge(0),
357fTriggerRestrictEta(-1),
358fEtaOrdering(kFALSE),
359fCutConversions(kFALSE),
360fCutResonances(kFALSE),
361fRejectResonanceDaughters(-1),
362 fOnlyOneEtaSide(0),
363fInjectedSignals(kFALSE),
364 fRemoveWeakDecays(kFALSE),
365fRemoveDuplicates(kFALSE),
366 fapplyTrigefficiency(kFALSE),
367 fapplyAssoefficiency(kFALSE),
368 ffillefficiency(kFALSE),
369 fmesoneffrequired(kFALSE),
370 fkaonprotoneffrequired(kFALSE),
371 //fAnalysisUtils(0x0),
372 fDCAXYCut(0)
373{
374
375 for ( Int_t i = 0; i < 16; i++) {
376 fHistQA[i] = NULL;
377 }
378
379for ( Int_t i = 0; i < 6; i++ ){
380 fTHnrecomatchedallPid[i] = NULL;
381 fTHngenprimPidTruth[i] = NULL;
382 effcorection[i]=NULL;
383 //effmap[i]=NULL;
384
385 }
386
387 for(Int_t ipart=0;ipart<NSpecies;ipart++)
388 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
389 fnsigmas[ipart][ipid]=999.;
390
391 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
392
393 // The last in the above list should not have a comma after it
394
395 // Constructor
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
399
400 DefineOutput(1, TList::Class()); // for output list
83cdcf92 401 DefineOutput(2, TList::Class());
c683985a 402
403}
404
405//________________________________________________________________________
406AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
407{
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()) {
411 delete fOutput;
412
413 }
83cdcf92 414
415if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
416 delete fOutputList;
417
418 }
419
c683985a 420 if (fPID) delete fPID;
421
422 }
423//________________________________________________________________________
83cdcf92 424
425//////////////////////////////////////////////////////////////////////////////////////////////////
426
427TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
428 // returns histo named name
429 return (TH2F*) fOutputList->FindObject(name);
430}
431
432//////////////////////////////////////////////////////////////////////////////////////////////////
433
c683985a 434Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
435
436{
437 //
438 // Puts the argument in the range [-pi/2,3 pi/2].
439 //
440
441 if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
442 if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
443
444 return DPhi;
445
446}
447//________________________________________________________________________
448void AliTwoParticlePIDCorr::UserCreateOutputObjects()
449{
450 // Create histograms
451 // Called once (on the worker node)
452 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
453 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
454 fPID = inputHandler->GetPIDResponse();
455
456 //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
457
458//get the efficiency correction map
459
83cdcf92 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);
c683985a 464
465 fOutput = new TList();
466 fOutput->SetOwner(); // IMPORTANT!
467
83cdcf92 468 fOutputList = new TList;
469 fOutputList->SetOwner();
470 fOutputList->SetName("PIDQAList");
471
472
c683985a 473 Int_t centmultbins=10;
474 Double_t centmultmin=0.0;
475 Double_t centmultmax=100.0;
476 if(fSampleType=="pPb" || fSampleType=="PbPb")
477 {
478 centmultbins=10;
479 centmultmin=0.0;
480 centmultmax=100.0;
481 }
482 if(fSampleType=="pp")
483 {
484 centmultbins=10;
485 centmultmin=0.0;
486 centmultmax=200.0;
487 }
488
489fhistcentrality=new TH1F("fhistcentrality","fhistcentrality",centmultbins*4,centmultmin,centmultmax);
490fOutput->Add(fhistcentrality);
491
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);
500
501fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
502fOutput->Add(fEtaSpectrasso);
503
504fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
505fOutput->Add(fphiSpectraasso);
506
c683985a 507 if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
508 fOutput->Add(fCentralityCorrelation);
509 }
510
511fHistoTPCdEdx = new TH2F("hHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
83cdcf92 512fOutputList->Add(fHistoTPCdEdx);
bb0bbd58 513fHistoTOFbeta = new TH2F(Form("hHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
83cdcf92 514 fOutputList->Add(fHistoTOFbeta);
c683985a 515
516 fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
83cdcf92 517 fOutputList->Add(fTPCTOFPion3d);
c683985a 518
519 fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
83cdcf92 520 fOutputList->Add(fTPCTOFKaon3d);
c683985a 521
522 fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
83cdcf92 523 fOutputList->Add(fTPCTOFProton3d);
c683985a 524
525if(ffillhistQAReco)
526 {
527 fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);
83cdcf92 528 fOutputList->Add(fPionPt);
c683985a 529 fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);
83cdcf92 530 fOutputList->Add(fPionEta);
c683985a 531 fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);
83cdcf92 532 fOutputList->Add(fPionPhi);
c683985a 533
534 fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);
83cdcf92 535 fOutputList->Add(fKaonPt);
c683985a 536 fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);
83cdcf92 537 fOutputList->Add(fKaonEta);
c683985a 538 fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);
83cdcf92 539 fOutputList->Add(fKaonPhi);
c683985a 540
541 fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);
83cdcf92 542 fOutputList->Add(fProtonPt);
c683985a 543 fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);
83cdcf92 544 fOutputList->Add(fProtonEta);
c683985a 545 fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);
83cdcf92 546 fOutputList->Add(fProtonPhi);
c683985a 547 }
548
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);
565
566for(Int_t i = 0; i < 16; i++)
567 {
568 fOutput->Add(fHistQA[i]);
569 }
570
571 kTrackVariablesPair=6 ;
572
573 if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7;
574
575 if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7;
576
577 if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8;
578
579
580// two particle histograms
83cdcf92 581 Int_t anaSteps = 1; // analysis steps
582 const char* title = "d^{2}N_{ch}/d#varphid#eta";
583
c683985a 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];
588
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"
83cdcf92 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"
c683985a 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";
598
599 if(fcontainPIDtrig){
600 defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
601 }
602 if(fcontainPIDasso){
603 defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
604 }
605 // =========================================================
606 // Customization (adopted from AliUEHistograms)
607 // =========================================================
608
609 TObjArray* lines = defaultBinningStr.Tokenize("\n");
610 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
611 {
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";
616 else
617 AliInfo(Form("Using custom binning for %s", tag.Data()));
618 }
619 delete lines;
620 fBinningString += fCustomBinning;
621
83cdcf92 622 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
c683985a 623
624 // =========================================================
625 // Now set the bins
626 // =========================================================
627
628 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
629 axisTitlePair[0] = " multiplicity";
630
631 dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
632 axisTitlePair[1] = "v_{Z} (cm)";
633
634 dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
635 axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
636
637 dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
638 axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
639
640 dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
641 axisTitlePair[4] = "#Delta#eta";
642
643 dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
644 axisTitlePair[5] = "#Delta#varphi (rad)";
645
646 if(fcontainPIDtrig && !fcontainPIDasso){
647 dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
648 axisTitlePair[6] = "PIDTrig";
649 }
650
651 if(!fcontainPIDtrig && fcontainPIDasso){
652 dBinsPair[6] = GetBinning(fBinningString, "PIDAsso", iBinPair[6]);
653 axisTitlePair[6] = "PIDAsso";
654 }
655
656if(fcontainPIDtrig && fcontainPIDasso){
657
658 dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
659 axisTitlePair[6] = "PIDTrig";
660
661 dBinsPair[7] = GetBinning(fBinningString, "PIDAsso", iBinPair[7]);
662 axisTitlePair[7] = "PIDAsso";
663 }
664
665 Int_t nEtaBin = -1;
666 Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
667
668 Int_t nPteffbin = -1;
669 Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
670
671
672 fminPtTrig=dBinsPair[2][0];
673 fmaxPtTrig=dBinsPair[2][iBinPair[2]];
674 fminPtAsso=dBinsPair[3][0];
675 fmaxPtAsso=dBinsPair[3][iBinPair[3]];
676
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
680
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]};
686
687TString Histrename;
688for(Int_t jj=0;jj<6;jj++)//PID type binning
689 {
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]);
702
703Histrename="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]);
715 }
716 }
717
83cdcf92 718//AliThns for Correlation plots(data & MC)
719
c683985a 720 if(ffilltrigassoUNID)
721 {
83cdcf92 722 fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
723for (Int_t j=0; j<kTrackVariablesPair; j++) {
724 fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
725 fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
c683985a 726 }
727 fOutput->Add(fTHnCorrUNID);
83cdcf92 728
729 fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
730for (Int_t j=0; j<kTrackVariablesPair; j++) {
731 fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
732 fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
c683985a 733 }
734 fOutput->Add(fTHnCorrUNIDmix);
735 }
736
737 if(ffilltrigIDassoID)
738 {
83cdcf92 739fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
740for (Int_t j=0; j<kTrackVariablesPair; j++) {
741 fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
742 fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
c683985a 743 }
744 fOutput->Add(fTHnCorrID);
745
83cdcf92 746fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
747for (Int_t j=0; j<kTrackVariablesPair; j++) {
748 fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
749 fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
c683985a 750 }
751 fOutput->Add(fTHnCorrIDmix);
752 }
753
754 if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
755 {
83cdcf92 756fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
757for (Int_t j=0; j<kTrackVariablesPair; j++) {
758 fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
759 fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
c683985a 760 }
761 fOutput->Add(fTHnCorrIDUNID);
762
83cdcf92 763
764fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
765for (Int_t j=0; j<kTrackVariablesPair; j++) {
766 fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
767 fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
c683985a 768 }
769 fOutput->Add(fTHnCorrIDUNIDmix);
770 }
771
772
773
774 //ThnSparse for Correlation plots(truth MC)
775 if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
83cdcf92 776
777fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
778for (Int_t j=0; j<kTrackVariablesPair; j++) {
779 fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
780 fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
c683985a 781 }
782 fOutput->Add(fCorrelatonTruthPrimary);
783
83cdcf92 784
785fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
786for (Int_t j=0; j<kTrackVariablesPair; j++) {
787 fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
788 fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
c683985a 789 }
83cdcf92 790 fOutput->Add(fCorrelatonTruthPrimarymix);
c683985a 791 }
792
c683985a 793 //binning for trigger no. counting
794
c683985a 795 Int_t* fBinst;
796 Int_t dims=3;
797 if(fcontainPIDtrig) dims=4;
c683985a 798 fBinst= new Int_t[dims];
799 for(Int_t i=0; i<3;i++)
800 {
801 fBinst[i]=iBinPair[i];
c683985a 802 }
803 if(fcontainPIDtrig){
83cdcf92 804 fBinst[3]=iBinPair[6];
c683985a 805 }
806
807 //ThSparse for trigger counting(data & reco MC)
808 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
809 {
83cdcf92 810 fTHnTrigcount = new AliTHn("fTHnTrigcount", "fTHnTrigcount", anaSteps, dims, fBinst);
c683985a 811for(Int_t i=0; i<3;i++){
83cdcf92 812 fTHnTrigcount->SetBinLimits(i, dBinsPair[i]);
813 fTHnTrigcount->SetVarTitle(i, axisTitlePair[i]);
c683985a 814 }
83cdcf92 815 if(fcontainPIDtrig)
816 {
817 fTHnTrigcount->SetBinLimits(3, dBinsPair[6]);
818 fTHnTrigcount->SetVarTitle(3, axisTitlePair[6]);
819 }
c683985a 820 fOutput->Add(fTHnTrigcount);
821 }
822
823 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
83cdcf92 824 //AliTHns for trigger counting(truth MC)
825 fTHnTrigcountMCTruthPrim = new AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", anaSteps, dims, fBinst);
c683985a 826for(Int_t i=0; i<3;i++){
83cdcf92 827 fTHnTrigcountMCTruthPrim->SetBinLimits(i, dBinsPair[i]);
828 fTHnTrigcountMCTruthPrim->SetVarTitle(i, axisTitlePair[i]);
c683985a 829 }
83cdcf92 830 if(fcontainPIDtrig)
831 {
832 fTHnTrigcountMCTruthPrim->SetBinLimits(3, dBinsPair[6]);
833 fTHnTrigcountMCTruthPrim->SetVarTitle(3, axisTitlePair[6]);
834 }
c683985a 835 fOutput->Add(fTHnTrigcountMCTruthPrim);
836 }
837
838if(fAnalysisType=="MCAOD"){
839 if(ffillhistQATruth)
840 {
841 MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
83cdcf92 842 fOutputList->Add(MCtruthpt);
c683985a 843
844 MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
83cdcf92 845 fOutputList->Add(MCtrutheta);
c683985a 846
847 MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
83cdcf92 848 fOutputList->Add(MCtruthphi);
c683985a 849
850 MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
83cdcf92 851 fOutputList->Add(MCtruthpionpt);
c683985a 852
853 MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
83cdcf92 854 fOutputList->Add(MCtruthpioneta);
c683985a 855
856 MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
83cdcf92 857 fOutputList->Add(MCtruthpionphi);
c683985a 858
859 MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
83cdcf92 860 fOutputList->Add(MCtruthkaonpt);
c683985a 861
862 MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
83cdcf92 863 fOutputList->Add(MCtruthkaoneta);
c683985a 864
865 MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
83cdcf92 866 fOutputList->Add(MCtruthkaonphi);
c683985a 867
868 MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
83cdcf92 869 fOutputList->Add(MCtruthprotonpt);
c683985a 870
871 MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
83cdcf92 872 fOutputList->Add(MCtruthprotoneta);
c683985a 873
874 MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
83cdcf92 875 fOutputList->Add(MCtruthprotonphi);
c683985a 876 }
877 fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
83cdcf92 878 fOutputList->Add(fPioncont);
c683985a 879
880 fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
83cdcf92 881 fOutputList->Add(fKaoncont);
c683985a 882
883 fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
83cdcf92 884 fOutputList->Add(fProtoncont);
c683985a 885 }
886
83cdcf92 887fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
888 fEventno->GetXaxis()->SetTitle("Centrality");
889 fEventno->GetYaxis()->SetTitle("Z_Vtx");
890fOutput->Add(fEventno);
891fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
892 fEventnobaryon->GetXaxis()->SetTitle("Centrality");
893 fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
894fOutput->Add(fEventnobaryon);
895fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
896 fEventnomeson->GetXaxis()->SetTitle("Centrality");
897 fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
898fOutput->Add(fEventnomeson);
899
0684fc6a 900fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
901fOutput->Add(fhistJetTrigestimate);
902
83cdcf92 903
c683985a 904//Mixing
905DefineEventPool();
906
907 if(fapplyTrigefficiency || fapplyAssoefficiency)
908 {
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]};
913
914 TString Histrexname;
915for(Int_t jj=0;jj<6;jj++)// PID type binning
916 {
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]);
929 }
930// TFile *fsifile = new TFile(fefffilename,"READ");
4a2a9605 931
932 if (TString(fefffilename).BeginsWith("alien:"))
933 TGrid::Connect("alien:");
c683985a 934 TFile *fileT=TFile::Open(fefffilename);
935 TString Nameg;
936for(Int_t jj=0;jj<6;jj++)//type binning
937 {
938Nameg="effmap";Nameg+=jj;
939//effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
940effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
941
942//effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
943 }
944//fsifile->Close();
945fileT->Close();
946
947 }
948
949//fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
950// fOutput->Add(fControlConvResoncances);
951
952
83cdcf92 953
954
955
956 //*****************************************************PIDQA histos*****************************************************//
957
958
959 //nsigma plot
960 for(Int_t ipart=0;ipart<NSpecies;ipart++){
961 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
962 Double_t miny=-30;
963 Double_t maxy=30;
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);
969 }
970 }
971
972 //nsigmaRec plot
973 for(Int_t ipart=0;ipart<NSpecies;ipart++){
974 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
975 Double_t miny=-10;
976 Double_t maxy=10;
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);
983 }
984 }
985
986 //nsigmaDC plot
987 if(fUseExclusiveNSigma) {
988 for(Int_t ipart=0;ipart<NSpecies;ipart++){
989 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
990 Double_t miny=-10;
991 Double_t maxy=10;
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);
998 }
999 }
1000 }
1001 //nsigmaMC plot
1002 if (fAnalysisType == "MCAOD"){
1003 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1004 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1005 Double_t miny=-30;
1006 Double_t maxy=30;
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);
1013 }
1014 }
1015 }
1016 //PID signal plot
1017 for(Int_t idet=0;idet<fNDetectors;idet++){
1018 for(Int_t ipart=0;ipart<NSpecies;ipart++){
1019 Double_t maxy=500;
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);
1025 }
1026 }
1027 //PID signal plot, before PID cut
1028 for(Int_t idet=0;idet<fNDetectors;idet++){
1029 Double_t maxy=500;
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);
1035 }
1036
c683985a 1037 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
83cdcf92 1038 PostData(2, fOutputList);
1039
1040 AliInfo("Finished setting up the Output");
1041
1042 TH1::AddDirectory(oldStatus);
1043
1044
1045
c683985a 1046}
1047//-------------------------------------------------------------------------------
1048void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1049
1050
1051 if(fAnalysisType == "AOD") {
1052
1053 doAODevent();
1054
1055 }//AOD--analysis-----
1056
1057 else if(fAnalysisType == "MCAOD") {
1058
1059 doMCAODevent();
1060
1061 }
1062
1063 else return;
1064
1065}
1066//-------------------------------------------------------------------------
1067void AliTwoParticlePIDCorr::doMCAODevent()
1068{
1069 AliVEvent *event = InputEvent();
1070 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1071 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1072 if (!aod) {
1073 AliError("Cannot get the AOD event");
1074 return;
1075 }
1076
1077// count all events(physics triggered)
1078 fEventCounter->Fill(1);
1079
1080 // get centrality object and check quality(valid for p-Pb and Pb-Pb)
1081 Double_t cent_v0=0.0;
1082
1083 if(fSampleType=="pPb" || fSampleType=="PbPb")
1084 {
1085 AliCentrality *centrality=0;
1086 if(aod)
1087 centrality = aod->GetHeader()->GetCentralityP();
1088 // if (centrality->GetQuality() != 0) return ;
1089
1090 if(centrality)
1091 {
1092 cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
1093 }
1094 else
1095 {
1096 cent_v0= -1;
1097 }
1098 }
1099
1100//check the PIDResponse handler
1101 if (!fPID) return;
1102
1103// get mag. field required for twotrack efficiency cut
1104 Float_t bSign = 0;
1105 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1106
1107 //check for TClonesArray(truth track MC information)
1108fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1109 if (!fArrayMC) {
1110 AliFatal("Error: MC particles branch not found!\n");
1111 return;
1112 }
1113
1114 //check for AliAODMCHeader(truth event MC information)
1115 AliAODMCHeader *header=NULL;
1116 header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1117 if(!header) {
1118 printf("MC header branch not found!\n");
1119 return;
1120 }
1121
1122//Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
1123Float_t zVtxmc =header->GetVtxZ();
1124 if(TMath::Abs(zVtxmc)>fzvtxcut) return;
1125
1126 // For productions with injected signals, figure out above which label to skip particles/tracks
1127 Int_t skipParticlesAbove = 0;
1128
1129 if (fInjectedSignals)
1130 {
1131 AliGenEventHeader* eventHeader = 0;
1132 Int_t headers = 0;
1133
1134// AOD
1135 if (!header)
1136 AliFatal("fInjectedSignals set but no MC header found");
1137
1138 headers = header->GetNCocktailHeaders();
1139 eventHeader = header->GetCocktailHeader(0);
1140
1141 if (!eventHeader)
1142 {
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);
1147 return;
1148 }
1149skipParticlesAbove = 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));
1151 }
1152
1153 //vertex selection(is it fine for PP?)
ef349d3c 1154 if ( fVertextype==1){
c683985a 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;
1168 }
ef349d3c 1169 else if(fVertextype==2) {//for pp and pb-pb case
c683985a 1170 Int_t nVertex = aod->GetNumberOfVertices();
1171 if( nVertex > 0 ) {
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;
1179
1180 // Select a quality vertex by number of tracks?
1181 if( nTracksPrim < fnTracksVertex ) {
1182 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1183 return ;
1184 }
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)
1187 // return kFALSE;
1188 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1189 }
1190 else return;
1191
1192 }
83cdcf92 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();
1197 Double32_t fCov[6];
1198 trkVtx->GetCovarianceMatrix(fCov);
1199 if(fCov[5] == 0) return;//proper vertex resolution
1200 }
1201 else {
1202 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
1203 return;//as there is no proper sample type
1204 }
c683985a 1205
1206
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
1208
1209 //count events having a proper vertex
1210 fEventCounter->Fill(2);
1211
1212 if (TMath::Abs(zvtx) > fzvtxcut) return;
1213
1214fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut for trkVtx only
1215
1216//now we have events passed physics trigger, centrality,eventzvtx cut
1217
1218//count events after vertex cut
1219 fEventCounter->Fill(5);
1220
1221if(!aod) return; //for safety
1222
1223if (fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
1224
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)
1226
1227 Double_t cent_v0_truth=0.0;
1228
1229 //TObjArray* tracksMCtruth=0;
1230TObjArray* 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!
1232
1233 eventno++;
1234
1235 //There is a small difference between zvtx and zVtxmc??????
1236 //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1237 //cout<<"***********************************************zvtx="<<zvtx<<endl;
1238
1239//now process the truth particles(for both efficiency & correlation function)
1240Int_t nMCTrack = fArrayMC->GetEntriesFast();
1241
1242for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
1243{ //MC truth track loop starts
1244
1245AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1246
1247 if(!partMC){
1248 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1249 continue;
1250 }
1251
1252//consider only charged particles
1253 if(partMC->Charge() == 0) continue;
1254
1255//consider only primary particles; neglect all secondary particles including from weak decays
1256 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1257
1258
1259//remove injected signals(primaries above <maxLabel>)
1260 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1261
1262//remove duplicates
1263 Bool_t isduplicate=kFALSE;
1264 if (fRemoveDuplicates)
1265 {
1266 for (Int_t j=iMC+1; j<nMCTrack; ++j)
1267 {//2nd trutuh loop starts
1268AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1269 if(!partMC2){
1270 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1271 continue;
1272 }
1273 if (partMC->GetLabel() == partMC2->GetLabel())
1274 {
1275isduplicate=kTRUE;
1276 break;
1277 }
1278 }//2nd truth loop ends
1279 }
1280 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1281
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;
1285
1286 if(!partMC) continue;//for safety
1287
1288 //To determine multiplicity in case of PP
1289 nooftrackstruth++;
1290 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1291//only physical primary(all/unidentified)
1292if(ffillhistQATruth)
1293 {
1294 MCtruthpt->Fill(partMC->Pt());
1295 MCtrutheta->Fill(partMC->Eta());
1296 MCtruthphi->Fill(partMC->Phi());
1297 }
1298 //get particle ID
1299Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1300Int_t particletypeTruth=-999;
1301 if (TMath::Abs(pdgtruth)==211)
1302 {
1303 particletypeTruth=SpPion;
1304if(ffillhistQATruth)
1305 {
1306 MCtruthpionpt->Fill(partMC->Pt());
1307 MCtruthpioneta->Fill(partMC->Eta());
1308 MCtruthpionphi->Fill(partMC->Phi());
1309 }
1310 }
1311 if (TMath::Abs(pdgtruth)==321)
1312 {
1313 particletypeTruth=SpKaon;
1314if(ffillhistQATruth)
1315 {
1316 MCtruthkaonpt->Fill(partMC->Pt());
1317 MCtruthkaoneta->Fill(partMC->Eta());
1318 MCtruthkaonphi->Fill(partMC->Phi());
1319 }
1320 }
1321if(TMath::Abs(pdgtruth)==2212)
1322 {
1323 particletypeTruth=SpProton;
1324if(ffillhistQATruth)
1325 {
1326 MCtruthprotonpt->Fill(partMC->Pt());
1327 MCtruthprotoneta->Fill(partMC->Eta());
1328 MCtruthprotonphi->Fill(partMC->Phi());
1329 }
1330 }
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)
1332
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
1335
1336 Double_t primmctruth[4] = {cent_v0, zVtxmc,partMC->Pt(), partMC->Eta()};
1337 if(ffillefficiency)
1338 {
1339 fTHngenprimPidTruth[5]->Fill(primmctruth);//for all primary truth particles(4)
1340if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[3]->Fill(primmctruth);//for primary truth mesons(3)
1341if (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
1345 }
1346
1347 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
ef349d3c 1348if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
c683985a 1349 {
1350LRCParticlePID* 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
1354 }
1355 }//MC truth track loop ends
1356
1357//*********************still in event loop
ef349d3c 1358 Float_t weghtval=1.0;
c683985a 1359
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
1362
1363 //now cent_v0_truth should be used for all correlation function calculation
1364
1365if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1366 {
1367 //Fill Correlations for MC truth particles(same event)
1368if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
ef349d3c 1369 Fillcorrelation(tracksMCtruth,0,cent_v0_truth,zVtxmc,weghtval,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
c683985a 1370
1371//start mixing
1372AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0_truth, zVtxmc+200);
1373if (pool2 && pool2->IsReady())
1374 {//start mixing only when pool->IsReady
1375if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1376 {//proceed only when no. of trigger particles >0 in current event
ef349d3c 1377 Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
c683985a 1378for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1379 { //pool event loop start
1380 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1381 if(!bgTracks6) continue;
ef349d3c 1382 Fillcorrelation(tracksMCtruth,bgTracks6,cent_v0_truth,zVtxmc,nmix,bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
c683985a 1383
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
1387
1388 //still in main event loop
1389
1390 if(tracksMCtruth){
1391if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1392 }
1393 }
1394
1395 //still in main event loop
1396
1397if(tracksMCtruth) delete tracksMCtruth;
1398
1399//now deal with reco tracks
1400 //TObjArray* tracksUNID=0;
1401 TObjArray* tracksUNID = new TObjArray;
1402 tracksUNID->SetOwner(kTRUE);
1403
1404 //TObjArray* tracksID=0;
1405 TObjArray* tracksID = new TObjArray;
1406 tracksID->SetOwner(kTRUE);
1407
1408
1409 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1410
1411 Double_t trackscount=0.0;
1412
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
1421
1422//remove injected signals
1423 if(fInjectedSignals)
1424 {
1425 AliAODMCParticle* mother = recomatched;
1426
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)
1430 {
1431 mother = 0;
1432 break;
1433 }
1434
1435 mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1436 if (!mother)
1437 break;
1438 }
1439 if (!mother)
1440 {
1441 Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
1442 continue;
1443 }
1444 if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
1445 }//remove injected signals
1446
1447 if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
1448
1449 Bool_t isduplicate2=kFALSE;
1450if (fRemoveDuplicates)
1451 {
1452 for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
1453 {//2nd loop starts
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())));
1457if(!recomatched2) continue;
1458
1459if (track->GetLabel() == track2->GetLabel())
1460 {
1461isduplicate2=kTRUE;
1462 break;
1463 }
1464 }//2nd loop ends
1465 }
1466 if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
1467
1468 fHistQA[11]->Fill(track->GetTPCNcls());
1469 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
1470
1471 if(tracktype==0) continue;
1472 if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
1473{
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)
1476 trackscount++;
1477
1478//check for eta , phi holes
1479 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1480 fphiSpectraasso->Fill(track->Phi(),track->Pt());
1481
1482 Int_t particletypeMC=-9999;
1483
1484//tag all particles as unidentified
1485 particletypeMC=unidentified;
1486
1487 Float_t effmatrix=1.;
1488
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)
1492
1493 //Clone & Reduce track list(TObjArray) for unidentified particles
ef349d3c 1494if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
c683985a 1495 {
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)
1501 }
1502 Double_t allrecomatchedpid[4] = {cent_v0, zVtxmc,recomatched->Pt(), recomatched->Eta()};
1503if(ffillefficiency) fTHnrecomatchedallPid[5]->Fill(allrecomatchedpid);//for all
1504
1505 //now start the particle identification process:)
1506
1507//get the pdg code of the corresponding truth particle
1508 Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
1509
83cdcf92 1510switch(TMath::Abs(pdgCode)){
1511 case 2212:
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]);
1517 }
1518 }
1519 break;
1520 case 321:
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]);
1526 }
1527 }
1528 break;
1529 case 211:
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]);
1535 }
1536 }
1537 break;
1538 }
1539
1540
c683985a 1541Float_t dEdx = track->GetTPCsignal();
1542 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1543
1544 if(HasTOFPID(track))
1545{
bb0bbd58 1546Double_t beta = GetBeta(track);
c683985a 1547fHistoTOFbeta->Fill(track->Pt(), beta);
1548 }
1549
1550 //do track identification(nsigma method)
83cdcf92 1551 particletypeMC=GetParticle(track,fFIllPIDQAHistos);//******************************problem is here
c683985a 1552
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]);
1558 }
1559 if(particletypeMC==SpUndefined) continue;
1560
1561 //Pt, Eta , Phi distribution of the reconstructed identified particles
1562if(ffillhistQAReco)
1563 {
1564if (particletypeMC==SpPion)
1565 {
1566 fPionPt->Fill(track->Pt());
1567 fPionEta->Fill(track->Eta());
1568 fPionPhi->Fill(track->Phi());
1569 }
1570if (particletypeMC==SpKaon)
1571 {
1572 fKaonPt->Fill(track->Pt());
1573 fKaonEta->Fill(track->Eta());
1574 fKaonPhi->Fill(track->Phi());
1575 }
1576if (particletypeMC==SpProton)
1577 {
1578 fProtonPt->Fill(track->Pt());
1579 fProtonEta->Fill(track->Eta());
1580 fProtonPhi->Fill(track->Phi());
1581 }
1582 }
1583
1584 //for misidentification fraction calculation(do it with tuneonPID)
1585 if(particletypeMC==SpPion )
1586 {
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());
1591 }
1592if(particletypeMC==SpKaon )
1593 {
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());
1598 }
1599 if(particletypeMC==SpProton )
1600 {
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());
1605 }
1606
1607 //fill tracking efficiency
1608 if(ffillefficiency)
1609 {
1610 if(particletypeMC==SpPion || particletypeMC==SpKaon)
1611 {
1612if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[3]->Fill(allrecomatchedpid);//for mesons
1613 }
1614 if(particletypeMC==SpKaon || particletypeMC==SpProton)
1615 {
1616if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[4]->Fill(allrecomatchedpid);//for kaons+protons
1617 }
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
1621 }
1622
ef349d3c 1623if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
c683985a 1624 {
1625if (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);
1630 }
1631 }// if(tracktype==1) condition structure ands
1632
1633}//reco track loop ends
1634
1635 //*************************************************************still in event loop
1636
1637//same event delta-eta-deltaphi plot
1638if(fSampleType=="pp") cent_v0=trackscount;//multiplicity
1639
1640if(trackscount>0.0)
1641 {
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
1644
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??????
1646
1647//count selected events having centrality betn 0-100%
1648 fEventCounter->Fill(6);
1649
83cdcf92 1650//***************************************event no. counting
1651Bool_t isbaryontrig=kFALSE;
1652Bool_t ismesontrig=kFALSE;
1653if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
1654
1655if(tracksID && tracksID->GetEntriesFast()>0)
1656 {
1657for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
1658 { //trigger loop starts
1659 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
1660 if(!trig) continue;
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;
1665 }//trig loop ends
1666 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
1667 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
1668 }
1669
c683985a 1670 //same event delte-eta, delta-phi plot
1671if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
1672 {//same event calculation starts
0684fc6a 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)
c683985a 1675 }
1676
1677if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
1678 {//same event calculation starts
0684fc6a 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)
c683985a 1681 }
1682
1683//still in main event loop
1684//start mixing
1685 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
1686AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
1687if (pool && pool->IsReady())
1688 {//start mixing only when pool->IsReady
ef349d3c 1689 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
c683985a 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
0684fc6a 1695 Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,nmix1,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
c683985a 1696 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
0684fc6a 1697 Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,nmix1,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
c683985a 1698 }// pool event loop ends mixing case
1699
1700} //if pool->IsReady() condition ends mixing case
1701 if(tracksUNID) {
1702if(pool)
1703 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
1704 }
1705 }//mixing with unidentified particles
1706
1707 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
1708AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
1709if (pool1 && pool1->IsReady())
1710 {//start mixing only when pool->IsReady
ef349d3c 1711 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
c683985a 1712for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
1713 { //pool event loop start
1714 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
1715 if(!bgTracks2) continue;
1716if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
0684fc6a 1717 Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
c683985a 1718if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
0684fc6a 1719 Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,nmix2,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
c683985a 1720
1721 }// pool event loop ends mixing case
1722} //if pool1->IsReady() condition ends mixing case
1723
1724if(tracksID) {
1725if(pool1)
1726 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
1727 }
1728 }//mixing with identified particles
1729
1730 //no. of events analyzed
1731fEventCounter->Fill(7);
1732 }
1733
1734if(tracksUNID) delete tracksUNID;
1735
1736if(tracksID) delete tracksID;
1737
1738
1739PostData(1, fOutput);
1740
1741}
1742//________________________________________________________________________
1743void AliTwoParticlePIDCorr::doAODevent()
1744{
1745
1746 //get AOD
1747 AliVEvent *event = InputEvent();
1748 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1749 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1750 if (!aod) {
1751 AliError("Cannot get the AOD event");
1752 return;
1753 }
1754
1755// count all events
1756 fEventCounter->Fill(1);
1757
1758 // get centrality object and check quality
1759 Double_t cent_v0=0;
1760
1761 if(fSampleType=="pPb" || fSampleType=="PbPb")
1762 {
1763 AliCentrality *centrality=0;
1764 if(aod)
1765 centrality = aod->GetHeader()->GetCentralityP();
1766 // if (centrality->GetQuality() != 0) return ;
1767
1768 if(centrality)
1769 {
1770 cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
1771 }
1772 else
1773 {
1774 cent_v0= -1;
1775 }
1776 }
1777
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
1780
1781// Pileup selection ************************************************
1782// if(frejectPileUp) //applicable only for TPC only tracks,not for hybrid tracks?.
1783 // {
1784 //if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(aod)) return;
1785 // }
1786
1787 if (!fPID) return;//this should be available with each event even if we don't do PID selection
1788
1789 //vertex selection(is it fine for PP?)
ef349d3c 1790 if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
c683985a 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;
1804 }
ef349d3c 1805 else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
c683985a 1806 Int_t nVertex = aod->GetNumberOfVertices();
1807 if( nVertex > 0 ) {
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;
1815
1816 // Select a quality vertex by number of tracks?
1817 if( nTracksPrim < fnTracksVertex ) {
1818 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1819 return ;
1820 }
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)
1823 // return kFALSE;
1824 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1825 }
1826 else return;
1827
1828 }
83cdcf92 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();
1833 Double32_t fCov[6];
1834 trkVtx->GetCovarianceMatrix(fCov);
1835 if(fCov[5] == 0) return;//proper vertex resolution
1836 }
1837 else {
1838 AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
1839 return;//as there is no proper sample type
1840 }
c683985a 1841
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
1843
1844//count events having a proper vertex
1845 fEventCounter->Fill(2);
1846
1847 if (TMath::Abs(zvtx) > fzvtxcut) return;
1848
1849//count events after vertex cut
1850 fEventCounter->Fill(5);
1851
1852
1853 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
1854
1855 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
1856
1857
1858 if(!aod) return; //for safety
1859
1860if(fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
1861
1862 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
1863 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
1864
1865 TObjArray* tracksID= new TObjArray;//only pions, kaons,protonsi.e. after doing the PID selection
1866 tracksID->SetOwner(kTRUE); // IMPORTANT!
1867
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)
1869
1870 eventno++;
1871
0684fc6a 1872 Bool_t fTrigPtmin1=kFALSE;
1873 Bool_t fTrigPtmin2=kFALSE;
1874 Bool_t fTrigPtJet=kFALSE;
1875
c683985a 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;
1884
1885 if(!track) continue;//for safety
1886
1887//check for eta , phi holes
1888 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1889 fphiSpectraasso->Fill(track->Phi(),track->Pt());
1890
1891 trackscount++;
1892
1893 //if no applyefficiency , set the eff factor=1.0
1894 Float_t effmatrix=1.0;
1895
1896 //tag all particles as unidentified that passed filterbit & kinematic cuts
1897 particletype=unidentified;
1898
1899
0684fc6a 1900 if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
1901 if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
1902 if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
1903
1904
c683985a 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
1906
1907
1908 //to reduce memory consumption in pool
ef349d3c 1909 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
c683985a 1910 {
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)
1917 }
1918
1919//now start the particle identificaion process:)
1920
1921//track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
1922
1923 Float_t dEdx = track->GetTPCsignal();
1924 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1925
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))
1928{
bb0bbd58 1929 Double_t beta = GetBeta(track);
c683985a 1930 fHistoTOFbeta->Fill(track->Pt(), beta);
1931 }
bb0bbd58 1932
c683985a 1933
1934//track identification(using nsigma method)
83cdcf92 1935 particletype=GetParticle(track,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
c683985a 1936
1937
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]);
1943 }
1944
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)!!!!!!!!!!!
1947
1948 //Pt, Eta , Phi distribution of the reconstructed identified particles
1949if(ffillhistQAReco)
1950 {
1951if (particletype==SpPion)
1952 {
1953 fPionPt->Fill(track->Pt());
1954 fPionEta->Fill(track->Eta());
1955 fPionPhi->Fill(track->Phi());
1956 }
1957if (particletype==SpKaon)
1958 {
1959 fKaonPt->Fill(track->Pt());
1960 fKaonEta->Fill(track->Eta());
1961 fKaonPhi->Fill(track->Phi());
1962 }
1963if (particletype==SpProton)
1964 {
1965 fProtonPt->Fill(track->Pt());
1966 fProtonEta->Fill(track->Eta());
1967 fProtonPhi->Fill(track->Phi());
1968 }
1969 }
1970
ef349d3c 1971 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
c683985a 1972 {
1973if (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);
1978 }
1979} //track loop ends but still in event loop
1980
1981if(trackscount<1.0){
1982 if(tracksUNID) delete tracksUNID;
1983 if(tracksID) delete tracksID;
1984 return;
1985 }
1986
0684fc6a 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);
1990
ef349d3c 1991 Float_t weightval=1.0;
1992
1993
c683985a 1994// cout<<fminPtAsso<<"***"<<fmaxPtAsso<<"*********************"<<fminPtTrig<<"***"<<fmaxPtTrig<<"*****************"<<kTrackVariablesPair<<endl;
1995if(fSampleType=="pp") cent_v0=trackscount;//multiplicity
1996
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
1999
2000if(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??????
2001
2002//count selected events having centrality betn 0-100%
2003 fEventCounter->Fill(6);
2004
83cdcf92 2005//***************************************event no. counting
2006Bool_t isbaryontrig=kFALSE;
2007Bool_t ismesontrig=kFALSE;
2008if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2009
2010if(tracksID && tracksID->GetEntriesFast()>0)
2011 {
2012for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2013 { //trigger loop starts
2014 LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2015 if(!trig) continue;
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;
2020 }//trig loop ends
2021 if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2022 if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2023 }
2024
2025
c683985a 2026//same event delta-eta-deltaphi plot
2027
2028if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2029 {//same event calculation starts
0684fc6a 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)
c683985a 2032 }
2033
2034if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2035 {//same event calculation starts
0684fc6a 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)
c683985a 2038 }
2039
2040//still in main event loop
2041
2042
2043//start mixing
2044 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2045AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
2046if (pool && pool->IsReady())
2047 {//start mixing only when pool->IsReady
ef349d3c 2048 Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
c683985a 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
0684fc6a 2054 Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,nmix1,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
c683985a 2055 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
0684fc6a 2056 Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,nmix1,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
c683985a 2057 }// pool event loop ends mixing case
2058
2059} //if pool->IsReady() condition ends mixing case
2060 if(tracksUNID) {
2061if(pool)
2062 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2063 }
2064 }//mixing with unidentified particles
2065
2066
2067 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2068AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
2069if (pool1 && pool1->IsReady())
2070 {//start mixing only when pool->IsReady
ef349d3c 2071 Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
c683985a 2072for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2073 { //pool event loop start
2074 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2075 if(!bgTracks2) continue;
2076if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
0684fc6a 2077 Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
c683985a 2078if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
0684fc6a 2079 Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,nmix2,bSign,fPtOrderDataReco,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
c683985a 2080
2081 }// pool event loop ends mixing case
2082} //if pool1->IsReady() condition ends mixing case
2083
2084if(tracksID) {
2085if(pool1)
2086 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2087 }
2088 }//mixing with identified particles
2089
2090
2091 //no. of events analyzed
2092fEventCounter->Fill(7);
2093
2094//still in main event loop
2095
2096
2097if(tracksUNID) delete tracksUNID;
2098
2099if(tracksID) delete tracksID;
2100
2101
2102PostData(1, fOutput);
2103
2104} // *************************event loop ends******************************************//_______________________________________________________________________
2105TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
2106{
2107 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
2108
2109 TObjArray* tracksClone = new TObjArray;
2110 tracksClone->SetOwner(kTRUE);
2111
2112 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2113 {
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);
2118 }
2119
2120 return tracksClone;
2121}
2122
2123//--------------------------------------------------------------------------------
ef349d3c 2124void 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)
c683985a 2125{
2126
2127 //before calling this function check that either trackstrig & tracksasso are available
2128
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();
2134
2135 //if(trackstrig)
2136 Int_t jmax=trackstrig->GetEntriesFast();
2137 if(tracksasso)
2138 jmax=tracksasso->GetEntriesFast();
2139
2140// identify K, Lambda candidates and flag those particles
2141 // a TObject bit is used for this
2142const UInt_t kResonanceDaughterFlag = 1 << 14;
2143 if (fRejectResonanceDaughters > 0)
2144 {
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)
2150 {
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));
2155 }
2156
2157for (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);
2161
2162 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2163 {
2164 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2165for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
2166 {
2167 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2168
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)
2170if (trig->IsEqual(asso)) continue;
2171
2172if (trig->Charge() * asso->Charge() > 0) continue;
2173
2174 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2175
2176if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
2177 {
2178 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2179
2180 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
2181 {
2182 trig->SetBit(kResonanceDaughterFlag);
2183 asso->SetBit(kResonanceDaughterFlag);
2184
2185// Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
2186 }
2187 }
2188 }
2189 }
2190 }
2191
ef349d3c 2192 //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
c683985a 2193
ef349d3c 2194 Float_t TriggerPtMin=fminPtTrig;
2195 Float_t TriggerPtMax=fmaxPtTrig;
2196 Int_t HighestPtTriggerIndx=-99999;
2197
2198 if(fSelectHighestPtTrig)//**************add this data member to the constructor
2199 {
2200for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2201 { //trigger loop starts(highest Pt trigger selection)
2202 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2203 if(!trig) continue;
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;}
2209
2210 Float_t trigeta=trig->Eta();
2211
2212 // some optimization
2213 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2214 continue;
2215
2216if (fOnlyOneEtaSide != 0)
2217 {
2218 if (fOnlyOneEtaSide * trigeta < 0)
2219 continue;
2220 }
2221 if (fTriggerSelectCharge != 0)
2222 if (trig->Charge() * fTriggerSelectCharge < 0)
2223 continue;
2224
2225 if (fRejectResonanceDaughters > 0)
2226 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2227
2228 if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
2229 {
2230 //TriggerPt=trigpt;//*********think about it
2231
2232 HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
2233 TriggerPtMin=trigpt;
2234 }
2235
2236 }//trigger loop ends(highest Pt trigger selection)
2237
2238 }//******************SelectHighestPtTrig condition ends
2239
2240
2241 //two particle correlation filling
c683985a 2242for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2243 { //trigger loop starts
2244 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2245 if(!trig) continue;
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;}
2251
2252 Float_t trigeta=trig->Eta();
2253
2254 // some optimization
2255 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2256 continue;
2257
2258if (fOnlyOneEtaSide != 0)
2259 {
2260 if (fOnlyOneEtaSide * trigeta < 0)
2261 continue;
2262 }
2263 if (fTriggerSelectCharge != 0)
2264 if (trig->Charge() * fTriggerSelectCharge < 0)
2265 continue;
2266
2267 if (fRejectResonanceDaughters > 0)
2268 if (trig->TestBit(kResonanceDaughterFlag)) continue;
2269
ef349d3c 2270 if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
2271 if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
2272 }
2273
c683985a 2274 Float_t trigphi=trig->Phi();
2275 Float_t trackefftrig=1.0;
2276 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
4a2a9605 2277 //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
c683985a 2278 Double_t* trigval;
2279 Int_t dim=3;
2280 if(fcontainPIDtrig) dim=4;
2281 trigval= new Double_t[dim];
2282 trigval[0] = cent;
2283 trigval[1] = vtx;
2284 trigval[2] = trigpt;
2285if(fcontainPIDtrig) trigval[3] = particlepidtrig;
2286
2287 //filling only for same event case(mixcase=kFALSE)
2288
2289 //trigger particle counting only when mixcase=kFALSE i.e for same event correlation function calculation
2290if(mixcase==kFALSE)
2291 {
2292 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
83cdcf92 2293 if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
c683985a 2294 }
2295 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
83cdcf92 2296 if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
c683985a 2297 }
2298if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
83cdcf92 2299 if(fillup=="trigUNIDassoID") fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
c683985a 2300 }
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
2302if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
83cdcf92 2303 if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
c683985a 2304 }
2305 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
83cdcf92 2306 if(fillup=="trigIDassoUNID" ) fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
c683985a 2307 }
2308if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
83cdcf92 2309 if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
c683985a 2310 }
2311
83cdcf92 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!!!!
c683985a 2313 }
2314
2315 //asso loop starts within trigger loop
2316 for(Int_t j=0;j<jmax;j++)
2317 {
2318 LRCParticlePID *asso=0;
2319 if(!tracksasso)
2320 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
2321 else
2322 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2323
2324 if(!asso) continue;
2325
2326 //to avoid overflow qnd underflow
2327 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
2328
0684fc6a 2329 //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
c683985a 2330
83cdcf92 2331 if(!tracksasso && j==i) continue;
c683985a 2332
bb0bbd58 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)
c683985a 2334 // if (tracksasso && trig->IsEqual(asso)) continue;
2335
2336 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
2337
2338 if (fPtOrder)
2339 if (asso->Pt() >= trig->Pt()) continue;
2340
2341 Int_t particlepidasso=asso->getparticle();
2342 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
2343
2344
2345if (fAssociatedSelectCharge != 0)
2346if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
2347
2348 if (fSelectCharge > 0)
2349 {
2350 // skip like sign
2351 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
2352 continue;
2353
2354 // skip unlike sign
2355 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
2356 continue;
2357 }
2358
2359if (fEtaOrdering)
2360 {
2361 if (trigeta < 0 && asso->Eta() < trigeta)
2362 continue;
2363 if (trigeta > 0 && asso->Eta() > trigeta)
2364 continue;
2365 }
2366
2367if (fRejectResonanceDaughters > 0)
2368 if (asso->TestBit(kResonanceDaughterFlag))
2369 {
2370// Printf("Skipped j=%d", j);
2371 continue;
2372 }
2373
2374 // conversions
2375 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
2376 {
2377 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2378
2379 if (mass < 0.1)
2380 {
2381 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2382
2383 //fControlConvResoncances->Fill(0.0, mass);
2384
2385 if (mass < 0.04*0.04)
2386 continue;
2387 }
2388 }
2389
2390 // K0s
2391 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2392 {
2393 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
2394
2395 const Float_t kK0smass = 0.4976;
2396
2397 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
2398 {
2399 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
2400
2401 //fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
2402
2403 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
2404 continue;
2405 }
2406 }
2407
2408 // Lambda
2409 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2410 {
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);
2413
2414 const Float_t kLambdaMass = 1.115;
2415
2416 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
2417 {
2418 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
2419
2420 //fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
2421
2422 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2423 continue;
2424 }
2425 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
2426 {
2427 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2428
2429 //fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
2430
2431 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2432 continue;
2433 }
2434 }
2435
2436 if (twoTrackEfficiencyCut)
2437 {
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();
2446
2447 Float_t deta= trigeta - eta[j];
2448
2449 // optimization
2450 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
2451 {
2452
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);
2456
2457 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
2458
2459 Float_t dphistarminabs = 1e5;
2460 Float_t dphistarmin = 1e5;
2461
2462 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
2463 {
2464 for (Double_t rad=0.8; rad<2.51; rad+=0.01)
2465 {
2466 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
2467
2468 Float_t dphistarabs = TMath::Abs(dphistar);
2469
2470 if (dphistarabs < dphistarminabs)
2471 {
2472 dphistarmin = dphistar;
2473 dphistarminabs = dphistarabs;
2474 }
2475 }
2476
2477if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
2478 {
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);
2480 continue;
2481 }
2482//fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
2483
2484 }
2485 }
2486 }
2487
ef349d3c 2488 Float_t weightperevent=weight;
c683985a 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());
2494
2495 //here get the two particle efficiency correction factor
ef349d3c 2496 Float_t effweight=trackefftrig*trackeffasso*weightperevent;
2497 // if(mixcase==kFALSE) cout<<"*******************effweight="<<effweight<<endl;
c683985a 2498 Double_t* vars;
2499 vars= new Double_t[kTrackVariablesPair];
2500 vars[0]=cent;
2501 vars[1]=vtx;
2502 vars[2]=trigpt;
2503 vars[3]=asso->Pt();
2504 vars[4]=deleta;
2505 vars[5]=delphi;
2506if(fcontainPIDtrig && !fcontainPIDasso) vars[6]=particlepidtrig;
2507if(!fcontainPIDtrig && fcontainPIDasso) vars[6]=particlepidasso;
2508 if(fcontainPIDtrig && fcontainPIDasso){
2509 vars[6]=particlepidtrig;
2510 vars[7]=particlepidasso;
2511 }
2512
2513 //Fill Correlation ThnSparses
2514 if(fillup=="trigassoUNID")
2515 {
83cdcf92 2516 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,0,1.0/effweight);
2517 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
c683985a 2518 }
2519 if(fillup=="trigIDassoID")
2520 {
83cdcf92 2521 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,0,1.0/effweight);
2522 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
c683985a 2523 }
2524 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
2525 {//in this case effweight should be 1 always
83cdcf92 2526 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
2527 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
c683985a 2528 }
2529 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
2530 {
83cdcf92 2531 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
2532 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
c683985a 2533 }
2534
2535delete[] vars;
2536 }//asso loop ends
2537delete[] trigval;
2538 }//trigger loop ends
2539
2540}
2541
2542//--------------------------------------------------------------------------------
2543Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
2544{
2545 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
2546 Int_t effVars[4];
2547 Float_t effvalue=1.;
2548
2549 if(parpid==unidentified)
2550 {
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);
2556 }
2557if(parpid==SpPion || parpid==SpKaon)
2558 {
2559 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2560 {
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);
2566 }
2567 else{
2568 if(parpid==SpPion)
2569 {
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);
2575 }
2576
2577 if(parpid==SpKaon)
2578 {
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);
2584 }
2585 }
2586 }
2587
2588 if(parpid==SpProton)
2589 {
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);
2595 }
2596
2597 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2598 {
2599 if(parpid==SpProton || parpid==SpKaon)
2600 {
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);
2606 }
2607 }
2608 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
2609 if(effvalue==0.) effvalue=1.;
2610
2611 return effvalue;
2612
2613}
2614//-----------------------------------------------------------------------
2615
2616Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield)
2617{
2618
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());
2626 Float_t dxy, dz;
2627 dxy = track->DCA();
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);
2639 }
2640//accepted tracks
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
2646// DCA XY
2647 if (fdcacut && fDCAXYCut)
2648 {
2649 if (!vertex)
2650 return 0;
2651
2652 Double_t pos[2];
2653 Double_t covar[3];
2654 AliAODTrack* clone =(AliAODTrack*) track->Clone();
2655 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
2656 delete clone;
2657 if (!success)
2658 return 0;
2659
2660// Printf("%f", ((AliAODTrack*)part)->DCA());
2661// Printf("%f", pos[0]);
2662 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
2663 return 0;
2664 }
2665
ef349d3c 2666 if (fSharedClusterCut >= 0)
2667 {
2668 Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
2669 if (frac > fSharedClusterCut)
2670 return 0;
2671 }
c683985a 2672 fHistQA[8]->Fill(pt);
2673 fHistQA[9]->Fill(track->Eta());
2674 fHistQA[10]->Fill(track->Phi());
2675 return 1;
2676 }
2677 //________________________________________________________________________________
83cdcf92 2678void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
c683985a 2679{
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 )
2681Float_t pt=track->Pt();
2682
2683//it is assumed that every track that passed the filterbit have proper TPC response(!!)
2684Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
2685Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
2686Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
2687
2688Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
2689Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
2690
2691 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
2692 {
2693
2694nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
2695nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
2696nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
2697//---combined
2698nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
2699nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
2700nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
2701
2702
2703 }
2704else{
2705 // --- combined
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);
2710
2711 }
2712
2713//set data member fnsigmas
2714 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
2715 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
2716 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
2717
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;
2722
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;
2727
83cdcf92 2728 if(FIllQAHistos){
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]);
2735 }
2736 }
2737 }
c683985a 2738
2739}
2740//----------------------------------------------------------------------------
83cdcf92 2741Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
c683985a 2742{
2743 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
2744if(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.;
2747 switch (fPIDType){
2748 case NSigmaTPC:
2749 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
2750 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
2751 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
2752 break;
2753 case NSigmaTOF:
2754 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
2755 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
2756 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
2757 break;
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]) ;
2762 break;
2763 }
2764
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
83cdcf92 2768 if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
2769if(FillQAHistos){
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]);
2774 }
2775 }
2776 return SpKaon;
2777 }
2778 if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
2779 if(FillQAHistos){
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]);
2784 }
2785 }
2786return SpPion;
2787 }
2788 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
2789if(FillQAHistos){
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]);
2794 }
2795 }
2796return SpProton;
2797 }
c683985a 2798
2799// else, return undefined
2800 return SpUndefined;
2801 }
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)
2806 return SpUndefined;
2807 }
2808
2809}
2810
2811//------------------------------------------------------------------------------------------
83cdcf92 2812Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
c683985a 2813 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
2814
2815 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
2816 //fill DC histos
2817 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
2818
83cdcf92 2819 Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
c683985a 2820
2821
2822 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
2823
2824 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
2825 switch (fPIDType) {
2826 case NSigmaTPC:
2827 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
2828 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
2829 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
2830 break;
2831 case NSigmaTOF:
2832 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
2833 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
2834 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
2835 break;
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]) ;
2840 break;
2841 }
2842
2843 //there is chance of overlapping only for particles having pt below 4.0 GEv
2844 if(trk->Pt()<=4.0){
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;
2848
2849 }
2850
83cdcf92 2851if(FIllQAHistos){
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]);
2859 }
2860 }
2861 }
2862 }
2863
c683985a 2864
2865 return fHasDoubleCounting;
2866}
2867
2868//////////////////////////////////////////////////////////////////////////////////////////////////
83cdcf92 2869Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
c683985a 2870 //return the specie according to the minimum nsigma value
2871 //no double counting, this has to be evaluated using CheckDoubleCounting()
c683985a 2872
83cdcf92 2873 Int_t ID=SpUndefined;
2874
2875 CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
c683985a 2876
83cdcf92 2877 ID=FindMinNSigma(trk,FIllQAHistos);
2878
2879if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
2880 Bool_t *HasDC;
2881 HasDC=GetDoubleCounting(trk,FIllQAHistos);
2882 for(Int_t ipart=0;ipart<NSpecies;ipart++){
2883 if(HasDC[ipart]==kTRUE) ID = SpUndefined;
2884 }
2885 }
2886
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());
c683985a 2894 }
c683985a 2895 }
83cdcf92 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());
2902 }
c683985a 2903
83cdcf92 2904 return ID;
c683985a 2905}
2906
2907//-------------------------------------------------------------------------------------
2908Bool_t
2909AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
2910{
2911 // check PID signal
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.
2918
2919 return kTRUE;
2920}
2921//___________________________________________________________
2922
2923Bool_t
2924AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
2925{
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;
2935if(fRemoveTracksT0Fill)
2936 {
2937Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
2938 if (startTimeMask < 0)return kFALSE;
2939 }
2940 return kTRUE;
2941}
2942
2943//________________________________________________________________________
bb0bbd58 2944Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
c683985a 2945{
2946 //it is called only when TOF PID is available
bb0bbd58 2947//TOF beta calculation
2948 Double_t tofTime=track->GetTOFsignal();
2949
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
2955 tof=tof*c;
2956 return length/tof;
2957
2958
2959 /*
c683985a 2960 Double_t p = track->P();
2961 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
2962 Double_t timei[5];
2963 track->GetIntegratedTimes(timei);
2964 return timei[0]/time;
bb0bbd58 2965 */
c683985a 2966}
2967//------------------------------------------------------------------------------------------------------
2968
2969Float_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)
2970{
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);
2976 photon = p1+p2;
2977 photon.M()*/
2978
2979 Float_t tantheta1 = 1e10;
2980
2981 if (eta1 < -1e-10 || eta1 > 1e-10)
2982 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
2983
2984 Float_t tantheta2 = 1e10;
2985 if (eta2 < -1e-10 || eta2 > 1e-10)
2986 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
2987
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);
2990
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 ) ) );
2992
2993 return mass2;
2994}
2995//---------------------------------------------------------------------------------
2996
2997Float_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)
2998{
2999 // calculate inv mass squared approximately
3000
3001 Float_t tantheta1 = 1e10;
3002
3003 if (eta1 < -1e-10 || eta1 > 1e-10)
3004 {
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);
3007 }
3008
3009 Float_t tantheta2 = 1e10;
3010 if (eta2 < -1e-10 || eta2 > 1e-10)
3011 {
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);
3014 }
3015
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);
3018
3019 // fold onto 0...pi
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;
3025
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);
3031 else
3032 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
3033
3034 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
3035
3036// Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
3037
3038 return mass2;
3039}
3040//--------------------------------------------------------------------------------
3041Float_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)
3042{
3043 //
3044 // calculates dphistar
3045 //
3046
3047 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
3048
3049 static const Double_t kPi = TMath::Pi();
3050
3051 // circularity
3052// if (dphistar > 2 * kPi)
3053// dphistar -= 2 * kPi;
3054// if (dphistar < -2 * kPi)
3055// dphistar += 2 * kPi;
3056
3057 if (dphistar > 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;
3063
3064 return dphistar;
3065}
3066//_________________________________________________________________________
3067void AliTwoParticlePIDCorr ::DefineEventPool()
3068{
ef349d3c 3069Int_t MaxNofEvents=1000;
c683985a 3070const Int_t NofVrtxBins=10+(1+10)*2;
3071Double_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
3074 };
3075 if (fSampleType=="pp"){
ef349d3c 3076const 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....
3078fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
c683985a 3079 }
3080if (fSampleType=="pPb" || fSampleType=="PbPb")
3081 {
3082const Int_t NofCentBins=15;
3083Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
ef349d3c 3084fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
c683985a 3085 }
ef349d3c 3086fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
c683985a 3087
3088//if(!fPoolMgr) return kFALSE;
3089//return kTRUE;
3090
3091}
3092//------------------------------------------------------------------------
3093Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
3094{
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
3099 // phi: .....
3100 //
3101 // returns bin edges which have to be deleted by the caller
3102
3103 TString config(configuration);
3104 TObjArray* lines = config.Tokenize("\n");
3105 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
3106 {
3107 TString line(lines->At(i)->GetName());
3108 if (line.BeginsWith(TString(tag) + ":"))
3109 {
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();
3116
3117 nBins = binning->GetEntriesFast() - 1;
3118
3119 delete binning;
3120 delete lines;
3121 return bins;
3122 }
3123 }
3124
3125 delete lines;
3126 AliFatal(Form("Tag %s not found in %s", tag, configuration));
3127 return 0;
3128}
c683985a 3129//________________________________________________________________________
3130void AliTwoParticlePIDCorr::Terminate(Option_t *)
3131{
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; }
3136
3137
3138}
3139//------------------------------------------------------------------
3140