AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskMultiparticleFemtoscopy.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2008, 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  /********************************
17  * femtoscopy with multiparticle *
18  * technology *
19  * *
20  * author: Ante Bilandzic *
21  * (abilandzic@gmail.com) *
22  ********************************/
23 
24 #include "Riostream.h"
25 #include "TDatabasePDG.h"
27 #include "AliLog.h"
28 #include "AliESDEvent.h"
29 #include "AliESDInputHandler.h"
30 #include "AliAODEvent.h"
31 #include "AliAODInputHandler.h"
32 #include "AliMCEvent.h"
33 #include "AliMCEventHandler.h"
34 #include "AliAnalysisManager.h"
35 #include "AliCentrality.h"
36 #include "AliStack.h"
37 #include "TFile.h"
38 
39 using std::cout;
40 using std::endl;
41 
43 
44 //================================================================================================================
45 
47  AliAnalysisTaskSE(name),
48  fHistList(NULL),
49  fAnalysisType(NULL),
50  fPIDResponse(NULL),
51  fMaxNoGlobalTracksAOD(5), // TBI this is landmine
52  fProcessBothKineAndReco(kFALSE),
53  fProcessOnlyKine(kFALSE),
54  fProcessOnlyReco(kFALSE),
55  fRejectFakeTracks(kTRUE),
56  fMC(NULL),
57  // 1.) Control histograms:
58  fControlHistogramsList(NULL),
59  fControlHistogramsFlagsPro(NULL),
60  fFillControlHistograms(kFALSE),
61  fControlHistogramsEventList(NULL),
62  fControlHistogramsEventFlagsPro(NULL),
63  fFillControlHistogramsEvent(kFALSE),
64  fGetNumberOfTracksHist(NULL),
65  fGetNumberOfGlobalTracksHist(NULL),
66  fGetNumberOfV0sHist(NULL),
67  fGetNumberOfCascadesHist(NULL),
68  fGetMagneticFieldHist(NULL),
69  fGetEventTypeHist(NULL),
70  fGetCentralityHist(NULL),
71  fGetNContributorsHist(NULL),
72  fGetChi2perNDFHist(NULL),
73  fGetNDaughtersHist(NULL),
74  fControlHistogramsNonIdentifiedParticlesList(NULL),
75  fControlHistogramsNonIdentifiedParticlesFlagsPro(NULL),
76  fFillControlHistogramsNonIdentifiedParticles(kFALSE),
77  fChargeHist(NULL),
78  fGetTPCNclsHist(NULL),
79  fGetTPCsignalNHist(NULL),
80  fGetITSNclsHist(NULL),
81  fdEdxVsPtHist(NULL),
82  fPtHist(NULL),
83  fEtaHist(NULL),
84  fPhiHist(NULL),
85  fMassHist(NULL),
86  fGetFilterMap(NULL),
87  fGetPdgCode(NULL),
88  fControlHistogramsNonIdentifiedParticlesFTSFList(NULL),
89  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro(NULL),
90  fFillControlHistogramsNonIdentifiedParticlesFTSF(kFALSE),
91  fFilterBitFTSF(128),
92  fChargeFTSFHist(NULL),
93  fGetTPCNclsFTSFHist(NULL),
94  fGetTPCsignalNFTSFHist(NULL),
95  fGetITSNclsFTSFHist(NULL),
96  fdEdxVsPtFTSFHist(NULL),
97  fPtFTSFHist(NULL),
98  fEtaFTSFHist(NULL),
99  fPhiFTSFHist(NULL),
100  fMassFTSFHist(NULL),
101  fGetFilterMapFTSF(NULL),
102  fGetPdgCodeFTSF(NULL),
103  fControlHistogramsIdentifiedParticlesList(NULL),
104  fControlHistogramsIdentifiedParticlesFlagsPro(NULL),
105  fFillControlHistogramsIdentifiedParticles(kFALSE),
106  fFillControlHistogramsWithGlobalTrackInfo(kFALSE),
107  fInclusiveSigmaCutsPro(NULL),
108  fExclusiveSigmaCutsPro(NULL),
109  fUseDefaultInclusiveSigmaCuts(kTRUE),
110  fUseDefaultExclusiveSigmaCuts(kTRUE),
111  fControlHistogramsV0sList(NULL),
112  fControlHistogramsV0sFlagsPro(NULL),
113  fFillControlHistogramsV0s(kFALSE),
114  fGetNProngsHist(NULL),
115  fMassK0ShortHist(NULL),
116  fMassLambdaHist(NULL),
117  fMassAntiLambdaHist(NULL),
118  fOpenAngleV0Hist(NULL),
119  fRadiusV0Hist(NULL),
120  fDcaV0ToPrimVertexHist(NULL),
121  fMomV0XHist(NULL),
122  fMomV0YHist(NULL),
123  fMomV0ZHist(NULL),
124  fPtV0Hist(NULL),
125  fPseudoRapV0Hist(NULL),
126  fPAHist(NULL),
127  // 2.) Event-by-event histograms:
128  fEBEHistogramsList(NULL),
129  fEBEObjectsFlagsPro(NULL),
130  //fFillEBEHistograms(kTRUE),
131  fUniqueIDHistEBE(NULL),
132  // 3.) Correlation functions:
133  fCorrelationFunctionsList(NULL),
134  fCorrelationFunctionsFlagsPro(NULL),
135  f2pCorrelationFunctionsFlagsPro(NULL),
136  f3pCorrelationFunctionsFlagsPro(NULL),
137  f4pCorrelationFunctionsFlagsPro(NULL),
138  fFillCorrelationFunctions(kFALSE),
139  fNormalizeCorrelationFunctions(kFALSE),
140  fCorrelationFunctionsIndices(NULL),
141  fFill3pCorrelationFunctions(kFALSE),
142  fFill4pCorrelationFunctions(kFALSE),
143  fNormalizationOption(0),
144  fnMergedBins(-44),
145  // 4.) Background:
146  fBackgroundList(NULL),
147  fBackgroundFlagsPro(NULL),
148  f2pBackgroundFlagsPro(NULL),
149  f3pBackgroundFlagsPro(NULL),
150  f4pBackgroundFlagsPro(NULL),
151  fBackgroundOption(0),
152  fEstimate2pBackground(kFALSE),
153  fEstimate3pBackground(kFALSE),
154  fEstimate4pBackground(kFALSE),
155  fMaxBufferSize1(10),
156  // 5.) Buffers:
157  fBuffersList(NULL),
158  fBuffersFlagsPro(NULL),
159  fFillBuffers(kFALSE),
160  fMaxBuffer(-44),
161  // 6.) QA:
162  fQAList(NULL),
163  fQAFlagsPro(NULL),
164  fFillQA(kFALSE),
165  fBailOutAfterQA(kFALSE),
166  fFillQAEvents(kFALSE),
167  fFillQAParticles(kFALSE),
168  fQAEventsList(NULL),
169  fQAParticlesList(NULL),
170  fQAFilterBitScan(NULL),
171  fQAIDvsFilterBit(NULL),
172  // 7.) Common event cuts:
173  fRejectEventsWithoutPrimaryVertex(kTRUE),
174  fMinMagneticField(0.001),
175  fCutOnNumberOfTracks(kFALSE),
176  fMinNumberOfTracks(-44),
177  fMaxNumberOfTracks(-44),
178  fCutOnNumberOfGlobalTracks(kFALSE),
179  fMinNumberOfGlobalTracks(-44),
180  fMaxNumberOfGlobalTracks(-44),
181  fCutOnNumberOfV0s(kFALSE),
182  fMinNumberOfV0s(-44),
183  fMaxNumberOfV0s(-44),
184  fCutOnNumberOfCascades(kFALSE),
185  fMinNumberOfCascades(-44),
186  fMaxNumberOfCascades(-44),
187  fCutOnVertexX(kFALSE),
188  fMinVertexX(-44.),
189  fMaxVertexX(-44.),
190  fCutOnVertexY(kFALSE),
191  fMinVertexY(-44.),
192  fMaxVertexY(-44.),
193  fCutOnVertexZ(kFALSE),
194  fMinVertexZ(-44.),
195  fMaxVertexZ(-44.),
196  fCutOnNContributors(kFALSE),
197  fMinNContributors(-44),
198  fMaxNContributors(-44),
199  // 8.) Global track cuts:
200  fGlobalTrackCutsList(NULL),
201  fGlobalTrackCutsFlagsPro(NULL),
202  fApplyGlobalTrackCuts(kFALSE),
203  // *.) Testing new ways to calculate correlation functions:
204  fCorrelationFunctionsTESTList(NULL),
205  fCorrelationFunctionsTESTFlagsPro(NULL),
206  fBoost(kFALSE),
207  fBoostVelocity(0.,0.,0.),
208  fnQ2bins(100),
209  fnQ2min(0.),
210  fnQ2max(10.),
211  fnQ3bins(100),
212  fnQ3min(0.),
213  fnQ3max(10.),
214  // *.) Testing new ways to calculate background:
215  fBackgroundTESTList(NULL),
216  fBackgroundTESTFlagsPro(NULL),
217  // *.) 'hybrid approach':
218  fHybridApproachList(NULL),
219  fHybridApproachFlagsPro(NULL),
220  fDoHybridApproach(kFALSE),
221  // *.) Online monitoring:
222  fOnlineMonitoring(kFALSE),
223  fUpdateOutputFile(kFALSE),
224  fUpdateFrequency(-44),
225  fUpdateWhichOutputFile(NULL),
226  fMaxNumberOfEvents(-44),
227  // *.) Debugging:
228  fDoSomeDebugging(kFALSE),
229  fWaitForSpecifiedEvent(kFALSE),
230  fRun(0),
231  fBunchCross(0),
232  fOrbit(0),
233  fPeriod(0)
234  {
235  // Constructor.
236 
237  AliDebug(2,"AliAnalysisTaskMultiparticleFemtoscopy::AliAnalysisTaskMultiparticleFemtoscopy(const char *name, Bool_t useParticleWeights)");
238 
239  // Base list:
240  fHistList = new TList();
241  fHistList->SetName("MPF"); // MultiParticle Femtoscopy
242  fHistList->SetOwner(kTRUE);
243 
244  // Initialize all arrays:
245  this->InitializeArrays();
246  this->InitializeArraysForControlHistograms();
247  this->InitializeArraysForEBEObjects();
248  this->InitializeArraysForCorrelationFunctions();
249  this->InitializeArraysForBackground();
250  this->InitializeArraysForBuffers();
251  this->InitializeArraysForQA();
252  this->InitializeArraysForGlobalTrackCuts();
253  this->InitializeArraysForCorrelationFunctionsTEST();
254  this->InitializeArraysForBackgroundTEST();
255 
256  // Define input and output slots here
257  // Input slot #0 works with an AliFlowEventSimple
258  //DefineInput(0, AliFlowEventSimple::Class());
259  // Input slot #1 is needed for the weights input file:
260  //if(useParticleWeights)
261  //{
262  // DefineInput(1, TList::Class());
263  //}
264  // Output slot #0 is reserved
265  // Output slot #1 writes into a TList container
266 
267  DefineOutput(1, TList::Class());
268 
269  if(useParticleWeights)
270  {
271  // TBI add support eventually for particle weights
272  }
273 
274 } // AliAnalysisTaskMultiparticleFemtoscopy::AliAnalysisTaskMultiparticleFemtoscopy(const char *name, Bool_t useParticleWeights):
275 
276 //================================================================================================================
277 
280  fHistList(NULL),
281  fAnalysisType(NULL),
282  fPIDResponse(NULL),
283  fMaxNoGlobalTracksAOD(5), // TBI this is landmine
284  fProcessBothKineAndReco(kFALSE),
285  fProcessOnlyKine(kFALSE),
286  fProcessOnlyReco(kFALSE),
287  fRejectFakeTracks(kTRUE),
288  fMC(NULL),
289  // 1.) Control histograms:
290  fControlHistogramsList(NULL),
291  fControlHistogramsFlagsPro(NULL),
292  fFillControlHistograms(kFALSE),
293  fControlHistogramsEventList(NULL),
294  fControlHistogramsEventFlagsPro(NULL),
295  fFillControlHistogramsEvent(kFALSE),
296  fGetNumberOfTracksHist(NULL),
297  fGetNumberOfGlobalTracksHist(NULL),
298  fGetNumberOfV0sHist(NULL),
299  fGetNumberOfCascadesHist(NULL),
300  fGetMagneticFieldHist(NULL),
301  fGetEventTypeHist(NULL),
302  fGetCentralityHist(NULL),
303  fGetNContributorsHist(NULL),
304  fGetChi2perNDFHist(NULL),
305  fGetNDaughtersHist(NULL),
306  fControlHistogramsNonIdentifiedParticlesList(NULL),
307  fControlHistogramsNonIdentifiedParticlesFlagsPro(NULL),
308  fFillControlHistogramsNonIdentifiedParticles(kFALSE),
309  fChargeHist(NULL),
310  fGetTPCNclsHist(NULL),
311  fGetTPCsignalNHist(NULL),
312  fGetITSNclsHist(NULL),
313  fdEdxVsPtHist(NULL),
314  fPtHist(NULL),
315  fEtaHist(NULL),
316  fPhiHist(NULL),
317  fMassHist(NULL),
318  fGetFilterMap(NULL),
319  fGetPdgCode(NULL),
320  fControlHistogramsNonIdentifiedParticlesFTSFList(NULL),
321  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro(NULL),
322  fFillControlHistogramsNonIdentifiedParticlesFTSF(kFALSE),
323  fFilterBitFTSF(128),
324  fChargeFTSFHist(NULL),
325  fGetTPCNclsFTSFHist(NULL),
326  fGetTPCsignalNFTSFHist(NULL),
327  fGetITSNclsFTSFHist(NULL),
328  fdEdxVsPtFTSFHist(NULL),
329  fPtFTSFHist(NULL),
330  fEtaFTSFHist(NULL),
331  fPhiFTSFHist(NULL),
332  fMassFTSFHist(NULL),
333  fGetFilterMapFTSF(NULL),
334  fGetPdgCodeFTSF(NULL),
335  fControlHistogramsIdentifiedParticlesList(NULL),
336  fControlHistogramsIdentifiedParticlesFlagsPro(NULL),
337  fFillControlHistogramsIdentifiedParticles(kFALSE),
338  fFillControlHistogramsWithGlobalTrackInfo(kFALSE),
339  fInclusiveSigmaCutsPro(NULL),
340  fExclusiveSigmaCutsPro(NULL),
341  fUseDefaultInclusiveSigmaCuts(kFALSE),
342  fUseDefaultExclusiveSigmaCuts(kFALSE),
343  fControlHistogramsV0sList(NULL),
344  fControlHistogramsV0sFlagsPro(NULL),
345  fFillControlHistogramsV0s(kFALSE),
346  fGetNProngsHist(NULL),
347  fMassK0ShortHist(NULL),
348  fMassLambdaHist(NULL),
349  fMassAntiLambdaHist(NULL),
350  fOpenAngleV0Hist(NULL),
351  fRadiusV0Hist(NULL),
352  fDcaV0ToPrimVertexHist(NULL),
353  fMomV0XHist(NULL),
354  fMomV0YHist(NULL),
355  fMomV0ZHist(NULL),
356  fPtV0Hist(NULL),
357  fPseudoRapV0Hist(NULL),
358  fPAHist(NULL),
359  // 2.) Event-by-event histograms:
360  fEBEHistogramsList(NULL),
361  fEBEObjectsFlagsPro(NULL),
362  //fFillEBEHistograms(kFALSE),
363  fUniqueIDHistEBE(NULL),
364  // 3.) Correlation functions:
365  fCorrelationFunctionsList(NULL),
366  fCorrelationFunctionsFlagsPro(NULL),
367  f2pCorrelationFunctionsFlagsPro(NULL),
368  f3pCorrelationFunctionsFlagsPro(NULL),
369  f4pCorrelationFunctionsFlagsPro(NULL),
370  fFillCorrelationFunctions(kFALSE),
371  fNormalizeCorrelationFunctions(kFALSE),
372  fCorrelationFunctionsIndices(NULL),
373  fFill3pCorrelationFunctions(kFALSE),
374  fFill4pCorrelationFunctions(kFALSE),
375  fNormalizationOption(0),
376  fnMergedBins(-44),
377  // 4.) Background:
378  fBackgroundList(NULL),
379  fBackgroundFlagsPro(NULL),
380  f2pBackgroundFlagsPro(NULL),
381  f3pBackgroundFlagsPro(NULL),
382  f4pBackgroundFlagsPro(NULL),
383  fBackgroundOption(0),
384  fEstimate2pBackground(kFALSE),
385  fEstimate3pBackground(kFALSE),
386  fEstimate4pBackground(kFALSE),
387  fMaxBufferSize1(0),
388  // 5.) Buffers:
389  fBuffersList(NULL),
390  fBuffersFlagsPro(NULL),
391  fFillBuffers(kFALSE),
392  fMaxBuffer(-44),
393  // 6.) QA:
394  fQAList(NULL),
395  fQAFlagsPro(NULL),
396  fFillQA(kFALSE),
397  fBailOutAfterQA(kFALSE),
398  fFillQAEvents(kFALSE),
399  fFillQAParticles(kFALSE),
400  fQAEventsList(NULL),
401  fQAParticlesList(NULL),
402  fQAFilterBitScan(NULL),
403  fQAIDvsFilterBit(NULL),
404  // 7.) Common event cuts:
405  fRejectEventsWithoutPrimaryVertex(kFALSE),
406  fMinMagneticField(0.001),
407  fCutOnNumberOfTracks(kFALSE),
408  fMinNumberOfTracks(-44),
409  fMaxNumberOfTracks(-44),
410  fCutOnNumberOfGlobalTracks(kFALSE),
411  fMinNumberOfGlobalTracks(-44),
412  fMaxNumberOfGlobalTracks(-44),
413  fCutOnNumberOfV0s(kFALSE),
414  fMinNumberOfV0s(-44),
415  fMaxNumberOfV0s(-44),
416  fCutOnNumberOfCascades(kFALSE),
417  fMinNumberOfCascades(-44),
418  fMaxNumberOfCascades(-44),
419  fCutOnVertexX(kFALSE),
420  fMinVertexX(-44.),
421  fMaxVertexX(-44.),
422  fCutOnVertexY(kFALSE),
423  fMinVertexY(-44.),
424  fMaxVertexY(-44.),
425  fCutOnVertexZ(kFALSE),
426  fMinVertexZ(-44.),
427  fMaxVertexZ(-44.),
428  fCutOnNContributors(kFALSE),
429  fMinNContributors(-44),
430  fMaxNContributors(-44),
431  // 8.) Global track cuts:
432  fGlobalTrackCutsList(NULL),
433  fGlobalTrackCutsFlagsPro(NULL),
434  fApplyGlobalTrackCuts(kFALSE),
435  // *.) Testing new ways to calculate correlation functions:
436  fCorrelationFunctionsTESTList(NULL),
437  fCorrelationFunctionsTESTFlagsPro(NULL),
438  fBoost(kFALSE),
439  fBoostVelocity(0.,0.,0.),
440  fnQ2bins(-44),
441  fnQ2min(-44.),
442  fnQ2max(-44.),
443  fnQ3bins(-44),
444  fnQ3min(-44.),
445  fnQ3max(-44.),
446  // *.) Testing new ways to calculate background:
447  fBackgroundTESTList(NULL),
448  fBackgroundTESTFlagsPro(NULL),
449  // *.) 'hybrid approach':
450  fHybridApproachList(NULL),
451  fHybridApproachFlagsPro(NULL),
452  fDoHybridApproach(kFALSE),
453  // *.) Online monitoring:
454  fOnlineMonitoring(kFALSE),
455  fUpdateOutputFile(kFALSE),
456  fUpdateFrequency(-44),
457  fUpdateWhichOutputFile(NULL),
458  fMaxNumberOfEvents(-44),
459  // *.) Debugging:
460  fDoSomeDebugging(kFALSE),
461  fWaitForSpecifiedEvent(kFALSE),
462  fRun(0),
463  fBunchCross(0),
464  fOrbit(0),
465  fPeriod(0)
466 {
467  // Dummy constructor.
468 
469  AliDebug(2,"AliAnalysisTaskMultiparticleFemtoscopy::AliAnalysisTaskMultiparticleFemtoscopy()");
470 
471 } // AliAnalysisTaskMultiparticleFemtoscopy::AliAnalysisTaskMultiparticleFemtoscopy():
472 
473 //================================================================================================================
474 
476 {
477  // Destructor.
478 
479  if(fHistList) delete fHistList;
480  if(fPIDResponse) delete fPIDResponse;
481 
482  for(Int_t index=0;index<fMaxNoGlobalTracksAOD;index++)
483  {
484  if(fGlobalTracksAOD[index]) delete fGlobalTracksAOD[index];
485  }
486 
487 } // AliAnalysisTaskMultiparticleFemtoscopy::~AliAnalysisTaskMultiparticleFemtoscopy()
488 
489 //================================================================================================================
490 
492 {
493  // Called at every worker node to initialize.
494 
495  // a) Insanity checks;
496  // b) Trick to avoid name clashes, part 1;
497  // c) Set all flags which are not directly set via setters;
498  // d) Book and nest all lists;
499  // e) Book all objects;
500  // f) Trick to avoid name clashes, part 2.
501 
502  // a) Insanity checks:
504 
505  // b) Trick to avoid name clashes, part 1:
506  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
507  TH1::AddDirectory(kFALSE);
508 
509  // c) Set all flags which are not directly set via setters:
512 
513  // d) Book and nest all lists:
514  this->BookAndNestAllLists();
515 
516  // e) Book all objects:
517  this->BookEverything();
522  this->BookEverythingForBuffers();
523  this->BookEverythingForQA();
528 
529  // f) Trick to avoid name clashes, part 2:
530  TH1::AddDirectory(oldHistAddStatus);
531 
532  PostData(1,fHistList);
533 
534 } // void AliAnalysisTaskMultiparticleFemtoscopy::UserCreateOutputObjects()
535 
536 //================================================================================================================
537 
539 {
540  // Main loop (called for each event).
541 
542  // a) Determine Ali{MC,ESD,AOD}Event and the analysis type;
543  // b) Insanity checks; // TBI disabled at the moment (see implementation)
544  // c) QA;
545  // d) Start analysis over MC, AOD, ESD, MC_AOD or MC_ESD;
546  // e) Reset event-by-event objects;
547  // f) PostData;
548  // g) Online monitoring.
549 
550  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::UserExec(Option_t *)";
551 
552  // a) Determine Ali{MC,ESD,AOD}Event and the analysis type:
553  AliMCEvent *aMC = MCEvent(); // from TaskSE
554  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
555  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
557  {
558  Fatal(sMethodName.Data(),"One (and only one!) of fProcessBothKineAndReco, fProcessOnlyKine and fProcessOnlyReco must be kTRUE in AddTask* macro!!!!");
559  }
560  else if(fProcessBothKineAndReco)
561  {
562  if(aMC && aAOD){*fAnalysisType="MC_AOD";}
563  else if(aMC && aESD){*fAnalysisType="MC_ESD";}
564  else{Fatal(sMethodName.Data(),"fProcessBothKineAndReco is kTRUE, but TBI ...");}
565  } // if(fProcessBothKineAndReco)
566  else if(fProcessOnlyKine)
567  {
568  if(aMC){*fAnalysisType="MC";}
569  else{Fatal(sMethodName.Data(),"fProcessOnlyKine is kTRUE, but TBI ...");}
570  }
571  else if(fProcessOnlyReco)
572  {
573  if(aAOD){*fAnalysisType="AOD";}
574  else if(aESD){*fAnalysisType="ESD";}
575  else{Fatal(sMethodName.Data(),"fProcessOnlyReco is kTRUE, but TBI ...");}
576  }
577  //cout<<Form("\n=> AATMF::UserExec(): fAnalysisType = \"%s\"\n",fAnalysisType->Data())<<endl;
578 
579  // b) Insanity checks:
581 
582  // c) QA:
583  if(fFillQA)
584  {
585  if(fAnalysisType->EqualTo("MC")&&aMC){QA(aMC);}
586  else if(fAnalysisType->EqualTo("AOD")&&aAOD){QA(aAOD);}
587  else if(fAnalysisType->EqualTo("ESD")&&aESD){QA(aESD);}
588  // TBI Do I have to support also cases "MC_AOD" and"MC_ESD"?
589  } // if(fFillQA)
590  if(fBailOutAfterQA){return;} // by default, this flag is set to kFALSE
591 
592  // d) Start analysis over MC, AOD, ESD, MC_AOD or MC_ESD:
593  if(fAnalysisType->Contains("MC"))
594  {
595  fMC = aMC;
596  }
597  if(fAnalysisType->EqualTo("MC"))
598  {
599  if(aMC){MC(aMC);}
600  }
601  else if(fAnalysisType->EqualTo("AOD"))
602  {
603  if(aAOD){AOD(aAOD);}
604  }
605  else if(fAnalysisType->EqualTo("ESD"))
606  {
607  if(aESD){ESD(aESD);}
608  }
609  else if(fAnalysisType->EqualTo("MC_AOD"))
610  {
611  if(aMC&&aAOD){AOD(aAOD);}
612  }
613  else if(fAnalysisType->EqualTo("MC_ESD"))
614  {
615  if(aMC&&aESD){ESD(aESD);}
616  }
617 
618  // e) Reset event-by-event objects:
619  this->ResetEBEObjects();
620 
621  // f) PostData:
622  PostData(1,fHistList);
623 
624  // g) Online monitoring:
626 
627 } // void AliAnalysisTaskMultiparticleFemtoscopy::UserExec(Option_t *)
628 
629 //================================================================================================================
630 
632 {
633  // At the end of the day...
634 
635  // a) Get pointer to the grandmother of all lists "fHistList";
636  // b) Get pointers to all objects in "fHistList";
637  // c) Normalize correlation functions;
638  // d) Dump the results.
639 
640  // a) Get pointer to the grandmother of all lists:
641  fHistList = (TList*)GetOutputData(1); // TBI doesn't work locally :'(
642  if(!fHistList){exit(1);}
643 
644  // b) Get pointers to all objects in "fHistList":
646 
647  // c) Normalize correlation functions:
649 
650  // d) Dump the results:
651  //TDirectoryFile *df = new TDirectoryFile("outputMPFanalysis","");
652  //df->Add(fHistList);
653  TFile *f = new TFile("AnalysisResults.root","RECREATE");
654  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
655 
656  delete f;
657 
658 } // end of void AliAnalysisTaskMultiparticleFemtoscopy::Terminate(Option_t *)
659 
660 //================================================================================================================
661 
663 {
664  // Local Quality Assurance. TBI
665 
666  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::QA(AliVEvent *ave)";
667 
668  // a) Determine Ali{MC,ESD,AOD}Event;
669 
670  // a) Determine Ali{MC,ESD,AOD}Event:
671  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
672  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
673  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
674 
675  if(aMC)
676  {
677  // TBI
678  }
679  else if(aESD)
680  {
681  // TBI
682  }
683  else if(aAOD)
684  {
685  // a) Loop over 'aTracks'
686  // b) Cut directly on aTracks:
687  // 1.) e.g. is it TPC-only => via filter bits
688  // 2.) kinematics => implement a function which from aTrack access its kinematics, and it cuts on it directly
689  // 3.) PID => implement a function which from aTrack access its gTracks, and from it the PID, and then it cuts on it directly
690  // 4.) global track parameters => implement a function which from aTrack access its gTracks, and then it cuts on its parameters directly
691  // c)
692 
693 
694  // a) Loop over 'aTracks':
695  Int_t nTracks = aAOD->GetNumberOfTracks();
696  for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
697  {
698  AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
699  Int_t nFilterBits = 14;
700  for(Int_t fb=0;fb<nFilterBits;fb++) // 'fb' is a 'left shifting opetator', i.e. 1<<fb = filterbit
701  {
702  if(atrack->TestFilterBit(1<<fb))
703  {
704  fQAFilterBitScan->Fill(fb);
705  fQAIDvsFilterBit->Fill(fb,atrack->GetID());
706  } // if(atrack->TestFilterBit(1<<fb))
707  } // for(Int_t fb=1;fb<nFilterBits;fb++)
708 
709 
710 
711  // t.b.c.
712 
713 
714 
715  } // for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
716 
717 
718 
719 
720 
721 
722 
723 /*
724 A) PARTICLE CUTS:
725 
726 - Assuming two sets of histograms for track cuts:
727  - fQAParticleHist[0][distribution_index][cut_index] => "before rain"
728  - fQAParticleHist[1][distribution_index][cut_index] => "after rain"
729 - Used in conjunction with four set of cuts:
730  - fParticleCuts[0] => direct cuts on aTrack
731  - fParticleCuts[1] => given aTrack, these are the cuts on corresponding gTrack
732  - fParticleCutsQA[0] => QA direct cuts on aTrack
733  - fParticleCutsQA[1] => given aTrack, these are the QA cuts on corresponding gTrack
734 - Example 1:
735  - Given aTrack
736  - fParticleCuts[0] && fParticleCuts[1] => Fill QA histos "before"
737  - fParticleCutsQA[0] && fParticleCutsQA[1] => Fill QA histos "after"
738 
739 
740 B) EVENT CUTS:
741 
742 - Assuming two sets of histograms for event cuts:
743  - fQAEventHist[0][distribution_index][cut_number] => "before rain"
744  - fQAEventHist[1][distribution_index][cut_number] => "after rain"
745 - Used in conjunction with two set of cuts:
746  - fEventCuts[0] => "default" cuts on event variables
747  - fEventCuts[1] => given the "default" cuts on event variables, cut event further
748 
749 - Example 1:
750  - Given aAOD
751  - if fEventCuts[0] => Fill QA histos "before"
752  - if fEventCuts[1] => Fill QA histos "after"
753 
754 */
755 
756 
757  }
758 
759  return;
760 
761 } // void AliAnalysisTaskMultiparticleFemtoscopy::QA(AliVEvent *ave)
762 
763 //================================================================================================================
764 
766 {
767  // Normalize correlation functions with the background.
768 
769  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::NormalizeCorrelationFunctions()";
770 
771  // Say something nice to the audience...
772  cout<<endl;
773  cout<<"=> Normalization option:"<<endl;
774  switch(fNormalizationOption)
775  {
776  case 0:
777  cout<<"\"just scale\""<<endl;
778  break;
779 
780  case 1:
781  cout<<Form(" \"use concrete interval\": min = %f, max = %f",fNormalizationInterval[0],fNormalizationInterval[1])<<endl;
782  cout<<"TBI: not implemented yet"<<endl;
783  exit(0);
784  break;
785 
786  default:
787  cout<<Form("And the fatal 'fNormalizationOption' value is... %d. Congratulations!!",fNormalizationOption)<<endl;
788  Fatal(sMethodName.Data(),"switch(fNormalizationOption)");
789  } // switch(fNormalizationOption)
790  if(fnMergedBins!=-44)
791  {
792  cout<<Form("=> Histograms will be rebinned: fnMergedBins = %d",fnMergedBins)<<endl;
793  }
794  cout<<endl;
795 
796  cout<<"TBI: not finalized yet"<<endl;
797 
798  for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
799  {
800  for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
801  {
802  if(fCorrelationFunctions[pid1][pid2] && fCorrelationFunctions[pid1][pid2]->GetEntries() > 0 &&
803  f2pBackground[pid1][pid2] && f2pBackground[pid1][pid2]->GetEntries() > 0)
804  {
805  fCorrelationFunctions[pid1][pid2]->Divide(f2pBackground[pid1][pid2]);
806  } // if(..)
807  } // for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
808  } // for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
809 
810 
811  for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
812  {
813  for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
814  {
815  for(Int_t pid3=0;pid3<10;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
816  {
817 
818  if(f3pCorrelationFunctions[pid1][pid2][pid3] && f3pCorrelationFunctions[pid1][pid2][pid3]->GetEntries() > 0 &&
819  f3pBackground[pid1][pid2][pid3] && f3pBackground[pid1][pid2][pid3]->GetEntries() > 0)
820  {
821  f3pCorrelationFunctions[pid1][pid2][pid3]->Divide(f3pBackground[pid1][pid2][pid3]);
822  } // if(..)
823 
824  } // for(Int_t pid3=0;pid3<10;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
825  } // for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
826  } // for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
827 
828 
829  for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
830  {
831  for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
832  {
833  for(Int_t pid3=0;pid3<10;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
834  {
835  for(Int_t pid4=0;pid4<10;pid4++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
836  {
837 
838  if(f4pCorrelationFunctions[pid1][pid2][pid3][pid4] && f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->GetEntries() > 0 &&
839  f4pBackground[pid1][pid2][pid3][pid4] && f4pBackground[pid1][pid2][pid3][pid4]->GetEntries() > 0)
840  {
841  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->Divide(f4pBackground[pid1][pid2][pid3][pid4]);
842  } // if(..)
843  } // for(Int_t pid4=0;pid4<10;pid4++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
844  } // for(Int_t pid3=0;pid3<10;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
845  } // for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
846  } // for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
847 
848 
849 } // void AliAnalysisTaskMultiparticleFemtoscopy::NormalizeCorrelationFunctions()
850 
851 //================================================================================================================
852 
854 {
855  // Fill control histograms for global event observables.
856 
857  // TBI: not really validated, synchronize all if statements eventually with the analogous ones in UserExec(...)
858 
859  // To do:
860  // 1) Add support for MC and ESD.
861 
862  // a) Determine Ali{MC,ESD,AOD}Event;
863  // b) ...
864 
865  // a) Determine Ali{MC,ESD,AOD}Event:
866  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
867  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
868  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
869 
870  if(aMC && fProcessOnlyKine) // TBI not validated
871  {
872  // MC event:
873  fGetNumberOfTracksHist->Fill(aMC->GetNumberOfTracks());
874  // ... TBI fill the rest
875  } // if(aMC)
876  else if(aESD)
877  {
878  // TBI
879  } // else if(aESD)
880  else if(aAOD) // TBI not validated, nevertheless works at the moment only for fProcessBothKineAndReco = kTRUE and fProcessOnlyReco = kTRUE
881  {
882  // AOD event:
883  fGetNumberOfTracksHist->Fill(aAOD->GetNumberOfTracks()); // not all tracks are unique, see comments in function GlobalTracksAOD(AliAODEvent *aAOD, Int_t index)
884  fGetNumberOfGlobalTracksHist->Fill(fGlobalTracksAOD[0]->GetSize()); // this is then my multiplicity. It shall be OK, since in fGlobalTracksAOD[0] I have filled only tracks with ID >= 0
885  fGetNumberOfV0sHist->Fill(aAOD->GetNumberOfV0s()); // some V0s share the daughter, in this histogram they are not separated
886  fGetNumberOfCascadesHist->Fill(aAOD->GetNumberOfCascades()); // TBI not validated
887  fGetMagneticFieldHist->Fill(aAOD->GetMagneticField()); // TBI not validated
888  fGetEventTypeHist->Fill(aAOD->GetEventType()); // TBI not validated
889  if(aAOD->GetCentrality())
890  {
891  fGetCentralityHist->Fill(aAOD->GetCentrality()->GetCentralityPercentile("V0M")); // TBI not validated
892  }
893  // AOD primary vertex:
894  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
895  fVertexXYZ[0]->Fill(avtx->GetX());
896  fVertexXYZ[1]->Fill(avtx->GetY());
897  fVertexXYZ[2]->Fill(avtx->GetZ());
898  fGetNContributorsHist->Fill(avtx->GetNContributors()); // TBI not validated
899  fGetChi2perNDFHist->Fill(avtx->GetChi2perNDF()); // TBI not validated
900  fGetNDaughtersHist->Fill(avtx->GetNDaughters()); // TBI not validated
901  } // else if(aAOD)
902 
903 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsEvent(AliVEvent *ave)
904 
905 //================================================================================================================
906 
908 {
909  // Fill control histograms for particles.
910  // Remark: The idea is to have one loop over the particles and to fill everything which needs to be filled.
911 
912  // To do:
913  // a) Not really validated for fProcessOnlyKine = kTRUE
914  // b) ESD not supported yet
915 
916  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsParticle(AliVEvent *ave)";
917 
918  // a) Determine Ali{MC,ESD,AOD}Event:
919  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
920  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
921  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
922 
923  if(aMC && fProcessOnlyKine) // TBI not validated
924  {
925  // MC tracks:
926  for(Int_t iTrack=0;iTrack<aMC->GetNumberOfTracks();iTrack++)
927  {
928  // a) Determine "amctrack";
929  // b) Insanity checks for "amctrack";
930  // c) Fill the control histograms for non-identified particles;
931  // d) Fill the control histograms for identified particles.
932 
933  // a) Determine "amctrack":
934  AliAODMCParticle *amcparticle = (AliAODMCParticle*)aMC->GetTrack(iTrack);
935 
936  // b) Insanity checks for "amctrack":
937  if(!amcparticle){Fatal(sMethodName.Data(),"!amctrack");}
938 
939  // c) Fill the control histograms for non-identified particles:
941 
942  // d) Fill the control histograms for identified particles:
944 
945  } // for(Int_t iTrack=0;iTrack<aMC->GetNumberOfTracks();iTrack++)
946  } // if(aMC)
947 
948  //=================================================================================================================
949 
950  else if(aESD)
951  {
952  // TBI
953  } // else if(aESD)
954 
955  //=================================================================================================================
956 
957  else if(aAOD) // TBI: Works for fProcessBothKineAndReco = kTRUE and fProcessOnlyReco = kTRUE
958  {
959  // AOD tracks:
960  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
961  {
962  // a) Determine "atrack" (a.k.a. "any track in AOD");
963  // b) Insanity checks for "atrack";
964  // c) Determine the corresponding "gtrack" (a.k.a. "normal global" track);
965  // d) Insanity checks for "gtrack";
966  // e) Fill the control histograms for non-identified particles; TBI: not validated
967  // f) Fill the control histograms for non-identified particles (f.t.s.f); // for the specific filterbit TBI: not validated
968  // g) Fill the control histograms for identified particles (validated, works both for fProcessBothKineAndReco = kTRUE and fProcessOnlyReco = kTRUE, shall be used for purities)
969 
970  // a) Determine "atrack" (a.k.a. "any track in AOD"):
971  AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
972 
973  // b) Insanity checks for "atrack":
974  if(!atrack){Fatal(sMethodName.Data(),"!atrack");}
975  Int_t n = InsanityChecksForTracks(atrack); // return values are documented in this function
976  if(0!=n){Fatal(sMethodName.Data(),"InsanityChecksForTracks(atrack), n = %d",n);}
977 
978  // c) Determine the corresponding "gtrack" (a.k.a. "normal global" track):
979  if(0 == fGlobalTracksAOD[0]->GetSize()){GlobalTracksAOD(aAOD,0);} // TBI most likely, this is an overkill
980  if(0 == fGlobalTracksAOD[0]->GetSize()){return;}
981  Int_t id = atrack->GetID();
982  AliAODTrack *gtrack = dynamic_cast<AliAODTrack*>(id>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id+1))));
983 
984  // d) Insanity checks for "gtrack":
985  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
986  Int_t ng = InsanityChecksForGlobalTracks(gtrack); // return values are documented in this function
987  if(0!=ng){Fatal(sMethodName.Data(),"InsanityChecksForGlobalTracks(gtrack), ng = %d",ng);}
988 
989  // e) Fill the control histograms for non-identified particles: TBI is there a more elegant solution to avoid double counting, besides checking for ID?
990  if(fFillControlHistogramsNonIdentifiedParticles && atrack->GetID()>=0){this->FillControlHistogramsNonIdentifiedParticles(gtrack);} // for AOD, these are 'normal global' tracks
991 
992  // f) Fill the control histograms for non-identified particles (f.t.s.f): // for the specific filterbit TBI validated at the moment only for TPC-only, due to h.w. atrack->GetID()<0
993  if(fFillControlHistogramsNonIdentifiedParticlesFTSF && atrack->GetID()<0 && atrack->TestFilterBit(fFilterBitFTSF)){this->FillControlHistogramsNonIdentifiedParticlesFTSF(atrack);}
994 
995  // g) Fill the control histograms for identified particles:
997 
998  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
999 
1000  } // else if(aAOD)
1001 
1002 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsParticle(AliVEvent *ave)
1003 
1004 //================================================================================================================
1005 
1007 {
1008  // Fill control histograms for non-identified particles.
1009 
1010  // a) Insanity checks;
1011  // b) Check out common selection criteria for 'normal global' tracks;
1012  // c) Fill control histograms.
1013 
1014  // a) Insanity checks:
1015  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticles(AliAODTrack *gtrack)";
1016  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");} // TBI keep this for some time, eventually just continue
1017 
1018  // b) Check out common selection criteria for 'normal global' tracks:
1019  if(!PassesGlobalTrackCuts(gtrack)){return;} // TBI in the method PassesGlobalTrackCuts track is hardwired to AliAODTrack. Try to generalize
1020 
1021  // c) Fill control histograms:
1022  fChargeHist->Fill(gtrack->Charge()+0.5); // see how this histogram was booked for convention used
1023  fGetTPCNclsHist->Fill(gtrack->GetTPCNcls());
1024  fGetTPCsignalNHist->Fill(gtrack->GetTPCsignalN());
1025  fGetITSNclsHist->Fill(gtrack->GetITSNcls());
1026  fdEdxVsPtHist->Fill(gtrack->GetTPCmomentum(),gtrack->GetTPCsignal());
1027  fPtHist->Fill(gtrack->Pt());
1028  fEtaHist->Fill(gtrack->Eta());
1029  fPhiHist->Fill(gtrack->Phi());
1030  fMassHist->Fill(gtrack->M());
1031 
1032 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticles(AliAODTrack *gtrack)
1033 
1034 //================================================================================================================
1035 
1037 {
1038  // Fill control histograms for all charged MC particles.
1039 
1040  // TBI: re-validate
1041 
1042  // a) Insanity checks;
1043  // b) Check cut selection criteria;
1044  // c) Fill control histograms.
1045 
1046  // a) Insanity checks:
1047  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticles(AliAODMCParticle *amcparticle)";
1048  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
1049 
1050  // b) Check cut selection criteria:
1051  //if(!PassesGlobalTrackCuts(gtrack)){return;} // TBI in the method PassesGlobalTrackCuts track is hardwired to AliAODTrack. Try to generalize
1052 
1053  // c) Fill control histograms:
1054  fChargeHist->Fill(amcparticle->Charge()+0.5); // see how this histogram was booked for convention used
1055  fPtHist->Fill(amcparticle->Pt());
1056  fEtaHist->Fill(amcparticle->Eta());
1057  fPhiHist->Fill(amcparticle->Phi());
1058  fMassHist->Fill(amcparticle->M());
1059  fGetPdgCode->Fill(amcparticle->GetPdgCode());
1060 
1061 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticles(AliAODMCParticle *amcparticle)
1062 
1063 //================================================================================================================
1064 
1066 {
1067  // Fill control histograms for non identified particles, for the specified filterbit.
1068 
1069  // a) Insanity checks;
1070  // b) Check out the selection criteria for atracks;
1071  // c) Fill control histograms.
1072 
1073  // a) Insanity checks:
1074  TString sMethodName = "AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticlesFTSF(AliAODTrack *atrack)";
1075  if(!atrack){Fatal(sMethodName.Data(),"!atrack");} // TBI keep this for some time, eventually just continue
1076 
1077  // b) Check out the selection criteria for atracks:
1078  if(!PassesCommonTrackCuts(atrack)){return;}
1079 
1080  // c) Fill control histograms.
1081  fChargeFTSFHist->Fill(atrack->Charge()+0.5); // see how this histogram was booked for convention used
1082  fGetTPCNclsFTSFHist->Fill(atrack->GetTPCNcls());
1083  fGetTPCsignalNFTSFHist->Fill(atrack->GetTPCsignalN());
1084  fGetITSNclsFTSFHist->Fill(atrack->GetITSNcls());
1085  fdEdxVsPtFTSFHist->Fill(atrack->GetTPCmomentum(),atrack->GetTPCsignal());
1086  fPtFTSFHist->Fill(atrack->Pt());
1087  fEtaFTSFHist->Fill(atrack->Eta());
1088  fPhiFTSFHist->Fill(atrack->Phi());
1089  fMassFTSFHist->Fill(atrack->M());
1090 
1091 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticlesFTSF(AliAODTrack *atrack)
1092 
1093 //================================================================================================================
1094 
1096 {
1097  // Fill control histograms for identified particles.
1098 
1099  // Remark 0: if fFillControlHistogramsWithGlobalTrackInfo = kTRUE track parameters are taken from 'gtrack', otherwise (and by default) from 'atrack'
1100  // Remark 1: if fFillControlHistogramsWithGlobalTrackInfo = kTRUE, only PassesGlobalTrackCuts(gtrack) is checked
1101  // if fFillControlHistogramsWithGlobalTrackInfo = kFALSE, only PassesCommonTrackCuts(atrack) is checked.
1102  // In both cases, PID functions Pion(...), etc. are called for 'gtrack', as only it has a PID info
1103 
1104  // a) Insanity checks;
1105  // b) Check cut selection criteria;
1106  // c) Fill control histograms.
1107 
1108  // a) Insanity checks:
1109  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsIdentifiedParticles(AliAODTrack *atrack, AliAODTrack *gtrack)";
1110  if(!atrack){Fatal(sMethodName.Data(),"!atrack");} // TBI keep this for some time, eventually just continue
1111  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");} // TBI keep this for some time, eventually just continue
1112 
1113  // b) Check cut selection criteria:
1115  {
1116  if(atrack->GetID()<0){return;} // this is needed to avoid double counting, I think. TBI well, re-think...
1117  if(!PassesGlobalTrackCuts(gtrack)){return;}
1118  }
1119  else
1120  {
1121  if(!PassesCommonTrackCuts(atrack)){return;}
1122  } // if(fFillControlHistogramsWithGlobalTrackInfo)
1123 
1124  // c) Fill control histograms:
1125  AliAODTrack *agtrack = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
1126  if(fFillControlHistogramsWithGlobalTrackInfo){agtrack = gtrack;}
1127  else {agtrack = atrack;}
1128  if(!agtrack){Fatal(sMethodName.Data(),"!agtrack");}
1129  Int_t nCharge = 1;
1130  Bool_t bPrimary = kTRUE;
1131  for(Int_t charge=0;charge<2;charge++) // 0 = +q, 1 = -q
1132  {
1133  if(1==charge){nCharge = -1;} // TBI this is just ugly
1134  for(Int_t ps=0;ps<2;ps++) // 0 = kPrimary, 1 = kFromDecayVtx
1135  {
1136  if(1==ps){bPrimary = kFALSE;}
1137  // c2) Pions:
1138  if(Pion(gtrack,nCharge,bPrimary))
1139  {
1140  fPtPIDHist[2][charge][ps]->Fill(agtrack->Pt());
1141  fPPIDHist[2][charge][ps][0]->Fill(agtrack->Px());
1142  fPPIDHist[2][charge][ps][1]->Fill(agtrack->Py());
1143  fPPIDHist[2][charge][ps][2]->Fill(agtrack->Pz());
1144  fMassPIDHist[2][charge][ps]->Fill(agtrack->M());
1145  fEtaPIDHist[2][charge][ps]->Fill(agtrack->Eta());
1146  fPhiPIDHist[2][charge][ps]->Fill(agtrack->Phi());
1147  }
1148  // c3) Kaons:
1149  if(Kaon(gtrack,nCharge,bPrimary))
1150  {
1151  fPtPIDHist[3][charge][ps]->Fill(agtrack->Pt());
1152  fPPIDHist[3][charge][ps][0]->Fill(agtrack->Px());
1153  fPPIDHist[3][charge][ps][1]->Fill(agtrack->Py());
1154  fPPIDHist[3][charge][ps][2]->Fill(agtrack->Pz());
1155  fMassPIDHist[3][charge][ps]->Fill(agtrack->M());
1156  fEtaPIDHist[3][charge][ps]->Fill(agtrack->Eta());
1157  fPhiPIDHist[3][charge][ps]->Fill(agtrack->Phi());
1158  }
1159  // c4) Protons:
1160  if(Proton(gtrack,nCharge,bPrimary))
1161  {
1162  fPtPIDHist[4][charge][ps]->Fill(agtrack->Pt());
1163  fPPIDHist[4][charge][ps][0]->Fill(agtrack->Px());
1164  fPPIDHist[4][charge][ps][1]->Fill(agtrack->Py());
1165  fPPIDHist[4][charge][ps][2]->Fill(agtrack->Pz());
1166  fMassPIDHist[4][charge][ps]->Fill(agtrack->M());
1167  fEtaPIDHist[4][charge][ps]->Fill(agtrack->Eta());
1168  fPhiPIDHist[4][charge][ps]->Fill(agtrack->Phi());
1169  }
1170  } // for(Int_t ps=0;ps<2;ps++) // 0 = kPrimary, 1 = kFromDecayVtx
1171  } // for(Int_t charge=0;charge<2;charge++) // 0 = +q, 1 = -q
1172 
1173 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsIdentifiedParticles(AliAODTrack *atrack, AliAODTrack *gtrack)
1174 
1175 //================================================================================================================
1176 
1178 {
1179  // Fill control histograms for identified Monte Carlo particles.
1180 
1181  // a) Insanity checks;
1182  // b) Check cut selection criteria;
1183  // c) Fill control histograms.
1184 
1185  // a) Insanity checks:
1186  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsIdentifiedParticles(AliAODMCParticle *amcparticle)";
1187  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
1188 
1189  // b) Check cut selection criteria:
1190  if(!PassesCommonTrackCuts(amcparticle)){return;} // TBI have a look at implementation
1191 
1192  // c) Fill control histograms:
1193  Int_t index = -44;
1194  // c2) Pions:
1195  if(211==TMath::Abs(amcparticle->GetPdgCode()))
1196  {
1197  index = 2;
1198  } // if(211==TMath::Abs(amcparticle->GetPdgCode()))
1199  // c3) Kaons:
1200  else if(321==TMath::Abs(amcparticle->GetPdgCode()))
1201  {
1202  index = 3;
1203  } // if(321==TMath::Abs(amcparticle->GetPdgCode()))
1204  // c4) Protons:
1205  else if(2212==TMath::Abs(amcparticle->GetPdgCode()))
1206  {
1207  index = 4;
1208  } // if(2212==TMath::Abs(amcparticle->GetPdgCode()))
1209  if(-44 != index)
1210  {
1211  Int_t charge = ((amcparticle->GetPdgCode()>0. ? 0 : 1)); // Okay... TBI
1212  Int_t isPhysicalPrimary = ((amcparticle->IsPhysicalPrimary() ? 0 : 1)); // Okay... TBI
1213  fPtPIDHist[index][charge][isPhysicalPrimary]->Fill(amcparticle->Pt());
1214  fPPIDHist[index][charge][isPhysicalPrimary][0]->Fill(amcparticle->Px());
1215  fPPIDHist[index][charge][isPhysicalPrimary][1]->Fill(amcparticle->Py());
1216  fPPIDHist[index][charge][isPhysicalPrimary][2]->Fill(amcparticle->Pz());
1217  fMassPIDHist[index][charge][isPhysicalPrimary]->Fill(amcparticle->M()); // in this context, clearly this is just a control histogram
1218  fEtaPIDHist[index][charge][isPhysicalPrimary]->Fill(amcparticle->Eta());
1219  fPhiPIDHist[index][charge][isPhysicalPrimary]->Fill(amcparticle->Phi());
1220  } // if(-44 != index)
1221 
1222 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsIdentifiedParticles(AliAODMCParticle *amcparticle)
1223 
1224 //================================================================================================================
1225 
1227 {
1228  // Estimate background.
1229 
1230  // To do:
1231  // 1) Add support for MC and ESD.
1232 
1233  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::EstimateBackground(AliVEvent *ave)";
1234 
1235  // a) Determine Ali{MC,ESD,AOD}Event;
1236  // b) ...
1237 
1238  // a) Determine Ali{MC,ESD,AOD}Event:
1239  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
1240  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
1241  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
1242 
1243  if(aMC)
1244  {
1245  TString pattern = ".11.13.211.321.2212."; // TBI this can be done programatically TBI only for these particles I need to estimate background
1246 
1247  return; // TBI re-validate the lines below, within if statement, by following the analogy with AOD case
1248 
1249  if(!PassesMixedEventCuts(aMC)){return;} // TBI this is empty at the moment for MC
1250  Int_t nTracks = aMC->GetNumberOfTracks();
1251  if(!fMixedEvents0[0])
1252  {
1253  TClonesArray ca0("AliAODMCParticle"); // temporary holder, its entries will be copied to fMixedEvents0[0]
1254  Int_t ca0Counter = 0;
1255  for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
1256  {
1257  AliAODMCParticle *amcparticle = dynamic_cast<AliAODMCParticle*>(aMC->GetTrack(iTrack));
1258  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
1259  if(!PassesCommonTrackCuts(amcparticle)){continue;} // TBI re-think, see implemntation of this method
1260  if(!(pattern.Contains(Form(".%d.",TMath::Abs(amcparticle->GetPdgCode()))))){continue;}
1261  // Fill fMixedEvents0[0] with tracks which survived all checks:
1262  ca0[ca0Counter++] = amcparticle;
1263  } // for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
1264  fMixedEvents0[0] = (TClonesArray*)ca0.Clone();
1265  if(!fMixedEvents0[0]){Fatal(sMethodName.Data(),"!fMixedEvents0[0]");}
1266  }
1267  else if(!fMixedEvents0[1])
1268  {
1269  TClonesArray ca1("AliAODMCParticle"); // temporary holder, its entries will be copied to fMixedEvents0[1]
1270  Int_t ca1Counter = 0;
1271  for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
1272  {
1273  AliAODMCParticle *amcparticle = dynamic_cast<AliAODMCParticle*>(aMC->GetTrack(iTrack));
1274  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
1275  if(!PassesCommonTrackCuts(amcparticle)){continue;} // TBI re-think, see implemntation of this method
1276  if(!(pattern.Contains(Form(".%d.",TMath::Abs(amcparticle->GetPdgCode()))))){continue;}
1277  // Fill fMixedEvents0[0] with tracks which survived all checks:
1278  ca1[ca1Counter++] = amcparticle;
1279  } // for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
1280  fMixedEvents0[1] = (TClonesArray*)ca1.Clone();
1281  if(!fMixedEvents0[1]){Fatal(sMethodName.Data(),"!fMixedEvents0[1]");}
1282  }
1283  // Shall we do something?
1284  if(fMixedEvents0[0] && fMixedEvents0[1])
1285  {
1287  }
1288  } // if(aMC)
1289  else if(aESD)
1290  {
1291  // TBI
1292  } // else if(aESD)
1293  else if(aAOD)
1294  {
1295  if(!PassesMixedEventCuts(aAOD)){return;} // TBI returns always true at the moment...
1296 
1297  switch(fBackgroundOption)
1298  {
1299  // 0 : shifting
1300  // 1 : permutations + binning in vertex z-position
1301  // Key objects are fMixedEvents1[10][50]. First index indicates the vertex-z range, when the all entries in the second index are filled, then:
1302  // a) Make all possible permutations of 2-events, 3-events, etc.
1303  // b) For each permutation call Calculate2pBackground(...), Calculate3pBackground(...), etc.
1304 
1305  case 0: // shifting
1306  if(0==fMixedEvents0[0]->GetEntries())
1307  {
1308  if(fGlobalTracksAOD[1]) fGlobalTracksAOD[1]->Delete();
1309  GlobalTracksAOD(aAOD,1);
1310  if(0 == fGlobalTracksAOD[1]->GetSize()){return;}
1311  *fMixedEvents0[0] = *(aAOD->GetTracks());
1312  if(0==fMixedEvents0[0]->GetEntries()){return;}
1313  }
1314  else if(0==fMixedEvents0[1]->GetEntries())
1315  {
1316  if(fGlobalTracksAOD[2]) fGlobalTracksAOD[2]->Delete();
1317  GlobalTracksAOD(aAOD,2);
1318  if(0 == fGlobalTracksAOD[2]->GetSize()){return;}
1319  *fMixedEvents0[1] = *(aAOD->GetTracks());
1320  if(0==fMixedEvents0[1]->GetEntries()){return;}
1321  }
1322  else if(0==fMixedEvents0[2]->GetEntries())
1323  {
1324  if(fGlobalTracksAOD[3]) fGlobalTracksAOD[3]->Delete();
1325  GlobalTracksAOD(aAOD,3);
1326  if(0 == fGlobalTracksAOD[3]->GetSize()){return;}
1327  *fMixedEvents0[2] = *(aAOD->GetTracks());
1328  if(0==fMixedEvents0[2]->GetEntries()){return;}
1329  }
1330  // Shall we do something?
1331  if(0!=fMixedEvents0[0]->GetEntries() && 0!=fMixedEvents0[1]->GetEntries() && 0!=fMixedEvents0[2]->GetEntries())
1332  {
1335  //if(fEstimate4pBackground) Calculate4pBackground(fMixedEvents0[0],fMixedEvents0[1],fMixedEvents0[2],fMixedEvents0[3]);
1336  // Shift mixed events:
1337  // [1] -> [0]
1338  fMixedEvents0[0]->Delete();
1339  *fMixedEvents0[0] = *fMixedEvents0[1];
1341  // [2] -> [1]
1342  fMixedEvents0[1]->Delete();
1343  *fMixedEvents0[1] = *fMixedEvents0[2];
1345  // Clean [2]:
1346  fMixedEvents0[2]->Delete();
1347  fGlobalTracksAOD[3]->Delete();
1348  } // if(0!=fMixedEvents0[0]->GetEntries() && 0!=fMixedEvents0[1]->GetEntries() && 0!=fMixedEvents0[2]->GetEntries())
1349  break; // case 0: // shifting
1350 
1351  //=============================================================
1352 
1353  case 1: // permutations + binning in vertex z-position
1354  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
1355  if(!avtx){return;}
1356  Int_t nBinsVertexZrange = sizeof(fMixedEvents1)/sizeof(fMixedEvents1[0]); // extracts numbers of rows in TClonesArray *fMixedEvents1[10][5];
1357  if(nBinsVertexZrange<=0){Fatal(sMethodName.Data(),"nBinsVertexZrange<=0");}
1358  Int_t indexX = -44; // this is eventually the first index in TClonesArray *fMixedEvents1[10][5];
1359  Float_t flBinWidthVertexZrange = (fMaxVertexZ-fMinVertexZ)/(1.*nBinsVertexZrange);
1360  for(Int_t b=0;b<nBinsVertexZrange;b++)
1361  {
1362  //cout<<Form("%f - %f : %f",fMinVertexZ+1.*b*flBinWidthVertexZrange,fMinVertexZ+1.*(b+1)*flBinWidthVertexZrange,avtx->GetZ())<<endl;
1363  if(avtx->GetZ()>=fMinVertexZ+1.*b*flBinWidthVertexZrange && avtx->GetZ()<fMinVertexZ+1.*(b+1)*flBinWidthVertexZrange)
1364  {
1365  indexX = b;
1366  break;
1367  }
1368  } // for(Int_t b=0;b<nBinsVertexZrange;b++)
1369  if(-44==indexX){Fatal(sMethodName.Data(),"-44==indexX");}
1370  Int_t indexY = -44; // this is eventually the second index in TClonesArray *fMixedEvents1[10][5];
1371  for(Int_t bufferNo=0;bufferNo<fMaxBufferSize1;bufferNo++) // determining the 1st empty buffer
1372  {
1373  if(0==fMixedEvents1[indexX][bufferNo]->GetEntries()) // empty buffer
1374  {
1375  indexY = bufferNo;
1376  break;
1377  }
1378  } // for(Int_t bufferNo=0;bufferNo<5;bufferNo++)
1379  // Save in the buffer this event:
1380  *fMixedEvents1[indexX][indexY] = *(aAOD->GetTracks());
1381  GlobalTracksAOD(aAOD,indexX,indexY);
1382  // If the buffer is full, calculate background:
1383  if(fMaxBufferSize1-1==indexY)
1384  {
1385  cout<<"Flushing the buffer for background..."<<endl;
1386  cout<<Form("fMaxBufferSize1 = %d",fMaxBufferSize1)<<endl;
1387  cout<<Form("fMixedEvents1[indexX][indexY] = fMixedEvents1[%d][%d]",indexX,indexY)<<endl;
1388  // Shall we do something for 2p background?
1389  if(fEstimate2pBackground) // this buffer is full TBI hardwired 4. Making all possible combinations of 2 events
1390  {
1391  for(Int_t bufferNo1=0;bufferNo1<fMaxBufferSize1;bufferNo1++) // 1st event in the mixed events
1392  {
1393  for(Int_t bufferNo2=bufferNo1+1;bufferNo2<fMaxBufferSize1;bufferNo2++) // 2nd event in the mixed events
1394  {
1395  Calculate2pBackground(fMixedEvents1[indexX][bufferNo1],fMixedEvents1[indexX][bufferNo2],fGlobalTracksAOD1[indexX][bufferNo1],fGlobalTracksAOD1[indexX][bufferNo2]);
1396  } // for(Int_t bufferNo2=bufferNo1+1;bufferNo2<5;bufferNo2++) // 2nd event in the mixed events
1397  } // for(Int_t bufferNo1=0;bufferNo1<5;bufferNo1++) // first event in the mixed events
1398  } // if(fEstimate2pBackground)
1399 
1400  // Shall we do something for 3p background?
1401  if(fEstimate3pBackground) // this buffer is full TBI hardwired 4. Making all possible combinations of 3 events
1402  {
1403  for(Int_t bufferNo1=0;bufferNo1<fMaxBufferSize1;bufferNo1++) // 1st event in the mixed events
1404  {
1405  for(Int_t bufferNo2=0;bufferNo2<fMaxBufferSize1;bufferNo2++) // 2nd event in the mixed events
1406  {
1407  if(bufferNo2<=bufferNo1){continue;}
1408  for(Int_t bufferNo3=0;bufferNo3<fMaxBufferSize1;bufferNo3++) // 3rd event in the mixed events
1409  {
1410  if(bufferNo2<=bufferNo1 || bufferNo3<=bufferNo1 || bufferNo3<=bufferNo2){continue;}
1411  //cout<<Form("%d %d %d",bufferNo1,bufferNo2,bufferNo3)<<endl;
1412  Calculate3pBackground(fMixedEvents1[indexX][bufferNo1],fMixedEvents1[indexX][bufferNo2],fMixedEvents1[indexX][bufferNo3],fGlobalTracksAOD1[indexX][bufferNo1],fGlobalTracksAOD1[indexX][bufferNo2],fGlobalTracksAOD1[indexX][bufferNo3]);
1413  } // for(Int_t bufferNo3=0;bufferNo3<fMaxBufferSize1;bufferNo3++) // 3rd event in the mixed events
1414  } // for(Int_t bufferNo2=bufferNo1+1;bufferNo2<5;bufferNo2++) // 2nd event in the mixed events
1415  } // for(Int_t bufferNo1=0;bufferNo1<5;bufferNo1++) // first event in the mixed events
1416  } // if(fEstimate3pBackground)
1417 
1418  // Shall we do something for 4p background?
1419  if(fEstimate4pBackground) // this buffer is full TBI hardwired 4. Making all possible combinations of 4 events
1420  {
1421  // ...
1422  } // if(fEstimate4pBackground)
1423 
1424  // Clean the buffer:
1425  for(Int_t bufferNo=0;bufferNo<fMaxBufferSize1;bufferNo++)
1426  {
1427  fMixedEvents1[indexX][bufferNo]->Delete();
1428  fGlobalTracksAOD1[indexX][bufferNo]->Delete();
1429  } // for(Int_t bufferNo=0;bufferNo<5;bufferNo++)
1430 
1431  } // if(fMaxBufferSize1-1==indexY)
1432 
1433  break; // permutations + binning in vertex z-position
1434 
1435  //=============================================================
1436 
1437 /*
1438  default:
1439  cout<<Form("And the fatal 'fBackgroundOption' value is... %d. Congratulations!!",fBackgroundOption)<<endl;
1440  Fatal(sMethodName.Data(),"switch(fBackgroundOption)");
1441 */
1442 
1443  } // switch(fBackgroundOption)
1444 
1445  } // else if(aAOD)
1446 
1447 } // void AliAnalysisTaskMultiparticleFemtoscopy::EstimateBackground(AliVEvent *ave)
1448 
1449 //================================================================================================================
1450 
1452 {
1453  // Estimate background for TEST, using only 'shifting' at the moment.
1454 
1455  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::EstimateBackgroundTEST(AliVEvent *ave)";
1456 
1457  // a) Determine Ali{MC,ESD,AOD}Event;
1458  // b) ...
1459 
1460  // a) Determine Ali{MC,ESD,AOD}Event:
1461  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
1462  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
1463  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
1464 
1465  if(aMC)
1466  {
1467  // TBI
1468  } // if(aMC)
1469  else if(aESD)
1470  {
1471  // TBI
1472  } // else if(aESD)
1473  else if(aAOD)
1474  {
1475  if(!PassesMixedEventCuts(aAOD)){return;} // TBI returns always true at the moment...
1476 
1477  if(0==fMixedEventsTEST[0]->GetEntries())
1478  {
1479  if(fGlobalTracksAODTEST[0]) fGlobalTracksAODTEST[0]->Delete();
1480  GlobalTracksAODTEST(aAOD,0);
1481  if(0 == fGlobalTracksAODTEST[0]->GetSize()){return;}
1482  *fMixedEventsTEST[0] = *(aAOD->GetTracks());
1483  if(0==fMixedEventsTEST[0]->GetEntries()){return;}
1484  }
1485  else if(0==fMixedEventsTEST[1]->GetEntries())
1486  {
1487  if(fGlobalTracksAODTEST[1]) fGlobalTracksAODTEST[1]->Delete();
1488  GlobalTracksAODTEST(aAOD,1);
1489  if(0 == fGlobalTracksAODTEST[1]->GetSize()){return;}
1490  *fMixedEventsTEST[1] = *(aAOD->GetTracks());
1491  if(0==fMixedEventsTEST[1]->GetEntries()){return;}
1492  }
1493  else if(0==fMixedEventsTEST[2]->GetEntries())
1494  {
1495  if(fGlobalTracksAODTEST[2]) fGlobalTracksAODTEST[2]->Delete();
1496  GlobalTracksAODTEST(aAOD,2);
1497  if(0 == fGlobalTracksAODTEST[2]->GetSize()){return;}
1498  *fMixedEventsTEST[2] = *(aAOD->GetTracks());
1499  if(0==fMixedEventsTEST[2]->GetEntries()){return;}
1500  }
1501  // Shall we do something?
1502  if(0!=fMixedEventsTEST[0]->GetEntries() && 0!=fMixedEventsTEST[1]->GetEntries() && 0!=fMixedEventsTEST[2]->GetEntries())
1503  {
1504 
1507 
1508  //if(fEstimate4pBackground) Calculate4pBackground(fMixedEventsTEST[0],fMixedEventsTEST[1],fMixedEventsTEST[2],fMixedEventsTEST[3]);
1509  // Shift mixed events:
1510  // [1] -> [0]
1511  fMixedEventsTEST[0]->Delete();
1514  // [2] -> [1]
1515  fMixedEventsTEST[1]->Delete();
1518  // Clean [2]:
1519  fMixedEventsTEST[2]->Delete();
1520  fGlobalTracksAODTEST[2]->Delete();
1521  } // if(0!=fMixedEventsTEST[0]->GetEntries() && 0!=fMixedEventsTEST[1]->GetEntries() && 0!=fMixedEventsTEST[2]->GetEntries())
1522 
1523  } // else if(aAOD)
1524 
1525 } // void AliAnalysisTaskMultiparticleFemtoscopy::EstimateBackgroundTEST(AliVEvent *ave)
1526 
1527 //================================================================================================================
1528 
1530 {
1531  // Estimate background for TEST, using only 'shifting' at the moment.
1532 
1533  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::DoHybridApproach(AliVEvent *ave)";
1534 
1535  // a) Determine Ali{MC,ESD,AOD}Event;
1536  // b) ...
1537 
1538  // a) Determine Ali{MC,ESD,AOD}Event:
1539  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
1540  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
1541  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
1542 
1543  if(aMC)
1544  {
1545  // TBI
1546  } // if(aMC)
1547  else if(aESD)
1548  {
1549  // TBI
1550  } // else if(aESD)
1551  else if(aAOD)
1552  {
1553  if(!PassesMixedEventCuts(aAOD)){return;} // TBI returns always true at the moment...
1554 
1555  // a) Calculate 1st term:
1556  HybridApproach1stTerm(aAOD); // calculate N3(p1,p2,p3)
1557 
1558  // b) Mixed events buffers:
1559  if(0==fMixedEventsHA[0]->GetEntries())
1560  {
1561  if(fGlobalTracksAODHA[0]) fGlobalTracksAODHA[0]->Delete();
1562  GlobalTracksAODHA(aAOD,0);
1563  if(0 == fGlobalTracksAODHA[0]->GetSize()){return;}
1564  *fMixedEventsHA[0] = *(aAOD->GetTracks());
1565  if(0==fMixedEventsHA[0]->GetEntries()){return;}
1566  }
1567  else if(0==fMixedEventsHA[1]->GetEntries())
1568  {
1569  if(fGlobalTracksAODHA[1]) fGlobalTracksAODHA[1]->Delete();
1570  GlobalTracksAODHA(aAOD,1);
1571  if(0 == fGlobalTracksAODHA[1]->GetSize()){return;}
1572  *fMixedEventsHA[1] = *(aAOD->GetTracks());
1573  if(0==fMixedEventsHA[1]->GetEntries()){return;}
1574  }
1575  else if(0==fMixedEventsHA[2]->GetEntries())
1576  {
1577  if(fGlobalTracksAODHA[2]) fGlobalTracksAODHA[2]->Delete();
1578  GlobalTracksAODHA(aAOD,2);
1579  if(0 == fGlobalTracksAODHA[2]->GetSize()){return;}
1580  *fMixedEventsHA[2] = *(aAOD->GetTracks());
1581  if(0==fMixedEventsHA[2]->GetEntries()){return;}
1582  }
1583 
1584  // c) If all three buffers for mixed events are full, do something for the remaining terms:
1585  if(0!=fMixedEventsHA[0]->GetEntries() && 0!=fMixedEventsHA[1]->GetEntries() && 0!=fMixedEventsHA[2]->GetEntries())
1586  {
1587  // d) Cross-terms, e.g. N_{2}^A(p1,p2)N^B_{1}(p3). TBI Since I am focused in this exercise only on 3 identical pions, because of symmetry I need to calculate only one cross-term
1588  // Remark: By design, the current event can go only to the last buffer ([2]), so it's safe to mix it both with [0] and [1], to gain statistics
1589  HybridApproach2ndTerm(aAOD,fMixedEventsHA[0],fGlobalTracksAODHA[0]); // calculate N_{2}^A(p1,p2)N^B_{1}(p3), where B is [0]
1590  HybridApproach2ndTerm(aAOD,fMixedEventsHA[1],fGlobalTracksAODHA[1]); // calculate N_{2}^A(p1,p2)N^B_{1}(p3), where B is [1]
1591 
1592  // e) Genuine mixed-event terms:
1593  HybridApproach5thTerm(fMixedEventsHA[0],fMixedEventsHA[1],fMixedEventsHA[2],fGlobalTracksAODHA[0],fGlobalTracksAODHA[1],fGlobalTracksAODHA[2]); // calculate N^A_{1}(p1)N^B_{1}(p2)N^C_{1}(p3)
1594 
1595  // Shift mixed events:
1596  // [1] -> [0]
1597  fMixedEventsHA[0]->Delete();
1598  *fMixedEventsHA[0] = *fMixedEventsHA[1];
1600  // [2] -> [1]
1601  fMixedEventsHA[1]->Delete();
1602  *fMixedEventsHA[1] = *fMixedEventsHA[2];
1604  // Clean [2]:
1605  fMixedEventsHA[2]->Delete();
1606  fGlobalTracksAODHA[2]->Delete();
1607 
1608  } // if(0!=fMixedEventsHA[0]->GetEntries() && 0!=fMixedEventsHA[1]->GetEntries() && 0!=fMixedEventsHA[2]->GetEntries())
1609 
1610  } // else if(aAOD)
1611 
1612 } // void AliAnalysisTaskMultiparticleFemtoscopy::DoHybridApproach(AliVEvent *ave)
1613 
1614 //================================================================================================================
1615 
1617 {
1618  // Do all debugging in this function.
1619 
1620  // TBI: add support for fProcessBothKineAndReco
1621 
1622  // a) Determine Ali{MC,ESD,AOD}Event;
1623  // b) Wait for specified event.
1624 
1625  // a) Determine Ali{MC,ESD,AOD}Event:
1626  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
1627  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
1628  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
1629 
1630  // b) Wait for specified event:
1632  {
1633  if(aMC)
1634  {
1635  // TBI
1636  } // if(aMC)
1637  else if(aESD)
1638  {
1639  // TBI
1640  } // else if(aESD)
1641  else if(aAOD)
1642  {
1643  cout<<Form("aAOD->GetRunNumber() = %d",aAOD->GetRunNumber())<<endl;
1644  cout<<Form("aAOD->GetBunchCrossNumber() = %d",aAOD->GetBunchCrossNumber())<<endl;
1645  cout<<Form("aAOD->GetOrbitNumber() = %d",aAOD->GetOrbitNumber())<<endl;
1646  cout<<Form("aAOD->GetPeriodNumber() = %d",aAOD->GetPeriodNumber())<<endl;
1647  if(!SpecifiedEvent(aAOD->GetRunNumber(),aAOD->GetBunchCrossNumber(),aAOD->GetOrbitNumber(),aAOD->GetPeriodNumber())){return;}
1648  } // else if(aAOD)
1649  } // if(fWaitForSpecifiedEvent)
1650 
1651 } // void AliAnalysisTaskMultiparticleFemtoscopy::DoSomeDebugging(AliVEvent *ave)
1652 
1653 //================================================================================================================
1654 
1656 {
1657  // Determine the current event number.
1658 
1659  Int_t currentEventNumber = -44;
1660 
1661  if(fAnalysisType->EqualTo("MC"))
1662  {
1663  currentEventNumber = (Int_t)fGetNumberOfTracksHist->GetEntries(); // TBI there is an issue clearly if fFillControlHistogramsEvent is not enabled. Yes, this is really shaky...
1664  }
1665 
1666  else if(fAnalysisType->EqualTo("ESD"))
1667  {
1668  // TBI
1669  }
1670 
1671  else if(fAnalysisType->Contains("AOD")) // TBI: Okay...?
1672  {
1673  currentEventNumber = (Int_t)fGetNumberOfTracksHist->GetEntries(); // TBI there is an issue clearly if fFillControlHistogramsEvent is not enabled. Yes, this is really shaky...
1674  }
1675 
1676  return currentEventNumber;
1677 
1678 } // Int_t AliAnalysisTaskMultiparticleFemtoscopy::CurrentEventNumber()
1679 
1680 //================================================================================================================
1681 
1683 {
1684  // Per request, do some online monitoring.
1685 
1686  // a) Update regularly the output file;
1687  // b) Bail out after specified number of events.
1688 
1689  // a) Update regularly the output file:
1690  if(fUpdateOutputFile)
1691  {
1692  Int_t currentEventNo = this->CurrentEventNumber(); // TBI not supported yet for MC and ESD
1693  if(0 == currentEventNo % fUpdateFrequency)
1694  {
1695  cout<<Form("nEvts: %d",currentEventNo)<<endl;
1696  TFile *f = new TFile(fUpdateWhichOutputFile->Data(),"recreate");
1697  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
1698  f->Close();
1699  }
1700  } // if(fUpdateOutputFile)
1701 
1702  // b) Bail out after specified number of events:
1703  if(fMaxNumberOfEvents > 0)
1704  {
1705  Int_t currentEventNo = this->CurrentEventNumber(); // TBI not supported yet for MC and ESD
1706  if(fMaxNumberOfEvents == currentEventNo)
1707  {
1708  cout<<Form("\nPer request, bailing out after %d events in the file %s .\n",fMaxNumberOfEvents,fUpdateWhichOutputFile->Data())<<endl;
1709  TFile *f = new TFile(fUpdateWhichOutputFile->Data(),"recreate");
1710  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
1711  f->Close();
1712  exit(0);
1713  }
1714  } // if(fMaxNumberOfEvents > 0)
1715 
1716 } // void AliAnalysisTaskMultiparticleFemtoscopy::OnlineMonitoring()
1717 
1718 //================================================================================================================
1719 
1721 {
1722  // Insanity checks for UserCreateOutputObjects() method.
1723 
1724  // a) Insanity checks for global track cuts;
1725  // b) Disabled temporarily 4-p stuff. TBI a) fMixedEvents0[3] needs to be upgraded to fMixedEvents0[4], and b) EstimateBackground(AliVEvent *ave) needs to be re-implemented
1726  // c) Insanity checks for background.
1727 
1728  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::InsanityChecksUserCreateOutputObjects()";
1729 
1730  // a) Insanity checks for global track cuts:
1731  Int_t returnValueICFGTC = this->InsanityChecksForGlobalTrackCuts();
1732  if(0!=returnValueICFGTC)
1733  {
1734  cout<<Form("\n\nSomething is fundamentally wrong with global track cuts!!!!")<<endl;
1735  cout<<Form("InsanityChecksForGlobalTrackCuts() returns %d, check its implementation for an explanation of return values.\n\n",returnValueICFGTC)<<endl;
1736  Fatal(sMethodName.Data(),"if(0!=returnValueICFGTC)");
1737  } // if(0!=returnValueICFGTC)
1738 
1739  // b) Disabled temporarily 4-p stuff. TBI
1741  {
1742  cout<<Form("\n\nCalculation of 4-p stuff is not validated yet!!!! \n\n")<<endl; // TBI
1743  Fatal(sMethodName.Data(),"fFill4pCorrelationFunctions || fEstimate4pBackground");
1744  }
1745 
1746  // c) Insanity checks for background:
1747  if(!(fBackgroundOption == 0 || fBackgroundOption == 1))
1748  {
1749  cout<<Form("\n\nThis option for calculating background is not supported yet!!!! \n\n")<<endl; // TBI
1750  Fatal(sMethodName.Data(),"!(fBackgroundOption == 0 || fBackgroundOption == 1)");
1751  } // if(fBackgroundOption !=0 || fBackgroundOption != 1)
1752  if(fMaxBufferSize1>=50)
1753  {
1754  cout<<Form("\n\nfMaxBufferSize1 can be at maximum 50 !!!! \n\n")<<endl; // TBI
1755  Fatal(sMethodName.Data(),"fMaxBufferSize1>=50");
1756  }
1757 
1758 } // void AliAnalysisTaskMultiparticleFemtoscopy::InsanityChecksUserCreateOutputObjects()
1759 
1760 //================================================================================================================
1761 
1763 {
1764  // Insanity checks for UserExec() method.
1765 
1766  // a) Insanity checks specific for all analyses types;
1767  // b) Insanity checks specific only for MC analyses;
1768  // c) Insanity checks specific only for ESD analyses;
1769  // d) Insanity checks specific only for AOD analyses.
1770 
1771  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::InsanityChecksUserExec()";
1772 
1773  // a) Insanity checks specific for all analyses types: // TBI is this really specific for all analyses types
1774  if(fBailOutAfterQA && !fFillQA){Fatal(sMethodName.Data(),"fBailOutAfterQA && !fFillQA");}
1775  if(fOnlineMonitoring && !fFillControlHistogramsEvent){Fatal(sMethodName.Data(),"fOnlineMonitoring && !fFillControlHistogramsEvent.\n\nAt the moment, fOnlineMonitoring can be used only if fFillControlHistogramsEvent is enabled.");}
1776  if(fFillControlHistogramsNonIdentifiedParticlesFTSF && !(fFilterBitFTSF==128)){Fatal(sMethodName.Data(),"FillControlHistogramsNonIdentifiedParticlesFTSF && !(fFilterBitFTSF==128)\n\nAt the moment, fFillControlHistogramsNonIdentifiedParticlesFTSF is validated only for TPC-only tracks.");}
1777 
1778 
1779 
1780  return; // TBI re-think and re-validate the rest
1781 
1782 
1783 
1784  // a2) TBI:
1785  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1786  {
1787  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
1788  {
1789  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
1790  {
1791  if(!fPIDCA[pid][pa][ps]) {Fatal(sMethodName.Data(),"!fPIDCA[pid][pa][ps]");}
1792  if(0 != fPIDCA[pid][pa][ps]->GetEntriesFast()){Fatal(sMethodName.Data(),"0 != fPIDCA[pid][pa][ps]->GetEntriesFast()"); }
1793  }
1794  }
1795  }
1796  for(Int_t pid=0;pid<1;pid++) // [0=Lambda,1=...] // TBI is this really specific for all analyses types
1797  {
1798  if(!fPIDV0sCA[pid]){Fatal(sMethodName.Data(),"!fPIDV0sCA[pid]");}
1799  if(0 != fPIDV0sCA[pid]->GetEntriesFast()){Fatal(sMethodName.Data(),"0 != fPIDV0sCA[pid]->GetEntriesFast()");}
1800  }
1801 
1802  // b) Insanity checks specific only for MC analyses:
1803  if(fAnalysisType->EqualTo("MC"))
1804  {
1805  // TBI
1806  }
1807 
1808  // c) Insanity checks specific only for ESD analyses:
1809  if(fAnalysisType->EqualTo("ESD")) // TBI 'else if' ?
1810  {
1811  // TBI
1812  }
1813 
1814  // d) Insanity checks specific only for AOD analyses:
1815  if(fAnalysisType->Contains("AOD")) // TBI 'else if' ?
1816  {
1817  if(!fGlobalTracksAOD[0]){Fatal(sMethodName.Data(),"!fGlobalTracksAOD[0]");}
1818  if(0 != fGlobalTracksAOD[0]->GetSize()){Fatal(sMethodName.Data(),"0 != fGlobalTracksAOD[0]->GetSize()");}
1819  //if(fProcessBothKineAndReco && fEstimateBackground){Fatal(sMethodName.Data(),"fProcessBothKineAndReco && fEstimateBackground");} TBI re-think and re-enable
1820  // Remark: Note that the similar check is not needed for instance for fGlobalTracksAOD[1] and fGlobalTracksAOD[2],
1821  // which are used to estimate background, since some events can be stored in a buffer, until a suitable co-event is found,
1822  // to calculate background.
1823  } // if(fAnalysisType->Contains("AOD")
1824 
1825 } // void AliAnalysisTaskMultiparticleFemtoscopy::InsanityChecksUserExec()
1826 
1827 //=======================================================================================================================
1828 
1830 {
1831  // Monte Carlo analysis is performed in this method.
1832 
1833  // a) Debugging;
1834  // b) Common event selection criteria;
1835  // c) Fill control histogram for global event observables;
1836  // d) Fill control histogram for particles;
1837  // e) Calculate correlation functions;
1838  // f) Calculate background;
1839  // g) V0s.
1840 
1841  // a) Debugging:
1842  if(fDoSomeDebugging){this->DoSomeDebugging(aMC);}
1843 
1844  // b) Common event selection criteria:
1845  if(!this->PassesCommonEventCuts(aMC)){return;}
1846 
1847  // c) Fill control histogram for global event observables:
1849 
1850  // d) Fill control histograms for particles:
1852 
1853  // e) Calculate correlation functions:
1855 
1856  // f) Calculate background:
1858 
1859  return;
1860 
1861  // g) V0s:
1862  // TBI
1863 
1864 } // void AliAnalysisTaskMultiparticleFemtoscopy::MC(AliMCEvent *aMC)
1865 
1866 //=======================================================================================================================
1867 
1869 {
1870  // ESD analysis is performed in this method.
1871 
1872  // ...
1873  if(aESD)
1874  {
1875  aESD = NULL; // TBI eliminating warnings temporarily
1876  }
1877 
1878 } // void AliAnalysisTaskMultiparticleFemtoscopy::ESD(AliESDEvent *aESD)
1879 
1880 //=======================================================================================================================
1881 
1883 {
1884  // AOD analysis is performed in this method.
1885 
1886  // a) Debugging;
1887  // b) Common event selection criteria;
1888  // c) Filter out "normal global" tracks for default analysis and cut on their number;
1889  // d) Fill control histogram for global event observables;
1890  // e) Fill control histogram for particles;
1891  // f) Calculate correlation functions;
1892  // g) Calculate background;
1893  // h) Calculate test correlation functions;
1894  // i) Calculate test background;
1895  // j) 'Hybrid appriach';
1896  // *) V0s.
1897 
1898  // a) Debugging:
1899  if(fDoSomeDebugging){this->DoSomeDebugging(aAOD);}
1900 
1901  // b) Common event selection criteria:
1902  if(!this->PassesCommonEventCuts(aAOD)){return;}
1903 
1904  // c) Filter out "normal global" tracks for default analysis and cut on their number:
1905  // Remark 1: 'TPC-only' tracks and 'global constrained to vertex' come with negative ID, and are therefore not stored in fGlobalTracksAOD[*]
1906  // Remark 2: Default analysis means that it is not used for mixed events
1907  this->GlobalTracksAOD(aAOD,0); // [0] stands for default analysis
1908  if(0 == fGlobalTracksAOD[0]->GetSize()) return; // yes, go to next event
1910  {
1911  if(fGlobalTracksAOD[0]->GetSize() < fMinNumberOfGlobalTracks) return;
1912  if(fGlobalTracksAOD[0]->GetSize() > fMaxNumberOfGlobalTracks) return;
1913  }
1914 
1915  // d) Fill control histogram for global event observables:
1917 
1918  // e) Fill control histograms for particles:
1920 
1921  // f) Calculate correlation functions:
1923 
1924  // g) Calculate background:
1926 
1927  // h) Calculate test correlation functions:
1928  Bool_t bFillCorrelationFunctionsTEST = kFALSE;
1929  for(Int_t t=0;t<10;t++)
1930  {
1931  bFillCorrelationFunctionsTEST = bFillCorrelationFunctionsTEST || fFillCorrelationFunctionsTEST[t];
1932  }
1933  if(bFillCorrelationFunctionsTEST){this->CalculateCorrelationFunctionsTEST(aAOD);}
1934 
1935  // i) Calculate test background:
1936  Bool_t bFillBackgroundTEST = kFALSE;
1937  for(Int_t t=0;t<10;t++)
1938  {
1939  bFillBackgroundTEST = bFillBackgroundTEST || fFillBackgroundTEST[t];
1940  }
1941  if(bFillBackgroundTEST){this->EstimateBackgroundTEST(aAOD);}
1942 
1943  // j) 'Hybrid appriach':
1944  if(fDoHybridApproach){this->DoHybridApproach(aAOD);}
1945 
1946  return;
1947 
1948 
1949  // *) V0s:
1950  V0s(aAOD); // TBI implement flag to enable/disable this call
1951 
1952 } // void AliAnalysisTaskMultiparticleFemtoscopy::AOD(AliAODEvent *aAOD)
1953 
1954 //=======================================================================================================================
1955 
1957 {
1958  // Analysis with V0s.
1959 
1960  // a) Determine Ali{MC,ESD,AOD}Event;
1961 
1962  // a) Determine Ali{MC,ESD,AOD}Event:
1963  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
1964  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
1965  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
1966 
1967  if(aMC)
1968  {
1969  // TBI
1970  }
1971  else if(aESD)
1972  {
1973  // TBI
1974  }
1975  else if(aAOD)
1976  {
1977  // AOD V0s:
1978  TClonesArray *caV0s = aAOD->GetV0s();
1979  Int_t index = 0;
1980  AliAODv0 *aAODv0 = NULL;
1981  Int_t nProngs = -44;
1982  while(caV0s->At(index))
1983  {
1984  aAODv0 = (AliAODv0*) caV0s->At(index++);
1985  if(!aAODv0) break;
1986 
1987 
1988  //AliAODv0 *temp = (AliAODv0*)fPIDV0sCA[0]->ConstructedAt(index-1);
1989  //temp = (AliAODv0*)aAODv0->Clone();
1990 
1991  //cout<<Form("indddex = %d, fPIDV0sCA[0]->GetEntries() = %d",index,fPIDV0sCA[0]->GetEntries())<<endl;
1992 
1993  // TBI...
1994  continue;
1995 
1996  nProngs = aAODv0->GetNProngs();
1997  fGetNProngsHist->Fill(nProngs);
1998  fMassK0ShortHist->Fill(aAODv0->MassK0Short());
1999  fMassLambdaHist->Fill(aAODv0->MassLambda());
2000  fMassAntiLambdaHist->Fill(aAODv0->MassAntiLambda());
2001  fOpenAngleV0Hist->Fill(aAODv0->OpenAngleV0());
2002  fRadiusV0Hist->Fill(aAODv0->RadiusV0());
2003  fDcaV0ToPrimVertexHist->Fill(aAODv0->DcaV0ToPrimVertex());
2004  fMomV0XHist->Fill(aAODv0->MomV0X());
2005  fMomV0YHist->Fill(aAODv0->MomV0Y());
2006  fMomV0ZHist->Fill(aAODv0->MomV0Z());
2007  fPtV0Hist->Fill(pow(aAODv0->Pt2V0(),0.5));
2008  fPseudoRapV0Hist->Fill(aAODv0->PseudoRapV0());
2009  fPAHist->Fill(aAODv0->Alpha(),aAODv0->PtArmV0());
2010  // Check sharing:
2011 
2012  fUniqueIDHistEBE->Fill(aAODv0->GetPosID());
2013  fUniqueIDHistEBE->Fill(aAODv0->GetNegID());
2014 
2015  /*
2016  //cout<<Form("V0PosID: %d , V0NegID: %d",aAODv0->GetPosID(),aAODv0->GetNegID())<<endl;
2017  Int_t trackPos = aAODv0->GetPosID()>=0 ? fGlobalTracksAOD[0]->GetValue(aAODv0->GetPosID()) : fGlobalTracksAOD[0]->GetValue(-(aAODv0->GetPosID()+1));
2018  Int_t trackNeg = aAODv0->GetNegID()>=0 ? fGlobalTracksAOD[0]->GetValue(aAODv0->GetNegID()) : fGlobalTracksAOD[0]->GetValue(-(aAODv0->GetNegID()+1));
2019  //cout<<Form("global : %d, global : %d", trackPos, trackNeg)<<endl;
2020  */
2021 
2022  if(-1 != fUniqueIDHistEBE->FindFirstBinAbove(1.44,1)) // TBI
2023  {
2024  cout<<Form("fUniqueIDHistEBE->FindFirstBinAbove(1.44) %d:",(Int_t)fUniqueIDHistEBE->FindFirstBinAbove(1.44,1))<<endl; //
2025  }
2026  } // while(caV0s->At(index))
2027 
2028 
2029  /*
2030  Int_t index2 = 0;
2031  while(caV0s->At(index2))
2032  {
2033  cout<<((AliAODv0*) caV0s->At(index2))->Pt()<<endl;
2034  cout<<((AliAODv0*) fPIDV0sCA[0]->At(index2))->Pt()<<endl;
2035  index2++;
2036  cout<<endl;
2037  }
2038  */
2039 
2040  // AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,globaltrack);
2041  // AliPIDResponse::EDetPidStatus statusTOF = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,globaltrack);
2042 
2043  }
2044 
2045 } // void AliAnalysisTaskMultiparticleFemtoscopy::V0s(AliVEvent *ave)
2046 
2047 //=======================================================================================================================
2048 
2050 {
2051  // Reset all event-by-event objects.
2052 
2053  // a) Reset event-by-event objects specific for all analyses types;
2054  // b) Reset event-by-event objects specific only for MC analyses;
2055  // c) Reset event-by-event objects specific only for ESD analyses;
2056  // d) Reset event-by-event objects specific only for AOD analyses.
2057 
2058  // a) Reset event-by-event objects specific for all analyses types:
2059  if(fUniqueIDHistEBE) fUniqueIDHistEBE->Reset();
2060  // TBI add comment
2061  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
2062  {
2063  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
2064  {
2065  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
2066  {
2067  if(fPIDCA[pid][pa][ps]) fPIDCA[pid][pa][ps]->Delete();
2068  }
2069  }
2070  }
2071  // TBI add comment
2072  for(Int_t pid=0;pid<1;pid++) // [0=Lambda,1=...]
2073  {
2074  if(fPIDV0sCA[pid]) fPIDV0sCA[pid]->Delete();
2075  }
2076 
2077  // b) Reset event-by-event objects specific only for MC analyses:
2078  if(fAnalysisType->EqualTo("MC"))
2079  {
2080  // TBI
2081  }
2082 
2083  // c) Reset event-by-event objects specific only for ESD analyses:
2084  if(fAnalysisType->EqualTo("ESD")) // TBI 'else if' ?
2085  {
2086  // TBI
2087  }
2088 
2089  // d) Reset event-by-event objects specific only for AOD analyses:
2090  if(fAnalysisType->Contains("AOD")) // TBI 'else if' ?
2091  {
2092  if(fGlobalTracksAOD[0]) fGlobalTracksAOD[0]->Delete();
2093  // Remark: Note that fGlobalTracksAOD[1] and fGlobalTracksAOD[2] used for event mixing do not have
2094  // to be reset e-b-e. Where then I reset them? TBI
2095  }
2096 
2097 } // void AliAnalysisTaskMultiparticleFemtoscopy::ResetEBEObjects()
2098 
2099 //=======================================================================================================================
2100 
2102 {
2103  // Is this primary/secondary (anti-)pion?
2104 
2105  // a) Insanity checks;
2106  // b) Trivial checks; // TBI
2107  // c) Track quality cuts;
2108  // d) PID;
2109  // e) PID MC (if available).
2110 
2111  // a) Insanity checks
2112  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Pion(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)";
2113  if(!(1 == charge || -1 == charge)){Fatal(sMethodName.Data(),"!(1 == charge || -1 == charge)");}
2114  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
2115 
2116  // b) Trivial checks:
2117  if(charge != (Int_t)gtrack->Charge()) return kFALSE;
2118  if(bPrimary && gtrack->GetType() != AliAODTrack::kPrimary)
2119  {
2120  return kFALSE;
2121  }
2122  else if(!bPrimary && gtrack->GetType() != AliAODTrack::kFromDecayVtx)
2123  {
2124  return kFALSE;
2125  }
2126 
2127  // c) Track quality cuts:
2128 /*
2129  if(bPrimary)
2130  {
2131  // TBI
2132  }
2133 */
2134 
2135  // d) PID:
2136  // For pT < 0.75 use only TPC
2137  if(gtrack->GetTPCmomentum() < 0.75) // TBI hardwired 0.75
2138  {
2139  AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,gtrack);
2140  if(!statusTPC) return kFALSE;
2141  // sigmas:
2142  Double_t dSigmaTPCPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kPion)));
2143  Double_t dSigmaTPCKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kKaon)));
2144  Double_t dSigmaTPCProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kProton)));
2145  // inclusive cuts and exclusive cuts:
2146  if(!(dSigmaTPCPion < fInclusiveSigmaCuts[2] && dSigmaTPCProton > fExclusiveSigmaCuts[2][4] && dSigmaTPCKaon > fExclusiveSigmaCuts[2][3])) return kFALSE;
2147  } // if(gtrack->GetTPCmomentum() <= 0.75) // TBI hardwired 0.75
2148  else if(gtrack->GetTPCmomentum() >= 0.75 )
2149  {
2150  // TBI use combined TPC and TOf
2151  return kFALSE; // TBI remove eventually
2152  }
2153 
2154  // e) PID MC (if available):
2156  {
2157  if(!fMC){Fatal(sMethodName.Data(),"!fMC");}
2158  AliAODMCParticle *mcParticle = dynamic_cast<AliAODMCParticle*>(fMC->GetTrack(TMath::Abs(gtrack->GetLabel())));
2159  if(!mcParticle){cout<<"WARNING: mcParticle is NULL in Pion(...)! gtrack = "<<gtrack<<endl; return kFALSE; } // TBI well, it remains undetermined, investigate further why in some very rare cases I cannot get this pointer. if I get too many warnings, that purity estimation will be affected
2160  if(charge < 0 && mcParticle->GetPdgCode() == -211) return kTRUE;
2161  else if(charge > 0 && mcParticle->GetPdgCode() == 211) return kTRUE;
2162  else return kFALSE;
2163  } // if(fProcessBothKineAndReco)
2164 
2165  return kTRUE;
2166 
2167 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Pion(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)
2168 
2169 
2170 //=======================================================================================================================
2171 
2173 {
2174  // Is this primary/secondary (anti-)kaon?
2175 
2176  // a) Insanity checks;
2177  // b) Trivial checks;
2178  // c) Track quality cuts;
2179  // d) PID;
2180  // e) PID MC (if available).
2181 
2182  // a) Insanity checks:
2183  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Kaon(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)";
2184  if(!(1 == charge || -1 == charge)){Fatal(sMethodName.Data(),"!(1 == charge || -1 == charge)");}
2185  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
2186 
2187  // b) Trivial checks:
2188  if(charge != (Int_t)gtrack->Charge()) return kFALSE;
2189  if(bPrimary && gtrack->GetType() != AliAODTrack::kPrimary)
2190  {
2191  return kFALSE;
2192  }
2193  else if(!bPrimary && gtrack->GetType() != AliAODTrack::kFromDecayVtx)
2194  {
2195  return kFALSE;
2196  }
2197 
2198  // c) Track quality cuts:
2199  /*
2200  if(bPrimary)
2201  {
2202  // TBI
2203  }
2204  */
2205 
2206  // d) PID:
2207  // For pT < 0.75 use only TPC
2208  if(gtrack->GetTPCmomentum() < 0.75) // TBI hardwired 0.75
2209  {
2210  AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,gtrack);
2211  if(!statusTPC) return kFALSE;
2212  // sigmas:
2213  Double_t dSigmaTPCKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kKaon)));
2214  Double_t dSigmaTPCProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kProton)));
2215  Double_t dSigmaTPCPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kPion)));
2216  // inclusive and exclusive cuts:
2217  if(!(dSigmaTPCKaon < fInclusiveSigmaCuts[3] && dSigmaTPCProton > fExclusiveSigmaCuts[3][4] && dSigmaTPCPion > fExclusiveSigmaCuts[3][2])) return kFALSE;
2218  } // if(gtrack->GetTPCmomentum() <= 0.75) // TBI hardwired 0.75
2219  else if(gtrack->GetTPCmomentum() >= 0.75 )
2220  {
2221  // TBI use combined TPC and TOf
2222  return kFALSE; // TBI remove eventually
2223  }
2224 
2225  // e) PID MC (if available):
2227  {
2228  if(!fMC){Fatal(sMethodName.Data(),"!fMC");}
2229  AliAODMCParticle *mcParticle = dynamic_cast<AliAODMCParticle*>(fMC->GetTrack(TMath::Abs(gtrack->GetLabel())));
2230  if(!mcParticle){cout<<"WARNING: mcParticle is NULL in Kaon(...)! gtrack = "<<gtrack<<endl; return kFALSE; } // TBI well, it remains undetermined, investigate further why in some very rare cases I cannot get this pointer. if I get too many warnings, that purity estimation will be affected
2231  if(charge < 0 && mcParticle->GetPdgCode() == -321) return kTRUE;
2232  else if(charge > 0 && mcParticle->GetPdgCode() == 321) return kTRUE;
2233  else return kFALSE;
2234  } // if(fProcessBothKineAndReco)
2235 
2236  return kTRUE;
2237 
2238 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Kaon(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)
2239 
2240 //=======================================================================================================================
2241 
2243 {
2244  // Is this primary/secondary (anti-)proton?
2245 
2246  // a) Insanity checks;
2247  // b) Trivial checks;
2248  // c) Track quality cuts;
2249  // d) PID;
2250  // e) PID MC (if available).
2251 
2252  // a) Insanity checks:
2253  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Proton(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)";
2254  if(!(1 == charge || -1 == charge)){Fatal(sMethodName.Data(),"!(1 == charge || -1 == charge)");}
2255  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
2256 
2257  // b) Trivial checks:
2258  if(charge != (Int_t)gtrack->Charge()) return kFALSE;
2259  if(bPrimary && gtrack->GetType() != AliAODTrack::kPrimary)
2260  {
2261  return kFALSE;
2262  }
2263  else if(!bPrimary && gtrack->GetType() != AliAODTrack::kFromDecayVtx)
2264  {
2265  return kFALSE;
2266  }
2267 
2268  // c) Track quality cuts:
2269 /*
2270  if(bPrimary)
2271  {
2272  // TBI
2273  }
2274 */
2275 
2276  // d) PID:
2277  // For pT < 0.75 use only TPC
2278  if(gtrack->GetTPCmomentum() < 0.75) // TBI hardwired 0.75
2279  {
2280  AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,gtrack);
2281  if(!statusTPC) return kFALSE;
2282  // sigmas:
2283  Double_t dSigmaTPCProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kProton)));
2284  Double_t dSigmaTPCPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kPion)));
2285  Double_t dSigmaTPCKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kKaon)));
2286  if(!(dSigmaTPCProton < fInclusiveSigmaCuts[4] && dSigmaTPCPion > fExclusiveSigmaCuts[4][2] && dSigmaTPCKaon > fExclusiveSigmaCuts[4][3])) return kFALSE;
2287  } // if(gtrack->GetTPCmomentum() <= 0.75) // TBI hardwired 0.75
2288  else if(gtrack->GetTPCmomentum() >= 0.75 )
2289  {
2290  // TBI use combined TPC and TOf
2291  return kFALSE; // TBI remove eventually
2292  }
2293 
2294  // e) PID MC (if available):
2296  {
2297  if(!fMC){Fatal(sMethodName.Data(),"!fMC");}
2298  AliAODMCParticle *mcParticle = dynamic_cast<AliAODMCParticle*>(fMC->GetTrack(TMath::Abs(gtrack->GetLabel())));
2299  if(!mcParticle){cout<<"WARNING: mcParticle is NULL in Proton(...)! gtrack = "<<gtrack<<endl; return kFALSE; } // TBI well, it remains undetermined, investigate further why in some very rare cases I cannot get this pointer. if I get too many warnings, that purity estimation will be affected
2300  if(charge < 0 && mcParticle->GetPdgCode() == -2212) return kTRUE;
2301  else if(charge > 0 && mcParticle->GetPdgCode() == 2212) return kTRUE;
2302  else return kFALSE;
2303  } // if(fProcessBothKineAndReco)
2304 
2305  return kTRUE;
2306 
2307 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Proton(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)
2308 
2309 //=======================================================================================================================
2310 
2312 {
2313  // Initialize arrays for all objects not classified yet.
2314 
2315  for(Int_t index=0;index<10;index++) // [0] is used in the default analysis, [1] and [2] for event mixing, etc.
2316  {
2317  fGlobalTracksAOD[index] = NULL;
2318  }
2319 
2320 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArrays()
2321 
2322 //=======================================================================================================================
2323 
2325 {
2326  // Initialize all arrays for control histograms.
2327 
2328  for(Int_t xyz=0;xyz<3;xyz++)
2329  {
2330  fVertexXYZ[xyz] = NULL;
2331  }
2332 
2333  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
2334  {
2335  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
2336  {
2337  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
2338  {
2339  fMassPIDHist[pid][pa][ps] = NULL;
2340  fPtPIDHist[pid][pa][ps] = NULL;
2341  for(Int_t xyz=0;xyz<3;xyz++) // kPrimary/kFromDecayVtx
2342  {
2343  fPPIDHist[pid][pa][ps][xyz] = NULL;
2344  }
2345  fEtaPIDHist[pid][pa][ps] = NULL;
2346  fPhiPIDHist[pid][pa][ps] = NULL;
2347  }
2348  }
2349  } // for(Int_t pid=0;pid<4;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
2350 
2351  // Inclusive sigma cuts:
2352  for(Int_t pf=0;pf<5;pf++) // PID function [0=Electron(...),1=Muon(...),2=Pion(...),3=Kaon(...),4=Proton(...)]
2353  {
2354  fInclusiveSigmaCuts[pf] = 0.;
2355  }
2356 
2357  // Exclusive sigma cuts:
2358  for(Int_t pf=0;pf<5;pf++) // PID function [0=Electron(...),1=Muon(...),2=Pion(...),3=Kaon(...),4=Proton(...)]
2359  {
2360  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
2361  {
2362  fExclusiveSigmaCuts[pf][pid] = 0.;
2363  }
2364  }
2365 
2366  // Normalization:
2367  fNormalizationInterval[0] = -0.44;
2368  fNormalizationInterval[1] = -0.44;
2369 
2370 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForControlHistograms()
2371 
2372 //=======================================================================================================================
2373 
2375 {
2376  // Initialize all arrays for e-b-e objects.
2377 
2378  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
2379  {
2380  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
2381  {
2382  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
2383  {
2384  fPIDCA[pid][pa][ps] = NULL;
2385  }
2386  }
2387  } // for(Int_t pid=0;pid<4;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
2388 
2389 
2390  for(Int_t pid=0;pid<1;pid++)
2391  {
2392  fPIDV0sCA[pid] = NULL;
2393  }
2394 
2395 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForEBEObjects()
2396 
2397 //=======================================================================================================================
2398 
2400 {
2401  // Initialize all arrays for correlation functions.
2402 
2403  // a) Initialize sublists;
2404  // b) Initialize 2p correlation functions;
2405  // c) Initialize 3p correlation functions;
2406  // d) Initialize 4p correlation functions.
2407 
2408  // a) Initialize sublists:
2409  for(Int_t cfs=0;cfs<3;cfs++)
2410  {
2411  fCorrelationFunctionsSublist[cfs] = NULL;
2412  }
2413 
2414  // b) Initialize 2p correlation functions:
2415  for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2416  {
2417  for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2418  {
2419  fCorrelationFunctions[pid1][pid2] = NULL;
2420  }
2421  }
2422 
2423  // c) Initialize 3p correlation functions:
2424  for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2425  {
2426  for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2427  {
2428  for(Int_t pid3=0;pid3<10;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2429  {
2430  f3pCorrelationFunctions[pid1][pid2][pid3] = NULL;
2431  }
2432  }
2433  }
2434 
2435  // c) Initialize 4p correlation functions:
2436  for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2437  {
2438  for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2439  {
2440  for(Int_t pid3=0;pid3<10;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2441  {
2442  for(Int_t pid4=0;pid4<10;pid4++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2443  {
2444  f4pCorrelationFunctions[pid1][pid2][pid3][pid4] = NULL;
2445  }
2446  }
2447  }
2448  }
2449 
2450 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForCorrelationFunctions()
2451 
2452 //=======================================================================================================================
2453 
2455 {
2456  // Initialize all arrays for background.
2457 
2458  for(Int_t bs=0;bs<3;bs++)
2459  {
2460  fBackgroundSublist[bs] = NULL; // lists to hold all background correlations, for 2p [0], 3p [1], 4p [2], etc., separately
2461  }
2462 
2463  for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2464  {
2465  for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2466  {
2467  f2pBackground[pid1][pid2] = NULL;
2468  }
2469  }
2470 
2471  for(Int_t pid1=0;pid1<10;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2472  {
2473  for(Int_t pid2=0;pid2<10;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2474  {
2475  for(Int_t pid3=0;pid3<10;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
2476  {
2477  f3pBackground[pid1][pid2][pid3] = NULL;
2478  }
2479  }
2480  }
2481 
2482  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
2483  {
2484  fMixedEvents0[me] = NULL;
2485  }
2486 
2487  for(Int_t vzr=0;vzr<10;vzr++) // vertez-z range
2488  {
2489  for(Int_t me=0;me<fMaxBufferSize1;me++) // mixed events
2490  {
2491  fMixedEvents1[vzr][me] = NULL;
2492  fGlobalTracksAOD1[vzr][me] = NULL;
2493  }
2494  }
2495 
2496 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForBackground()
2497 
2498 //=======================================================================================================================
2499 
2501 {
2502  // Initialize all arrays for buffers.
2503 
2504  for(Int_t am=0;am<2;am++) // analysis method [0=AOD||ESD,1=MC]
2505  {
2506  for(Int_t e=0;e<10;e++) // number of events to buffer, max = apparently 10
2507  {
2508  for(Int_t p=0;p<10000;p++) // charged particles, max = apparently 10000
2509  {
2510  fChargedParticlesCA[am][e][p] = NULL;
2511  // Example: fChargedParticlesCA[0][4][6] is 6th charged particles in 5th buffered event, in AOD||ESD analysis
2512  } // for(Int_t p=0;p<10000;p++) // charged particles, max = apparently 10000
2513  } // for(Int_t e=0;e<10;e++) // number of events to buffer, max = apparently 10
2514  } // for(Int_t am=0;am<2;am++) // analysis method [0=AOD||ESD,1=MC]
2515 
2516  for(Int_t e=0;e<10;e++) // number of events to buffer, max = apparently 10
2517  {
2518  fChargedParticlesEM[e] = NULL;
2519  } // for(Int_t e=0;e<10;e++) // number of events to buffer, max = apparently 10
2520 
2521 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForBuffers()
2522 
2523 //=======================================================================================================================
2524 
2526 {
2527  // Initialize all arrays for QA.
2528 
2529  for(Int_t ba=0;ba<2;ba++) // 0="before rain",1="after rain"
2530  {
2531  for(Int_t di=0;di<10;di++) // number of distinct distributions, maximum is 10 (at the moment, as it seems...)
2532  {
2533  for(Int_t ci=0;ci<10;ci++) // number of distinct track cuts, maximum is 10 (at the moment, as it seems...)
2534  {
2535  fQAParticleHist[ba][di][ci] = NULL;
2536  // Example: fQAParticleHist[0][2][4] is distribution of 3rd observable BEFORE the 5th cut was applied
2537  // fQAParticleHist[1][2][4] is distribution of 3rd observable AFTER the 5th cut was applied
2538  } // for(Int_t ci=0;ci<10;ci++) // number of distinct track cuts, maximum is 10 (at the moment)
2539  } // for(Int_t di=0;di<10;di++) // number of distinct distributions, maximum is 10 (at the moment)
2540  } // for(Int_t ba=0;ba<2;ba++) // 0="before rain",1="after rain"
2541 
2542 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForQA()
2543 
2544 //=======================================================================================================================
2545 
2547 {
2548  // Initialize all arrays for common global tracks cuts. In essence, these are the default cuts for "normal" global tracks in AOD.
2549 
2550  // The default values, can be overwritten with the setters SetPtRange(...), etc.
2551  fPtRange[0] = 0.2; // pt min
2552  fPtRange[1] = 5.0; // pt max
2553  fEtaRange[0] = -0.8; // eta min
2554  fEtaRange[1] = 0.8; // eta max
2555  fPhiRange[0] = 0.; // phi min
2556  fPhiRange[1] = TMath::TwoPi(); // phi max
2557 
2558 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForGlobalTrackCuts()
2559 
2560 //=======================================================================================================================
2561 
2563 {
2564  // Isanity checks for global track cuts.
2565 
2566  // Return values:
2567  // 1) pt range;
2568  // 2) eta range;
2569  // 3) phi range;
2570 
2571  if(fPtRange[0]<0.||fPtRange[1]<0.){cout<<Form("fPtRange[0] = %f, fPtRange[1] = %f",fPtRange[0],fPtRange[1])<<endl; return 1;} // pt is positive
2572  if(fPtRange[0]>fPtRange[1]){cout<<Form("fPtRange[0] = %f, fPtRange[1] = %f",fPtRange[0],fPtRange[1])<<endl; return 1;} // ptmin > ptmax
2573 
2574  if(fEtaRange[0]>fEtaRange[1]){cout<<Form("fEtaRange[0] = %f, fEtaRange[1] = %f",fEtaRange[0],fEtaRange[1])<<endl; return 2;} // etamin > etamax
2575 
2576  if(fPhiRange[0]<0.||fPhiRange[1]<0.){cout<<Form("fPhiRange[0] = %f, fPhiRange[1] = %f",fPhiRange[0],fPhiRange[1])<<endl; return 3;} // phi is positive
2577  //if(fPhiRange[0]>TMath::TwoPi()||fPhiRange[1]>TMath::TwoPi()){return 3;} // phi shall be smaller than 2\pi TBI
2578  if(fPhiRange[0]>fPhiRange[1]){cout<<Form("fPhiRange[0] = %f, fPhiRange[1] = %f",fPhiRange[0],fPhiRange[1])<<endl; return 3;} // phimin > phimax
2579 
2580  return 0; // Life is like a wheel...
2581 
2582 } // Int_t AliAnalysisTaskMultiparticleFemtoscopy::InsanityChecksForGlobalTrackCuts()
2583 
2584 //=======================================================================================================================
2585 
2587 {
2588  // Initialize all arrays for test correlation functions.
2589 
2590  const Int_t nTestsMax = 10; // see what is hardwired in .h file
2591 
2592  for(Int_t t=0;t<nTestsMax;t++)
2593  {
2595  fFillCorrelationFunctionsTEST[t] = kFALSE;
2596  }
2597 
2598  for(Int_t t=0;t<nTestsMax;t++)
2599  {
2600  for(Int_t q=0;q<2;q++) // Q2 and Q3
2601  {
2602  for(Int_t tt=0;tt<7;tt++) // 7 distinct terms in the expression of 3p cumulant. TBI for 2p, this is an overkill, but for 4p this will need to be enlarged
2603  {
2604  for(Int_t d=0;d<10;d++) // differential options, see what is hardwired in .h file
2605  {
2606  fCorrelationFunctionsTEST[t][q][tt][d] = NULL;
2607  }
2608  } // for(Int_t t=0;t<nTests;t++)
2609  } // for(Int_t q=0;q<2;q++) // Q2 and Q3
2610  } // for(Int_t t=0;t<nTests;t++)
2611 
2612  for(Int_t t=0;t<nTestsMax;t++)
2613  {
2614  for(Int_t q=0;q<2;q++) // Q2 and Q3
2615  {
2616  for(Int_t tt=0;tt<4;tt++) // 4 distinct cumulants, three 2p and one 3p. TBI for 2p, this is an overkill, but for 4p this will need to be enlarged
2617  {
2618  for(Int_t d=0;d<10;d++) // differential options, see what is hardwired in .h file
2619  {
2620  fSignalCumulantsTEST[t][q][tt][d] = NULL;
2621  }
2622  } // for(Int_t t=0;t<nTests;t++)
2623  } // for(Int_t q=0;q<2;q++) // Q2 and Q3
2624  } // for(Int_t t=0;t<nTests;t++)
2625 
2626  for(Int_t t=0;t<2;t++)
2627  {
2628  fSignalYieldTEST[t] = NULL;
2629  }
2630 
2631 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForCorrelationFunctionsTEST()
2632 
2633 //=======================================================================================================================
2634 
2636 {
2637  // Initialize all arrays for test background.
2638 
2639  const Int_t nTestsMax = 10; // see what is hardwired in .h file
2640 
2641  for(Int_t t=0;t<nTestsMax;t++)
2642  {
2643  fBackgroundTESTSublist[t] = NULL;
2644  fFillBackgroundTEST[t] = kFALSE;
2645  }
2646 
2647  for(Int_t t=0;t<nTestsMax;t++)
2648  {
2649  for(Int_t q=0;q<2;q++) // Q2 and Q3
2650  {
2651  for(Int_t tt=0;tt<7;tt++) // 7 distinct terms in the expression of 3p cumulant. TBI for 2p, this is an overkill, but for 4p this will need to be enlarged
2652  {
2653  for(Int_t d=0;d<10;d++) // differential options, see what is hardwired in .h file
2654  {
2655  fBackgroundTEST[t][q][tt][d] = NULL;
2656  }
2657  } // for(Int_t tt=0;tt<7;tt++)
2658  } // for(Int_t q=0;q<2;q++) // Q2 and Q3
2659  } // for(Int_t t=0;t<nTestsMax;t++)
2660 
2661  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
2662  {
2663  fMixedEventsTEST[me] = NULL;
2664  fGlobalTracksAODTEST[me] = NULL;
2665  }
2666 
2667  for(Int_t t=0;t<2;t++)
2668  {
2669  fBackgroundYieldTEST[t] = NULL;
2670  }
2671 
2672 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForBackgroundTEST()
2673 
2674 //=======================================================================================================================
2675 
2677 {
2678  // Initialize all arrays for 'hybrid approach'.
2679 
2680  for(Int_t t=0;t<5;t++) // five distinct terms in the numerator of Eq. (18)
2681  {
2682  fDistributionsVsQ3[t] = NULL;
2683  }
2684 
2685  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
2686  {
2687  fMixedEventsHA[me] = NULL;
2688  fGlobalTracksAODHA[me] = NULL;
2689  }
2690 
2691 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForHybridApproach()
2692 
2693 //=======================================================================================================================
2694 
2696 {
2697  // Book everything for test correlation functions.
2698 
2699  // a) Book the profile holding all the flags for test correlation functions;
2700  // b) Book correlation functions for all tests;
2701  // c) Store temporarely also the yield.
2702 
2703  // a) Book the profile holding all the flags for test correlation functions:
2704  const Int_t nTests = 10;
2705  fCorrelationFunctionsTESTFlagsPro = new TProfile("fCorrelationFunctionsTESTFlagsPro","Flags and settings for test correlation functions",nTests,0,nTests);
2706  fCorrelationFunctionsTESTFlagsPro->SetTickLength(-0.01,"Y");
2707  fCorrelationFunctionsTESTFlagsPro->SetMarkerStyle(25);
2708  fCorrelationFunctionsTESTFlagsPro->SetLabelSize(0.04);
2709  fCorrelationFunctionsTESTFlagsPro->SetLabelOffset(0.02,"Y");
2710  fCorrelationFunctionsTESTFlagsPro->SetStats(kFALSE);
2711  fCorrelationFunctionsTESTFlagsPro->SetFillColor(kGray);
2712  fCorrelationFunctionsTESTFlagsPro->SetLineColor(kBlack);
2713  for(Int_t t=0;t<nTests;t++)
2714  {
2715  fCorrelationFunctionsTESTFlagsPro->GetXaxis()->SetBinLabel(t+1,Form("fFillCorrelationFunctionsTEST[%d]",t)); fCorrelationFunctionsTESTFlagsPro->Fill(t+0.5,fFillCorrelationFunctionsTEST[t]);
2716  }
2718 
2719  // b) Book correlation functions for all tests:
2720  const Int_t nCumulants = 4; // TBI this applies only to vs. Q3
2721  const Int_t n3pCumulantTerms = 7;
2722  const Int_t nQs = 2; // Q2 and Q3 at the moment, if you change it here, change it also in exprs. like for(Int_t q=0;q<2;q++)
2723  const Int_t n2pCumulantTerms = 3;
2724  TString s2pCumulantTerms[n2pCumulantTerms] = {"#LTX_{1}#GT","#LTX_{2}#GT","#LTX_{1}X_{2}#GT"};
2725  TString s3pCumulantTerms[n3pCumulantTerms] = {"#LTX_{1}#GT","#LTX_{2}#GT","#LTX_{3}#GT","#LTX_{1}X_{2}#GT","#LTX_{1}X_{3}#GT","#LTX_{2}X_{3}#GT","#LTX_{1}X_{2}X_{3}#GT"};
2726  TString sCumulants[nCumulants] = {"#LTX_{1}X_{2}#GT_{c}","#LTX_{1}X_{3}#GT_{c}","#LTX_{2}X_{3}#GT_{c}","#LTX_{1}X_{2}X_{3}#GT_{c}"}; // TBI this applies only to vs. Q3
2727  TString sXYZ[3] = {"x","y","z"};
2728  TString sQs[nQs] = {"Q_{2}","Q_{3}"};
2729  for(Int_t t=0;t<nTests;t++)
2730  {
2732  {
2733 
2734  // Correlations vs. Q2:
2735  for(Int_t ct=0;ct<n2pCumulantTerms;ct++)
2736  {
2737  for(Int_t xyz=0;xyz<3;xyz++)
2738  {
2739  fCorrelationFunctionsTEST[t][0][ct][xyz] = new TProfile(Form("fCorrelationFunctionsTEST[%d][0][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ2bins,fnQ2min,fnQ2max);
2740  fCorrelationFunctionsTEST[t][0][ct][xyz]->SetStats(kFALSE);
2741  fCorrelationFunctionsTEST[t][0][ct][xyz]->SetMarkerStyle(kFullSquare);
2742  fCorrelationFunctionsTEST[t][0][ct][xyz]->SetMarkerColor(kBlue);
2743  fCorrelationFunctionsTEST[t][0][ct][xyz]->SetLineColor(kBlue);
2744  fCorrelationFunctionsTEST[t][0][ct][xyz]->GetXaxis()->SetTitle(sQs[0].Data());
2745  fCorrelationFunctionsTEST[t][0][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",s2pCumulantTerms[ct].Data(),sXYZ[xyz].Data()));
2747  } // for(Int_t xyz=0;xyz<3;xyz++)
2748  } // for(Int_t ct=0;ct<n2pCumulantTerms;ct++)
2749 
2750  // Cumulants vs. Q2:
2751  for(Int_t xyz=0;xyz<3;xyz++)
2752  {
2753  fSignalCumulantsTEST[t][0][0][xyz] = new TProfile(Form("fSignalCumulantsTEST[%d][0][%d][%d]",t,0,xyz),Form("TEST: %d",t),fnQ2bins,fnQ2min,fnQ2max);
2754  fSignalCumulantsTEST[t][0][0][xyz]->SetStats(kFALSE);
2755  fSignalCumulantsTEST[t][0][0][xyz]->SetMarkerStyle(kFullSquare);
2756  fSignalCumulantsTEST[t][0][0][xyz]->SetMarkerColor(kGreen+2);
2757  fSignalCumulantsTEST[t][0][0][xyz]->SetLineColor(kGreen+2);
2758  fSignalCumulantsTEST[t][0][0][xyz]->GetXaxis()->SetTitle(sQs[0].Data());
2759  fSignalCumulantsTEST[t][0][0][xyz]->GetYaxis()->SetTitle(Form("#LTX_{1}X_{2}#GT_{c,%s}",sXYZ[xyz].Data()));
2761  } // for(Int_t xyz=0;xyz<3;xyz++)
2762 
2763  // Correlations vs. Q3: TBI this clearly can be unified with the above.
2764  for(Int_t ct=0;ct<n3pCumulantTerms;ct++)
2765  {
2766  for(Int_t xyz=0;xyz<3;xyz++)
2767  {
2768  fCorrelationFunctionsTEST[t][1][ct][xyz] = new TProfile(Form("fCorrelationFunctionsTEST[%d][1][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ3bins,fnQ3min,fnQ3max);
2769  fCorrelationFunctionsTEST[t][1][ct][xyz]->SetStats(kFALSE);
2770  fCorrelationFunctionsTEST[t][1][ct][xyz]->SetMarkerStyle(kFullSquare);
2771  fCorrelationFunctionsTEST[t][1][ct][xyz]->SetMarkerColor(kBlue);
2772  fCorrelationFunctionsTEST[t][1][ct][xyz]->SetLineColor(kBlue);
2773  fCorrelationFunctionsTEST[t][1][ct][xyz]->GetXaxis()->SetTitle(sQs[1].Data());
2774  fCorrelationFunctionsTEST[t][1][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",s3pCumulantTerms[ct].Data(),sXYZ[xyz].Data()));
2776  } // for(Int_t xyz=0;xyz<3;xyz++)
2777  } // for(Int_t ct=0;ct<n3pCumulantTerms;ct++)
2778 
2779  // Cumulants vs. Q3:
2780  for(Int_t ct=0;ct<nCumulants;ct++) // TBI this applies only to vs. Q3
2781  {
2782  for(Int_t xyz=0;xyz<3;xyz++)
2783  {
2784  fSignalCumulantsTEST[t][1][ct][xyz] = new TProfile(Form("fSignalCumulantsTEST[%d][1][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ3bins,fnQ3min,fnQ3max);
2785  fSignalCumulantsTEST[t][1][ct][xyz]->SetStats(kFALSE);
2786  fSignalCumulantsTEST[t][1][ct][xyz]->SetMarkerStyle(kFullSquare);
2787  fSignalCumulantsTEST[t][1][ct][xyz]->SetMarkerColor(kGreen+2);
2788  fSignalCumulantsTEST[t][1][ct][xyz]->SetLineColor(kGreen+2);
2789  fSignalCumulantsTEST[t][1][ct][xyz]->GetXaxis()->SetTitle(sQs[1].Data());
2790  fSignalCumulantsTEST[t][1][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",sCumulants[ct].Data(),sXYZ[xyz].Data()));
2792  } // for(Int_t xyz=0;xyz<3;xyz++)
2793  } // for(Int_t ct=0;ct<nCumulants;ct++)
2794 
2795  } // if(fFillCorrelationFunctionsTEST[t])
2796 
2797  } // for(Int_t t=0;t<nTests;t++)
2798 
2799  // c) Store temporarely also the yield:
2800  // vs. Q2:
2801  fSignalYieldTEST[0] = new TH1F("fSignalYieldTEST[0]","dN(X_{1}X_{2})/Q_{2} for signal (TEST_0)",fnQ2bins,fnQ2min,fnQ2max);
2802  fSignalYieldTEST[0]->SetStats(kFALSE);
2803  fSignalYieldTEST[0]->SetLineColor(kBlue);
2804  fSignalYieldTEST[0]->SetFillColor(kBlue-10);
2805  fSignalYieldTEST[0]->GetXaxis()->SetTitle("Q_{2}");
2806  fSignalYieldTEST[0]->GetYaxis()->SetTitle("counts");
2808  // vs. Q3:
2809  fSignalYieldTEST[1] = new TH1F("fSignalYieldTEST[1]","dN(X_{1}X_{2}X_{3})/Q_{3} for signal (TEST_1)",fnQ3bins,fnQ3min,fnQ3max);
2810  fSignalYieldTEST[1]->SetStats(kFALSE);
2811  fSignalYieldTEST[1]->SetLineColor(kBlue);
2812  fSignalYieldTEST[1]->SetFillColor(kBlue-10);
2813  fSignalYieldTEST[1]->GetXaxis()->SetTitle("Q_{3}");
2814  fSignalYieldTEST[1]->GetYaxis()->SetTitle("counts");
2816 
2817 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForCorrelationFunctionsTEST()
2818 
2819 //=======================================================================================================================
2820 
2822 {
2823  // Book everything for test background.
2824 
2825  // a) Book the profile holding all the flags for test background;
2826  // b) Book background objects for all tests;
2827  // c) Book fMixedEventsTEST[3] and fGlobalTracksAODTEST[3];
2828  // d) Store temporary the yield as well.
2829 
2830  // a) Book the profile holding all the flags for test background:
2831  const Int_t nTests = 10;
2832  fBackgroundTESTFlagsPro = new TProfile("fBackgroundTESTFlagsPro","Flags and settings for test background",nTests,0,nTests);
2833  fBackgroundTESTFlagsPro->SetTickLength(-0.01,"Y");
2834  fBackgroundTESTFlagsPro->SetMarkerStyle(25);
2835  fBackgroundTESTFlagsPro->SetLabelSize(0.04);
2836  fBackgroundTESTFlagsPro->SetLabelOffset(0.02,"Y");
2837  fBackgroundTESTFlagsPro->SetStats(kFALSE);
2838  fBackgroundTESTFlagsPro->SetFillColor(kGray);
2839  fBackgroundTESTFlagsPro->SetLineColor(kBlack);
2840  for(Int_t t=0;t<nTests;t++)
2841  {
2842  fBackgroundTESTFlagsPro->GetXaxis()->SetBinLabel(t+1,Form("fFillBackgroundTEST[%d]",t)); fBackgroundTESTFlagsPro->Fill(t+0.5,fFillBackgroundTEST[t]);
2843  }
2845 
2846  // b) Book background objects for all tests:
2847  const Int_t nCumulants = 4; // TBI this applies only to vs. Q3
2848  const Int_t n3pCumulantTerms = 7;
2849  const Int_t nQs = 2; // Q2 and Q3 at the moment, if you change it here, change it also in exprs. like for(Int_t q=0;q<2;q++)
2850  const Int_t n2pCumulantTerms = 3;
2851  TString s2pCumulantTerms[n2pCumulantTerms] = {"#LTX_{1}#GT","#LTX_{2}#GT","#LTX_{1}X_{2}#GT"};
2852  TString s3pCumulantTerms[n3pCumulantTerms] = {"#LTX_{1}#GT","#LTX_{2}#GT","#LTX_{3}#GT","#LTX_{1}X_{2}#GT","#LTX_{1}X_{3}#GT","#LTX_{2}X_{3}#GT","#LTX_{1}X_{2}X_{3}#GT"};
2853  TString sCumulants[nCumulants] = {"#LTX_{1}X_{2}#GT_{c}","#LTX_{1}X_{3}#GT_{c}","#LTX_{2}X_{3}#GT_{c}","#LTX_{1}X_{2}X_{3}#GT_{c}"}; // TBI this applies only to vs. Q3
2854  TString sXYZ[3] = {"x","y","z"};
2855  TString sQs[nQs] = {"Q_{2}","Q_{3}"};
2856  for(Int_t t=0;t<nTests;t++)
2857  {
2858  if(fFillBackgroundTEST[t])
2859  {
2860 
2861  // Background correlations vs. Q2:
2862  for(Int_t ct=0;ct<n2pCumulantTerms;ct++)
2863  {
2864  for(Int_t xyz=0;xyz<3;xyz++)
2865  {
2866  fBackgroundTEST[t][0][ct][xyz] = new TProfile(Form("fBackgroundTEST[%d][0][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ2bins,fnQ2min,fnQ2max);
2867  fBackgroundTEST[t][0][ct][xyz]->SetStats(kFALSE);
2868  fBackgroundTEST[t][0][ct][xyz]->SetMarkerStyle(kFullSquare);
2869  fBackgroundTEST[t][0][ct][xyz]->SetMarkerColor(kBlue);
2870  fBackgroundTEST[t][0][ct][xyz]->SetLineColor(kBlue);
2871  fBackgroundTEST[t][0][ct][xyz]->GetXaxis()->SetTitle(sQs[0].Data());
2872  fBackgroundTEST[t][0][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",s2pCumulantTerms[ct].Data(),sXYZ[xyz].Data()));
2873  fBackgroundTESTSublist[t]->Add(fBackgroundTEST[t][0][ct][xyz]);
2874  } // for(Int_t xyz=0;xyz<3;xyz++)
2875  } // for(Int_t ct=0;ct<n2pCumulantTerms;ct++)
2876 
2877  // Background cumulants vs. Q2:
2878  for(Int_t xyz=0;xyz<3;xyz++)
2879  {
2880  fBackgroundCumulantsTEST[t][0][0][xyz] = new TProfile(Form("fBackgroundCumulantsTEST[%d][0][%d][%d]",t,0,xyz),Form("TEST: %d",t),fnQ2bins,fnQ2min,fnQ2max);
2881  fBackgroundCumulantsTEST[t][0][0][xyz]->SetStats(kFALSE);
2882  fBackgroundCumulantsTEST[t][0][0][xyz]->SetMarkerStyle(kFullSquare);
2883  fBackgroundCumulantsTEST[t][0][0][xyz]->SetMarkerColor(kGreen+2);
2884  fBackgroundCumulantsTEST[t][0][0][xyz]->SetLineColor(kGreen+2);
2885  fBackgroundCumulantsTEST[t][0][0][xyz]->GetXaxis()->SetTitle(sQs[0].Data());
2886  fBackgroundCumulantsTEST[t][0][0][xyz]->GetYaxis()->SetTitle(Form("#LTX_{1}X_{2}#GT_{c,%s}",sXYZ[xyz].Data()));
2887  fBackgroundTESTSublist[t]->Add(fBackgroundCumulantsTEST[t][0][0][xyz]);
2888  } // for(Int_t xyz=0;xyz<3;xyz++)
2889 
2890  // Background correlations vs. Q3: TBI this clearly can be unified with the above.
2891  for(Int_t ct=0;ct<n3pCumulantTerms;ct++)
2892  {
2893  for(Int_t xyz=0;xyz<3;xyz++)
2894  {
2895  fBackgroundTEST[t][1][ct][xyz] = new TProfile(Form("fBackgroundTEST[%d][1][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ3bins,fnQ3min,fnQ3max);
2896  fBackgroundTEST[t][1][ct][xyz]->SetStats(kFALSE);
2897  fBackgroundTEST[t][1][ct][xyz]->SetMarkerStyle(kFullSquare);
2898  fBackgroundTEST[t][1][ct][xyz]->SetMarkerColor(kBlue);
2899  fBackgroundTEST[t][1][ct][xyz]->SetLineColor(kBlue);
2900  fBackgroundTEST[t][1][ct][xyz]->GetXaxis()->SetTitle(sQs[1].Data());
2901  fBackgroundTEST[t][1][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",s3pCumulantTerms[ct].Data(),sXYZ[xyz].Data()));
2902  fBackgroundTESTSublist[t]->Add(fBackgroundTEST[t][1][ct][xyz]);
2903  } // for(Int_t xyz=0;xyz<3;xyz++)
2904  } // for(Int_t ct=0;ct<n3pCumulantTerms;ct++)
2905 
2906  // Background cumulants vs. Q3:
2907  for(Int_t ct=0;ct<nCumulants;ct++) // TBI this applies only to vs. Q3
2908  {
2909  for(Int_t xyz=0;xyz<3;xyz++)
2910  {
2911  fBackgroundCumulantsTEST[t][1][ct][xyz] = new TProfile(Form("fBackgroundCumulantsTEST[%d][1][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ3bins,fnQ3min,fnQ3max);
2912  fBackgroundCumulantsTEST[t][1][ct][xyz]->SetStats(kFALSE);
2913  fBackgroundCumulantsTEST[t][1][ct][xyz]->SetMarkerStyle(kFullSquare);
2914  fBackgroundCumulantsTEST[t][1][ct][xyz]->SetMarkerColor(kGreen+2);
2915  fBackgroundCumulantsTEST[t][1][ct][xyz]->SetLineColor(kGreen+2);
2916  fBackgroundCumulantsTEST[t][1][ct][xyz]->GetXaxis()->SetTitle(sQs[1].Data());
2917  fBackgroundCumulantsTEST[t][1][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",sCumulants[ct].Data(),sXYZ[xyz].Data()));
2918  fBackgroundTESTSublist[t]->Add(fBackgroundCumulantsTEST[t][1][ct][xyz]);
2919  } // for(Int_t xyz=0;xyz<3;xyz++)
2920  } // for(Int_t ct=0;ct<nCumulants;ct++)
2921  } // if(fFillCorrelationFunctionsTEST[t])
2922 
2923  } // for(Int_t t=0;t<nTests;t++)
2924 
2925 
2926  // c) Book fMixedEventsTEST[3] and fGlobalTracksAODTEST[3]:
2927  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
2928  {
2929  fMixedEventsTEST[me] = new TClonesArray("AliAODTrack",10000);
2930  fGlobalTracksAODTEST[me] = new TExMap();
2931  }
2932 
2933  // d) Store temporary the yield as well:
2934  // vs. Q2:
2935  fBackgroundYieldTEST[0] = new TH1F("fBackgroundYieldTEST[0]","dN(X_{1}X_{2})/Q_{2} for mixed-events (TEST_0)",fnQ2bins,fnQ2min,fnQ2max);
2936  fBackgroundYieldTEST[0]->SetStats(kFALSE);
2937  fBackgroundYieldTEST[0]->SetLineColor(kRed);
2938  fBackgroundYieldTEST[0]->SetFillColor(kRed-10);
2939  fBackgroundYieldTEST[0]->GetXaxis()->SetTitle("Q_{2}");
2940  fBackgroundYieldTEST[0]->GetYaxis()->SetTitle("counts");
2942  // vs. Q3:
2943  fBackgroundYieldTEST[1] = new TH1F("fBackgroundYieldTEST[1]","dN(X_{1}X_{2}X_{3})/Q_{3} for mixed-events (TEST_1)",fnQ3bins,fnQ3min,fnQ3max);
2944  fBackgroundYieldTEST[1]->SetStats(kFALSE);
2945  fBackgroundYieldTEST[1]->SetLineColor(kRed);
2946  fBackgroundYieldTEST[1]->SetFillColor(kRed-10);
2947  fBackgroundYieldTEST[1]->GetXaxis()->SetTitle("Q_{3}");
2948  fBackgroundYieldTEST[1]->GetYaxis()->SetTitle("counts");
2950 
2951 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForBackgroundTEST()
2952 
2953 //=======================================================================================================================
2954 
2956 {
2957  // Book everything for 'hybrid approach'.
2958 
2959  // a) Book the profile holding all the flags;
2960  // b) Book fMixedEventsHA[3] and fGlobalTracksAODHA[3];
2961  // c) ...
2962 
2963  // a) Book the profile holding all the flags:
2964  fHybridApproachFlagsPro = new TProfile("fHybridApproachFlagsPro","Flags and settings for 'hybrid approach'",1,0,1);
2965  fHybridApproachFlagsPro->SetTickLength(-0.01,"Y");
2966  fHybridApproachFlagsPro->SetMarkerStyle(25);
2967  fHybridApproachFlagsPro->SetLabelSize(0.04);
2968  fHybridApproachFlagsPro->SetLabelOffset(0.02,"Y");
2969  fHybridApproachFlagsPro->SetStats(kFALSE);
2970  fHybridApproachFlagsPro->SetFillColor(kGray);
2971  fHybridApproachFlagsPro->SetLineColor(kBlack);
2972  fHybridApproachFlagsPro->GetXaxis()->SetBinLabel(1,"fDoHybridApproach"); fHybridApproachFlagsPro->Fill(0.5,fDoHybridApproach);
2974 
2975  if(!fDoHybridApproach){return;}
2976 
2977  TString sTerm[5] = {"N^{A}_{3}(p_{1},p_{2},p_{3})","N^{A}_{2}(p_{1},p_{2})N^{B}_{1}(p_{3})","N^{A}_{2}(p_{2},p_{3})N^{B}_{1}(p_{1})","N^{A}_{2}(p_{3},p_{1})N^{B}_{1}(p_{2})","N^{A}_{1}(p_{1})N^{B}_{1}(p_{2})N^{C}_{1}(p_{3})"};
2978  for(Int_t t=0;t<5;t++) // five distinct terms in the numerator of Eq. (18)
2979  {
2980  fDistributionsVsQ3[t] = new TH1F(Form("fDistributionsVsQ3[%d]",t),sTerm[t].Data(),5*fnQ3bins,fnQ3min,fnQ3max); // TBI hardwired refined binning, remove eventually
2981  //fDistributionsVsQ3[t]->SetStats(kFALSE);
2982  fDistributionsVsQ3[t]->SetMarkerStyle(kFullSquare);
2983  fDistributionsVsQ3[t]->SetMarkerColor(kBlue);
2984  fDistributionsVsQ3[t]->SetLineColor(kBlue);
2985  fDistributionsVsQ3[t]->SetFillColor(kBlue-10);
2986  fDistributionsVsQ3[t]->GetXaxis()->SetTitle("Q_{3}");
2987  fDistributionsVsQ3[t]->GetYaxis()->SetTitle("counts");
2989  }
2990 
2991  // b) Book fMixedEventsHA[3] and fGlobalTracksAODHA[3]:
2992  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
2993  {
2994  fMixedEventsHA[me] = new TClonesArray("AliAODTrack",10000);
2995  fGlobalTracksAODHA[me] = new TExMap();
2996  }
2997 
2998  // c) ...
2999 
3000 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForHybridApproach()
3001 
3002 //=======================================================================================================================
3003 
3005 {
3006  // Book and nest all lists nested in the base list fHistList.
3007 
3008  // a) Book and nest lists for control histograms;
3009  // b) Book and nest lists for eveny-by-event histograms;
3010  // c) Book and nest lists for correlation functions;
3011  // d) Book and nest lists for background;
3012  // e) Book and nest lists for buffers;
3013  // f) Book and nest lists for QA;
3014  // g) Book and nest lists for common global track cuts;
3015  // h) Book and nest lists for test correlation functions;
3016  // i) Book and nest lists for test background;
3017  // j) Book and nest lists for 'hybrid approach'.
3018 
3019  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::BookAndNestAllLists()";
3020  if(!fHistList){Fatal(sMethodName.Data(),"fHistList is NULL");}
3021 
3022  // a) Book and nest lists for control histograms:
3023  fControlHistogramsList = new TList();
3024  fControlHistogramsList->SetName("Control_histograms");
3025  fControlHistogramsList->SetOwner(kTRUE);
3028  fControlHistogramsEventList->SetName("Event");
3029  fControlHistogramsEventList->SetOwner(kTRUE);
3032  fControlHistogramsNonIdentifiedParticlesList->SetName("Non-identified_particles");
3036  fControlHistogramsNonIdentifiedParticlesFTSFList->SetName("Non-identified_particles (f.t.s.f.)"); // for the specified filterbit
3040  fControlHistogramsIdentifiedParticlesList->SetName("Identified_particles");
3044  fControlHistogramsV0sList->SetName("V0s");
3045  fControlHistogramsV0sList->SetOwner(kTRUE);
3047 
3048  // b) Book and nest lists for eveny-by-event histograms:
3049  fEBEHistogramsList = new TList();
3050  fEBEHistogramsList->SetName("Event-by-event_histograms");
3051  fEBEHistogramsList->SetOwner(kTRUE);
3053 
3054  // c) Book and nest lists for correlation functions:
3056  fCorrelationFunctionsList->SetName("CorrelationFunctions");
3057  fCorrelationFunctionsList->SetOwner(kTRUE);
3059  // Correlation functions sublists:
3060  for(Int_t cfs=0;cfs<3;cfs++)
3061  {
3062  fCorrelationFunctionsSublist[cfs] = new TList();
3063  fCorrelationFunctionsSublist[cfs]->SetName(Form("f%dpCorrelationFunctions",cfs+2));
3064  fCorrelationFunctionsSublist[cfs]->SetOwner(kTRUE);
3066  } // for(Int_t cfs=0;cfs<3;cfs++)
3067 
3068  // d) Book and nest lists for background:
3069  fBackgroundList = new TList();
3070  fBackgroundList->SetName("Background");
3071  fBackgroundList->SetOwner(kTRUE);
3072  fHistList->Add(fBackgroundList);
3073  // Background sublists:
3074  for(Int_t bs=0;bs<3;bs++)
3075  {
3076  fBackgroundSublist[bs] = new TList();
3077  fBackgroundSublist[bs]->SetName(Form("f%dpBackground",bs+2));
3078  fBackgroundSublist[bs]->SetOwner(kTRUE);
3080  } // for(Int_t bs=0;bs<3;bs++)
3081 
3082  // e) Book and nest lists for buffers:
3083  fBuffersList = new TList();
3084  fBuffersList->SetName("Buffers");
3085  fBuffersList->SetOwner(kTRUE);
3086  fHistList->Add(fBuffersList);
3087 
3088  // f) Book and nest lists for QA:
3089  fQAList = new TList();
3090  fQAList->SetName("QA");
3091  fQAList->SetOwner(kTRUE);
3092  fHistList->Add(fQAList);
3093  if(fFillQAEvents)
3094  {
3095  fQAEventsList = new TList();
3096  fQAEventsList->SetName("Events");
3097  fQAEventsList->SetOwner(kTRUE);
3098  fQAList->Add(fQAEventsList);
3099  }
3100  if(fFillQAParticles)
3101  {
3102  fQAParticlesList = new TList();
3103  fQAParticlesList->SetName("Particles");
3104  fQAParticlesList->SetOwner(kTRUE);
3105  fQAList->Add(fQAParticlesList);
3106  }
3107 
3108  // g) Book and nest lists for common global track cuts:
3109  fGlobalTrackCutsList = new TList();
3110  fGlobalTrackCutsList->SetName("Global_track_cuts");
3111  fGlobalTrackCutsList->SetOwner(kTRUE);
3113 
3114  // h) Book and nest lists for test correlation functions:
3116  fCorrelationFunctionsTESTList->SetName("CorrelationFunctionsTEST");
3117  fCorrelationFunctionsTESTList->SetOwner(kTRUE);
3119  const Int_t nTests = 10;
3120  for(Int_t t=0;t<nTests;t++)
3121  {
3122  if(!fFillCorrelationFunctionsTEST[t]){continue;}
3124  fCorrelationFunctionsTESTSublist[t]->SetName(Form("TEST_%d",t));
3125  fCorrelationFunctionsTESTSublist[t]->SetOwner(kTRUE);
3127  } // for(Int_t t=0;t<nTests;t++)
3128 
3129  // i) Book and nest lists for test background:
3130  fBackgroundTESTList = new TList();
3131  fBackgroundTESTList->SetName("BackgroundTEST");
3132  fBackgroundTESTList->SetOwner(kTRUE);
3134  for(Int_t t=0;t<nTests;t++)
3135  {
3136  if(!fFillBackgroundTEST[t]){continue;}
3137  fBackgroundTESTSublist[t] = new TList();
3138  fBackgroundTESTSublist[t]->SetName(Form("TEST_%d",t));
3139  fBackgroundTESTSublist[t]->SetOwner(kTRUE);
3141  } // for(Int_t t=0;t<nTests;t++)
3142 
3143  // j) Book and nest lists for 'hybrid approach':
3144  fHybridApproachList = new TList();
3145  fHybridApproachList->SetName("HybridApproach");
3146  fHybridApproachList->SetOwner(kTRUE);
3148 
3149 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookAndNestAllLists()
3150 
3151 //=======================================================================================================================
3152 
3154 {
3155  // Book all unclassified objects temporary here. TBI
3156 
3157  fAnalysisType = new TString();
3158 
3159  fPIDResponse = new AliPIDResponse();
3160 
3161  for(Int_t index=0;index<fMaxNoGlobalTracksAOD;index++)
3162  {
3163  fGlobalTracksAOD[index] = new TExMap();
3164  }
3165 
3166 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverything()
3167 
3168 //=======================================================================================================================
3169 
3171 {
3172  // Book all the stuff for control histograms.
3173 
3174  // a) Book the profile holding all the flags for control histograms;
3175  // b) Common variables;
3176  // c) Book all control histograms...
3177  // c0) Event;
3178  // c1) Non-identified particles (for AOD these are "normal global" tracks);
3179  // c2) Non-identified particles for the specified filterbit (f.t.s.f.) (by default TPC-only);
3180  // c3) Identified particles;
3181  // c4) V0s.
3182 
3183  // a) Book the profile holding all the flags for control histograms: TBI stil incomplete
3184  fControlHistogramsFlagsPro = new TProfile("fControlHistogramsFlagsPro","Flags and settings for control histograms",1,0,1);
3185  fControlHistogramsFlagsPro->SetTickLength(-0.01,"Y");
3186  fControlHistogramsFlagsPro->SetMarkerStyle(25);
3187  fControlHistogramsFlagsPro->SetLabelSize(0.04);
3188  fControlHistogramsFlagsPro->SetLabelOffset(0.02,"Y");
3189  fControlHistogramsFlagsPro->SetStats(kFALSE);
3190  fControlHistogramsFlagsPro->SetFillColor(kGray);
3191  fControlHistogramsFlagsPro->SetLineColor(kBlack);
3192  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistograms"); fControlHistogramsFlagsPro->Fill(0.5,fFillControlHistograms);
3194 
3195  if(!fFillControlHistograms){return;}
3196 
3197  // b) Common variables:
3198  TString sParticleLabel[5] = {"e","#mu","#pi","K","p"};
3199  Double_t dNominalMass[5] = {TDatabasePDG::Instance()->GetParticle(11)->Mass(),TDatabasePDG::Instance()->GetParticle(13)->Mass(),TDatabasePDG::Instance()->GetParticle(211)->Mass(),TDatabasePDG::Instance()->GetParticle(321)->Mass(),TDatabasePDG::Instance()->GetParticle(2212)->Mass()};
3200  // ...
3201 
3202  // c0) Event:
3203  // Book the profile holding all the flags for control histograms for global event observables:
3204  fControlHistogramsEventFlagsPro = new TProfile("fControlHistogramsEventFlagsPro","Flags and settings for TBI",1,0,1);
3205  fControlHistogramsEventFlagsPro->SetTickLength(-0.01,"Y");
3206  fControlHistogramsEventFlagsPro->SetMarkerStyle(25);
3207  fControlHistogramsEventFlagsPro->SetLabelSize(0.04);
3208  fControlHistogramsEventFlagsPro->SetLabelOffset(0.02,"Y");
3209  fControlHistogramsEventFlagsPro->SetStats(kFALSE);
3210  fControlHistogramsEventFlagsPro->SetFillColor(kGray);
3211  fControlHistogramsEventFlagsPro->SetLineColor(kBlack);
3212  fControlHistogramsEventFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsEvent"); fControlHistogramsEventFlagsPro->Fill(0.5,fFillControlHistogramsEvent);
3215  {
3216  fGetNumberOfTracksHist = new TH1I("fGetNumberOfTracksHist","aAOD->GetNumberOfTracks() (Remark: Not all of tracks are unique.)",10000,0,10000);
3217  //fGetNumberOfTracksHist->SetStats(kFALSE);
3218  fGetNumberOfTracksHist->SetFillColor(kBlue-10);
3220  fGetNumberOfGlobalTracksHist = new TH1I("fGetNumberOfGlobalTracksHist","fGlobalTracksAOD[0]->GetSize()",10000,0,10000);
3221  //fGetNumberOfGlobalTracksHist->SetStats(kFALSE);
3222  fGetNumberOfGlobalTracksHist->SetFillColor(kBlue-10);
3224  fGetNumberOfV0sHist = new TH1I("fGetNumberOfV0sHist","aAOD->GetNumberOfV0s() (Remark: Some V0s share the daughter.)",10000,0,10000);
3225  fGetNumberOfV0sHist->SetStats(kFALSE);
3226  fGetNumberOfV0sHist->SetFillColor(kBlue-10);
3228  fGetNumberOfCascadesHist = new TH1I("fGetNumberOfCascadesHist","aAOD->GetNumberOfCascades() (TBI: Not validated.)",10000,0,10000);
3229  fGetNumberOfCascadesHist->SetStats(kFALSE);
3230  fGetNumberOfCascadesHist->SetFillColor(kBlue-10);
3232  fGetMagneticFieldHist = new TH1D("fGetMagneticFieldHist","aAOD->GetMagneticField()",200,-10.,10.);
3233  fGetMagneticFieldHist->SetFillColor(kBlue-10);
3235  fGetEventTypeHist = new TH1I("fGetEventTypeHist","aAOD->GetEventType()",1000,0,10);
3236  fGetEventTypeHist->SetFillColor(kBlue-10);
3238  fGetCentralityHist = new TH1D("fGetCentralityHist","aAOD->GetCentrality()",100,0.,100.);
3239  fGetCentralityHist->SetFillColor(kBlue-10);
3241  TString sxyz[3] = {"X","Y","Z"};
3242  for(Int_t xyz=0;xyz<3;xyz++)
3243  {
3244  fVertexXYZ[xyz] = new TH1F(Form("fVertex%s",sxyz[xyz].Data()),Form("avtz->Get%s()",sxyz[xyz].Data()),100000,-50.,50);
3245  fVertexXYZ[xyz]->SetStats(kFALSE);
3247  }
3248  fGetNContributorsHist = new TH1I("fGetNContributorsHist","avtx->GetNContributors()",10000,0,10000);
3249  fGetNContributorsHist->SetStats(kFALSE);
3250  fGetNContributorsHist->SetFillColor(kBlue-10);
3252  fGetChi2perNDFHist = new TH1F("fGetChi2perNDFHist","avtx->GetChi2perNDF()",5000,0.,50.);
3253  fGetChi2perNDFHist->SetStats(kFALSE);
3254  fGetChi2perNDFHist->SetFillColor(kBlue-10);
3256  fGetNDaughtersHist = new TH1I("GetNDaughtersHist","avtx->GetNDaughters()",10000,0,10000);
3257  fGetNDaughtersHist->SetStats(kFALSE);
3258  fGetNDaughtersHist->SetFillColor(kBlue-10);
3260  // ...
3261  } // if(fFillControlHistogramsEvent)
3262 
3263  // c1) Non-identified particles:
3264  // Book the profile holding all the flags for control histograms for non-identified particles:
3265  fControlHistogramsNonIdentifiedParticlesFlagsPro = new TProfile("fControlHistogramsNonIdentifiedParticlesFlagsPro","Flags and settings for TBI",1,0,1);
3266  fControlHistogramsNonIdentifiedParticlesFlagsPro->SetTickLength(-0.01,"Y");
3269  fControlHistogramsNonIdentifiedParticlesFlagsPro->SetLabelOffset(0.02,"Y");
3273  fControlHistogramsNonIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsNonIdentifiedParticles"); fControlHistogramsNonIdentifiedParticlesFlagsPro->Fill(0.5,fFillControlHistogramsNonIdentifiedParticles);
3276  {
3277  fChargeHist = new TH1I("fChargeHist","atrack->Charge()",5,-2,3);
3278  fChargeHist->SetStats(kFALSE);
3279  fChargeHist->SetFillColor(kBlue-10);
3280  fChargeHist->GetXaxis()->SetBinLabel(1,"-2");
3281  fChargeHist->GetXaxis()->SetBinLabel(2,"-1");
3282  fChargeHist->GetXaxis()->SetBinLabel(3,"0");
3283  fChargeHist->GetXaxis()->SetBinLabel(4,"1");
3284  fChargeHist->GetXaxis()->SetBinLabel(5,"2");
3286  fGetTPCNclsHist = new TH1I("fGetTPCNclsHist","atrack->fGetTPCNclsHist()",200,0,200);
3287  fGetTPCNclsHist->SetStats(kFALSE);
3288  fGetTPCNclsHist->SetFillColor(kBlue-10);
3289  fGetTPCNclsHist->GetXaxis()->SetTitle("TPCNcls");
3291  fGetTPCsignalNHist = new TH1I("fGetTPCsignalNHist","atrack->fGetTPCsignalNHist()",200,0,200);
3292  fGetTPCsignalNHist->SetStats(kFALSE);
3293  fGetTPCsignalNHist->SetFillColor(kBlue-10);
3294  fGetTPCsignalNHist->GetXaxis()->SetTitle("TPCsignalN");
3296  fGetITSNclsHist = new TH1I("fGetITSNclsHist","atrack->fGetITSNclsHist()",200,0,200);
3297  fGetITSNclsHist->SetStats(kFALSE);
3298  fGetITSNclsHist->SetFillColor(kBlue-10);
3299  fGetITSNclsHist->GetXaxis()->SetTitle("ITSNcls");
3301  fdEdxVsPtHist = new TH2F("fdEdxVsPtHist","atrack->GetTPCmomentum(),atrack->GetTPCsignal()",1000,0.,20.,1000,-500.,500.);
3302  fdEdxVsPtHist->SetStats(kFALSE);
3304  fPtHist = new TH1F("fPtHist","atrack->Pt()",1000,0.,20.);
3305  fPtHist->SetStats(kFALSE);
3306  fPtHist->SetFillColor(kBlue-10);
3307  fPtHist->SetMinimum(0.);
3309  fEtaHist = new TH1F("fEtaHist","atrack->Eta()",200,-2.,2.);
3310  fEtaHist->SetStats(kFALSE);
3311  fEtaHist->SetFillColor(kBlue-10);
3312  fEtaHist->SetMinimum(0.);
3314  fPhiHist = new TH1F("fPhiHist","atrack->Phi()",360,0.,TMath::TwoPi());
3315  fPhiHist->SetStats(kFALSE);
3316  fPhiHist->SetFillColor(kBlue-10);
3317  fPhiHist->SetMinimum(0.);
3319  // TBI
3320  fMassHist = new TH1F("fMassHist","atrack->M()",10000,0.,10.);
3321  fMassHist->SetStats(kFALSE);
3322  fMassHist->SetFillColor(kBlue-10);
3323  fMassHist->SetMinimum(0.);
3324  for(Int_t nm=0;nm<5;nm++) // nominal masses
3325  {
3326  fMassHist->GetXaxis()->SetBinLabel(fMassHist->FindBin(dNominalMass[nm]),Form("m_{%s}",sParticleLabel[nm].Data()));
3327  }
3329  fGetFilterMap = new TH1I("fGetFilterMap","atrack->fGetFilterMap()",10000,0,10000);
3330  fGetFilterMap->SetStats(kFALSE);
3331  fGetFilterMap->SetFillColor(kBlue-10);
3332  fGetFilterMap->SetMinimum(0.);
3334  fGetPdgCode = new TH1I("fGetPdgCode","atrack->fGetPdgCode()",20000,-10000,10000);
3335  fGetPdgCode->SetStats(kFALSE);
3336  fGetPdgCode->SetFillColor(kBlue-10);
3337  fGetPdgCode->SetMinimum(0.);
3339  } // if(fFillControlHistogramsNonIdentifiedParticles)
3340 
3341  // c2) Non-identified particles for the specified filterbit (f.t.s.f.) (by default TPC-only);
3342  // Book the profile holding all the flags for control histograms for non-identified particles for the specified filterbit (f.t.s.f.) (by default TPC-only);
3343  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro = new TProfile("fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro","Flags and settings for TBI",2,0,2);
3344  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro->SetTickLength(-0.01,"Y");
3347  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro->SetLabelOffset(0.02,"Y");
3351  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsNonIdentifiedParticlesFTSF"); fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro->Fill(0.5,fFillControlHistogramsNonIdentifiedParticlesFTSF);
3355  {
3356  fChargeFTSFHist = new TH1I("fChargeFTSFHist",Form("atrack->Charge(), fb = %d",fFilterBitFTSF),5,-2,3);
3357  fChargeFTSFHist->SetStats(kFALSE);
3358  fChargeFTSFHist->SetFillColor(kBlue-10);
3359  fChargeFTSFHist->GetXaxis()->SetBinLabel(1,"-2");
3360  fChargeFTSFHist->GetXaxis()->SetBinLabel(2,"-1");
3361  fChargeFTSFHist->GetXaxis()->SetBinLabel(3,"0");
3362  fChargeFTSFHist->GetXaxis()->SetBinLabel(4,"1");
3363  fChargeFTSFHist->GetXaxis()->SetBinLabel(5,"2");
3365  fGetTPCNclsFTSFHist = new TH1I("fGetTPCNclsFTSFHist",Form("atrack->fGetTPCNclsFTSFHist(), fb = %d",fFilterBitFTSF),200,0,200);
3366  fGetTPCNclsFTSFHist->SetStats(kFALSE);
3367  fGetTPCNclsFTSFHist->SetFillColor(kBlue-10);
3368  fGetTPCNclsFTSFHist->GetXaxis()->SetTitle("TPCNcls");
3370  fGetTPCsignalNFTSFHist = new TH1I("fGetTPCsignalNFTSFHist",Form("atrack->fGetTPCsignalNFTSFHist(), fb = %d",fFilterBitFTSF),200,0,200);
3371  fGetTPCsignalNFTSFHist->SetStats(kFALSE);
3372  fGetTPCsignalNFTSFHist->SetFillColor(kBlue-10);
3373  fGetTPCsignalNFTSFHist->GetXaxis()->SetTitle("TPCsignalN");
3375  fGetITSNclsFTSFHist = new TH1I("fGetITSNclsFTSFHist",Form("atrack->fGetITSNclsFTSFHist(), fb = %d",fFilterBitFTSF),200,0,200);
3376  fGetITSNclsFTSFHist->SetStats(kFALSE);
3377  fGetITSNclsFTSFHist->SetFillColor(kBlue-10);
3378  fGetITSNclsFTSFHist->GetXaxis()->SetTitle("ITSNcls");
3380  fdEdxVsPtFTSFHist = new TH2F("fdEdxVsPtFTSFHist",Form("atrack->GetTPCmomentum(),atrack->GetTPCsignal(), fb = %d",fFilterBitFTSF),1000,0.,20.,1000,-500.,500.);
3381  fdEdxVsPtFTSFHist->SetStats(kFALSE);
3383  fPtFTSFHist = new TH1F("fPtFTSFHist",Form("atrack->Pt(), fb = %d",fFilterBitFTSF),1000,0.,20.);
3384  fPtFTSFHist->SetStats(kFALSE);
3385  fPtFTSFHist->SetFillColor(kBlue-10);
3386  fPtFTSFHist->SetMinimum(0.);
3388  fEtaFTSFHist = new TH1F("fEtaFTSFHist",Form("atrack->Eta(), fb = %d",fFilterBitFTSF),200,-2.,2.);
3389  fEtaFTSFHist->SetStats(kFALSE);
3390  fEtaFTSFHist->SetFillColor(kBlue-10);
3391  fEtaFTSFHist->SetMinimum(0.);
3393  fPhiFTSFHist = new TH1F("fPhiFTSFHist",Form("atrack->Phi(), fb = %d",fFilterBitFTSF),360,0.,TMath::TwoPi());
3394  fPhiFTSFHist->SetStats(kFALSE);
3395  fPhiFTSFHist->SetFillColor(kBlue-10);
3396  fPhiFTSFHist->SetMinimum(0.);
3398  // TBI
3399  fMassFTSFHist = new TH1F("fMassFTSFHist",Form("atrack->M(), fb = %d",fFilterBitFTSF),10000,0.,10.);
3400  fMassFTSFHist->SetStats(kFALSE);
3401  fMassFTSFHist->SetFillColor(kBlue-10);
3402  fMassFTSFHist->SetMinimum(0.);
3403  for(Int_t nm=0;nm<5;nm++) // nominal masses
3404  {
3405  fMassFTSFHist->GetXaxis()->SetBinLabel(fMassFTSFHist->FindBin(dNominalMass[nm]),Form("m_{%s}",sParticleLabel[nm].Data()));
3406  }
3408  fGetFilterMap = new TH1I("fGetFilterMap",Form("atrack->fGetFilterMap(), fb = %d",fFilterBitFTSF),10000,0,10000);
3409  fGetFilterMap->SetStats(kFALSE);
3410  fGetFilterMap->SetFillColor(kBlue-10);
3411  fGetFilterMap->SetMinimum(0.);
3413  fGetPdgCode = new TH1I("fGetPdgCode",Form("atrack->fGetPdgCode(), fb = %d",fFilterBitFTSF),20000,-10000,10000);
3414  fGetPdgCode->SetStats(kFALSE);
3415  fGetPdgCode->SetFillColor(kBlue-10);
3416  fGetPdgCode->SetMinimum(0.);
3418  } // if(fFillControlFTSFHistogramsNonIdentifiedParticlesFTSF)
3419 
3420  // c3) Identified particles:
3421  fControlHistogramsIdentifiedParticlesFlagsPro = new TProfile("fControlHistogramsIdentifiedParticlesFlagsPro","Flags and settings for identified particles",4,0,4);
3422  fControlHistogramsIdentifiedParticlesFlagsPro->SetTickLength(-0.01,"Y");
3425  fControlHistogramsIdentifiedParticlesFlagsPro->SetLabelOffset(0.02,"Y");
3427  fControlHistogramsIdentifiedParticlesFlagsPro->SetFillColor(kGray);
3428  fControlHistogramsIdentifiedParticlesFlagsPro->SetLineColor(kBlack);
3429  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsIdentifiedParticles"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(0.5,fFillControlHistogramsIdentifiedParticles);
3430  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(2,"fUseDefaultInclusiveSigmaCuts"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(1.5,(Int_t)fUseDefaultInclusiveSigmaCuts);
3431  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(3,"fUseDefaultExclusiveSigmaCuts"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(2.5,(Int_t)fUseDefaultExclusiveSigmaCuts);
3432  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(4,"fFillControlHistogramsWithGlobalTrackInfo"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(3.5,(Int_t)fFillControlHistogramsWithGlobalTrackInfo);
3434 
3436  {
3437  // Default inclusive sigma cuts:
3438  if(fUseDefaultInclusiveSigmaCuts)
3439  {
3440  fInclusiveSigmaCuts[2] = 3.; // i.e. in function Pion(...) the inclusive cut for pions is 2 sigma
3441  fInclusiveSigmaCuts[3] = 3.; // i.e. in function Kaon(...) the inclusive cut for kaons is 2 sigma
3442  fInclusiveSigmaCuts[4] = 3.; // i.e. in function Proton(...) the inclusive cut for protons is 2 sigma
3443  }
3444  const Int_t nPidFunctions = 5; //
3445  TString sPidFunctions[nPidFunctions] = {"Electron(...)","Muon(...)","Pion(...)","Kaon(...)","Proton(...)"};
3446  const Int_t nParticleSpecies = 5;
3447  TString sParticleSpecies[nParticleSpecies] = {"e","#mu","#pi","K","p"};
3448  fInclusiveSigmaCutsPro = new TProfile("fInclusiveSigmaCutsPro","Inclusive sigma cuts",nPidFunctions,0.,nPidFunctions);
3449  fInclusiveSigmaCutsPro->GetYaxis()->SetTitle("n#sigma");
3450  for(Int_t pidFunction=0;pidFunction<5;pidFunction++)
3451  {
3452  fInclusiveSigmaCutsPro->SetTickLength(-0.01,"Y");
3453  fInclusiveSigmaCutsPro->SetMarkerStyle(25);
3454  fInclusiveSigmaCutsPro->SetLabelSize(0.04);
3455  fInclusiveSigmaCutsPro->SetLabelOffset(0.02,"Y");
3456  fInclusiveSigmaCutsPro->SetStats(kFALSE);
3457  fInclusiveSigmaCutsPro->SetFillColor(kGray);
3458  fInclusiveSigmaCutsPro->SetLineColor(kBlack);
3459  fInclusiveSigmaCutsPro->GetXaxis()->SetBinLabel(pidFunction+1,sPidFunctions[pidFunction].Data());
3460  fInclusiveSigmaCutsPro->Fill(pidFunction+0.5,fInclusiveSigmaCuts[pidFunction]);
3461  } // for(Int_t pidFunction=0;pidFunction<5;pidFunction++)
3463 
3464  // Default exclusive sigma cuts:
3465  if(fUseDefaultExclusiveSigmaCuts)
3466  {
3467  // Pion(...)
3468  fExclusiveSigmaCuts[2][3] = 4.; // i.e. in function Pion(...) the exclusive cut for kaons is 4 sigma
3469  fExclusiveSigmaCuts[2][4] = 4.; // i.e. in function Pion(...) the exclusive cut for protons is 4 sigma
3470  // Kaon(...)
3471  fExclusiveSigmaCuts[3][2] = 4.; // i.e. in function Kaon(...) the exclusive cut for pions is 4 sigma
3472  fExclusiveSigmaCuts[3][4] = 4.; // i.e. in function Kaon(...) the exclusive cut for protons is 4 sigma
3473  // Proton(...)
3474  fExclusiveSigmaCuts[4][2] = 4.; // i.e. in function Proton(...) the exclusive cut for pions is 4 sigma
3475  fExclusiveSigmaCuts[4][3] = 4.; // i.e. in function Proton(...) the exclusive cut for kaons is 4 sigma
3476  } // if(fUseDefaultExclusiveSigmaCuts)
3477  fExclusiveSigmaCutsPro = new TProfile2D("fExclusiveSigmaCutsPro","Exclusive sigma cuts",nPidFunctions,0.,nPidFunctions,nParticleSpecies,0.,nParticleSpecies);
3478  fExclusiveSigmaCutsPro->SetTickLength(-0.01,"Y");
3479  fExclusiveSigmaCutsPro->SetMarkerStyle(25);
3480  fExclusiveSigmaCutsPro->SetLabelSize(0.04);
3481  //fExclusiveSigmaCutsPro->SetLabelOffset(0.01,"X");
3482  //fExclusiveSigmaCutsPro->SetLabelOffset(0.01,"Y");
3483  fExclusiveSigmaCutsPro->SetTitleOffset(0.9,"Z");
3484  fExclusiveSigmaCutsPro->SetStats(kFALSE);
3485  fExclusiveSigmaCutsPro->SetFillColor(kGray);
3486  fExclusiveSigmaCutsPro->SetLineColor(kBlack);
3487  fExclusiveSigmaCutsPro->GetZaxis()->SetTitle("n#sigma");
3488  for(Int_t pidFunction=0;pidFunction<nPidFunctions;pidFunction++)
3489  {
3490  fExclusiveSigmaCutsPro->GetXaxis()->SetBinLabel(pidFunction+1,sPidFunctions[pidFunction].Data());
3491  for(Int_t particleSpecies=0;particleSpecies<nParticleSpecies;particleSpecies++)
3492  {
3493  if(0==pidFunction){fExclusiveSigmaCutsPro->GetYaxis()->SetBinLabel(particleSpecies+1,sParticleSpecies[particleSpecies].Data());}
3494  fExclusiveSigmaCutsPro->Fill(pidFunction+0.5,particleSpecies+0.5,fExclusiveSigmaCuts[pidFunction][particleSpecies]);
3495  } // for(Int_t pidFunction=0;pidFunction<nPidFunctions;pidFunction++)
3496  } // for(Int_t pidFunction=0;pidFunction<nPidFunctions;pidFunction++)
3498 
3499  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
3500  {
3501  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
3502  {
3503  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
3504  {
3505  if(fFillControlHistogramsWithGlobalTrackInfo)
3506  {
3507  fMassPIDHist[pid][pa][ps] = new TH1F(Form("fMassPIDHist[%d][%d][%d]",pid,pa,ps),Form("fMassPIDHist[%d][%d][%d] (%s) ('gtrack' parameters)",pid,pa,ps,sParticleLabel[pid].Data()),10000,0.,10.);
3508  }
3509  else
3510  {
3511  fMassPIDHist[pid][pa][ps] = new TH1F(Form("fMassPIDHist[%d][%d][%d]",pid,pa,ps),Form("fMassPIDHist[%d][%d][%d] (%s) ('atrack' parameters)",pid,pa,ps,sParticleLabel[pid].Data()),10000,0.,10.);
3512  }
3513  fMassPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("m [GeV/c^{2}]");
3514  for(Int_t nm=0;nm<5;nm++) // nominal masses
3515  {
3516  fMassPIDHist[pid][pa][ps]->GetXaxis()->SetBinLabel(fMassPIDHist[pid][pa][ps]->FindBin(dNominalMass[nm]),Form("m_{%s}",sParticleLabel[nm].Data()));
3517  }
3518  fMassPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
3520  if(fFillControlHistogramsWithGlobalTrackInfo)
3521  {
3522  fPtPIDHist[pid][pa][ps] = new TH1F(Form("fPtPIDHist[%d][%d][%d]",pid,pa,ps),Form("fPtPIDHist[%d][%d][%d] ('gtrack' parameters)",pid,pa,ps),1000,0.,10.);
3523  fPPIDHist[pid][pa][ps][0] = new TH1F(Form("fPPIDHist[%d][%d][%d][0]",pid,pa,ps),Form("fPPIDHist[%d][%d][%d][0] ('gtrack' parameters)",pid,pa,ps),2000,-10.,10.);
3524  fPPIDHist[pid][pa][ps][1] = new TH1F(Form("fPPIDHist[%d][%d][%d][1]",pid,pa,ps),Form("fPPIDHist[%d][%d][%d][1] ('gtrack' parameters)",pid,pa,ps),2000,-10.,10.);
3525  fPPIDHist[pid][pa][ps][2] = new TH1F(Form("fPPIDHist[%d][%d][%d][2]",pid,pa,ps),Form("fPPIDHist[%d][%d][%d][2] ('gtrack' parameters)",pid,pa,ps),2000,-10.,10.);
3526  }
3527  else
3528  {
3529  fPtPIDHist[pid][pa][ps] = new TH1F(Form("fPtPIDHist[%d][%d][%d]",pid,pa,ps),Form("fPtPIDHist[%d][%d][%d] ('atrack' parameters)",pid,pa,ps),1000,0.,10.);
3530  fPPIDHist[pid][pa][ps][0] = new TH1F(Form("fPPIDHist[%d][%d][%d][0]",pid,pa,ps),Form("fPPIDHist[%d][%d][%d][0] ('atrack' parameters)",pid,pa,ps),2000,-10.,10.);
3531  fPPIDHist[pid][pa][ps][1] = new TH1F(Form("fPPIDHist[%d][%d][%d][1]",pid,pa,ps),Form("fPPIDHist[%d][%d][%d][1] ('atrack' parameters)",pid,pa,ps),2000,-10.,10.);
3532  fPPIDHist[pid][pa][ps][2] = new TH1F(Form("fPPIDHist[%d][%d][%d][2]",pid,pa,ps),Form("fPPIDHist[%d][%d][%d][2] ('atrack' parameters)",pid,pa,ps),2000,-10.,10.);
3533  }
3534  fPtPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
3535  fPtPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
3537  fPPIDHist[pid][pa][ps][0]->GetXaxis()->SetTitle("p_{x} [GeV/c]");
3538  fPPIDHist[pid][pa][ps][0]->SetFillColor(kBlue-10);
3540  fPPIDHist[pid][pa][ps][1]->GetXaxis()->SetTitle("p_{y} [GeV/c]");
3541  fPPIDHist[pid][pa][ps][1]->SetFillColor(kBlue-10);
3543  fPPIDHist[pid][pa][ps][2]->GetXaxis()->SetTitle("p_{z} [GeV/c]");
3544  fPPIDHist[pid][pa][ps][2]->SetFillColor(kBlue-10);
3546 
3547  if(fFillControlHistogramsWithGlobalTrackInfo)
3548  {
3549  fEtaPIDHist[pid][pa][ps] = new TH1F(Form("fEtaPIDHist[%d][%d][%d]",pid,pa,ps),Form("fEtaPIDHist[%d][%d][%d] ('gtrack' parameters)",pid,pa,ps),200000,-2.,2.);
3550  }
3551  else
3552  {
3553  fEtaPIDHist[pid][pa][ps] = new TH1F(Form("fEtaPIDHist[%d][%d][%d]",pid,pa,ps),Form("fEtaPIDHist[%d][%d][%d] ('atrack' parameters)",pid,pa,ps),200000,-2.,2.);
3554  }
3555  fEtaPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
3557  if(fFillControlHistogramsWithGlobalTrackInfo)
3558  {
3559  fPhiPIDHist[pid][pa][ps] = new TH1F(Form("fPhiPIDHist[%d][%d][%d]",pid,pa,ps),Form("fPhiPIDHist[%d][%d][%d] ('gtrack' parameters)",pid,pa,ps),360,0.,TMath::TwoPi());
3560  }
3561  else
3562  {
3563  fPhiPIDHist[pid][pa][ps] = new TH1F(Form("fPhiPIDHist[%d][%d][%d]",pid,pa,ps),Form("fPhiPIDHist[%d][%d][%d] ('atrack' parameters)",pid,pa,ps),360,0.,TMath::TwoPi());
3564  }
3565  fPhiPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("#phi");
3566  fPhiPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
3568  } // for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
3569  } // for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
3570  } // for(Int_t pid=0;pid<4;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
3571  } // if(fFillControlHistogramsIdentifiedParticles)
3572 
3573  // c4) V0s:
3574  // Book the profile holding all the flags for V0s:
3575  fControlHistogramsV0sFlagsPro = new TProfile("fControlHistogramsV0sFlagsPro","Flags and settings for V0s",1,0,1);
3576  fControlHistogramsV0sFlagsPro->SetTickLength(-0.01,"Y");
3577  fControlHistogramsV0sFlagsPro->SetMarkerStyle(25);
3578  fControlHistogramsV0sFlagsPro->SetLabelSize(0.04);
3579  fControlHistogramsV0sFlagsPro->SetLabelOffset(0.02,"Y");
3580  fControlHistogramsV0sFlagsPro->SetStats(kFALSE);
3581  fControlHistogramsV0sFlagsPro->SetFillColor(kGray);
3582  fControlHistogramsV0sFlagsPro->SetLineColor(kBlack);
3583  fControlHistogramsV0sFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsV0s"); fControlHistogramsV0sFlagsPro->Fill(0.5,fFillControlHistogramsV0s);
3586  {
3587  fGetNProngsHist = new TH1I("fGetNProngsHist","aAODv0->GetNProngs()",10,0,10);
3588  fGetNProngsHist->SetStats(kFALSE);
3589  fGetNProngsHist->SetFillColor(kBlue-10);
3591  // TBI
3592  fMassK0ShortHist = new TH1F("fMassK0ShortHist","aAODv0->MassK0Short()",1000000,0.,100.);
3593  //fMassK0ShortHist->SetStats(kFALSE);
3594  fMassK0ShortHist->SetFillColor(kBlue-10);
3595  Double_t dMassK0Short = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // nominal mass
3596  //fMassK0ShortHist->GetXaxis()->SetBinLabel(fMassK0ShortHist->FindBin(dMassK0Short),Form("m_{K_{S}^{0}} = %f",dMassK0Short));
3597  fMassK0ShortHist->SetBinContent(fMassK0ShortHist->FindBin(dMassK0Short),1e6);
3599  // TBI
3600  fMassLambdaHist = new TH1F("fMassLambdaHist","aAODv0->MassLambda()",1000000,0.,100.);
3601  //fMassLambdaHist->SetStats(kFALSE);
3602  fMassLambdaHist->SetFillColor(kBlue-10);
3603  Double_t dMassLambda = TDatabasePDG::Instance()->GetParticle(3122)->Mass(); // nominal mass
3604  //fMassLambdaHist->GetXaxis()->SetBinLabel(fMassLambdaHist->FindBin(dMassLambda),Form("m_{#Lambda^{0}} = %f",dMassLambda));
3605  fMassLambdaHist->SetBinContent(fMassLambdaHist->FindBin(dMassLambda),1e6);
3607  // TBI
3608  fMassAntiLambdaHist = new TH1F("fMassAntiLambdaHist","aAODv0->MassAntiLambda()",1000000,0.,100.);
3609  //fMassAntiLambdaHist->SetStats(kFALSE);
3610  fMassAntiLambdaHist->SetFillColor(kBlue-10);
3611  Double_t dMassAntiLambda = TDatabasePDG::Instance()->GetParticle(-3122)->Mass(); // nominal mass
3612  //fMassAntiLambdaHist->GetXaxis()->SetBinLabel(fMassAntiLambdaHist->FindBin(dMassAntiLambda),Form("m_{#bar{Lambda}^{0}} = %f",dMassAntiLambda));
3613  fMassAntiLambdaHist->SetBinContent(fMassAntiLambdaHist->FindBin(dMassAntiLambda),1e6);
3615  // TBI
3616  fOpenAngleV0Hist = new TH1F("fOpenAngleV0Hist","aAODv0->fOpenAngleV0()",10000,-0.044,TMath::Pi()+0.044);
3617  fOpenAngleV0Hist->SetStats(kFALSE);
3618  fOpenAngleV0Hist->SetFillColor(kBlue-10);
3620  // TBI
3621  fRadiusV0Hist = new TH1F("fRadiusV0Hist","aAODv0->fRadiusV0()",10000,0.,1000.);
3622  fRadiusV0Hist->SetStats(kFALSE);
3623  fRadiusV0Hist->SetFillColor(kBlue-10);
3625  // TBI
3626  fDcaV0ToPrimVertexHist = new TH1F("fDcaV0ToPrimVertexHist","aAODv0->fDcaV0ToPrimVertex()",10000,0.,1000.);
3627  fDcaV0ToPrimVertexHist->SetStats(kFALSE);
3628  fDcaV0ToPrimVertexHist->SetFillColor(kBlue-10);
3630  // TBI
3631  fMomV0XHist = new TH1F("fMomV0XHist","aAODv0->fMomV0X() = px(+) + px(-)",10000,-1000.,1000.);
3632  fMomV0XHist->SetStats(kFALSE);
3633  fMomV0XHist->SetFillColor(kBlue-10);
3635  // TBI
3636  fMomV0YHist = new TH1F("fMomV0YHist","aAODv0->fMomV0Y() = py(+) + py(-)",10000,-1000.,1000.);
3637  fMomV0YHist->SetStats(kFALSE);
3638  fMomV0YHist->SetFillColor(kBlue-10);
3640  // TBI
3641  fMomV0ZHist = new TH1F("fMomV0ZHist","aAODv0->fMomV0Z() = pz(+) + pz(-)",10000,-1000.,1000.);
3642  fMomV0ZHist->SetStats(kFALSE);
3643  fMomV0ZHist->SetFillColor(kBlue-10);
3645  // TBI
3646  fPtV0Hist = new TH1F("fPtV0Hist","pow(aAODv0->fPt2V0(),0.5)",10000,0.,100.);
3647  fPtV0Hist->SetStats(kFALSE);
3648  fPtV0Hist->SetFillColor(kBlue-10);
3650  // TBI
3651  fPseudoRapV0Hist = new TH1F("fPseudoRapV0Hist","aAODv0->PseudoRapV0()",1000,-10.,10.);
3652  fPseudoRapV0Hist->SetStats(kFALSE);
3653  fPseudoRapV0Hist->SetFillColor(kBlue-10);
3655  // TBI
3656  fPAHist = new TH2F("fPAHist","TBI",100,-2.,2.,100,0.,1.);
3657  fPAHist->SetStats(kFALSE);
3658  fPAHist->GetXaxis()->SetTitle("#alpha");
3659  fPAHist->GetYaxis()->SetTitle("p_{T}");
3661  } // if(fFillControlHistogramsV0s)
3662 
3663 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForControlHistograms()
3664 
3665 //=======================================================================================================================
3666 
3668 {
3669  // Book all the stuff for event-by-event objects.
3670 
3671  // a) Book the profile holding all the flags for event-by-event objects;
3672  // b) Book all event-by-event objects.
3673 
3674  // a) Book the profile holding all the flags for EBE objects:
3675  fEBEObjectsFlagsPro = new TProfile("fEBEObjectsFlagsPro","Flags and settings for event-by-event histograms",1,0,1);
3676  fEBEObjectsFlagsPro->SetTickLength(-0.01,"Y");
3677  fEBEObjectsFlagsPro->SetMarkerStyle(25);
3678  fEBEObjectsFlagsPro->SetLabelSize(0.04);
3679  fEBEObjectsFlagsPro->SetLabelOffset(0.02,"Y");
3680  fEBEObjectsFlagsPro->SetStats(kFALSE);
3681  fEBEObjectsFlagsPro->SetFillColor(kGray);
3682  fEBEObjectsFlagsPro->SetLineColor(kBlack);
3683  //fEBEObjectsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillEBEHistograms"); fEBEObjectsFlagsPro->Fill(0.5,fFillEBEHistograms);
3685 
3686  //if(!fFillEBEHistograms){return;} // TBI rethink
3687 
3688  // TBI
3689  fUniqueIDHistEBE = new TH1I("fUniqueIDHistEBE","TBI",40000,-20000,20000);
3690  fUniqueIDHistEBE->SetStats(kFALSE);
3691  fUniqueIDHistEBE->SetFillColor(kBlue-10);
3692  // TBI I do not want to store this histogram, right?
3693 
3694  // TBI add comment
3695  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
3696  {
3697  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
3698  {
3699  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
3700  {
3701  fPIDCA[pid][pa][ps] = new TClonesArray("AliAODTrack",10000);
3702  }
3703  }
3704  }
3705 
3706  // TBI add comment
3707  for(Int_t pid=0;pid<1;pid++) // [0=Lambda,1=...]
3708  {
3709  fPIDV0sCA[pid] = new TClonesArray("AliAODv0",10000);
3710  }
3711 
3712 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForEBEObjects()
3713 
3714 //=======================================================================================================================
3715 
3717 {
3718  // Book all the stuff for correlation functions objects.
3719 
3720  // a) Book the profile holding all the flags for correlation functions objects;
3721  // b) Book all 2p correlation functions;
3722  // c) Book TExMap *fCorrelationFunctionsIndices;
3723  // d) Book all 3p correlation functions;
3724  // e) Book all 4p correlation functions.
3725 
3726  // a) Book the profile holding all the flags for correlation functions objects:
3727  fCorrelationFunctionsFlagsPro = new TProfile("fCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",8,0,8);
3728  fCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
3729  fCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
3730  fCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
3731  fCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
3732  fCorrelationFunctionsFlagsPro->SetStats(kFALSE);
3733  fCorrelationFunctionsFlagsPro->SetFillColor(kGray);
3734  fCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
3735  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(0.5,fFillCorrelationFunctions);
3736  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(2,"fNormalizeCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(1.5,fNormalizeCorrelationFunctions);
3737  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(3,"fFill3pCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(2.5,fFill3pCorrelationFunctions);
3738  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(4,"fFill4pCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(3.5,fFill4pCorrelationFunctions);
3739  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(5,"fNormalizationOption"); fCorrelationFunctionsFlagsPro->Fill(4.5,fNormalizationOption);
3740  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(6,"fNormalizationInterval[0]"); fCorrelationFunctionsFlagsPro->Fill(5.5,fNormalizationInterval[0]);
3741  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(7,"fNormalizationInterval[1]"); fCorrelationFunctionsFlagsPro->Fill(6.5,fNormalizationInterval[1]);
3742  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(8,"fnMergedBins"); fCorrelationFunctionsFlagsPro->Fill(7.5,fnMergedBins);
3744  // 2p:
3745  f2pCorrelationFunctionsFlagsPro = new TProfile("f2pCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",1,0,1);
3746  f2pCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
3747  f2pCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
3748  f2pCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
3749  f2pCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
3750  f2pCorrelationFunctionsFlagsPro->SetStats(kFALSE);
3751  f2pCorrelationFunctionsFlagsPro->SetFillColor(kGray);
3752  f2pCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
3753  f2pCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f2pCorrelationFunctionsFlagsPro->Fill(0.5,0.); // TBI dummy at the moment
3755  // 3p:
3756  f3pCorrelationFunctionsFlagsPro = new TProfile("f3pCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",1,0,1);
3757  f3pCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
3758  f3pCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
3759  f3pCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
3760  f3pCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
3761  f3pCorrelationFunctionsFlagsPro->SetStats(kFALSE);
3762  f3pCorrelationFunctionsFlagsPro->SetFillColor(kGray);
3763  f3pCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
3764  f3pCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f3pCorrelationFunctionsFlagsPro->Fill(0.5,0.); // TBI dummy at the moment
3766  // 4p:
3767  f4pCorrelationFunctionsFlagsPro = new TProfile("f4pCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",1,0,1);
3768  f4pCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
3769  f4pCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
3770  f4pCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
3771  f4pCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
3772  f4pCorrelationFunctionsFlagsPro->SetStats(kFALSE);
3773  f4pCorrelationFunctionsFlagsPro->SetFillColor(kGray);
3774  f4pCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
3775  f4pCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f4pCorrelationFunctionsFlagsPro->Fill(0.5,0.); // TBI dummy at the moment
3777 
3778  if(!fFillCorrelationFunctions){return;} // TBI is this safe? It is not, because now if I want to fill 3-p, I always have to fill also 2-p
3779 
3780  const Int_t nParticleSpecies = 5; // Supported at the moment: 0=e,1=mu,2=pi,3=K,4=p
3781  TString sParticles[2*nParticleSpecies] = {"e^{+}","#mu^{+}","#pi^{+}","K^{+}","p^{+}","e^{-}","#mu^{-}","#pi^{-}","K^{-}","p^{-}"};
3782 
3783  // b) Book all 2p correlation functions:
3784  // Remark 0: First particle in the pair is always the one with positive charge.
3785  // Remark 1: Diagonal elements are particle/antiparticle pairs.
3786  for(Int_t pid1=0;pid1<2*nParticleSpecies;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3787  {
3788  for(Int_t pid2=pid1;pid2<2*nParticleSpecies;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3789  {
3790  // Correlation functions:
3791  fCorrelationFunctions[pid1][pid2] = new TH1F(Form("fCorrelationFunctions[%d][%d]",pid1,pid2),Form("fCorrelationFunctions[%d][%d] = (%s,%s)",pid1,pid2,sParticles[pid1].Data(),sParticles[pid2].Data()),10000,0.,10.); // TBI rethink the boundaries and nbins
3792  fCorrelationFunctions[pid1][pid2]->SetStats(kTRUE);
3793  fCorrelationFunctions[pid1][pid2]->SetFillColor(kBlue-10);
3794  fCorrelationFunctions[pid1][pid2]->SetXTitle("Q_{2}");
3795  fCorrelationFunctions[pid1][pid2]->SetYTitle("C(Q_{2})");
3797  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3798  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3799 
3800  // c) Book TExMap *fCorrelationFunctionsIndices:
3801  // Chosen convention: [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 5=e,6=mu,7=pi,8=K,9=p]
3802  fCorrelationFunctionsIndices = new TExMap();
3803  fCorrelationFunctionsIndices->Add(11,0);
3804  fCorrelationFunctionsIndices->Add(13,1);
3805  fCorrelationFunctionsIndices->Add(211,2);
3806  fCorrelationFunctionsIndices->Add(321,3);
3807  fCorrelationFunctionsIndices->Add(2212,4);
3808  fCorrelationFunctionsIndices->Add(-11,5);
3809  fCorrelationFunctionsIndices->Add(-13,6);
3810  fCorrelationFunctionsIndices->Add(-211,7);
3811  fCorrelationFunctionsIndices->Add(-321,8);
3812  fCorrelationFunctionsIndices->Add(-2212,9);
3813 
3814  // d) Book all 3p correlation functions:
3816  {
3817  for(Int_t pid1=0;pid1<2*nParticleSpecies;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3818  {
3819  for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3820  {
3821  for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3822  {
3823  // 3-p correlation functions:
3824  f3pCorrelationFunctions[pid1][pid2][pid3] = new TH1F(Form("f3pCorrelationFunctions[%d][%d][%d]",pid1,pid2,pid3),Form("f3pCorrelationFunctions[%d][%d][%d] = (%s,%s,%s)",pid1,pid2,pid3,sParticles[pid1].Data(),sParticles[pid2].Data(),sParticles[pid3].Data()),10000,0.,10.); // TBI rethink the boundaries and nbins
3825  f3pCorrelationFunctions[pid1][pid2][pid3]->SetStats(kTRUE);
3826  f3pCorrelationFunctions[pid1][pid2][pid3]->SetFillColor(kBlue-10);
3827  f3pCorrelationFunctions[pid1][pid2][pid3]->SetXTitle("Q_{3}");
3828  f3pCorrelationFunctions[pid1][pid2][pid3]->SetYTitle("C(Q_{3})");
3829  fCorrelationFunctionsSublist[1]->Add(f3pCorrelationFunctions[pid1][pid2][pid3]);
3830  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3831  } // for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3832  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3833  } // if(fFill3pCorrelationFunctions)
3834 
3835  // e) Book all 4p correlation functions:
3837  {
3838  for(Int_t pid1=0;pid1<2*nParticleSpecies;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3839  {
3840  for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3841  {
3842  for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3843  {
3844  for(Int_t pid4=0;pid4<2*nParticleSpecies;pid4++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3845  {
3846  // 4-p correlation functions:
3847  f4pCorrelationFunctions[pid1][pid2][pid3][pid4] = new TH1F(Form("f4pCorrelationFunctions[%d][%d][%d][%d]",pid1,pid2,pid3,pid4),Form("f4pCorrelationFunctions[%d][%d][%d][%d] = (%s,%s,%s,%s)",pid1,pid2,pid3,pid4,sParticles[pid1].Data(),sParticles[pid2].Data(),sParticles[pid3].Data(),sParticles[pid4].Data()),10000,0.,10.); // TBI rethink the boundaries and nbins
3848  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->SetStats(kTRUE);
3849  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->SetFillColor(kBlue-10);
3850  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->SetXTitle("Q_{4}");
3851  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->SetYTitle("C(Q_{4})");
3852  fCorrelationFunctionsSublist[2]->Add(f4pCorrelationFunctions[pid1][pid2][pid3][pid4]);
3853  } // for(Int_t pid4=0;pid4<2*nParticleSpecies;pid4++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3854  } // for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3855  } // for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3856  } // for(Int_t pid1=0;pid1<2*nParticleSpecies;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3857  } // if(fFill4pCorrelationFunctions)
3858 
3859 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForCorrelationFunctions()
3860 
3861 //=======================================================================================================================
3862 
3864 {
3865  // Book all the stuff for background objects.
3866 
3867  // a) Book the profile holding all the flags for background objects;
3868  // b) Book all histograms for 2p background;
3869  // c) Book all histograms for 3p background;
3870  // d) Book all histograms for 4p background;
3871  // e) Book buffer objects.
3872 
3873  // a) Book the profile holding all the flags for background objects:
3874  fBackgroundFlagsPro = new TProfile("fBackgroundFlagsPro","Flags and settings for background histograms",4,0,4);
3875  fBackgroundFlagsPro->SetTickLength(-0.01,"Y");
3876  fBackgroundFlagsPro->SetMarkerStyle(25);
3877  fBackgroundFlagsPro->SetLabelSize(0.04);
3878  fBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
3879  fBackgroundFlagsPro->SetStats(kFALSE);
3880  fBackgroundFlagsPro->SetFillColor(kGray);
3881  fBackgroundFlagsPro->SetLineColor(kBlack);
3882  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"fBackgroundOption"); fBackgroundFlagsPro->Fill(0.5,fBackgroundOption);
3883  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(2,"fEstimate2pBackground"); fBackgroundFlagsPro->Fill(1.5,fEstimate2pBackground);
3884  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(3,"fEstimate3pBackground"); fBackgroundFlagsPro->Fill(2.5,fEstimate3pBackground);
3885  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(4,"fEstimate4pBackground"); fBackgroundFlagsPro->Fill(3.5,fEstimate4pBackground);
3887  // 2p:
3888  f2pBackgroundFlagsPro = new TProfile("f2pBackgroundFlagsPro","Flags and settings for 2p background histograms",1,0,1);
3889  f2pBackgroundFlagsPro->SetTickLength(-0.01,"Y");
3890  f2pBackgroundFlagsPro->SetMarkerStyle(25);
3891  f2pBackgroundFlagsPro->SetLabelSize(0.04);
3892  f2pBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
3893  f2pBackgroundFlagsPro->SetStats(kFALSE);
3894  f2pBackgroundFlagsPro->SetFillColor(kGray);
3895  f2pBackgroundFlagsPro->SetLineColor(kBlack);
3896  f2pBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f2pBackgroundFlagsPro->Fill(0.5,0); // TBI dummy at the moment
3898  // 3p:
3899  f3pBackgroundFlagsPro = new TProfile("f3pBackgroundFlagsPro","Flags and settings for 3p background histograms",1,0,1);
3900  f3pBackgroundFlagsPro->SetTickLength(-0.01,"Y");
3901  f3pBackgroundFlagsPro->SetMarkerStyle(25);
3902  f3pBackgroundFlagsPro->SetLabelSize(0.04);
3903  f3pBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
3904  f3pBackgroundFlagsPro->SetStats(kFALSE);
3905  f3pBackgroundFlagsPro->SetFillColor(kGray);
3906  f3pBackgroundFlagsPro->SetLineColor(kBlack);
3907  f3pBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f3pBackgroundFlagsPro->Fill(0.5,0); // TBI dummy at the moment
3909  // 4p:
3910  f4pBackgroundFlagsPro = new TProfile("f4pBackgroundFlagsPro","Flags and settings for 4p background histograms",1,0,1);
3911  f4pBackgroundFlagsPro->SetTickLength(-0.01,"Y");
3912  f4pBackgroundFlagsPro->SetMarkerStyle(25);
3913  f4pBackgroundFlagsPro->SetLabelSize(0.04);
3914  f4pBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
3915  f4pBackgroundFlagsPro->SetStats(kFALSE);
3916  f4pBackgroundFlagsPro->SetFillColor(kGray);
3917  f4pBackgroundFlagsPro->SetLineColor(kBlack);
3918  f4pBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f4pBackgroundFlagsPro->Fill(0.5,0); // TBI dummy at the moment
3920 
3921  //if(!fEstimateBackground){return;} // TBI is this safe?
3922 
3923  const Int_t nParticleSpecies = 5; // Supported at the moment: 0=e,1=mu,2=pi,3=K,4=p
3924  TString sParticles[2*nParticleSpecies] = {"e^{+}","#mu^{+}","#pi^{+}","K^{+}","p^{+}","e^{-}","#mu^{-}","#pi^{-}","K^{-}","p^{-}"};
3925 
3926  // b) Book all histograms for 2p background:
3928  {
3929  // Remark 0: First particle in the pair is always the one with positive charge.
3930  // Remark 1: Diagonal elements are particle/antiparticle pairs.
3931  for(Int_t pid1=0;pid1<2*nParticleSpecies;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3932  {
3933  for(Int_t pid2=pid1;pid2<2*nParticleSpecies;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3934  {
3935  // Background:
3936  f2pBackground[pid1][pid2] = new TH1F(Form("f2pBackground[%d][%d]",pid1,pid2),Form("f2pBackground[%d][%d] = (%s,%s)",pid1,pid2,sParticles[pid1].Data(),sParticles[pid2].Data()),10000,0.,10.); // TBI rethink the boundaries and nbins
3937  f2pBackground[pid1][pid2]->SetStats(kTRUE);
3938  f2pBackground[pid1][pid2]->SetFillColor(kRed-10);
3939  f2pBackground[pid1][pid2]->SetLineColor(kRed);
3940  f2pBackground[pid1][pid2]->SetXTitle("k");
3941  f2pBackground[pid1][pid2]->SetYTitle("B(k)");
3942  fBackgroundSublist[0]->Add(f2pBackground[pid1][pid2]);
3943  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3944  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3945  } // if(fEstimate2pBackground)
3946 
3947  // c) Book all histograms for 3p background:
3949  {
3950  for(Int_t pid1=0;pid1<2*nParticleSpecies;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3951  {
3952  for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3953  {
3954  for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3955  {
3956  // Background:
3957  f3pBackground[pid1][pid2][pid3] = new TH1F(Form("f3pBackground[%d][%d][%d]",pid1,pid2,pid3),Form("f3pBackground[%d][%d][%d] = (%s,%s,%s)",pid1,pid2,pid3,sParticles[pid1].Data(),sParticles[pid2].Data(),sParticles[pid3].Data()),10000,0.,10.); // TBI rethink the boundaries and nbins
3958  f3pBackground[pid1][pid2][pid3]->SetStats(kTRUE);
3959  f3pBackground[pid1][pid2][pid3]->SetFillColor(kRed-10);
3960  f3pBackground[pid1][pid2][pid3]->SetLineColor(kRed);
3961  f3pBackground[pid1][pid2][pid3]->SetXTitle("Q_{3}");
3962  f3pBackground[pid1][pid2][pid3]->SetYTitle("B(Q_{3})");
3963  fBackgroundSublist[1]->Add(f3pBackground[pid1][pid2][pid3]);
3964  } // for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3965  } // for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3966  } // for(Int_t pid=0;pid<2*nParticleSpecies;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3967  } // if(fEstimate2pBackground)
3968 
3969  // d) Book all histograms for 4p background:
3971  {
3972  for(Int_t pid1=0;pid1<2*nParticleSpecies;pid1++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3973  {
3974  for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3975  {
3976  for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3977  {
3978  for(Int_t pid4=0;pid4<2*nParticleSpecies;pid4++) // [particle(+q): 0=e,1=mu,2=pi,3=K,4=p, anti-particle(-q): 0=e,1=mu,2=pi,3=K,4=p]
3979  {
3980  // Background:
3981  f4pBackground[pid1][pid2][pid3][pid4] = new TH1F(Form("f4pBackground[%d][%d][%d][%d]",pid1,pid2,pid3,pid4),Form("f4pBackground[%d][%d][%d][%d] = (%s,%s,%s,%s)",pid1,pid2,pid3,pid4,sParticles[pid1].Data(),sParticles[pid2].Data(),sParticles[pid3].Data(),sParticles[pid4].Data()),10000,0.,10.); // TBI rethink the boundaries and nbins
3982  f4pBackground[pid1][pid2][pid3][pid4]->SetStats(kTRUE);
3983  f4pBackground[pid1][pid2][pid3][pid4]->SetFillColor(kRed-10);
3984  f4pBackground[pid1][pid2][pid3][pid4]->SetLineColor(kRed);
3985  f4pBackground[pid1][pid2][pid3][pid4]->SetXTitle("Q_{4}");
3986  f4pBackground[pid1][pid2][pid3][pid4]->SetYTitle("B(Q_{4})");
3987  fBackgroundSublist[2]->Add(f4pBackground[pid1][pid2][pid3][pid4]);
3988  } // for(Int_t pid4=0;pid4<2*nParticleSpecies;pid4++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3989  } // for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3990  } // for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3991  } // for(Int_t pid=0;pid<2*nParticleSpecies;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3992  } // if(fEstimate2pBackground)
3993 
3994  // e) Book buffer objects:
3995  switch(fBackgroundOption)
3996  {
3997  case 0: // shifting (see documentation in EstimateBackground(...))
3998  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
3999  {
4000  fMixedEvents0[me] = new TClonesArray("AliAODTrack",10000);
4001  }
4002  break; // case 0: // shifting
4003 
4004  case 1: // permutations + binning in vertex z-position
4005  for(Int_t vzr=0;vzr<10;vzr++) // vertez-z range
4006  {
4007  for(Int_t me=0;me<fMaxBufferSize1;me++) // buffer for mixed events
4008  {
4009  fMixedEvents1[vzr][me] = new TClonesArray("AliAODTrack",10000);
4010  fGlobalTracksAOD1[vzr][me] = new TExMap();
4011  }
4012  }
4013  break; // permutations + binning in vertex z-position
4014 
4015  /*
4016  default:
4017  cout<<Form("And the fatal 'fBackgroundOption' value is... %d. Congratulations!!",fBackgroundOption)<<endl;
4018  Fatal(sMethodName.Data(),"switch(fBackgroundOption)");
4019  */
4020  } // switch(fBackgroundOption)
4021 
4022 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForBackground()
4023 
4024 //=======================================================================================================================
4025 
4027 {
4028  // Book all the stuff for buffers.
4029 
4030  // a) Book the profile holding all the flags for buffers;
4031  // b) Book TClonesArray *fChargedParticlesCA[2][10][10000];
4032  // c) Book TExMap *fChargedParticlesEM[10];
4033 
4034  // a) Book the profile holding all the flags for buffers:
4035  fBuffersFlagsPro = new TProfile("fBuffersFlagsPro","Flags and settings for buffers",2,0,2);
4036  fBuffersFlagsPro->SetTickLength(-0.01,"Y");
4037  fBuffersFlagsPro->SetMarkerStyle(25);
4038  fBuffersFlagsPro->SetLabelSize(0.04);
4039  fBuffersFlagsPro->SetLabelOffset(0.02,"Y");
4040  fBuffersFlagsPro->SetStats(kFALSE);
4041  fBuffersFlagsPro->SetFillColor(kGray);
4042  fBuffersFlagsPro->SetLineColor(kBlack);
4043  fBuffersFlagsPro->GetXaxis()->SetBinLabel(1,"fFillBuffers"); fBuffersFlagsPro->Fill(0.5,(Int_t)fFillBuffers);
4044  fBuffersFlagsPro->GetXaxis()->SetBinLabel(2,"fMaxBuffer"); fBuffersFlagsPro->Fill(1.5,fMaxBuffer);
4046 
4047  if(!fFillBuffers){return;} // TBI is this safe? Well, let's hope so for the time being...
4048 
4049  // b) Book TClonesArray *fChargedParticlesCA[2][10][10000]:
4050  for(Int_t am=0;am<1+(Int_t)fProcessBothKineAndReco;am++) // analysis method [0=AOD||ESD,1=MC]
4051  {
4052  for(Int_t e=0;e<fMaxBuffer;e++) // number of events to buffer, max = apparently 10
4053  {
4054  for(Int_t p=0;p<10000;p++) // charged particles, max = apparently 10000
4055  {
4056  fChargedParticlesCA[am][e][p] = NULL; // TBI not sure if I need this really here, this way...
4057  // Example: fChargedParticlesCA[0][4][6] is 6th charged particles in 5th buffered event, in AOD||ESD analysis
4058  } // for(Int_t p=0;p<10000;p++) // charged particles, max = apparently 10000
4059  } // for(Int_t e=0;e<10;e++) // number of events to buffer, max = apparently 10
4060  } // for(Int_t am=0;am<2;am++) // analysis method [0=AOD||ESD,1=MC]
4061 
4062  // c) Book TExMap *fChargedParticlesEM[10]:
4063  for(Int_t e=0;e<fMaxBuffer;e++) // number of events to buffer, max = apparently 10
4064  {
4065  fChargedParticlesEM[e] = new TExMap();
4066  } // for(Int_t e=0;e<fMaxBuffer;e++) // number of events to buffer, max = apparently 10
4067 
4068 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForBuffers()
4069 
4070 //=======================================================================================================================
4071 
4073 {
4074  // Book all the stuff for QA.
4075 
4076  // a) Book the profile holding all the flags for QA;
4077  // b) Book all objects for "QA events";
4078  // c) Book all objects for "QA particles".
4079 
4080  // a) Book the profile holding all the flags for QA:
4081  fQAFlagsPro = new TProfile("fQAFlagsPro","Flags and settings for QA",4,0.,4.);
4082  fQAFlagsPro->SetTickLength(-0.01,"Y");
4083  fQAFlagsPro->SetMarkerStyle(25);
4084  fQAFlagsPro->SetLabelSize(0.04);
4085  fQAFlagsPro->SetLabelOffset(0.02,"Y");
4086  fQAFlagsPro->SetStats(kFALSE);
4087  fQAFlagsPro->SetFillColor(kGray);
4088  fQAFlagsPro->SetLineColor(kBlack);
4089  fQAFlagsPro->GetXaxis()->SetBinLabel(1,"fFillQA"); fQAFlagsPro->Fill(0.5,(Int_t)fFillQA);
4090  fQAFlagsPro->GetXaxis()->SetBinLabel(2,"fBailOutAfterQA"); fQAFlagsPro->Fill(1.5,(Int_t)fBailOutAfterQA);
4091  fQAFlagsPro->GetXaxis()->SetBinLabel(3,"fFillQAEvents"); fQAFlagsPro->Fill(2.5,(Int_t)fFillQAEvents);
4092  fQAFlagsPro->GetXaxis()->SetBinLabel(4,"fFillQAParticles"); fQAFlagsPro->Fill(3.5,(Int_t)fFillQAParticles);
4093  fQAList->Add(fQAFlagsPro);
4094 
4095  if(!fFillQA){return;}
4096 
4097  // b) Book all objects for "QA events":
4098  if(fFillQAEvents)
4099  {
4100  // ...
4101  } // if(fFillQAEvents)
4102 
4103  // c) Book all objects for "QA particles":
4104  if(fFillQAParticles)
4105  {
4106  // c0) Book fQAFilterBitScan:
4107  Int_t nFilterBits = 14;
4108  fQAFilterBitScan = new TH1I("fQAFilterBitScan","fQAFilterBitScan",nFilterBits,0,nFilterBits);
4109  fQAFilterBitScan->SetStats(kFALSE);
4110  fQAFilterBitScan->SetFillColor(kBlue-10);
4111  fQAFilterBitScan->SetXTitle("FilterBit");
4112  fQAFilterBitScan->SetYTitle("# particles");
4113  for(Int_t fb=0;fb<nFilterBits;fb++) // 'fb' is a 'left shifting opetator', i.e. 1<<fb = filterbit
4114  {
4115  fQAFilterBitScan->GetXaxis()->SetBinLabel(fb+1,Form("%d",1<<fb));
4116  }
4118 
4119  // c1) Book fQAIDvsFilterBit:
4120  Int_t nIDsMax = 10000;
4121  fQAIDvsFilterBit = new TH2I("fQAIDvsFilterBit","fQAIDvsFilterBit",nFilterBits,0,nFilterBits,2*nIDsMax,-nIDsMax,nIDsMax);
4122  //fQAIDvsFilterBit->SetStats(kFALSE);
4123  //fQAIDvsFilterBit->SetFillColor(kBlue-10);
4124  fQAIDvsFilterBit->SetXTitle("FilterBit");
4125  fQAIDvsFilterBit->SetYTitle("atrack->GetID()");
4126  for(Int_t fb=0;fb<nFilterBits;fb++) // 'fb' is a 'left shifting opetator', i.e. 1<<fb = filterbit
4127  {
4128  fQAIDvsFilterBit->GetXaxis()->SetBinLabel(fb+1,Form("%d",1<<fb));
4129  }
4131 
4132  // c2) Book fQAParticleHist[2][10][10];
4133  TString sBeforeAfter[2] = {"before","after"};
4134  TString sDistributions[10] = {"p_{T}","#eta","","","","","","","",""}; // Only the ones with non-empty entries are booked and filled, see below if(sDistributions[di].EqualTo("")){continue;}
4135  TString sCuts[10] = {"PID","","","","","","","","",""}; // Only the ones with non-empty entries are booked and filled, see below if(sCuts[ci].EqualTo("")){continue;}
4136  for(Int_t ba=0;ba<2;ba++) // 0="before rain",1="after rain"
4137  {
4138  for(Int_t di=0;di<10;di++) // number of distinct distributions, maximum is 10 (at the moment, as it seems...)
4139  {
4140  if(sDistributions[di].EqualTo("")){continue;}
4141  for(Int_t ci=0;ci<10;ci++) // number of distinct track cuts, maximum is 10 (at the moment, as it seems...)
4142  {
4143  if(sCuts[ci].EqualTo("")){continue;}
4144  fQAParticleHist[ba][di][ci] = new TH1F(Form("fQAParticleHist[%d][%d][%d]",ba,di,ci),Form("fQAParticleHist[%d][%d][%d] = (%s,%s,%s)",ba,di,ci,sBeforeAfter[ba].Data(),sDistributions[di].Data(),sCuts[ci].Data() ),10000,0.,10.); // TBI rethink the boundaries and nbins
4145  fQAParticleHist[ba][di][ci]->SetStats(kFALSE);
4146  fQAParticleHist[ba][di][ci]->SetFillColor(kBlue-10);
4147  fQAParticleHist[ba][di][ci]->SetXTitle("TBI");
4148  fQAParticleHist[ba][di][ci]->SetYTitle("TBI");
4149  fQAParticlesList->Add(fQAParticleHist[ba][di][ci]);
4150  } // for(Int_t ci=0;ci<10;ci++) // number of distinct track cuts, maximum is 10 (at the moment)
4151  } // for(Int_t di=0;di<10;di++) // number of distinct distributions, maximum is 10 (at the moment)
4152  } // for(Int_t ba=0;ba<2;ba++) // 0="before rain",1="after rain"
4153  } // if(fFillQAParticles)
4154 
4155 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForQA()
4156 
4157 //=======================================================================================================================
4158 
4160 {
4161  // Book all objects for global track cuts (applied only on "normal" global tracks in AOD).
4162 
4163  // a) Book the profile holding all the flags;
4164  // b) ...
4165 
4166  // a) Book the profile holding all the flags:
4167  fGlobalTrackCutsFlagsPro = new TProfile("fGlobalTrackCutsFlagsPro","Flags and settings for global track cuts",7,0.,7.);
4168  fGlobalTrackCutsFlagsPro->SetTickLength(-0.01,"Y");
4169  fGlobalTrackCutsFlagsPro->SetMarkerStyle(25);
4170  fGlobalTrackCutsFlagsPro->SetLabelSize(0.04);
4171  fGlobalTrackCutsFlagsPro->SetLabelOffset(0.02,"Y");
4172  fGlobalTrackCutsFlagsPro->SetStats(kFALSE);
4173  fGlobalTrackCutsFlagsPro->SetFillColor(kGray);
4174  fGlobalTrackCutsFlagsPro->SetLineColor(kBlack);
4175  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(1,"fApplyGlobalTrackCuts"); fGlobalTrackCutsFlagsPro->Fill(0.5,(Int_t)fApplyGlobalTrackCuts);
4176  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(2,"fPtRange[0]"); fGlobalTrackCutsFlagsPro->Fill(1.5,fPtRange[0]);
4177  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(3,"fPtRange[1]"); fGlobalTrackCutsFlagsPro->Fill(2.5,fPtRange[1]);
4178  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(4,"fEtaRange[0]"); fGlobalTrackCutsFlagsPro->Fill(3.5,fEtaRange[0]);
4179  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(5,"fEtaRange[1]"); fGlobalTrackCutsFlagsPro->Fill(4.5,fEtaRange[1]);
4180  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(6,"fPhiRange[0]"); fGlobalTrackCutsFlagsPro->Fill(5.5,fPhiRange[0]);
4181  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(7,"fPhiRange[1]"); fGlobalTrackCutsFlagsPro->Fill(6.5,fPhiRange[1]);
4183 
4184 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForGlobalTrackCuts()
4185 
4186 //=======================================================================================================================
4187 
4189 {
4190  // Filter out unique global tracks in AOD and store them in fGlobalTracksAOD[index].
4191 
4192  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
4193  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
4194  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
4195  // Global constrained tracks have the same negative ID as the TPC-only tracks, both associated with the same "normal global" tracks. E.g. we can have
4196  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
4197  // 1: 0: 2.073798 -0.503640 2.935432
4198  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
4199  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
4200  // In fact, this is important, otherwise there is double or even triple counting in some cases.
4201  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
4202  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
4203  // b) All of them apparently have positive ID, i.e. these are global tracks
4204  // c) Clearly, we cannot use TestFilterBit() on them
4205  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
4206  // Remark 3: There is a performance penalty when fGlobalTracksAOD[1] and fGlobalTracksAOD[2] needed for mixed events are calculated.
4207  // Yes, I can get them directly from fGlobalTracksAOD[0], without calling this method for them again. TBI today
4208 
4209  // a) Insanity checks;
4210  // b) Determine the map.
4211 
4212  // a) Insanity checks:
4213  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD, Int_t index)";
4214  if(!fGlobalTracksAOD[index]){Fatal(sMethodName.Data(),"fGlobalTracksAOD[%d]",index);}
4215  if(0 != fGlobalTracksAOD[index]->GetSize()){fGlobalTracksAOD[index]->Delete();} // yes, this method determines mapping from scratch each time
4216 
4217  // b) Determine the map:
4218  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4219  {
4220  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
4221  if(aodTrack)
4222  {
4223  Int_t id = aodTrack->GetID();
4224  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
4225  if(id>=0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this, it seems that id>=0 is just enough, the second constraint is most likely just an overkill
4226  {
4227  fGlobalTracksAOD[index]->Add(id,iTrack); // "key" = id, "value" = iTrack
4228  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
4229  } // if(aodTrack)
4230  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4231 
4232 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD)
4233 
4234 //=======================================================================================================================
4235 
4237 {
4238  // Filter out unique global tracks in AOD and store them in fGlobalTracksAODTEST[index].
4239 
4240  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
4241  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
4242  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
4243  // Global constrained tracks have the same negative ID as the TPC-only tracks, both associated with the same "normal global" tracks. E.g. we can have
4244  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
4245  // 1: 0: 2.073798 -0.503640 2.935432
4246  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
4247  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
4248  // In fact, this is important, otherwise there is double or even triple counting in some cases.
4249  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
4250  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
4251  // b) All of them apparently have positive ID, i.e. these are global tracks
4252  // c) Clearly, we cannot use TestFilterBit() on them
4253  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
4254  // Remark 3: There is a performance penalty when fGlobalTracksAODTEST[1] and fGlobalTracksAODTEST[2] needed for mixed events are calculated.
4255  // Yes, I can get them directly from fGlobalTracksAODTEST[0], without calling this method for them again. TBI today
4256 
4257  // a) Insanity checks;
4258  // b) Determine the map.
4259 
4260  // a) Insanity checks:
4261  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAODTEST(AliAODEvent *aAOD, Int_t index)";
4262  if(!fGlobalTracksAODTEST[index]){Fatal(sMethodName.Data(),"fGlobalTracksAODTEST[%d]",index);}
4263  if(0 != fGlobalTracksAODTEST[index]->GetSize()){fGlobalTracksAODTEST[index]->Delete();} // yes, this method determines mapping from scratch each time
4264 
4265  // b) Determine the map:
4266  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4267  {
4268  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
4269  if(aodTrack)
4270  {
4271  Int_t id = aodTrack->GetID();
4272  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
4273  if(id>=0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this, it seems that id>=0 is just enough, the second constraint is most likely just an overkill
4274  {
4275  fGlobalTracksAODTEST[index]->Add(id,iTrack); // "key" = id, "value" = iTrack
4276  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
4277  } // if(aodTrack)
4278  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4279 
4280 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAODTEST(AliAODEvent *aAOD)
4281 
4282 //=======================================================================================================================
4283 
4285 {
4286  // Filter out unique global tracks in AOD and store them in fGlobalTracksAODHA[index].
4287 
4288  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
4289  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
4290  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
4291  // Global constrained tracks have the same negative ID as the TPC-only tracks, both associated with the same "normal global" tracks. E.g. we can have
4292  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
4293  // 1: 0: 2.073798 -0.503640 2.935432
4294  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
4295  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
4296  // In fact, this is important, otherwise there is double or even triple counting in some cases.
4297  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
4298  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
4299  // b) All of them apparently have positive ID, i.e. these are global tracks
4300  // c) Clearly, we cannot use TestFilterBit() on them
4301  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
4302  // Remark 3: There is a performance penalty when fGlobalTracksAODHA[1] and fGlobalTracksAODHA[2] needed for mixed events are calculated.
4303  // Yes, I can get them directly from fGlobalTracksAODHA[0], without calling this method for them again. TBI today
4304 
4305  // a) Insanity checks;
4306  // b) Determine the map.
4307 
4308  // a) Insanity checks:
4309  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAODHA(AliAODEvent *aAOD, Int_t index)";
4310  if(!fGlobalTracksAODHA[index]){Fatal(sMethodName.Data(),"fGlobalTracksAODHA[%d]",index);}
4311  if(0 != fGlobalTracksAODHA[index]->GetSize()){fGlobalTracksAODHA[index]->Delete();} // yes, this method determines mapping from scratch each time
4312 
4313  // b) Determine the map:
4314  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4315  {
4316  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
4317  if(aodTrack)
4318  {
4319  Int_t id = aodTrack->GetID();
4320  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
4321  if(id>=0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this, it seems that id>=0 is just enough, the second constraint is most likely just an overkill
4322  {
4323  fGlobalTracksAODHA[index]->Add(id,iTrack); // "key" = id, "value" = iTrack
4324  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
4325  } // if(aodTrack)
4326  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4327 
4328 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAODHA(AliAODEvent *aAOD, Int_t index)
4329 
4330 //=======================================================================================================================
4331 
4333 {
4334  // Filter out unique global tracks in AOD and store them in fGlobalTracksAOD1[10][5]. TBI hardwired "1" in an array, not in the function name, so this is a bit incosistent...
4335 
4336  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
4337  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
4338  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
4339  // Global constrained tracks have the same negative ID as the TPC-only tracks, both associated with the same "normal global" tracks. E.g. we can have
4340  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
4341  // 1: 0: 2.073798 -0.503640 2.935432
4342  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
4343  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
4344  // In fact, this is important, otherwise there is double or even triple counting in some cases.
4345  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
4346  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
4347  // b) All of them apparently have positive ID, i.e. these are global tracks
4348  // c) Clearly, we cannot use TestFilterBit() on them
4349  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
4350  // Remark 3: Used only for fBackgroundOption = 1, i.e. permutations + binning in vertex z-position
4351 
4352  // a) Insanity checks;
4353  // b) Determine the map.
4354 
4355  // a) Insanity checks:
4356  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD, Int_t indexX, Int_t indexY)";
4357  if(!fGlobalTracksAOD1[indexX][indexY]){Fatal(sMethodName.Data(),"fGlobalTracksAOD[%d][%d]",indexX,indexY);}
4358  if(0 != fGlobalTracksAOD1[indexX][indexY]->GetSize()){fGlobalTracksAOD1[indexX][indexY]->Delete();} // yes, this method determines mapping from scratch each time
4359 
4360  // b) Determine the map:
4361  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4362  {
4363  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
4364  if(aodTrack)
4365  {
4366  Int_t id = aodTrack->GetID();
4367  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
4368  if(id>=0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this, it seems that id>=0 is just enough, the second constraint is most likely just an overkill
4369  {
4370  fGlobalTracksAOD1[indexX][indexY]->Add(id,iTrack); // "key" = id, "value" = iTrack
4371  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
4372  } // if(aodTrack)
4373  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4374 
4375 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD, Int_t indexX, Int_t indexY)
4376 
4377 //=======================================================================================================================
4378 
4380 {
4381  // Check if this is event specified in a steering macro via the setter void SetWaitForSpecifiedEvent(UInt_t run, UShort_t bunchCross, UInt_t orbit, UInt_t period).
4382 
4383  if(run != fRun) return kFALSE;
4384  else if(bunchCross != fBunchCross) return kFALSE;
4385  else if(orbit != fOrbit) return kFALSE;
4386  else if(period != fPeriod) return kFALSE;
4387 
4388  return kTRUE;
4389 
4390 } // void AliAnalysisTaskMultiparticleFemtoscopy::SpecifiedEvent(UInt_t run, UShort_t bunchCross, UInt_t orbit, UInt_t period)
4391 
4392 //=======================================================================================================================
4393 
4395 {
4396  // Check if the track passes common global track cuts (irrespectively of PID).
4397 
4398  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesGlobalTrackCuts(AliAODTrack *gtrack)";
4399  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
4400 
4401  // To do: add data members and corresponding setters:
4402  // fTPCNclsMin, fTPCNclsMax
4403  // fTPCsignalNMin, fTPCsignalNMax
4404  if(gtrack->Pt()<fPtRange[0]) return kFALSE;
4405  if(gtrack->Pt()>=fPtRange[1]) return kFALSE;
4406  if(gtrack->Eta()<fEtaRange[0]) return kFALSE;
4407  if(gtrack->Eta()>=fEtaRange[1]) return kFALSE;
4408  if(gtrack->Phi()<fPhiRange[0]) return kFALSE;
4409  if(gtrack->Phi()>=fPhiRange[1]) return kFALSE;
4410  //if(gtrack->GetTPCNcls()<70) return kFALSE;
4411  //if(gtrack->GetTPCsignalN()<70) return kFALSE;
4412 
4413  if(fRejectFakeTracks && gtrack->GetLabel()<0) return kFALSE; // TBI is it more precise <=0 ?
4414 
4415  return kTRUE;
4416 
4417 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesGlobalTrackCuts(AliAODTrack *gtrack)
4418 
4419 //=======================================================================================================================
4420 
4422 {
4423  // Check if the track passes common analysis specific track (e.g. TPC-only) cuts.
4424  // Therefore we can cut independetly on global track parameters, and on TPC-only cut parameters.
4425 
4426  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODTrack *atrack)";
4427  if(!atrack){Fatal(sMethodName.Data(),"!atrack");}
4428 
4429  if(!atrack->TestFilterBit(128)) return kFALSE; // TPC-only TBI setter TBI#2 There might be some conflict with FillControlHistogramsNonIdentifiedParticlesFTSF
4430 
4431  // TBI: implement eventually separate cuts for 'atracks' and 'gtracks'
4432  if(atrack->Pt()<fPtRange[0]) return kFALSE;
4433  if(atrack->Pt()>=fPtRange[1]) return kFALSE;
4434  if(atrack->Eta()<fEtaRange[0]) return kFALSE;
4435  if(atrack->Eta()>=fEtaRange[1]) return kFALSE;
4436  if(atrack->Phi()<fPhiRange[0]) return kFALSE;
4437  if(atrack->Phi()>=fPhiRange[1]) return kFALSE;
4438  //if(gtrack->GetTPCNcls()<70) return kFALSE;
4439  //if(gtrack->GetTPCsignalN()<70) return kFALSE;
4440 
4441  if(fRejectFakeTracks && atrack->GetLabel()<0) return kFALSE; // TBI is it more precise <=0 ?
4442 
4443  return kTRUE;
4444 
4445 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODTrack *atrack)
4446 
4447 //=======================================================================================================================
4448 
4450 {
4451  // TBI not validated. this method applies only to MC, make it uniform with AOD.
4452 
4453  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODMCParticle *amcparticle)";
4454  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
4455 
4456  if(TMath::Abs(amcparticle->Charge())<1.e-10) return kFALSE; // skipping neutral particles TBI
4457  if(!amcparticle->IsPhysicalPrimary()) return kFALSE; // skipping secondaries TBI
4458  if(amcparticle->Pt()<0.2) return kFALSE;
4459  if(amcparticle->Pt()>=5.0) return kFALSE;
4460  if(amcparticle->Eta()<-0.8) return kFALSE;
4461  if(amcparticle->Eta()>=0.8) return kFALSE;
4462 
4463  return kTRUE;
4464 
4465 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODMCParticle *amcparticle)
4466 
4467 //=======================================================================================================================
4468 
4470 {
4471  // Check if the event passes common event cuts.
4472 
4473  // TBI: add support for fProcessBothKineAndReco
4474 
4475  // a) Determine Ali{MC,ESD,AOD}Event;
4476 
4477  // a) Determine Ali{MC,ESD,AOD}Event:
4478  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
4479  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
4480  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
4481 
4482  if(aMC && fProcessOnlyKine)
4483  {
4484  // TBI
4485  }
4486  else if(aESD)
4487  {
4488  // TBI
4489  }
4490  else if(aAOD) // shall work both for fProcessBothKineAndReco = kTRUE and fProcessOnlyReco = kTRUE
4491  {
4492  // a) Cuts on AliAODEvent:
4493  if(fRejectEventsWithoutPrimaryVertex && !aAOD->GetPrimaryVertex()) return kFALSE;
4494  if(TMath::Abs(aAOD->GetMagneticField())<fMinMagneticField) return kFALSE;
4496  {
4497  if(aAOD->GetNumberOfTracks() < fMinNumberOfTracks) return kFALSE;
4498  if(aAOD->GetNumberOfTracks() > fMaxNumberOfTracks) return kFALSE;
4499  }
4500  if(fCutOnNumberOfV0s)
4501  {
4502  if(aAOD->GetNumberOfV0s() < fMinNumberOfV0s) return kFALSE;
4503  if(aAOD->GetNumberOfV0s() > fMaxNumberOfV0s) return kFALSE;
4504  }
4506  {
4507  if(aAOD->GetNumberOfCascades() < fMinNumberOfCascades) return kFALSE;
4508  if(aAOD->GetNumberOfCascades() > fMaxNumberOfCascades) return kFALSE;
4509  }
4510 
4511  // b) Cuts on AliAODVertex:
4512  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
4513  if(fCutOnVertexX)
4514  {
4515  if(avtx->GetX() < fMinVertexX) return kFALSE;
4516  if(avtx->GetX() > fMaxVertexX) return kFALSE;
4517  }
4518  if(fCutOnVertexY)
4519  {
4520  if(avtx->GetY() < fMinVertexY) return kFALSE;
4521  if(avtx->GetY() > fMaxVertexY) return kFALSE;
4522  }
4523  if(fCutOnVertexZ)
4524  {
4525  if(avtx->GetZ() < fMinVertexZ) return kFALSE;
4526  if(avtx->GetZ() > fMaxVertexZ) return kFALSE;
4527  }
4529  {
4530  if(avtx->GetNContributors() < fMinNContributors) return kFALSE;
4531  if(avtx->GetNContributors() > fMaxNContributors) return kFALSE;
4532  }
4533  } // else if(aAOD)
4534 
4535  return kTRUE;
4536 
4537 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonEventCuts(AliVEvent *ave)
4538 
4539 //=======================================================================================================================
4540 
4542 {
4543  // Check if the event passes mixed event cuts.
4544 
4545  // Remark: PassesCommonEventCuts() is already checked beforehand in UserExec(), therefore here only some new cuts characteristic only for mixed events shall be implemented.
4546 
4547  // a) Determine Ali{MC,ESD,AOD}Event;
4548 
4549  // a) Determine Ali{MC,ESD,AOD}Event:
4550  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
4551  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
4552  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
4553 
4554  if(aMC)
4555  {
4556  // TBI
4557  }
4558  else if(aESD)
4559  {
4560  // TBI
4561  }
4562  else if(aAOD)
4563  {
4564 
4565  // TBI
4566  /* Check whether the cuts below are already being applied in PassesCommonEventCuts()
4567  if(!aAOD->GetPrimaryVertex()) return kFALSE;
4568  if(TMath::Abs(aAOD->GetMagneticField())<0.001) return kFALSE;
4569  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
4570  //if(TMath::Abs(avtx->GetX())>10.0) return kFALSE;
4571  //if(TMath::Abs(avtx->GetY())>10.0) return kFALSE;
4572  if(TMath::Abs(avtx->GetZ())>10.0) return kFALSE; // TBI setter
4573  if(avtx->GetNContributors()<=2) return kFALSE; // TBI setter
4574  */
4575  }
4576 
4577  return kTRUE;
4578 
4579 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesMixedEventCuts(AliVEvent *ave)
4580 
4581 //=======================================================================================================================
4582 
4584 {
4585  // Calculate correlation functions.
4586 
4587  // a) Insanity checks;
4588  // b) Two nested loops to calculate C(k), just an example; TBI
4589  // c) Three nested loops to calculate C(Q3), just an example; TBI
4590  // d) Four nested loops to calculate C(Q4), just an example. TBI
4591 
4592  // a) Insanity checks:
4593  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::CalculateCorrelationFunctions(AliAODEvent *aAOD)";
4594  if(!aAOD){Fatal(sMethodName.Data(),"!aAOD");}
4595  if(0 == fGlobalTracksAOD[0]->GetSize()){Fatal(sMethodName.Data(),"0 == fGlobalTracksAOD[0]->GetSize()");} // this case shall be already treated in UserExec
4596 
4597  // b) Two nested loops to calculate C(k), just an example:
4598  Int_t nTracks = aAOD->GetNumberOfTracks();
4599  for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
4600  {
4601  AliAODTrack *atrack1 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack1));
4602  // TBI Temporary track insanity checks:
4603  if(!atrack1){Fatal(sMethodName.Data(),"!atrack1");} // TBI keep this for some time, eventually just continue
4604  if(atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4605  if(atrack1->TestFilterBit(128) && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->TestFiletrBit(128) && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4606  if(!PassesCommonTrackCuts(atrack1)){continue;} // TBI re-think
4607  // Corresponding AOD global track:
4608  Int_t id1 = atrack1->GetID();
4609  AliAODTrack *gtrack1 = dynamic_cast<AliAODTrack*>(id1>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id1)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id1+1))));
4610  if(!gtrack1){Fatal(sMethodName.Data(),"!gtrack1");} // TBI keep this for some time, eventually just continue
4611  Int_t gid1 = (id1>=0 ? id1 : -(id1+1)); // ID of corresponding global track
4612  // Common track selection criteria for all "normal" global tracks:
4613  if(!PassesGlobalTrackCuts(gtrack1)){continue;}
4614 
4615  // Loop over the 2nd particle:
4616  for(Int_t iTrack2=iTrack1+1;iTrack2<nTracks;iTrack2++)
4617  {
4618  AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack2));
4619  // TBI Temporary track insanity checks:
4620  if(!atrack2){Fatal(sMethodName.Data(),"!atrack2");} // TBI keep this for some time, eventually just continue
4621  if(atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4622  if(atrack2->TestFilterBit(128) && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->TestFiletrBit(128) && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4623  if(!PassesCommonTrackCuts(atrack2)){continue;} // TBI re-think
4624  // Corresponding AOD global track:
4625  Int_t id2 = atrack2->GetID();
4626  AliAODTrack *gtrack2 = dynamic_cast<AliAODTrack*>(id2>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id2)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id2+1))));
4627  if(!gtrack2){Fatal(sMethodName.Data(),"!gtrack2");} // TBI keep this for some time, eventually just continue
4628  Int_t gid2 = (id2>=0 ? id2 : -(id2+1)); // ID of corresponding global track
4629  if(gid1==gid2){continue;} // Eliminate self-correlations:
4630 
4631  // Common track selection criteria for all "normal" global tracks:
4632  if(!PassesGlobalTrackCuts(gtrack2)){continue;}
4633 
4634  // Okay, so we have two tracks, let's decide whether kinemetics from 'atrack' or 'gtrack' is taken, let's check PID, and fill the correlation functions:
4635  AliAODTrack *agtrack1 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4636  AliAODTrack *agtrack2 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4638  {
4639  agtrack1 = gtrack1;
4640  agtrack2 = gtrack2;
4641  }
4642  else
4643  {
4644  agtrack1 = atrack1;
4645  agtrack2 = atrack2;
4646  } // if(fFillControlHistogramsWithGlobalTrackInfo)
4647  if(!agtrack1){Fatal(sMethodName.Data(),"!agtrack1");}
4648  if(!agtrack2){Fatal(sMethodName.Data(),"!agtrack2");}
4649 
4650  // 1.) Same particle species:
4651 
4652  // a) pion-pion:
4653  // a1) pi+pi+ [2][2]:
4654  if(Pion(gtrack1,1,kTRUE) && Pion(gtrack2,1,kTRUE))
4655  {
4656  fCorrelationFunctions[2][2]->Fill(Q2(agtrack1,agtrack2));
4657  }
4658  // a2) pi-pi- [7][7]:
4659  if(Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,-1,kTRUE))
4660  {
4661  fCorrelationFunctions[7][7]->Fill(Q2(agtrack1,agtrack2));
4662  }
4663  // a3) pi+pi- || pi-pi+ [2][7]:
4664  if((Pion(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE)) || (Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,1,kTRUE)))
4665  {
4666  fCorrelationFunctions[2][7]->Fill(Q2(agtrack1,agtrack2));
4667  }
4668 
4669  // b) kaon-kaon:
4670  // b1) K+K+ [3][3]:
4671  if(Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE))
4672  {
4673  fCorrelationFunctions[3][3]->Fill(Q2(agtrack1,agtrack2));
4674  }
4675  // b2) K-K- [8][8]:
4676  if(Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
4677  {
4678  fCorrelationFunctions[8][8]->Fill(Q2(agtrack1,agtrack2));
4679  }
4680  // b3) K+K- || K-K+ [3][8]:
4681  if((Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE)) || (Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,1,kTRUE)))
4682  {
4683  fCorrelationFunctions[3][8]->Fill(Q2(agtrack1,agtrack2));
4684  }
4685 
4686  // c) proton-proton:
4687  // c1) p+p+ [4][4]:
4688  if(Proton(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
4689  {
4690  fCorrelationFunctions[4][4]->Fill(Q2(agtrack1,agtrack2));
4691  }
4692  // c2) p-p- [9][9]:
4693  if(Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4694  {
4695  fCorrelationFunctions[9][9]->Fill(Q2(agtrack1,agtrack2));
4696  }
4697  // c3) p+p- || p-p+ [4][9]:
4698  if((Proton(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE)) || (Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,1,kTRUE)))
4699  {
4700  fCorrelationFunctions[4][9]->Fill(Q2(agtrack1,agtrack2));
4701  }
4702 
4703  // 2.) Mixed particle species:
4704  // a) pion-kaon
4705  // a1) pi+K+ [2][3]:
4706  if(Pion(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE))
4707  {
4708  fCorrelationFunctions[2][3]->Fill(Q2(agtrack1,agtrack2));
4709  }
4710  // a2) pi+K- [2][8]
4711  if(Pion(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
4712  {
4713  fCorrelationFunctions[2][8]->Fill(Q2(agtrack1,agtrack2));
4714  }
4715  // a3) K+pi- [3][7]
4716  if(Kaon(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE))
4717  {
4718  fCorrelationFunctions[3][7]->Fill(Q2(agtrack1,agtrack2));
4719  }
4720  // a4) pi-K- [7][8]
4721  if(Pion(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
4722  {
4723  fCorrelationFunctions[7][8]->Fill(Q2(agtrack1,agtrack2));
4724  }
4725  // b) pion-proton
4726  // b1) pi+p+ [2][4]:
4727  if(Pion(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
4728  {
4729  fCorrelationFunctions[2][4]->Fill(Q2(agtrack1,agtrack2));
4730  }
4731  // b2) pi+p- [2][9]
4732  if(Pion(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4733  {
4734  fCorrelationFunctions[2][9]->Fill(Q2(agtrack1,agtrack2));
4735  }
4736  // b3) p+pi- [4][7]
4737  if(Proton(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE))
4738  {
4739  fCorrelationFunctions[4][7]->Fill(Q2(agtrack1,agtrack2));
4740  }
4741  // b4) pi-p- [7][9]
4742  if(Pion(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4743  {
4744  fCorrelationFunctions[7][9]->Fill(Q2(agtrack1,agtrack2));
4745  }
4746  // c) kaon-proton
4747  // c1) K+p+ [3][4]:
4748  if(Kaon(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
4749  {
4750  fCorrelationFunctions[3][4]->Fill(Q2(agtrack1,agtrack2));
4751  }
4752  // c2) K+p- [3][9]
4753  if(Kaon(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4754  {
4755  fCorrelationFunctions[3][9]->Fill(Q2(agtrack1,agtrack2));
4756  }
4757  // c3) p+K- [4][8]
4758  if(Proton(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
4759  {
4760  fCorrelationFunctions[4][8]->Fill(Q2(agtrack1,agtrack2));
4761  }
4762  // c4) K-p- [8][9]
4763  if(Kaon(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4764  {
4765  fCorrelationFunctions[8][9]->Fill(Q2(agtrack1,agtrack2));
4766  }
4767  } // for(Int_t iTrack2=0;iTrack2<nTracks;iTrack2++)
4768  } // for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
4769 
4770  // c) Three nested loops to calculate C(Q3), just an example; TBI
4772 
4773  // d) Four nested loops to calculate C(Q4), just an example. TBI
4775 
4776 } // void AliAnalysisTaskMultiparticleFemtoscopy::CalculateCorrelationFunctions(AliAODEvent *aAOD)
4777 
4778 //=======================================================================================================================
4779 
4781 {
4782  // Calculate 3-particle correlation functions.
4783 
4784  // a) Insanity checks;
4785  // b) Three nested loops to calculate C(Q3), just an example.
4786 
4787  // a) Insanity checks:
4788  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::Calculate3pCorrelationFunctions(AliAODEvent *aAOD)";
4789  if(!aAOD){Fatal(sMethodName.Data(),"!aAOD");}
4790  if(0 == fGlobalTracksAOD[0]->GetSize()){Fatal(sMethodName.Data(),"0 == fGlobalTracksAOD[0]->GetSize()");} // this case shall be already treated in UserExec
4791 
4792  // b) Three nested loops to calculate C(Q3), just an example:
4793  Int_t nTracks = aAOD->GetNumberOfTracks();
4794  for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
4795  {
4796  AliAODTrack *atrack1 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack1));
4797  // TBI Temporary track insanity checks:
4798  if(!atrack1){Fatal(sMethodName.Data(),"!atrack1");} // TBI keep this for some time, eventually just continue
4799  if(atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4800  if(atrack1->TestFilterBit(128) && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->TestFiletrBit(128) && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4801  if(!PassesCommonTrackCuts(atrack1)){continue;} // TBI re-think
4802  // Corresponding AOD global track:
4803  Int_t id1 = atrack1->GetID();
4804  AliAODTrack *gtrack1 = dynamic_cast<AliAODTrack*>(id1>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id1)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id1+1))));
4805  if(!gtrack1){Fatal(sMethodName.Data(),"!gtrack1");} // TBI keep this for some time, eventually just continue
4806  Int_t gid1 = (id1>=0 ? id1 : -(id1+1)); // ID of corresponding global track
4807  // Common track selection criteria for all "normal" global tracks:
4808  if(!PassesGlobalTrackCuts(gtrack1)){continue;}
4809 
4810  // Loop over the 2nd particle:
4811  for(Int_t iTrack2=0;iTrack2<nTracks;iTrack2++)
4812  {
4813  if(iTrack2<=iTrack1){continue;} // Eliminate self-evident self-correlations, and permutations as well
4814  AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack2));
4815  // TBI Temporary track insanity checks:
4816  if(!atrack2){Fatal(sMethodName.Data(),"!atrack2");} // TBI keep this for some time, eventually just continue
4817  if(atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4818  if(atrack2->TestFilterBit(128) && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->TestFiletrBit(128) && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4819  if(!PassesCommonTrackCuts(atrack2)){continue;} // TBI re-think
4820  // Corresponding AOD global track:
4821  Int_t id2 = atrack2->GetID();
4822  AliAODTrack *gtrack2 = dynamic_cast<AliAODTrack*>(id2>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id2)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id2+1))));
4823  if(!gtrack2){Fatal(sMethodName.Data(),"!gtrack2");} // TBI keep this for some time, eventually just continue
4824  Int_t gid2 = (id2>=0 ? id2 : -(id2+1)); // ID of corresponding global track
4825  if(gid1==gid2){continue;} // Eliminate not-so-evident self-correlations:
4826  // Common track selection criteria for all "normal" global tracks:
4827  if(!PassesGlobalTrackCuts(gtrack2)){continue;}
4828 
4829  // Loop over the 3rd particle:
4830  for(Int_t iTrack3=0;iTrack3<nTracks;iTrack3++)
4831  {
4832  if(iTrack3<=iTrack2 || iTrack3<=iTrack1){continue;} // Eliminate self-evident self-correlations, and permutations as well
4833  AliAODTrack *atrack3 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack3));
4834  // TBI Temporary track insanity checks:
4835  if(!atrack3){Fatal(sMethodName.Data(),"!atrack3");} // TBI keep this for some time, eventually just continue
4836  if(atrack3->GetID()>=0 && atrack3->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack3->GetID()>=0 && atrack3->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4837  if(atrack3->TestFilterBit(128) && atrack3->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack3->TestFiletrBit(128) && atrack3->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4838  if(!PassesCommonTrackCuts(atrack3)){continue;} // TBI re-think
4839  // Corresponding AOD global track:
4840  Int_t id3 = atrack3->GetID();
4841  AliAODTrack *gtrack3 = dynamic_cast<AliAODTrack*>(id3>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id3)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id3+1))));
4842  if(!gtrack3){Fatal(sMethodName.Data(),"!gtrack3");} // TBI keep this for some time, eventually just continue
4843  Int_t gid3 = (id3>=0 ? id3 : -(id3+1)); // ID of corresponding global track
4844  if(gid3==gid2 || gid3==gid1){continue;} // Eliminate not-so-evident self-correlations
4845  // Common track selection criteria for all "normal" global tracks:
4846  if(!PassesGlobalTrackCuts(gtrack3)){continue;}
4847 
4848  // Okay, so we have three tracks, let's decide whether kinemetics from 'atrack' or 'gtrack' is taken, let's check PID, and fill the correlation functions:
4849  AliAODTrack *agtrack1 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4850  AliAODTrack *agtrack2 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4851  AliAODTrack *agtrack3 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4853  {
4854  agtrack1 = gtrack1;
4855  agtrack2 = gtrack2;
4856  agtrack3 = gtrack3;
4857  }
4858  else
4859  {
4860  agtrack1 = atrack1;
4861  agtrack2 = atrack2;
4862  agtrack3 = atrack3;
4863  } // if(fFillControlHistogramsWithGlobalTrackInfo)
4864  if(!agtrack1){Fatal(sMethodName.Data(),"!agtrack1");}
4865  if(!agtrack2){Fatal(sMethodName.Data(),"!agtrack2");}
4866  if(!agtrack3){Fatal(sMethodName.Data(),"!agtrack3");}
4867 
4868  // Cases of interest:
4869  // a) Same species and same charge;
4870  // b) Same species but different charge combinations (modulo permutations);
4871  // c) Two pions + something else (modulo permutations);
4872  // d) Two nucleons + something else (modulo permutations).
4873 
4874  // a) Same species:
4875  // pi+pi+pi+
4876  if(Pion(gtrack1,1,kTRUE) && Pion(gtrack2,1,kTRUE) && Pion(gtrack3,1,kTRUE))
4877  {
4878  f3pCorrelationFunctions[2][2][2]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4879  }
4880  // pi-pi-pi-
4881  if(Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,-1,kTRUE) && Pion(gtrack3,-1,kTRUE))
4882  {
4883  f3pCorrelationFunctions[7][7][7]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4884  }
4885  // K+K+K+
4886  if(Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE) && Kaon(gtrack3,1,kTRUE))
4887  {
4888  f3pCorrelationFunctions[3][3][3]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4889  }
4890  // K-K-K-
4891  if(Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE) && Kaon(gtrack3,-1,kTRUE))
4892  {
4893  f3pCorrelationFunctions[8][8][8]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4894  }
4895  // p+p+p+
4896  if(Proton(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE) && Proton(gtrack3,1,kTRUE))
4897  {
4898  f3pCorrelationFunctions[4][4][4]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4899  }
4900  // p-p-p-
4901  if(Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE) && Proton(gtrack3,-1,kTRUE))
4902  {
4903  f3pCorrelationFunctions[9][9][9]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4904  }
4905 
4906  // b) Same species but different charge combinations (modulo permutations):
4907  // pi+pi+pi-
4908  if(Pion(gtrack1,1,kTRUE) && Pion(gtrack2,1,kTRUE) && Pion(gtrack3,-1,kTRUE))
4909  {
4910  f3pCorrelationFunctions[2][2][7]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4911  }
4912  // pi+pi-pi-
4913  if(Pion(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE) && Pion(gtrack3,-1,kTRUE))
4914  {
4915  f3pCorrelationFunctions[2][7][7]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4916  }
4917  // K+K+K-
4918  if(Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE) && Kaon(gtrack3,-1,kTRUE))
4919  {
4920  f3pCorrelationFunctions[3][3][8]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4921  }
4922  // K+K-K-
4923  if(Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE) && Kaon(gtrack3,-1,kTRUE))
4924  {
4925  f3pCorrelationFunctions[3][8][8]->Fill(<