]> 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
Corrected end-of-line behavior
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / TriggerPID / AliTwoParticlePIDCorr.cxx
CommitLineData
896d3200 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"
27
28#include "AliCentrality.h"
29#include "Riostream.h"
30
31#include <TSpline.h>
32#include <AliPID.h>
33#include "AliESDpid.h"
34#include "AliAODpidUtil.h"
35#include <AliPIDResponse.h>
36
37#include <AliAnalysisManager.h>
38#include <AliInputEventHandler.h>
39#include "AliAODInputHandler.h"
40
41#include "AliAnalysisTaskSE.h"
42#include "AliAnalysisManager.h"
43#include "AliCentrality.h"
44
45#include "AliVEvent.h"
46#include "AliAODEvent.h"
47#include "AliAODTrack.h"
48#include "AliVTrack.h"
49
50#include "THnSparse.h"
51
52#include "AliAODMCHeader.h"
53#include "AliAODMCParticle.h"
54#include "AliMCEventHandler.h"
55#include "AliMCEvent.h"
56#include "AliMCParticle.h"
57#include "TParticle.h"
58#include <TDatabasePDG.h>
59#include <TParticlePDG.h>
60
61#include "AliGenCocktailEventHeader.h"
62#include "AliGenEventHeader.h"
63
64#include "AliEventPoolManager.h"
65//#include "AliAnalysisUtils.h"
66using namespace AliPIDNameSpace;
67using namespace std;
68
69ClassImp(AliTwoParticlePIDCorr)
70ClassImp(LRCParticlePID)
71//________________________________________________________________________
72AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
73:AliAnalysisTaskSE(),
74 fOutput(0),
75 fCentralityMethod("V0A"),
76 fSampleType("pPb"),
77fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
78 trkVtx(0),
79 zvtx(0),
80 fFilterBit(768),
81 fzvtxcut(10.0),
82 ffilltrigassoUNID(kFALSE),
83 ffilltrigUNIDassoID(kFALSE),
84 ffilltrigIDassoUNID(kTRUE),
85 ffilltrigIDassoID(kFALSE),
86 ffilltrigIDassoIDMCTRUTH(kFALSE),
87 fPtOrderMCTruth(kFALSE),
88 fTriggerSpeciesSelection(kFALSE),
89 fAssociatedSpeciesSelection(kFALSE),
90 fTriggerSpecies(SpPion),
91 fAssociatedSpecies(SpPion),
92 fCustomBinning(""),
93 fBinningString(""),
94 fcontainPIDtrig(kTRUE),
95 fcontainPIDasso(kFALSE),
96//frejectPileUp(kFALSE),
97 fminPt(0.2),
98 fmaxPt(10.0),
99 fmineta(-0.8),
100 fmaxeta(0.8),
101 fminprotonsigmacut(-6.0),
102 fmaxprotonsigmacut(-3.0),
103 fminpionsigmacut(0.0),
104 fmaxpionsigmacut(4.0),
105 fselectprimaryTruth(kTRUE),
106 fonlyprimarydatareco(kFALSE),
107 fdcacut(kFALSE),
108 fdcacutvalue(3.0),
109 ffillhistQAReco(kFALSE),
110 ffillhistQATruth(kFALSE),
111 kTrackVariablesPair(0),
112 fminPtTrig(0),
113 fmaxPtTrig(0),
114 fminPtComboeff(2.0),
115 fmaxPtComboeff(4.0),
116 fminPtAsso(0),
117 fmaxPtAsso(0),
118 fhistcentrality(0),
119 fEventCounter(0),
120 fEtaSpectrasso(0),
121 fphiSpectraasso(0),
122 MCtruthpt(0),
123 MCtrutheta(0),
124 MCtruthphi(0),
125 MCtruthpionpt(0),
126 MCtruthpioneta(0),
127 MCtruthpionphi(0),
128 MCtruthkaonpt(0),
129 MCtruthkaoneta(0),
130 MCtruthkaonphi(0),
131 MCtruthprotonpt(0),
132 MCtruthprotoneta(0),
133 MCtruthprotonphi(0),
134 fPioncont(0),
135 fKaoncont(0),
136 fProtoncont(0),
137 fCentralityCorrelation(0x0),
138 fHistoTPCdEdx(0x0),
139 fHistoTOFbeta(0x0),
140 fTPCTOFPion3d(0),
141 fTPCTOFKaon3d(0),
142 fTPCTOFProton3d(0),
143 fPionPt(0),
144 fPionEta(0),
145 fPionPhi(0),
146 fKaonPt(0),
147 fKaonEta(0),
148 fKaonPhi(0),
149 fProtonPt(0),
150 fProtonEta(0),
151 fProtonPhi(0),
152 fCorrelatonTruthPrimary(0),
153 fCorrelatonTruthPrimarymix(0),
154 fTHnCorrUNID(0),
155 fTHnCorrUNIDmix(0),
156 fTHnCorrID(0),
157 fTHnCorrIDmix(0),
158 fTHnCorrIDUNID(0),
159 fTHnCorrIDUNIDmix(0),
160 fTHnTrigcount(0),
161 fTHnTrigcountMCTruthPrim(0),
162 fPoolMgr(0x0),
163 fArrayMC(0),
164 fAnalysisType("MCAOD"),
165 fefffilename(""),
166 twoTrackEfficiencyCutValue(0.02),
167//fControlConvResoncances(0),
168 fPID(NULL),
169 eventno(0),
170 fPtTOFPIDmin(0.6),
171 fPtTOFPIDmax(4.0),
172 fRequestTOFPID(kTRUE),
173 fPIDType(NSigmaTPCTOF),
174 fNSigmaPID(3),
175 fUseExclusiveNSigma(kFALSE),
176 fRemoveTracksT0Fill(kFALSE),
177fSelectCharge(0),
178fTriggerSelectCharge(0),
179fAssociatedSelectCharge(0),
180fTriggerRestrictEta(-1),
181fEtaOrdering(kFALSE),
182fCutConversions(kFALSE),
183fCutResonances(kFALSE),
184fRejectResonanceDaughters(-1),
185 fOnlyOneEtaSide(0),
186fInjectedSignals(kFALSE),
187 fRemoveWeakDecays(kFALSE),
188fRemoveDuplicates(kFALSE),
189 fapplyTrigefficiency(kFALSE),
190 fapplyAssoefficiency(kFALSE),
191 ffillefficiency(kFALSE),
192 fmesoneffrequired(kFALSE),
193 fkaonprotoneffrequired(kFALSE),
194//fAnalysisUtils(0x0),
195 fDCAXYCut(0)
196
197{
198 for ( Int_t i = 0; i < 16; i++) {
199 fHistQA[i] = NULL;
200 }
201
202 for ( Int_t i = 0; i < 6; i++ ){
203 fTHnrecomatchedallPid[i] = NULL;
204 fTHngenprimPidTruth[i] = NULL;
205 effcorection[i]=NULL;
206 //effmap[i]=NULL;
207
208 }
209
210 for(Int_t ipart=0;ipart<NSpecies;ipart++)
211 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
212 fnsigmas[ipart][ipid]=999.;
213
214 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
215
216 }
217//________________________________________________________________________
218AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
219 :AliAnalysisTaskSE(name),
220 fOutput(0),
221 fCentralityMethod("V0A"),
222 fSampleType("pPb"),
223 fnTracksVertex(1), // QA tracks pointing to principal vertex (= 3 default)
224 trkVtx(0),
225 zvtx(0),
226 fFilterBit(768),
227 fzvtxcut(10.0),
228 ffilltrigassoUNID(kFALSE),
229 ffilltrigUNIDassoID(kFALSE),
230 ffilltrigIDassoUNID(kTRUE),
231 ffilltrigIDassoID(kFALSE),
232 ffilltrigIDassoIDMCTRUTH(kFALSE),
233 fPtOrderMCTruth(kFALSE),
234 fTriggerSpeciesSelection(kFALSE),
235 fAssociatedSpeciesSelection(kFALSE),
236 fTriggerSpecies(SpPion),
237 fAssociatedSpecies(SpPion),
238 fCustomBinning(""),
239 fBinningString(""),
240 fcontainPIDtrig(kTRUE),
241 fcontainPIDasso(kFALSE),
242 // frejectPileUp(kFALSE),
243 fminPt(0.2),
244 fmaxPt(10.0),
245 fmineta(-0.8),
246 fmaxeta(0.8),
247 fminprotonsigmacut(-6.0),
248 fmaxprotonsigmacut(-3.0),
249 fminpionsigmacut(0.0),
250 fmaxpionsigmacut(4.0),
251 fselectprimaryTruth(kTRUE),
252 fonlyprimarydatareco(kFALSE),
253 fdcacut(kFALSE),
254 fdcacutvalue(3.0),
255 ffillhistQAReco(kFALSE),
256 ffillhistQATruth(kFALSE),
257 kTrackVariablesPair(0),
258 fminPtTrig(0),
259 fmaxPtTrig(0),
260 fminPtComboeff(2.0),
261 fmaxPtComboeff(4.0),
262 fminPtAsso(0),
263 fmaxPtAsso(0),
264 fhistcentrality(0),
265 fEventCounter(0),
266 fEtaSpectrasso(0),
267 fphiSpectraasso(0),
268 MCtruthpt(0),
269 MCtrutheta(0),
270 MCtruthphi(0),
271 MCtruthpionpt(0),
272 MCtruthpioneta(0),
273 MCtruthpionphi(0),
274 MCtruthkaonpt(0),
275 MCtruthkaoneta(0),
276 MCtruthkaonphi(0),
277 MCtruthprotonpt(0),
278 MCtruthprotoneta(0),
279 MCtruthprotonphi(0),
280 fPioncont(0),
281 fKaoncont(0),
282 fProtoncont(0),
283 fCentralityCorrelation(0x0),
284 fHistoTPCdEdx(0x0),
285 fHistoTOFbeta(0x0),
286 fTPCTOFPion3d(0),
287 fTPCTOFKaon3d(0),
288 fTPCTOFProton3d(0),
289 fPionPt(0),
290 fPionEta(0),
291 fPionPhi(0),
292 fKaonPt(0),
293 fKaonEta(0),
294 fKaonPhi(0),
295 fProtonPt(0),
296 fProtonEta(0),
297 fProtonPhi(0),
298 fCorrelatonTruthPrimary(0),
299 fCorrelatonTruthPrimarymix(0),
300 fTHnCorrUNID(0),
301 fTHnCorrUNIDmix(0),
302 fTHnCorrID(0),
303 fTHnCorrIDmix(0),
304 fTHnCorrIDUNID(0),
305 fTHnCorrIDUNIDmix(0),
306 fTHnTrigcount(0),
307 fTHnTrigcountMCTruthPrim(0),
308 fPoolMgr(0x0),
309 fArrayMC(0),
310 fAnalysisType("MCAOD"),
311 fefffilename(""),
312 twoTrackEfficiencyCutValue(0.02),
313//fControlConvResoncances(0),
314 fPID(NULL),
315 eventno(0),
316 fPtTOFPIDmin(0.6),
317 fPtTOFPIDmax(4.0),
318 fRequestTOFPID(kTRUE),
319 fPIDType(NSigmaTPCTOF),
320 fNSigmaPID(3),
321 fUseExclusiveNSigma(kFALSE),
322 fRemoveTracksT0Fill(kFALSE),
323fSelectCharge(0),
324fTriggerSelectCharge(0),
325fAssociatedSelectCharge(0),
326fTriggerRestrictEta(-1),
327fEtaOrdering(kFALSE),
328fCutConversions(kFALSE),
329fCutResonances(kFALSE),
330fRejectResonanceDaughters(-1),
331 fOnlyOneEtaSide(0),
332fInjectedSignals(kFALSE),
333 fRemoveWeakDecays(kFALSE),
334fRemoveDuplicates(kFALSE),
335 fapplyTrigefficiency(kFALSE),
336 fapplyAssoefficiency(kFALSE),
337 ffillefficiency(kFALSE),
338 fmesoneffrequired(kFALSE),
339 fkaonprotoneffrequired(kFALSE),
340 //fAnalysisUtils(0x0),
341 fDCAXYCut(0)
342{
343
344 for ( Int_t i = 0; i < 16; i++) {
345 fHistQA[i] = NULL;
346 }
347
348for ( Int_t i = 0; i < 6; i++ ){
349 fTHnrecomatchedallPid[i] = NULL;
350 fTHngenprimPidTruth[i] = NULL;
351 effcorection[i]=NULL;
352 //effmap[i]=NULL;
353
354 }
355
356 for(Int_t ipart=0;ipart<NSpecies;ipart++)
357 for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
358 fnsigmas[ipart][ipid]=999.;
359
360 for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
361
362 // The last in the above list should not have a comma after it
363
364 // Constructor
365 // Define input and output slots here (never in the dummy constructor)
366 // Input slot #0 works with a TChain - it is connected to the default input container
367 // Output slot #1 writes into a TH1 container
368
369 DefineOutput(1, TList::Class()); // for output list
370
371}
372
373//________________________________________________________________________
374AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
375{
376 // Destructor. Clean-up the output list, but not the histograms that are put inside
377 // (the list is owner and will clean-up these histograms). Protect in PROOF case.
378 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
379 delete fOutput;
380
381 }
382 if (fPID) delete fPID;
383
384 }
385//________________________________________________________________________
386Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
387
388{
389 //
390 // Puts the argument in the range [-pi/2,3 pi/2].
391 //
392
393 if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
394 if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
395
396 return DPhi;
397
398}
399//________________________________________________________________________
400void AliTwoParticlePIDCorr::UserCreateOutputObjects()
401{
402 // Create histograms
403 // Called once (on the worker node)
404 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
405 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
406 fPID = inputHandler->GetPIDResponse();
407
408 //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
409
410//get the efficiency correction map
411
412
413 fOutput = new TList();
414 fOutput->SetOwner(); // IMPORTANT!
415
416 Int_t centmultbins=10;
417 Double_t centmultmin=0.0;
418 Double_t centmultmax=100.0;
419 if(fSampleType=="pPb" || fSampleType=="PbPb")
420 {
421 centmultbins=10;
422 centmultmin=0.0;
423 centmultmax=100.0;
424 }
425 if(fSampleType=="pp")
426 {
427 centmultbins=10;
428 centmultmin=0.0;
429 centmultmax=200.0;
430 }
431
432fhistcentrality=new TH1F("fhistcentrality","fhistcentrality",centmultbins*4,centmultmin,centmultmax);
433fOutput->Add(fhistcentrality);
434
435 fEventCounter = new TH1F("fEventCounter","EventCounter", 10, 0.5,10.5);
436 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
437 fEventCounter->GetXaxis()->SetBinLabel(2,"Have a vertex");
438 fEventCounter->GetXaxis()->SetBinLabel(5,"After vertex Cut");
439 fEventCounter->GetXaxis()->SetBinLabel(6,"Within 0-100% centrality");
440 fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
441 //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
442 fOutput->Add(fEventCounter);
443
444fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
445fOutput->Add(fEtaSpectrasso);
446
447fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
448fOutput->Add(fphiSpectraasso);
449
450
451 if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
452 fOutput->Add(fCentralityCorrelation);
453 }
454
455fHistoTPCdEdx = new TH2F("hHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
456fOutput->Add(fHistoTPCdEdx);
457fHistoTOFbeta = new TH2F(Form("hHistoTOFbeta"), ";p_{T} (GeV/c);v/c",70, 0., 7., 500, 0.1, 1.1);
458 fOutput->Add(fHistoTOFbeta);
459
460 fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
461 fOutput->Add(fTPCTOFPion3d);
462
463 fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
464 fOutput->Add(fTPCTOFKaon3d);
465
466 fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
467 fOutput->Add(fTPCTOFProton3d);
468
469if(ffillhistQAReco)
470 {
471 fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);
472 fOutput->Add(fPionPt);
473 fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);
474 fOutput->Add(fPionEta);
475 fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);
476 fOutput->Add(fPionPhi);
477
478 fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);
479 fOutput->Add(fKaonPt);
480 fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);
481 fOutput->Add(fKaonEta);
482 fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);
483 fOutput->Add(fKaonPhi);
484
485 fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);
486 fOutput->Add(fProtonPt);
487 fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);
488 fOutput->Add(fProtonEta);
489 fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);
490 fOutput->Add(fProtonPhi);
491 }
492
493 fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
494 fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
495 fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
496 fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx After Cut ", 50, -5., 5.);
497 fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
498 fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
499 fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
500 fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
501 fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
502 fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
503 fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
504 fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
505 fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
506 fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
507 fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
508 fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
509
510for(Int_t i = 0; i < 16; i++)
511 {
512 fOutput->Add(fHistQA[i]);
513 }
514
515 kTrackVariablesPair=6 ;
516
517 if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7;
518
519 if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7;
520
521 if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8;
522
523
524// two particle histograms
525 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
526 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
527 TString* axisTitlePair; // axis titles for track variables
528 axisTitlePair=new TString[kTrackVariablesPair];
529
530 TString defaultBinningStr;
531 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"
532 "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"
533 "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
534 "p_t_eff: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"
535 "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
536 "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
537 "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"
538 "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
539
540 if(fcontainPIDtrig){
541 defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
542 }
543 if(fcontainPIDasso){
544 defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
545 }
546 // =========================================================
547 // Customization (adopted from AliUEHistograms)
548 // =========================================================
549
550 TObjArray* lines = defaultBinningStr.Tokenize("\n");
551 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
552 {
553 TString line(lines->At(i)->GetName());
554 TString tag = line(0, line.Index(":")+1);
555 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
556 fBinningString += line + "\n";
557 else
558 AliInfo(Form("Using custom binning for %s", tag.Data()));
559 }
560 delete lines;
561 fBinningString += fCustomBinning;
562
563 AliInfo(Form("Used THn Binning:\n%s",fBinningString.Data()));
564
565 // =========================================================
566 // Now set the bins
567 // =========================================================
568
569 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
570 axisTitlePair[0] = " multiplicity";
571
572 dBinsPair[1] = GetBinning(fBinningString, "vertex", iBinPair[1]);
573 axisTitlePair[1] = "v_{Z} (cm)";
574
575 dBinsPair[2] = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
576 axisTitlePair[2] = "p_{T,trig.} (GeV/c)";
577
578 dBinsPair[3] = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
579 axisTitlePair[3] = "p_{T,assoc.} (GeV/c)";
580
581 dBinsPair[4] = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
582 axisTitlePair[4] = "#Delta#eta";
583
584 dBinsPair[5] = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
585 axisTitlePair[5] = "#Delta#varphi (rad)";
586
587 if(fcontainPIDtrig && !fcontainPIDasso){
588 dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
589 axisTitlePair[6] = "PIDTrig";
590 }
591
592 if(!fcontainPIDtrig && fcontainPIDasso){
593 dBinsPair[6] = GetBinning(fBinningString, "PIDAsso", iBinPair[6]);
594 axisTitlePair[6] = "PIDAsso";
595 }
596
597if(fcontainPIDtrig && fcontainPIDasso){
598
599 dBinsPair[6] = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
600 axisTitlePair[6] = "PIDTrig";
601
602 dBinsPair[7] = GetBinning(fBinningString, "PIDAsso", iBinPair[7]);
603 axisTitlePair[7] = "PIDAsso";
604 }
605
606 Int_t nEtaBin = -1;
607 Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
608
609 Int_t nPteffbin = -1;
610 Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
611
612
613 fminPtTrig=dBinsPair[2][0];
614 fmaxPtTrig=dBinsPair[2][iBinPair[2]];
615 fminPtAsso=dBinsPair[3][0];
616 fmaxPtAsso=dBinsPair[3][iBinPair[3]];
617
618 //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
619 //fmaxPtComboeff=fmaxPtTrig;
620//THnSparses for calculation of efficiency
621
622 if((fAnalysisType =="MCAOD") && ffillefficiency) {
623 const Int_t nDim = 4;// cent zvtx pt eta
624 Int_t fBinsCh[nDim] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//**********************change it
625 Double_t fMinCh[nDim] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
626 Double_t fMaxCh[nDim] = { dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
627
628TString Histrename;
629for(Int_t jj=0;jj<6;jj++)//PID type binning
630 {
631 Histrename="fTHnrecomatchedallPid";Histrename+=jj;
632 fTHnrecomatchedallPid[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);
633 fTHnrecomatchedallPid[jj]->Sumw2();
634 fTHnrecomatchedallPid[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
635 fTHnrecomatchedallPid[jj]->GetAxis(0)->SetTitle("Centrality");
636 fTHnrecomatchedallPid[jj]->GetAxis(1)->Set(iBinPair[1],dBinsPair[1]);
637 fTHnrecomatchedallPid[jj]->GetAxis(1)->SetTitle("zvtx");
638 fTHnrecomatchedallPid[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
639 fTHnrecomatchedallPid[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
640 fTHnrecomatchedallPid[jj]->GetAxis(3)->Set(nEtaBin,EtaBin);
641 fTHnrecomatchedallPid[jj]->GetAxis(3)->SetTitle("#eta");
642 fOutput->Add(fTHnrecomatchedallPid[jj]);
643
644Histrename="fTHngenprimPidTruth";Histrename+=jj;
645 fTHngenprimPidTruth[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);
646 fTHngenprimPidTruth[jj]->Sumw2();
647 fTHngenprimPidTruth[jj]->GetAxis(0)->Set(iBinPair[0],dBinsPair[0]);
648 fTHngenprimPidTruth[jj]->GetAxis(0)->SetTitle("Centrality");
649 fTHngenprimPidTruth[jj]->GetAxis(1)->Set(iBinPair[1], dBinsPair[1]);
650 fTHngenprimPidTruth[jj]->GetAxis(1)->SetTitle("zvtx");
651 fTHngenprimPidTruth[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
652 fTHngenprimPidTruth[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
653 fTHngenprimPidTruth[jj]->GetAxis(3)->Set(nEtaBin, EtaBin);
654 fTHngenprimPidTruth[jj]->GetAxis(3)->SetTitle("#eta");
655 fOutput->Add(fTHngenprimPidTruth[jj]);
656 }
657 }
658
659 Int_t fBins[kTrackVariablesPair];
660 Double_t fMin[kTrackVariablesPair];
661 Double_t fMax[kTrackVariablesPair];
662
663//ThnSparses for Correlation plots(data & MC)
664 for(Int_t i=0;i<kTrackVariablesPair;i++)
665 {
666 fBins[i] =iBinPair[i];
667 fMin[i]= dBinsPair[i][0];
668 fMax[i] = dBinsPair[i][iBinPair[i]];
669 }
670 if(ffilltrigassoUNID)
671 {
672 fTHnCorrUNID = new THnSparseF("fTHnCorrUNID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
673 fTHnCorrUNID->Sumw2();
674 for(Int_t i=0; i<kTrackVariablesPair;i++){
675 SetAsymmetricBin(fTHnCorrUNID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
676 }
677 fOutput->Add(fTHnCorrUNID);
678
679fTHnCorrUNIDmix = new THnSparseF("fTHnCorrUNIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
680 fTHnCorrUNIDmix->Sumw2();
681for(Int_t i=0; i<kTrackVariablesPair;i++){
682 SetAsymmetricBin(fTHnCorrUNIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
683 }
684 fOutput->Add(fTHnCorrUNIDmix);
685 }
686
687 if(ffilltrigIDassoID)
688 {
689 fTHnCorrID = new THnSparseF("fTHnCorrID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
690 fTHnCorrID->Sumw2();
691for(Int_t i=0; i<kTrackVariablesPair;i++){
692 SetAsymmetricBin(fTHnCorrID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
693 }
694 fOutput->Add(fTHnCorrID);
695
696fTHnCorrIDmix = new THnSparseF("fTHnCorrIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
697 fTHnCorrIDmix->Sumw2();
698for(Int_t i=0; i<kTrackVariablesPair;i++){
699 SetAsymmetricBin(fTHnCorrIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
700 }
701 fOutput->Add(fTHnCorrIDmix);
702 }
703
704 if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
705 {
706 fTHnCorrIDUNID = new THnSparseF("fTHnCorrIDUNID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
707 fTHnCorrIDUNID->Sumw2();
708for(Int_t i=0; i<kTrackVariablesPair;i++){
709 SetAsymmetricBin(fTHnCorrIDUNID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
710 }
711 fOutput->Add(fTHnCorrIDUNID);
712
713 fTHnCorrIDUNIDmix = new THnSparseF("fTHnCorrIDUNIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
714 fTHnCorrIDUNIDmix->Sumw2();
715for(Int_t i=0; i<kTrackVariablesPair;i++){
716 SetAsymmetricBin(fTHnCorrIDUNIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
717 }
718 fOutput->Add(fTHnCorrIDUNIDmix);
719 }
720
721
722
723 //ThnSparse for Correlation plots(truth MC)
724 if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
725 fCorrelatonTruthPrimary = new THnSparseF("fCorrelatonTruthPrimary","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
726 fCorrelatonTruthPrimary->Sumw2();
727for(Int_t i=0; i<kTrackVariablesPair;i++){
728 SetAsymmetricBin(fCorrelatonTruthPrimary,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
729 }
730 fOutput->Add(fCorrelatonTruthPrimary);
731
732 fCorrelatonTruthPrimarymix = new THnSparseF("fCorrelatonTruthPrimarymix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
733 fCorrelatonTruthPrimarymix->Sumw2();
734for(Int_t i=0; i<kTrackVariablesPair;i++){
735 SetAsymmetricBin(fCorrelatonTruthPrimarymix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
736 }
737 fOutput->Add(fCorrelatonTruthPrimarymix);
738 }
739
740
741 //binning for trigger no. counting
742
743 Double_t* fMint;
744 Double_t* fMaxt;
745 Int_t* fBinst;
746 Int_t dims=3;
747 if(fcontainPIDtrig) dims=4;
748 fMint= new Double_t[dims];
749 fMaxt= new Double_t[dims];
750 fBinst= new Int_t[dims];
751 for(Int_t i=0; i<3;i++)
752 {
753 fBinst[i]=iBinPair[i];
754 fMint[i]=dBinsPair[i][0];
755 fMaxt[i]=dBinsPair[i][iBinPair[i]];
756 }
757 if(fcontainPIDtrig){
758 fBinst[3]=iBinPair[6];
759 fMint[3]=dBinsPair[6][0];
760 fMaxt[3]=dBinsPair[6][iBinPair[6]];
761 }
762
763 //ThSparse for trigger counting(data & reco MC)
764 if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
765 {
766 fTHnTrigcount = new THnSparseF("fTHnTrigcount","cent:zvtx:pt", dims, fBinst, fMint, fMaxt);
767 fTHnTrigcount->Sumw2();
768for(Int_t i=0; i<3;i++){
769 SetAsymmetricBin(fTHnTrigcount,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
770 }
771 if(fcontainPIDtrig) SetAsymmetricBin(fTHnTrigcount,3,dBinsPair[6],iBinPair[6],axisTitlePair[6]);
772 fOutput->Add(fTHnTrigcount);
773 }
774
775 if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
776 //ThSparse for trigger counting(truth MC)
777fTHnTrigcountMCTruthPrim = new THnSparseF("fTHnTrigcountMCTruthPrim","cent:zvtx:pt", dims, fBinst, fMint, fMaxt);
778 fTHnTrigcountMCTruthPrim->Sumw2();
779for(Int_t i=0; i<3;i++){
780 SetAsymmetricBin(fTHnTrigcountMCTruthPrim,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
781 }
782if(fcontainPIDtrig) SetAsymmetricBin(fTHnTrigcountMCTruthPrim,3,dBinsPair[6],iBinPair[6],axisTitlePair[6]);
783 fOutput->Add(fTHnTrigcountMCTruthPrim);
784 }
785
786if(fAnalysisType=="MCAOD"){
787 if(ffillhistQATruth)
788 {
789 MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
790 fOutput->Add(MCtruthpt);
791
792 MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
793 fOutput->Add(MCtrutheta);
794
795 MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
796 fOutput->Add(MCtruthphi);
797
798 MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
799 fOutput->Add(MCtruthpionpt);
800
801 MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
802 fOutput->Add(MCtruthpioneta);
803
804 MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
805 fOutput->Add(MCtruthpionphi);
806
807 MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
808 fOutput->Add(MCtruthkaonpt);
809
810 MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
811 fOutput->Add(MCtruthkaoneta);
812
813 MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
814 fOutput->Add(MCtruthkaonphi);
815
816 MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
817 fOutput->Add(MCtruthprotonpt);
818
819 MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
820 fOutput->Add(MCtruthprotoneta);
821
822 MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
823 fOutput->Add(MCtruthprotonphi);
824 }
825 fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
826 fOutput->Add(fPioncont);
827
828 fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
829 fOutput->Add(fKaoncont);
830
831 fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
832 fOutput->Add(fProtoncont);
833 }
834
835//Mixing
836DefineEventPool();
837
838 if(fapplyTrigefficiency || fapplyAssoefficiency)
839 {
840 const Int_t nDimt = 4;// cent zvtx pt eta
841 Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
842 Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
843 Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
844
845 TString Histrexname;
846for(Int_t jj=0;jj<6;jj++)// PID type binning
847 {
848 Histrexname="effcorection";Histrexname+=jj;
849 effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
850 effcorection[jj]->Sumw2();
851 effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
852 effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
853 effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
854 effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
855 effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
856 effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
857 effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
858 effcorection[jj]->GetAxis(3)->SetTitle("#eta");
859 fOutput->Add(effcorection[jj]);
860 }
861// TFile *fsifile = new TFile(fefffilename,"READ");
862 TFile *fileT=TFile::Open(fefffilename);
863 TString Nameg;
864for(Int_t jj=0;jj<6;jj++)//type binning
865 {
866Nameg="effmap";Nameg+=jj;
867//effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
868effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
869
870//effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
871 }
872//fsifile->Close();
873fileT->Close();
874
875 }
876
877//fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
878// fOutput->Add(fControlConvResoncances);
879
880
881 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
882}
883//-------------------------------------------------------------------------------
884void AliTwoParticlePIDCorr::UserExec( Option_t * ){
885
886
887 if(fAnalysisType == "AOD") {
888
889 doAODevent();
890
891 }//AOD--analysis-----
892
893 else if(fAnalysisType == "MCAOD") {
894
895 doMCAODevent();
896
897 }
898
899 else return;
900
901}
902//-------------------------------------------------------------------------
903void AliTwoParticlePIDCorr::doMCAODevent()
904{
905 AliVEvent *event = InputEvent();
906 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
907 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
908 if (!aod) {
909 AliError("Cannot get the AOD event");
910 return;
911 }
912
913// count all events(physics triggered)
914 fEventCounter->Fill(1);
915
916 // get centrality object and check quality(valid for p-Pb and Pb-Pb)
917 Double_t cent_v0=0.0;
918
919 if(fSampleType=="pPb" || fSampleType=="PbPb")
920 {
921 AliCentrality *centrality=0;
922 if(aod)
923 centrality = aod->GetHeader()->GetCentralityP();
924 // if (centrality->GetQuality() != 0) return ;
925
926 if(centrality)
927 {
928 cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
929 }
930 else
931 {
932 cent_v0= -1;
933 }
934 }
935
936//check the PIDResponse handler
937 if (!fPID) return;
938
939// get mag. field required for twotrack efficiency cut
940 Float_t bSign = 0;
941 bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
942
943 //check for TClonesArray(truth track MC information)
944fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
945 if (!fArrayMC) {
946 AliFatal("Error: MC particles branch not found!\n");
947 return;
948 }
949
950 //check for AliAODMCHeader(truth event MC information)
951 AliAODMCHeader *header=NULL;
952 header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
953 if(!header) {
954 printf("MC header branch not found!\n");
955 return;
956 }
957
958//Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
959Float_t zVtxmc =header->GetVtxZ();
960 if(TMath::Abs(zVtxmc)>fzvtxcut) return;
961
962 // For productions with injected signals, figure out above which label to skip particles/tracks
963 Int_t skipParticlesAbove = 0;
964
965 if (fInjectedSignals)
966 {
967 AliGenEventHeader* eventHeader = 0;
968 Int_t headers = 0;
969
970// AOD
971 if (!header)
972 AliFatal("fInjectedSignals set but no MC header found");
973
974 headers = header->GetNCocktailHeaders();
975 eventHeader = header->GetCocktailHeader(0);
976
977 if (!eventHeader)
978 {
979 // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
980 // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
981 AliError("First event header not found. Skipping this event.");
982 //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
983 return;
984 }
985skipParticlesAbove = eventHeader->NProduced();
986 AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
987 }
988
989 //vertex selection(is it fine for PP?)
990 if ( fSampleType=="pPb"){
991 trkVtx = aod->GetPrimaryVertex();
992 if (!trkVtx || trkVtx->GetNContributors()<=0) return;
993 TString vtxTtl = trkVtx->GetTitle();
994 if (!vtxTtl.Contains("VertexerTracks")) return;
995 zvtx = trkVtx->GetZ();
996 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
997 if (!spdVtx || spdVtx->GetNContributors()<=0) return;
998 TString vtxTyp = spdVtx->GetTitle();
999 Double_t cov[6]={0};
1000 spdVtx->GetCovarianceMatrix(cov);
1001 Double_t zRes = TMath::Sqrt(cov[5]);
1002 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1003 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1004 }
1005 else if(fSampleType=="PbPb" || fSampleType=="pp") {//for pp and pb-pb case
1006 Int_t nVertex = aod->GetNumberOfVertices();
1007 if( nVertex > 0 ) {
1008 trkVtx = aod->GetPrimaryVertex();
1009 Int_t nTracksPrim = trkVtx->GetNContributors();
1010 zvtx = trkVtx->GetZ();
1011 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1012 // Reject TPC only vertex
1013 TString name(trkVtx->GetName());
1014 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
1015
1016 // Select a quality vertex by number of tracks?
1017 if( nTracksPrim < fnTracksVertex ) {
1018 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1019 return ;
1020 }
1021 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1022 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1023 // return kFALSE;
1024 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1025 }
1026 else return;
1027
1028 }
1029 else return;//as there is no proper sample type
1030
1031
1032 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
1033
1034 //count events having a proper vertex
1035 fEventCounter->Fill(2);
1036
1037 if (TMath::Abs(zvtx) > fzvtxcut) return;
1038
1039fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut for trkVtx only
1040
1041//now we have events passed physics trigger, centrality,eventzvtx cut
1042
1043//count events after vertex cut
1044 fEventCounter->Fill(5);
1045
1046if(!aod) return; //for safety
1047
1048if (fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
1049
1050 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)
1051
1052 Double_t cent_v0_truth=0.0;
1053
1054 //TObjArray* tracksMCtruth=0;
1055TObjArray* 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
1056 tracksMCtruth->SetOwner(kTRUE); //***********************************IMPORTANT!
1057
1058 eventno++;
1059
1060 //There is a small difference between zvtx and zVtxmc??????
1061 //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1062 //cout<<"***********************************************zvtx="<<zvtx<<endl;
1063
1064//now process the truth particles(for both efficiency & correlation function)
1065Int_t nMCTrack = fArrayMC->GetEntriesFast();
1066
1067for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
1068{ //MC truth track loop starts
1069
1070AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1071
1072 if(!partMC){
1073 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1074 continue;
1075 }
1076
1077//consider only charged particles
1078 if(partMC->Charge() == 0) continue;
1079
1080//consider only primary particles; neglect all secondary particles including from weak decays
1081 if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1082
1083
1084//remove injected signals(primaries above <maxLabel>)
1085 if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1086
1087//remove duplicates
1088 Bool_t isduplicate=kFALSE;
1089 if (fRemoveDuplicates)
1090 {
1091 for (Int_t j=iMC+1; j<nMCTrack; ++j)
1092 {//2nd trutuh loop starts
1093AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1094 if(!partMC2){
1095 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1096 continue;
1097 }
1098 if (partMC->GetLabel() == partMC2->GetLabel())
1099 {
1100isduplicate=kTRUE;
1101 break;
1102 }
1103 }//2nd truth loop ends
1104 }
1105 if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1106
1107//give only kinematic cuts at the truth level
1108 if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1109 if (partMC->Pt() < fminPt || partMC->Pt() > fmaxPt) continue;
1110
1111 if(!partMC) continue;//for safety
1112
1113 //To determine multiplicity in case of PP
1114 nooftrackstruth++;
1115 //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1116//only physical primary(all/unidentified)
1117if(ffillhistQATruth)
1118 {
1119 MCtruthpt->Fill(partMC->Pt());
1120 MCtrutheta->Fill(partMC->Eta());
1121 MCtruthphi->Fill(partMC->Phi());
1122 }
1123 //get particle ID
1124Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1125Int_t particletypeTruth=-999;
1126 if (TMath::Abs(pdgtruth)==211)
1127 {
1128 particletypeTruth=SpPion;
1129if(ffillhistQATruth)
1130 {
1131 MCtruthpionpt->Fill(partMC->Pt());
1132 MCtruthpioneta->Fill(partMC->Eta());
1133 MCtruthpionphi->Fill(partMC->Phi());
1134 }
1135 }
1136 if (TMath::Abs(pdgtruth)==321)
1137 {
1138 particletypeTruth=SpKaon;
1139if(ffillhistQATruth)
1140 {
1141 MCtruthkaonpt->Fill(partMC->Pt());
1142 MCtruthkaoneta->Fill(partMC->Eta());
1143 MCtruthkaonphi->Fill(partMC->Phi());
1144 }
1145 }
1146if(TMath::Abs(pdgtruth)==2212)
1147 {
1148 particletypeTruth=SpProton;
1149if(ffillhistQATruth)
1150 {
1151 MCtruthprotonpt->Fill(partMC->Pt());
1152 MCtruthprotoneta->Fill(partMC->Eta());
1153 MCtruthprotonphi->Fill(partMC->Phi());
1154 }
1155 }
1156 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)
1157
1158 // -- Fill THnSparse for efficiency and contamination calculation
1159 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
1160
1161 Double_t primmctruth[4] = {cent_v0, zVtxmc,partMC->Pt(), partMC->Eta()};
1162 if(ffillefficiency)
1163 {
1164 fTHngenprimPidTruth[5]->Fill(primmctruth);//for all primary truth particles(4)
1165if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[3]->Fill(primmctruth);//for primary truth mesons(3)
1166if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[4]->Fill(primmctruth);//for primary truth kaons+protons(4)
1167 if (TMath::Abs(pdgtruth)==211) fTHngenprimPidTruth[0]->Fill(primmctruth);//for pions
1168 if (TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[1]->Fill(primmctruth);//for kaons
1169 if(TMath::Abs(pdgtruth)==2212) fTHngenprimPidTruth[2]->Fill(primmctruth);//for protons
1170 }
1171
1172 Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1173if(partMC->Pt()>=fminPtAsso || partMC->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
1174 {
1175LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,partMC->Charge(),partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);
1176//copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1177 copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1178 tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1179 }
1180 }//MC truth track loop ends
1181
1182//*********************still in event loop
1183
1184 if (fSampleType=="pp") cent_v0_truth=nooftrackstruth;
1185 else cent_v0_truth=cent_v0;//the notation cent_v0 is reserved for reco/data case
1186
1187 //now cent_v0_truth should be used for all correlation function calculation
1188
1189if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1190 {
1191 //Fill Correlations for MC truth particles(same event)
1192if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1193 Fillcorrelation(tracksMCtruth,0,cent_v0_truth,zVtxmc,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1194
1195//start mixing
1196AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0_truth, zVtxmc+200);
1197if (pool2 && pool2->IsReady())
1198 {//start mixing only when pool->IsReady
1199if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1200 {//proceed only when no. of trigger particles >0 in current event
1201for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1202 { //pool event loop start
1203 TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1204 if(!bgTracks6) continue;
1205 Fillcorrelation(tracksMCtruth,bgTracks6,cent_v0_truth,zVtxmc,bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1206
1207 }// pool event loop ends mixing case
1208 }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1209} //if pool->IsReady() condition ends mixing case
1210
1211 //still in main event loop
1212
1213 if(tracksMCtruth){
1214if(pool2) pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1215 }
1216 }
1217
1218 //still in main event loop
1219
1220if(tracksMCtruth) delete tracksMCtruth;
1221
1222//now deal with reco tracks
1223 //TObjArray* tracksUNID=0;
1224 TObjArray* tracksUNID = new TObjArray;
1225 tracksUNID->SetOwner(kTRUE);
1226
1227 //TObjArray* tracksID=0;
1228 TObjArray* tracksID = new TObjArray;
1229 tracksID->SetOwner(kTRUE);
1230
1231
1232 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1233
1234 Double_t trackscount=0.0;
1235
1236// loop over reconstructed tracks
1237 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1238{ // reconstructed track loop starts
1239 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1240 if (!track) continue;
1241 //get the corresponding MC track at the truth level (doing reco matching)
1242 AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
1243 if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
1244
1245//remove injected signals
1246 if(fInjectedSignals)
1247 {
1248 AliAODMCParticle* mother = recomatched;
1249
1250 while (!mother->IsPhysicalPrimary())
1251 {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
1252 if (mother->GetMother() < 0)
1253 {
1254 mother = 0;
1255 break;
1256 }
1257
1258 mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1259 if (!mother)
1260 break;
1261 }
1262 if (!mother)
1263 {
1264 Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
1265 continue;
1266 }
1267 if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
1268 }//remove injected signals
1269
1270 if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
1271
1272 Bool_t isduplicate2=kFALSE;
1273if (fRemoveDuplicates)
1274 {
1275 for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
1276 {//2nd loop starts
1277 AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
1278 if (!track2) continue;
1279 AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
1280if(!recomatched2) continue;
1281
1282if (track->GetLabel() == track2->GetLabel())
1283 {
1284isduplicate2=kTRUE;
1285 break;
1286 }
1287 }//2nd loop ends
1288 }
1289 if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
1290
1291 fHistQA[11]->Fill(track->GetTPCNcls());
1292 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
1293
1294 if(tracktype==0) continue;
1295 if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
1296{
1297 if(!track) continue;//for safety
1298 //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
1299 trackscount++;
1300
1301//check for eta , phi holes
1302 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1303 fphiSpectraasso->Fill(track->Phi(),track->Pt());
1304
1305 Int_t particletypeMC=-9999;
1306
1307//tag all particles as unidentified
1308 particletypeMC=unidentified;
1309
1310 Float_t effmatrix=1.;
1311
1312// -- Fill THnSparse for efficiency calculation
1313 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
1314 //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)
1315
1316 //Clone & Reduce track list(TObjArray) for unidentified particles
1317if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
1318 {
1319 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
1320 effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);
1321 LRCParticlePID* copy = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1322 copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
1323 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
1324 }
1325 Double_t allrecomatchedpid[4] = {cent_v0, zVtxmc,recomatched->Pt(), recomatched->Eta()};
1326if(ffillefficiency) fTHnrecomatchedallPid[5]->Fill(allrecomatchedpid);//for all
1327
1328 //now start the particle identification process:)
1329
1330//get the pdg code of the corresponding truth particle
1331 Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
1332
1333Float_t dEdx = track->GetTPCsignal();
1334 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1335
1336 if(HasTOFPID(track))
1337{
1338Float_t beta = GetBeta(track);
1339fHistoTOFbeta->Fill(track->Pt(), beta);
1340 }
1341
1342 //do track identification(nsigma method)
1343 particletypeMC=GetParticle(track);//******************************problem is here
1344
1345//2-d TPCTOF map(for each Pt interval)
1346 if(HasTOFPID(track)){
1347 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
1348 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
1349 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
1350 }
1351 if(particletypeMC==SpUndefined) continue;
1352
1353 //Pt, Eta , Phi distribution of the reconstructed identified particles
1354if(ffillhistQAReco)
1355 {
1356if (particletypeMC==SpPion)
1357 {
1358 fPionPt->Fill(track->Pt());
1359 fPionEta->Fill(track->Eta());
1360 fPionPhi->Fill(track->Phi());
1361 }
1362if (particletypeMC==SpKaon)
1363 {
1364 fKaonPt->Fill(track->Pt());
1365 fKaonEta->Fill(track->Eta());
1366 fKaonPhi->Fill(track->Phi());
1367 }
1368if (particletypeMC==SpProton)
1369 {
1370 fProtonPt->Fill(track->Pt());
1371 fProtonEta->Fill(track->Eta());
1372 fProtonPhi->Fill(track->Phi());
1373 }
1374 }
1375
1376 //for misidentification fraction calculation(do it with tuneonPID)
1377 if(particletypeMC==SpPion )
1378 {
1379 if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
1380 if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
1381 if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
1382 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
1383 }
1384if(particletypeMC==SpKaon )
1385 {
1386 if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
1387 if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
1388 if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
1389 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
1390 }
1391 if(particletypeMC==SpProton )
1392 {
1393 if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
1394 if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
1395 if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
1396 if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
1397 }
1398
1399 //fill tracking efficiency
1400 if(ffillefficiency)
1401 {
1402 if(particletypeMC==SpPion || particletypeMC==SpKaon)
1403 {
1404if(TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[3]->Fill(allrecomatchedpid);//for mesons
1405 }
1406 if(particletypeMC==SpKaon || particletypeMC==SpProton)
1407 {
1408if(TMath::Abs(pdgCode)==321 || TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[4]->Fill(allrecomatchedpid);//for kaons+protons
1409 }
1410 if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211) fTHnrecomatchedallPid[0]->Fill(allrecomatchedpid);//for pions
1411 if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) fTHnrecomatchedallPid[1]->Fill(allrecomatchedpid);//for kaons
1412 if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[2]->Fill(allrecomatchedpid);//for protons
1413 }
1414
1415if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
1416 {
1417if (fapplyTrigefficiency || fapplyAssoefficiency)
1418 effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
1419 LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1420 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1421 tracksID->Add(copy1);
1422 }
1423 }// if(tracktype==1) condition structure ands
1424
1425}//reco track loop ends
1426
1427 //*************************************************************still in event loop
1428
1429//same event delta-eta-deltaphi plot
1430if(fSampleType=="pp") cent_v0=trackscount;//multiplicity
1431
1432if(trackscount>0.0)
1433 {
1434//fill the centrality/multiplicity distribution of the selected events
1435 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
1436
1437 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??????
1438
1439//count selected events having centrality betn 0-100%
1440 fEventCounter->Fill(6);
1441
1442 //same event delte-eta, delta-phi plot
1443if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
1444 {//same event calculation starts
1445 if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,bSign,kTRUE,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
1446 if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
1447 }
1448
1449if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
1450 {//same event calculation starts
1451 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
1452 if(ffilltrigIDassoID) Fillcorrelation(tracksID,0,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
1453 }
1454
1455//still in main event loop
1456//start mixing
1457 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
1458AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
1459if (pool && pool->IsReady())
1460 {//start mixing only when pool->IsReady
1461 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
1462 { //pool event loop start
1463 TObjArray* bgTracks = pool->GetEvent(jMix);
1464 if(!bgTracks) continue;
1465 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
1466 Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,bSign,kTRUE,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
1467 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
1468 Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
1469 }// pool event loop ends mixing case
1470
1471} //if pool->IsReady() condition ends mixing case
1472 if(tracksUNID) {
1473if(pool)
1474 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
1475 }
1476 }//mixing with unidentified particles
1477
1478 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
1479AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
1480if (pool1 && pool1->IsReady())
1481 {//start mixing only when pool->IsReady
1482for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
1483 { //pool event loop start
1484 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
1485 if(!bgTracks2) continue;
1486if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
1487 Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
1488if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
1489 Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
1490
1491 }// pool event loop ends mixing case
1492} //if pool1->IsReady() condition ends mixing case
1493
1494if(tracksID) {
1495if(pool1)
1496 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
1497 }
1498 }//mixing with identified particles
1499
1500 //no. of events analyzed
1501fEventCounter->Fill(7);
1502 }
1503
1504if(tracksUNID) delete tracksUNID;
1505
1506if(tracksID) delete tracksID;
1507
1508
1509PostData(1, fOutput);
1510
1511}
1512//________________________________________________________________________
1513void AliTwoParticlePIDCorr::doAODevent()
1514{
1515
1516 //get AOD
1517 AliVEvent *event = InputEvent();
1518 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1519 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1520 if (!aod) {
1521 AliError("Cannot get the AOD event");
1522 return;
1523 }
1524
1525// count all events
1526 fEventCounter->Fill(1);
1527
1528 // get centrality object and check quality
1529 Double_t cent_v0=0;
1530
1531 if(fSampleType=="pPb" || fSampleType=="PbPb")
1532 {
1533 AliCentrality *centrality=0;
1534 if(aod)
1535 centrality = aod->GetHeader()->GetCentralityP();
1536 // if (centrality->GetQuality() != 0) return ;
1537
1538 if(centrality)
1539 {
1540 cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
1541 }
1542 else
1543 {
1544 cent_v0= -1;
1545 }
1546 }
1547
1548 Float_t bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
1549 Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
1550
1551// Pileup selection ************************************************
1552// if(frejectPileUp) //applicable only for TPC only tracks,not for hybrid tracks?.
1553 // {
1554 //if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(aod)) return;
1555 // }
1556
1557 if (!fPID) return;//this should be available with each event even if we don't do PID selection
1558
1559 //vertex selection(is it fine for PP?)
1560 if ( fSampleType=="pPb"){
1561 trkVtx = aod->GetPrimaryVertex();
1562 if (!trkVtx || trkVtx->GetNContributors()<=0) return;
1563 TString vtxTtl = trkVtx->GetTitle();
1564 if (!vtxTtl.Contains("VertexerTracks")) return;
1565 zvtx = trkVtx->GetZ();
1566 const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
1567 if (!spdVtx || spdVtx->GetNContributors()<=0) return;
1568 TString vtxTyp = spdVtx->GetTitle();
1569 Double_t cov[6]={0};
1570 spdVtx->GetCovarianceMatrix(cov);
1571 Double_t zRes = TMath::Sqrt(cov[5]);
1572 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1573 if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1574 }
1575 else if(fSampleType=="PbPb" || fSampleType=="pp") {//for pp and pb-pb case
1576 Int_t nVertex = aod->GetNumberOfVertices();
1577 if( nVertex > 0 ) {
1578 trkVtx = aod->GetPrimaryVertex();
1579 Int_t nTracksPrim = trkVtx->GetNContributors();
1580 zvtx = trkVtx->GetZ();
1581 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1582 // Reject TPC only vertex
1583 TString name(trkVtx->GetName());
1584 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
1585
1586 // Select a quality vertex by number of tracks?
1587 if( nTracksPrim < fnTracksVertex ) {
1588 //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1589 return ;
1590 }
1591 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1592 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1593 // return kFALSE;
1594 // if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1595 }
1596 else return;
1597
1598 }
1599 else return;//as there is no proper sample type
1600
1601
1602 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
1603
1604//count events having a proper vertex
1605 fEventCounter->Fill(2);
1606
1607 if (TMath::Abs(zvtx) > fzvtxcut) return;
1608
1609//count events after vertex cut
1610 fEventCounter->Fill(5);
1611
1612
1613 //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
1614
1615 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
1616
1617
1618 if(!aod) return; //for safety
1619
1620if(fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0) return;//for pp case it is the multiplicity
1621
1622 TObjArray* tracksUNID= new TObjArray;//track info before doing PID
1623 tracksUNID->SetOwner(kTRUE); // IMPORTANT!
1624
1625 TObjArray* tracksID= new TObjArray;//only pions, kaons,protonsi.e. after doing the PID selection
1626 tracksID->SetOwner(kTRUE); // IMPORTANT!
1627
1628 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)
1629
1630 eventno++;
1631
1632 for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1633{ //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
1634 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1635 if (!track) continue;
1636 fHistQA[11]->Fill(track->GetTPCNcls());
1637 Int_t particletype=-9999;//required for PID filling
1638 Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
1639 if(tracktype!=1) continue;
1640
1641 if(!track) continue;//for safety
1642
1643//check for eta , phi holes
1644 fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1645 fphiSpectraasso->Fill(track->Phi(),track->Pt());
1646
1647 trackscount++;
1648
1649 //if no applyefficiency , set the eff factor=1.0
1650 Float_t effmatrix=1.0;
1651
1652 //tag all particles as unidentified that passed filterbit & kinematic cuts
1653 particletype=unidentified;
1654
1655
1656 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
1657
1658
1659 //to reduce memory consumption in pool
1660 if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)
1661 {
1662 //Clone & Reduce track list(TObjArray) for unidentified particles
1663 if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
1664 effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);
1665 LRCParticlePID* copy = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1666 copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1667 tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
1668 }
1669
1670//now start the particle identificaion process:)
1671
1672//track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
1673
1674 Float_t dEdx = track->GetTPCsignal();
1675 fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1676
1677 //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)
1678 if(HasTOFPID(track))
1679{
1680 Float_t beta = GetBeta(track);
1681 fHistoTOFbeta->Fill(track->Pt(), beta);
1682 }
1683
1684
1685//track identification(using nsigma method)
1686 particletype=GetParticle(track);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
1687
1688
1689//2-d TPCTOF map(for each Pt interval)
1690 if(HasTOFPID(track)){
1691 fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
1692 fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
1693 fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
1694 }
1695
1696//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!!!!!
1697 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)!!!!!!!!!!!
1698
1699 //Pt, Eta , Phi distribution of the reconstructed identified particles
1700if(ffillhistQAReco)
1701 {
1702if (particletype==SpPion)
1703 {
1704 fPionPt->Fill(track->Pt());
1705 fPionEta->Fill(track->Eta());
1706 fPionPhi->Fill(track->Phi());
1707 }
1708if (particletype==SpKaon)
1709 {
1710 fKaonPt->Fill(track->Pt());
1711 fKaonEta->Fill(track->Eta());
1712 fKaonPhi->Fill(track->Phi());
1713 }
1714if (particletype==SpProton)
1715 {
1716 fProtonPt->Fill(track->Pt());
1717 fProtonEta->Fill(track->Eta());
1718 fProtonPhi->Fill(track->Phi());
1719 }
1720 }
1721
1722if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)
1723 {
1724if (fapplyTrigefficiency || fapplyAssoefficiency)
1725 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
1726 LRCParticlePID* copy1 = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1727 copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1728 tracksID->Add(copy1);
1729 }
1730} //track loop ends but still in event loop
1731
1732if(trackscount<1.0){
1733 if(tracksUNID) delete tracksUNID;
1734 if(tracksID) delete tracksID;
1735 return;
1736 }
1737
1738// cout<<fminPtAsso<<"***"<<fmaxPtAsso<<"*********************"<<fminPtTrig<<"***"<<fmaxPtTrig<<"*****************"<<kTrackVariablesPair<<endl;
1739if(fSampleType=="pp") cent_v0=trackscount;//multiplicity
1740
1741//fill the centrality/multiplicity distribution of the selected events
1742 fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
1743
1744if(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??????
1745
1746//count selected events having centrality betn 0-100%
1747 fEventCounter->Fill(6);
1748
1749//same event delta-eta-deltaphi plot
1750
1751if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
1752 {//same event calculation starts
1753 if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,bSign,kTRUE,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
1754 if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID) Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
1755 }
1756
1757if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
1758 {//same event calculation starts
1759 if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID) Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
1760 if(ffilltrigIDassoID) Fillcorrelation(tracksID,0,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
1761 }
1762
1763//still in main event loop
1764
1765
1766//start mixing
1767 if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
1768AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
1769if (pool && pool->IsReady())
1770 {//start mixing only when pool->IsReady
1771 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
1772 { //pool event loop start
1773 TObjArray* bgTracks = pool->GetEvent(jMix);
1774 if(!bgTracks) continue;
1775 if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
1776 Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,bSign,kTRUE,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
1777 if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
1778 Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
1779 }// pool event loop ends mixing case
1780
1781} //if pool->IsReady() condition ends mixing case
1782 if(tracksUNID) {
1783if(pool)
1784 pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
1785 }
1786 }//mixing with unidentified particles
1787
1788
1789 if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
1790AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
1791if (pool1 && pool1->IsReady())
1792 {//start mixing only when pool->IsReady
1793for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
1794 { //pool event loop start
1795 TObjArray* bgTracks2 = pool1->GetEvent(jMix);
1796 if(!bgTracks2) continue;
1797if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
1798 Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
1799if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
1800 Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
1801
1802 }// pool event loop ends mixing case
1803} //if pool1->IsReady() condition ends mixing case
1804
1805if(tracksID) {
1806if(pool1)
1807 pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
1808 }
1809 }//mixing with identified particles
1810
1811
1812 //no. of events analyzed
1813fEventCounter->Fill(7);
1814
1815//still in main event loop
1816
1817
1818if(tracksUNID) delete tracksUNID;
1819
1820if(tracksID) delete tracksID;
1821
1822
1823PostData(1, fOutput);
1824
1825} // *************************event loop ends******************************************//_______________________________________________________________________
1826TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
1827{
1828 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
1829
1830 TObjArray* tracksClone = new TObjArray;
1831 tracksClone->SetOwner(kTRUE);
1832
1833 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1834 {
1835 LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
1836 LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval());
1837 copy100->SetUniqueID(particle->GetUniqueID());
1838 tracksClone->Add(copy100);
1839 }
1840
1841 return tracksClone;
1842}
1843
1844//--------------------------------------------------------------------------------
1845void AliTwoParticlePIDCorr::Fillcorrelation(TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)
1846{
1847
1848 //before calling this function check that either trackstrig & tracksasso are available
1849
1850 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
1851 TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
1852 TArrayF eta(input->GetEntriesFast());
1853 for (Int_t i=0; i<input->GetEntriesFast(); i++)
1854 eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
1855
1856 //if(trackstrig)
1857 Int_t jmax=trackstrig->GetEntriesFast();
1858 if(tracksasso)
1859 jmax=tracksasso->GetEntriesFast();
1860
1861// identify K, Lambda candidates and flag those particles
1862 // a TObject bit is used for this
1863const UInt_t kResonanceDaughterFlag = 1 << 14;
1864 if (fRejectResonanceDaughters > 0)
1865 {
1866 Double_t resonanceMass = -1;
1867 Double_t massDaughter1 = -1;
1868 Double_t massDaughter2 = -1;
1869 const Double_t interval = 0.02;
1870 switch (fRejectResonanceDaughters)
1871 {
1872 case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
1873 case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
1874 case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
1875 default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
1876 }
1877
1878for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
1879 trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
1880 for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
1881 tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
1882
1883 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
1884 {
1885 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
1886for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
1887 {
1888 LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
1889
1890 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
1891if (trig->IsEqual(asso)) continue;
1892
1893if (trig->Charge() * asso->Charge() > 0) continue;
1894
1895 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
1896
1897if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
1898 {
1899 mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
1900
1901 if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
1902 {
1903 trig->SetBit(kResonanceDaughterFlag);
1904 asso->SetBit(kResonanceDaughterFlag);
1905
1906// Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
1907 }
1908 }
1909 }
1910 }
1911 }
1912
1913 //two particle correlation filling
1914
1915for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
1916 { //trigger loop starts
1917 LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
1918 if(!trig) continue;
1919 Float_t trigpt=trig->Pt();
1920 //to avoid overflow qnd underflow
1921 if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
1922 Int_t particlepidtrig=trig->getparticle();
1923 if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
1924
1925 Float_t trigeta=trig->Eta();
1926
1927 // some optimization
1928 if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
1929 continue;
1930
1931if (fOnlyOneEtaSide != 0)
1932 {
1933 if (fOnlyOneEtaSide * trigeta < 0)
1934 continue;
1935 }
1936 if (fTriggerSelectCharge != 0)
1937 if (trig->Charge() * fTriggerSelectCharge < 0)
1938 continue;
1939
1940 if (fRejectResonanceDaughters > 0)
1941 if (trig->TestBit(kResonanceDaughterFlag)) continue;
1942
1943 Float_t trigphi=trig->Phi();
1944 Float_t trackefftrig=1.0;
1945 if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
1946 // cout<<"*******************trackefftrig="<<trackefftrig<<endl;
1947 Double_t* trigval;
1948 Int_t dim=3;
1949 if(fcontainPIDtrig) dim=4;
1950 trigval= new Double_t[dim];
1951 trigval[0] = cent;
1952 trigval[1] = vtx;
1953 trigval[2] = trigpt;
1954if(fcontainPIDtrig) trigval[3] = particlepidtrig;
1955
1956 //filling only for same event case(mixcase=kFALSE)
1957
1958 //trigger particle counting only when mixcase=kFALSE i.e for same event correlation function calculation
1959if(mixcase==kFALSE)
1960 {
1961 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
1962 if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1963 }
1964 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
1965 if(fillup=="trigassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1966 }
1967if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
1968 if(fillup=="trigUNIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1969 }
1970 //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
1971if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
1972 if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1973 }
1974 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
1975 if(fillup=="trigIDassoUNID" ) fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1976 }
1977if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
1978 if(fillup=="trigIDassoID") fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1979 }
1980
1981 if(fillup=="trigIDassoIDMCTRUTH") fTHnTrigcountMCTruthPrim->Fill(trigval,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!!!!
1982 }
1983
1984 //asso loop starts within trigger loop
1985 for(Int_t j=0;j<jmax;j++)
1986 {
1987 LRCParticlePID *asso=0;
1988 if(!tracksasso)
1989 asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
1990 else
1991 asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
1992
1993 if(!asso) continue;
1994
1995 //to avoid overflow qnd underflow
1996 if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
1997
1998 if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
1999
2000 if(!tracksasso && i==j) continue;
2001
2002 // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
2003 // if (tracksasso && trig->IsEqual(asso)) continue;
2004
2005 if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
2006
2007 if (fPtOrder)
2008 if (asso->Pt() >= trig->Pt()) continue;
2009
2010 Int_t particlepidasso=asso->getparticle();
2011 if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
2012
2013
2014if (fAssociatedSelectCharge != 0)
2015if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
2016
2017 if (fSelectCharge > 0)
2018 {
2019 // skip like sign
2020 if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
2021 continue;
2022
2023 // skip unlike sign
2024 if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
2025 continue;
2026 }
2027
2028if (fEtaOrdering)
2029 {
2030 if (trigeta < 0 && asso->Eta() < trigeta)
2031 continue;
2032 if (trigeta > 0 && asso->Eta() > trigeta)
2033 continue;
2034 }
2035
2036if (fRejectResonanceDaughters > 0)
2037 if (asso->TestBit(kResonanceDaughterFlag))
2038 {
2039// Printf("Skipped j=%d", j);
2040 continue;
2041 }
2042
2043 // conversions
2044 if (fCutConversions && asso->Charge() * trig->Charge() < 0)
2045 {
2046 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2047
2048 if (mass < 0.1)
2049 {
2050 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2051
2052 //fControlConvResoncances->Fill(0.0, mass);
2053
2054 if (mass < 0.04*0.04)
2055 continue;
2056 }
2057 }
2058
2059 // K0s
2060 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2061 {
2062 Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
2063
2064 const Float_t kK0smass = 0.4976;
2065
2066 if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
2067 {
2068 mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
2069
2070 //fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
2071
2072 if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
2073 continue;
2074 }
2075 }
2076
2077 // Lambda
2078 if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2079 {
2080 Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
2081 Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2082
2083 const Float_t kLambdaMass = 1.115;
2084
2085 if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
2086 {
2087 mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
2088
2089 //fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
2090
2091 if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2092 continue;
2093 }
2094 if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
2095 {
2096 mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2097
2098 //fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
2099
2100 if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2101 continue;
2102 }
2103 }
2104
2105 if (twoTrackEfficiencyCut)
2106 {
2107 // the variables & cuthave been developed by the HBT group
2108 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
2109 Float_t phi1 = trig->Phi();
2110 Float_t pt1 = trig->Pt();
2111 Float_t charge1 = trig->Charge();
2112 Float_t phi2 = asso->Phi();
2113 Float_t pt2 = asso->Pt();
2114 Float_t charge2 = asso->Charge();
2115
2116 Float_t deta= trigeta - eta[j];
2117
2118 // optimization
2119 if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
2120 {
2121
2122 // check first boundaries to see if is worth to loop and find the minimum
2123 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
2124 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
2125
2126 const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
2127
2128 Float_t dphistarminabs = 1e5;
2129 Float_t dphistarmin = 1e5;
2130
2131 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
2132 {
2133 for (Double_t rad=0.8; rad<2.51; rad+=0.01)
2134 {
2135 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
2136
2137 Float_t dphistarabs = TMath::Abs(dphistar);
2138
2139 if (dphistarabs < dphistarminabs)
2140 {
2141 dphistarmin = dphistar;
2142 dphistarminabs = dphistarabs;
2143 }
2144 }
2145
2146if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
2147 {
2148// 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);
2149 continue;
2150 }
2151//fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
2152
2153 }
2154 }
2155 }
2156
2157 Float_t trackeffasso=1.0;
2158 if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
2159 //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
2160 Float_t deleta=trigeta-eta[j];
2161 Float_t delphi=PhiRange(trigphi-asso->Phi());
2162
2163 //here get the two particle efficiency correction factor
2164 Float_t effweight=trackefftrig*trackeffasso;
2165 //cout<<"*******************effweight="<<effweight<<endl;
2166 Double_t* vars;
2167 vars= new Double_t[kTrackVariablesPair];
2168 vars[0]=cent;
2169 vars[1]=vtx;
2170 vars[2]=trigpt;
2171 vars[3]=asso->Pt();
2172 vars[4]=deleta;
2173 vars[5]=delphi;
2174if(fcontainPIDtrig && !fcontainPIDasso) vars[6]=particlepidtrig;
2175if(!fcontainPIDtrig && fcontainPIDasso) vars[6]=particlepidasso;
2176 if(fcontainPIDtrig && fcontainPIDasso){
2177 vars[6]=particlepidtrig;
2178 vars[7]=particlepidasso;
2179 }
2180
2181 //Fill Correlation ThnSparses
2182 if(fillup=="trigassoUNID")
2183 {
2184 if(mixcase==kFALSE) fTHnCorrUNID->Fill(vars,1.0/effweight);
2185 if(mixcase==kTRUE) fTHnCorrUNIDmix->Fill(vars,1.0/effweight);
2186 }
2187 if(fillup=="trigIDassoID")
2188 {
2189 if(mixcase==kFALSE) fTHnCorrID->Fill(vars,1.0/effweight);
2190 if(mixcase==kTRUE) fTHnCorrIDmix->Fill(vars,1.0/effweight);
2191 }
2192 if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
2193 {//in this case effweight should be 1 always
2194 if(mixcase==kFALSE) fCorrelatonTruthPrimary->Fill(vars,1.0/effweight);
2195 if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,1.0/effweight);
2196 }
2197 if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
2198 {
2199 if(mixcase==kFALSE) fTHnCorrIDUNID->Fill(vars,1.0/effweight);
2200 if(mixcase==kTRUE) fTHnCorrIDUNIDmix->Fill(vars,1.0/effweight);
2201 }
2202
2203delete[] vars;
2204 }//asso loop ends
2205delete[] trigval;
2206 }//trigger loop ends
2207
2208}
2209
2210//--------------------------------------------------------------------------------
2211Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
2212{
2213 //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
2214 Int_t effVars[4];
2215 Float_t effvalue=1.;
2216
2217 if(parpid==unidentified)
2218 {
2219 effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
2220 effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
2221 effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
2222 effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
2223 effvalue=effcorection[5]->GetBinContent(effVars);
2224 }
2225if(parpid==SpPion || parpid==SpKaon)
2226 {
2227 if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2228 {
2229 effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
2230 effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
2231 effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
2232 effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
2233 effvalue=effcorection[3]->GetBinContent(effVars);
2234 }
2235 else{
2236 if(parpid==SpPion)
2237 {
2238 effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
2239 effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
2240 effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
2241 effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
2242 effvalue=effcorection[0]->GetBinContent(effVars);
2243 }
2244
2245 if(parpid==SpKaon)
2246 {
2247 effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
2248 effVars[1] = effcorection[1]->GetAxis(1)->FindBin(evzvtx);
2249 effVars[2] = effcorection[1]->GetAxis(2)->FindBin(track->Pt());
2250 effVars[3] = effcorection[1]->GetAxis(3)->FindBin(track->Eta());
2251 effvalue=effcorection[1]->GetBinContent(effVars);
2252 }
2253 }
2254 }
2255
2256 if(parpid==SpProton)
2257 {
2258 effVars[0] = effcorection[2]->GetAxis(0)->FindBin(cent);
2259 effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
2260 effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
2261 effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
2262 effvalue=effcorection[2]->GetBinContent(effVars);
2263 }
2264
2265 if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2266 {
2267 if(parpid==SpProton || parpid==SpKaon)
2268 {
2269 effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
2270 effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
2271 effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
2272 effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
2273 effvalue=effcorection[4]->GetBinContent(effVars);
2274 }
2275 }
2276 // Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
2277 if(effvalue==0.) effvalue=1.;
2278
2279 return effvalue;
2280
2281}
2282//-----------------------------------------------------------------------
2283
2284Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield)
2285{
2286
2287 if(!track) return 0;
2288 Bool_t trackOK = track->TestFilterBit(fFilterBit);
2289 if(!trackOK) return 0;
2290 //select only primary traks(for data & reco MC tracks)
2291 if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
2292 if(track->Charge()==0) return 0;
2293 fHistQA[12]->Fill(track->GetTPCNcls());
2294 Float_t dxy, dz;
2295 dxy = track->DCA();
2296 dz = track->ZAtDCA();
2297 fHistQA[6]->Fill(dxy);
2298 fHistQA[7]->Fill(dz);
2299 Float_t chi2ndf = track->Chi2perNDF();
2300 fHistQA[13]->Fill(chi2ndf);
2301 Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
2302 fHistQA[14]->Fill(nCrossedRowsTPC);
2303 //Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
2304 if (track->GetTPCNclsF()>0) {
2305 Float_t ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
2306 fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
2307 }
2308//accepted tracks
2309 Float_t pt=track->Pt();
2310 if(pt< fminPt || pt> fmaxPt) return 0;
2311 if(TMath::Abs(track->Eta())> fmaxeta) return 0;
2312 if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
2313 //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required
2314// DCA XY
2315 if (fdcacut && fDCAXYCut)
2316 {
2317 if (!vertex)
2318 return 0;
2319
2320 Double_t pos[2];
2321 Double_t covar[3];
2322 AliAODTrack* clone =(AliAODTrack*) track->Clone();
2323 Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
2324 delete clone;
2325 if (!success)
2326 return 0;
2327
2328// Printf("%f", ((AliAODTrack*)part)->DCA());
2329// Printf("%f", pos[0]);
2330 if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
2331 return 0;
2332 }
2333
2334 fHistQA[8]->Fill(pt);
2335 fHistQA[9]->Fill(track->Eta());
2336 fHistQA[10]->Fill(track->Phi());
2337 return 1;
2338 }
2339 //________________________________________________________________________________
2340void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track)
2341{
2342//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 )
2343Float_t pt=track->Pt();
2344
2345//it is assumed that every track that passed the filterbit have proper TPC response(!!)
2346Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
2347Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
2348Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
2349
2350Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
2351Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
2352
2353 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
2354 {
2355
2356nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
2357nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
2358nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
2359//---combined
2360nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
2361nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
2362nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
2363
2364
2365 }
2366else{
2367 // --- combined
2368 // if TOF is missing and below fPtTOFPID only the TPC information is used
2369 nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
2370 nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
2371 nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
2372
2373 }
2374
2375//set data member fnsigmas
2376 fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
2377 fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
2378 fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
2379
2380 //for all tracks below fPtTOFPIDmin and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
2381 fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
2382 fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
2383 fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
2384
2385 //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)
2386 fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
2387 fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
2388 fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
2389
2390
2391}
2392//----------------------------------------------------------------------------
2393Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track)
2394{
2395 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
2396if(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
2397//get the identity of the particle with the minimum Nsigma
2398 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
2399 switch (fPIDType){
2400 case NSigmaTPC:
2401 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
2402 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
2403 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
2404 break;
2405 case NSigmaTOF:
2406 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
2407 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
2408 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
2409 break;
2410 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
2411 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
2412 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
2413 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
2414 break;
2415 }
2416
2417 if(track->Pt()<=fPtTOFPIDmax) {
2418 // guess the particle based on the smaller nsigma (within fNSigmaPID)
2419 if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
2420if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)) return SpKaon;
2421if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) return SpPion;
2422if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) return SpProton;
2423
2424// else, return undefined
2425 return SpUndefined;
2426 }
2427 else {//asymmetric nsigma cut around pion's mean position for tracks having Pt>4.0 Gev
2428 if(fminprotonsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxprotonsigmacut) return SpProton;
2429 if(fminpionsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxpionsigmacut) return SpPion;
2430// else, return undefined(here we don't detect kaons)
2431 return SpUndefined;
2432 }
2433
2434}
2435
2436//------------------------------------------------------------------------------------------
2437Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk){
2438 //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
2439
2440 //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
2441 //fill DC histos
2442 for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
2443
2444 Int_t MinNSigma=FindMinNSigma(trk);//not filling the NSigmaRec histos
2445
2446
2447 if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
2448
2449 Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
2450 switch (fPIDType) {
2451 case NSigmaTPC:
2452 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
2453 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPC]) ;
2454 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPC]) ;
2455 break;
2456 case NSigmaTOF:
2457 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
2458 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTOF]) ;
2459 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTOF]) ;
2460 break;
2461 case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
2462 nsigmaProton = TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
2463 nsigmaKaon = TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF]) ;
2464 nsigmaPion = TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF]) ;
2465 break;
2466 }
2467
2468 //there is chance of overlapping only for particles having pt below 4.0 GEv
2469 if(trk->Pt()<=4.0){
2470 if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
2471 if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
2472 if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
2473
2474 }
2475
2476
2477 return fHasDoubleCounting;
2478}
2479
2480//////////////////////////////////////////////////////////////////////////////////////////////////
2481Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk){
2482 //return the specie according to the minimum nsigma value
2483 //no double counting, this has to be evaluated using CheckDoubleCounting()
2484 //Printf("fPtTOFPID %.1f, fRequestTOFPID %d, fNSigmaPID %.1f, fPIDType %d",fPtTOFPID,fRequestTOFPID,fNSigmaPID,fPIDType);
2485
2486 CalculateNSigmas(trk);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
2487
2488 if(fUseExclusiveNSigma){
2489 Bool_t *HasDC;
2490 HasDC=GetDoubleCounting(trk);
2491 for(Int_t ipart=0;ipart<NSpecies;ipart++){
2492 if(HasDC[ipart]==kTRUE) return SpUndefined;
2493 }
2494 return FindMinNSigma(trk);//NSigmaRec distr filled here
2495 }
2496 else return FindMinNSigma(trk);//NSigmaRec distr filled here
2497
2498}
2499
2500//-------------------------------------------------------------------------------------
2501Bool_t
2502AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
2503{
2504 // check PID signal
2505 AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
2506 if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
2507 //ULong_t status=track->GetStatus();
2508 //if (!( (status & AliAODTrack::kTPCpid ) == AliAODTrack::kTPCpid )) return kFALSE;//remove light nuclei
2509 //if (track->GetTPCsignal() <= 0.) return kFALSE;
2510 // 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.
2511
2512 return kTRUE;
2513}
2514//___________________________________________________________
2515
2516Bool_t
2517AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
2518{
2519 // check TOF matched track
2520 //ULong_t status=track->GetStatus();
2521 //if (!( (status & AliAODTrack::kITSin ) == AliAODTrack::kITSin )) return kFALSE;
2522 AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
2523 if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
2524 if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
2525 //if(!((status & AliAODTrack::kTOFpid ) == AliAODTrack::kTOFpid )) return kFALSE;
2526 //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
2527 // if (probMis > 0.01) return kFALSE;
2528if(fRemoveTracksT0Fill)
2529 {
2530Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
2531 if (startTimeMask < 0)return kFALSE;
2532 }
2533 return kTRUE;
2534}
2535
2536//________________________________________________________________________
2537Float_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
2538{
2539 //it is called only when TOF PID is available
2540 Double_t p = track->P();
2541 Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
2542 Double_t timei[5];
2543 track->GetIntegratedTimes(timei);
2544 return timei[0]/time;
2545}
2546//------------------------------------------------------------------------------------------------------
2547
2548Float_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)
2549{
2550 // calculate inv mass squared
2551 // same can be achieved, but with more computing time with
2552 /*TLorentzVector photon, p1, p2;
2553 p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
2554 p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
2555 photon = p1+p2;
2556 photon.M()*/
2557
2558 Float_t tantheta1 = 1e10;
2559
2560 if (eta1 < -1e-10 || eta1 > 1e-10)
2561 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
2562
2563 Float_t tantheta2 = 1e10;
2564 if (eta2 < -1e-10 || eta2 > 1e-10)
2565 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
2566
2567 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
2568 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
2569
2570 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 ) ) );
2571
2572 return mass2;
2573}
2574//---------------------------------------------------------------------------------
2575
2576Float_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)
2577{
2578 // calculate inv mass squared approximately
2579
2580 Float_t tantheta1 = 1e10;
2581
2582 if (eta1 < -1e-10 || eta1 > 1e-10)
2583 {
2584 Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
2585 tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
2586 }
2587
2588 Float_t tantheta2 = 1e10;
2589 if (eta2 < -1e-10 || eta2 > 1e-10)
2590 {
2591 Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
2592 tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
2593 }
2594
2595 Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
2596 Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
2597
2598 // fold onto 0...pi
2599 Float_t deltaPhi = TMath::Abs(phi1 - phi2);
2600 while (deltaPhi > TMath::TwoPi())
2601 deltaPhi -= TMath::TwoPi();
2602 if (deltaPhi > TMath::Pi())
2603 deltaPhi = TMath::TwoPi() - deltaPhi;
2604
2605 Float_t cosDeltaPhi = 0;
2606 if (deltaPhi < TMath::Pi()/3)
2607 cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
2608 else if (deltaPhi < 2*TMath::Pi()/3)
2609 cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
2610 else
2611 cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
2612
2613 Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
2614
2615// Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
2616
2617 return mass2;
2618}
2619//--------------------------------------------------------------------------------
2620Float_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)
2621{
2622 //
2623 // calculates dphistar
2624 //
2625
2626 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2627
2628 static const Double_t kPi = TMath::Pi();
2629
2630 // circularity
2631// if (dphistar > 2 * kPi)
2632// dphistar -= 2 * kPi;
2633// if (dphistar < -2 * kPi)
2634// dphistar += 2 * kPi;
2635
2636 if (dphistar > kPi)
2637 dphistar = kPi * 2 - dphistar;
2638 if (dphistar < -kPi)
2639 dphistar = -kPi * 2 - dphistar;
2640 if (dphistar > kPi) // might look funny but is needed
2641 dphistar = kPi * 2 - dphistar;
2642
2643 return dphistar;
2644}
2645//_________________________________________________________________________
2646void AliTwoParticlePIDCorr ::DefineEventPool()
2647{
2648const Int_t MaxNofEvents=1000;
2649const Int_t MaxNofTracks=50000;
2650const Int_t NofVrtxBins=10+(1+10)*2;
2651Double_t ZvrtxBins[NofVrtxBins+1]={ -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10,
2652 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110,
2653 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210
2654 };
2655 if (fSampleType=="pp"){
2656const Int_t NofCentBins=5;
2657Double_t CentralityBins[NofCentBins+1]={0.,10., 20., 40.,80.,200.1};
2658fPoolMgr = new AliEventPoolManager(MaxNofEvents,MaxNofTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
2659 }
2660if (fSampleType=="pPb" || fSampleType=="PbPb")
2661 {
2662const Int_t NofCentBins=15;
2663Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
2664fPoolMgr = new AliEventPoolManager(MaxNofEvents,MaxNofTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
2665 }
2666fPoolMgr->SetTargetValues(MaxNofTracks, 0.1, 5);
2667
2668//if(!fPoolMgr) return kFALSE;
2669//return kTRUE;
2670
2671}
2672//------------------------------------------------------------------------
2673Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
2674{
2675 // This method is a copy from AliUEHist::GetBinning
2676 // takes the binning from <configuration> identified by <tag>
2677 // configuration syntax example:
2678 // 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
2679 // phi: .....
2680 //
2681 // returns bin edges which have to be deleted by the caller
2682
2683 TString config(configuration);
2684 TObjArray* lines = config.Tokenize("\n");
2685 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2686 {
2687 TString line(lines->At(i)->GetName());
2688 if (line.BeginsWith(TString(tag) + ":"))
2689 {
2690 line.Remove(0, strlen(tag) + 1);
2691 line.ReplaceAll(" ", "");
2692 TObjArray* binning = line.Tokenize(",");
2693 Double_t* bins = new Double_t[binning->GetEntriesFast()];
2694 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
2695 bins[j] = TString(binning->At(j)->GetName()).Atof();
2696
2697 nBins = binning->GetEntriesFast() - 1;
2698
2699 delete binning;
2700 delete lines;
2701 return bins;
2702 }
2703 }
2704
2705 delete lines;
2706 AliFatal(Form("Tag %s not found in %s", tag, configuration));
2707 return 0;
2708}
2709//_______________________________________________________________________________
2710void AliTwoParticlePIDCorr::SetAsymmetricBin(THnSparse *h,Int_t dim,Double_t *arraybin,Int_t arraybinsize,TString axisTitle)
2711{
2712 TAxis *axis = 0x0;
2713 axis = h->GetAxis(dim);
2714 axis->Set(arraybinsize,arraybin);
2715 axis->SetTitle(axisTitle);
2716}
2717
2718//________________________________________________________________________
2719void AliTwoParticlePIDCorr::Terminate(Option_t *)
2720{
2721 // Draw result to screen, or perform fitting, normalizations
2722 // Called once at the end of the query
2723 fOutput = dynamic_cast<TList*> (GetOutputData(1));
2724 if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
2725
2726
2727}
2728//------------------------------------------------------------------
2729