AliPhysics  d497547 (d497547)
 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  fEab_TEST6[t] = NULL;
2630  }
2631 
2632 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForCorrelationFunctionsTEST()
2633 
2634 //=======================================================================================================================
2635 
2637 {
2638  // Initialize all arrays for test background.
2639 
2640  const Int_t nTestsMax = 10; // see what is hardwired in .h file
2641 
2642  for(Int_t t=0;t<nTestsMax;t++)
2643  {
2644  fBackgroundTESTSublist[t] = NULL;
2645  fFillBackgroundTEST[t] = kFALSE;
2646  }
2647 
2648  for(Int_t t=0;t<nTestsMax;t++)
2649  {
2650  for(Int_t q=0;q<2;q++) // Q2 and Q3
2651  {
2652  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
2653  {
2654  for(Int_t d=0;d<10;d++) // differential options, see what is hardwired in .h file
2655  {
2656  fBackgroundTEST[t][q][tt][d] = NULL;
2657  }
2658  } // for(Int_t tt=0;tt<7;tt++)
2659  } // for(Int_t q=0;q<2;q++) // Q2 and Q3
2660  } // for(Int_t t=0;t<nTestsMax;t++)
2661 
2662  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
2663  {
2664  fMixedEventsTEST[me] = NULL;
2665  fGlobalTracksAODTEST[me] = NULL;
2666  }
2667 
2668  for(Int_t t=0;t<2;t++)
2669  {
2670  fBackgroundYieldTEST[t] = NULL;
2671  }
2672 
2673 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForBackgroundTEST()
2674 
2675 //=======================================================================================================================
2676 
2678 {
2679  // Initialize all arrays for 'hybrid approach'.
2680 
2681  for(Int_t t=0;t<5;t++) // five distinct terms in the numerator of Eq. (18)
2682  {
2683  fDistributionsVsQ3[t] = NULL;
2684  }
2685 
2686  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
2687  {
2688  fMixedEventsHA[me] = NULL;
2689  fGlobalTracksAODHA[me] = NULL;
2690  }
2691 
2692 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForHybridApproach()
2693 
2694 //=======================================================================================================================
2695 
2697 {
2698  // Book everything for test correlation functions.
2699 
2700  // a) Book the profile holding all the flags for test correlation functions;
2701  // b) Book correlation functions for all tests;
2702  // c) Store temporarely also the yield.
2703 
2704  // a) Book the profile holding all the flags for test correlation functions:
2705  const Int_t nTests = 10;
2706  fCorrelationFunctionsTESTFlagsPro = new TProfile("fCorrelationFunctionsTESTFlagsPro","Flags and settings for test correlation functions",nTests,0,nTests);
2707  fCorrelationFunctionsTESTFlagsPro->SetTickLength(-0.01,"Y");
2708  fCorrelationFunctionsTESTFlagsPro->SetMarkerStyle(25);
2709  fCorrelationFunctionsTESTFlagsPro->SetLabelSize(0.04);
2710  fCorrelationFunctionsTESTFlagsPro->SetLabelOffset(0.02,"Y");
2711  fCorrelationFunctionsTESTFlagsPro->SetStats(kFALSE);
2712  fCorrelationFunctionsTESTFlagsPro->SetFillColor(kGray);
2713  fCorrelationFunctionsTESTFlagsPro->SetLineColor(kBlack);
2714  for(Int_t t=0;t<nTests;t++)
2715  {
2716  fCorrelationFunctionsTESTFlagsPro->GetXaxis()->SetBinLabel(t+1,Form("fFillCorrelationFunctionsTEST[%d]",t)); fCorrelationFunctionsTESTFlagsPro->Fill(t+0.5,fFillCorrelationFunctionsTEST[t]);
2717  }
2719 
2720  // b) Book correlation functions for all tests:
2721  const Int_t nCumulants = 4; // TBI this applies only to vs. Q3
2722  const Int_t n3pCumulantTerms = 7;
2723  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++)
2724  const Int_t n2pCumulantTerms = 3;
2725  TString s2pCumulantTerms[n2pCumulantTerms] = {"#LTX_{1}#GT","#LTX_{2}#GT","#LTX_{1}X_{2}#GT"};
2726  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"};
2727  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
2728  TString sXYZ[3] = {"x","y","z"};
2729  TString sQs[nQs] = {"Q_{2}","Q_{3}"};
2730  for(Int_t t=0;t<nTests;t++)
2731  {
2733  {
2734 
2735  // Correlations vs. Q2:
2736  for(Int_t ct=0;ct<n2pCumulantTerms;ct++)
2737  {
2738  for(Int_t xyz=0;xyz<3;xyz++)
2739  {
2740  fCorrelationFunctionsTEST[t][0][ct][xyz] = new TProfile(Form("fCorrelationFunctionsTEST[%d][0][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ2bins,fnQ2min,fnQ2max);
2741  fCorrelationFunctionsTEST[t][0][ct][xyz]->SetStats(kFALSE);
2742  fCorrelationFunctionsTEST[t][0][ct][xyz]->SetMarkerStyle(kFullSquare);
2743  fCorrelationFunctionsTEST[t][0][ct][xyz]->SetMarkerColor(kBlue);
2744  fCorrelationFunctionsTEST[t][0][ct][xyz]->SetLineColor(kBlue);
2745  fCorrelationFunctionsTEST[t][0][ct][xyz]->GetXaxis()->SetTitle(sQs[0].Data());
2746  fCorrelationFunctionsTEST[t][0][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",s2pCumulantTerms[ct].Data(),sXYZ[xyz].Data()));
2748  } // for(Int_t xyz=0;xyz<3;xyz++)
2749  } // for(Int_t ct=0;ct<n2pCumulantTerms;ct++)
2750 
2751  // Cumulants vs. Q2:
2752  for(Int_t xyz=0;xyz<3;xyz++)
2753  {
2754  fSignalCumulantsTEST[t][0][0][xyz] = new TProfile(Form("fSignalCumulantsTEST[%d][0][%d][%d]",t,0,xyz),Form("TEST: %d",t),fnQ2bins,fnQ2min,fnQ2max);
2755  fSignalCumulantsTEST[t][0][0][xyz]->SetStats(kFALSE);
2756  fSignalCumulantsTEST[t][0][0][xyz]->SetMarkerStyle(kFullSquare);
2757  fSignalCumulantsTEST[t][0][0][xyz]->SetMarkerColor(kGreen+2);
2758  fSignalCumulantsTEST[t][0][0][xyz]->SetLineColor(kGreen+2);
2759  fSignalCumulantsTEST[t][0][0][xyz]->GetXaxis()->SetTitle(sQs[0].Data());
2760  fSignalCumulantsTEST[t][0][0][xyz]->GetYaxis()->SetTitle(Form("#LTX_{1}X_{2}#GT_{c,%s}",sXYZ[xyz].Data()));
2762  } // for(Int_t xyz=0;xyz<3;xyz++)
2763 
2764  // Correlations vs. Q3: TBI this clearly can be unified with the above.
2765  for(Int_t ct=0;ct<n3pCumulantTerms;ct++)
2766  {
2767  for(Int_t xyz=0;xyz<3;xyz++)
2768  {
2769  fCorrelationFunctionsTEST[t][1][ct][xyz] = new TProfile(Form("fCorrelationFunctionsTEST[%d][1][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ3bins,fnQ3min,fnQ3max);
2770  fCorrelationFunctionsTEST[t][1][ct][xyz]->SetStats(kFALSE);
2771  fCorrelationFunctionsTEST[t][1][ct][xyz]->SetMarkerStyle(kFullSquare);
2772  fCorrelationFunctionsTEST[t][1][ct][xyz]->SetMarkerColor(kBlue);
2773  fCorrelationFunctionsTEST[t][1][ct][xyz]->SetLineColor(kBlue);
2774  fCorrelationFunctionsTEST[t][1][ct][xyz]->GetXaxis()->SetTitle(sQs[1].Data());
2775  fCorrelationFunctionsTEST[t][1][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",s3pCumulantTerms[ct].Data(),sXYZ[xyz].Data()));
2777  } // for(Int_t xyz=0;xyz<3;xyz++)
2778  } // for(Int_t ct=0;ct<n3pCumulantTerms;ct++)
2779 
2780  // Cumulants vs. Q3:
2781  for(Int_t ct=0;ct<nCumulants;ct++) // TBI this applies only to vs. Q3
2782  {
2783  for(Int_t xyz=0;xyz<3;xyz++)
2784  {
2785  fSignalCumulantsTEST[t][1][ct][xyz] = new TProfile(Form("fSignalCumulantsTEST[%d][1][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ3bins,fnQ3min,fnQ3max);
2786  fSignalCumulantsTEST[t][1][ct][xyz]->SetStats(kFALSE);
2787  fSignalCumulantsTEST[t][1][ct][xyz]->SetMarkerStyle(kFullSquare);
2788  fSignalCumulantsTEST[t][1][ct][xyz]->SetMarkerColor(kGreen+2);
2789  fSignalCumulantsTEST[t][1][ct][xyz]->SetLineColor(kGreen+2);
2790  fSignalCumulantsTEST[t][1][ct][xyz]->GetXaxis()->SetTitle(sQs[1].Data());
2791  fSignalCumulantsTEST[t][1][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",sCumulants[ct].Data(),sXYZ[xyz].Data()));
2793  } // for(Int_t xyz=0;xyz<3;xyz++)
2794  } // for(Int_t ct=0;ct<nCumulants;ct++)
2795 
2796  } // if(fFillCorrelationFunctionsTEST[t])
2797 
2798  } // for(Int_t t=0;t<nTests;t++)
2799 
2800  // c) Store temporarely also the yield:
2801  // vs. Q2:
2802  fSignalYieldTEST[0] = new TH1F("fSignalYieldTEST[0]","dN(X_{1}X_{2})/Q_{2} for signal (TEST_0)",fnQ2bins,fnQ2min,fnQ2max);
2803  fSignalYieldTEST[0]->SetStats(kFALSE);
2804  fSignalYieldTEST[0]->SetLineColor(kBlue);
2805  fSignalYieldTEST[0]->SetFillColor(kBlue-10);
2806  fSignalYieldTEST[0]->GetXaxis()->SetTitle("Q_{2}");
2807  fSignalYieldTEST[0]->GetYaxis()->SetTitle("counts");
2809  // vs. Q3:
2810  fSignalYieldTEST[1] = new TH1F("fSignalYieldTEST[1]","dN(X_{1}X_{2}X_{3})/Q_{3} for signal (TEST_1)",fnQ3bins,fnQ3min,fnQ3max);
2811  fSignalYieldTEST[1]->SetStats(kFALSE);
2812  fSignalYieldTEST[1]->SetLineColor(kBlue);
2813  fSignalYieldTEST[1]->SetFillColor(kBlue-10);
2814  fSignalYieldTEST[1]->GetXaxis()->SetTitle("Q_{3}");
2815  fSignalYieldTEST[1]->GetYaxis()->SetTitle("counts");
2817 
2818  // TBI temporary here for "Test 6" (checking E_a^b as L.I. observable):
2819  // vs. Q2:
2820  fEab_TEST6[0] = new TH1F("fEab_TEST6[0]","",10*fnQ2bins,fnQ2min,fnQ2max);
2821  fEab_TEST6[0]->SetStats(kFALSE);
2822  fEab_TEST6[0]->SetLineColor(kBlue);
2823  fEab_TEST6[0]->SetFillColor(kBlue-10);
2824  fEab_TEST6[0]->GetXaxis()->SetTitle("Q_{2}");
2825  fEab_TEST6[0]->GetYaxis()->SetTitle("dN/dE_{a}^{b}");
2827 
2828  fEab_TEST6[1] = new TH1F("fEab_TEST6[1]","",10*fnQ2bins,fnQ2min,fnQ2max);
2829  fEab_TEST6[1]->SetStats(kFALSE);
2830  fEab_TEST6[1]->SetLineColor(kRed);
2831  fEab_TEST6[1]->SetFillColor(kRed-10);
2832  fEab_TEST6[1]->GetXaxis()->SetTitle("Q_{2}");
2833  fEab_TEST6[1]->GetYaxis()->SetTitle("dN/dE_{a}^{b}");
2835 
2836 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForCorrelationFunctionsTEST()
2837 
2838 //=======================================================================================================================
2839 
2841 {
2842  // Book everything for test background.
2843 
2844  // a) Book the profile holding all the flags for test background;
2845  // b) Book background objects for all tests;
2846  // c) Book fMixedEventsTEST[3] and fGlobalTracksAODTEST[3];
2847  // d) Store temporary the yield as well.
2848 
2849  // a) Book the profile holding all the flags for test background:
2850  const Int_t nTests = 10;
2851  fBackgroundTESTFlagsPro = new TProfile("fBackgroundTESTFlagsPro","Flags and settings for test background",nTests,0,nTests);
2852  fBackgroundTESTFlagsPro->SetTickLength(-0.01,"Y");
2853  fBackgroundTESTFlagsPro->SetMarkerStyle(25);
2854  fBackgroundTESTFlagsPro->SetLabelSize(0.04);
2855  fBackgroundTESTFlagsPro->SetLabelOffset(0.02,"Y");
2856  fBackgroundTESTFlagsPro->SetStats(kFALSE);
2857  fBackgroundTESTFlagsPro->SetFillColor(kGray);
2858  fBackgroundTESTFlagsPro->SetLineColor(kBlack);
2859  for(Int_t t=0;t<nTests;t++)
2860  {
2861  fBackgroundTESTFlagsPro->GetXaxis()->SetBinLabel(t+1,Form("fFillBackgroundTEST[%d]",t)); fBackgroundTESTFlagsPro->Fill(t+0.5,fFillBackgroundTEST[t]);
2862  }
2864 
2865  // b) Book background objects for all tests:
2866  const Int_t nCumulants = 4; // TBI this applies only to vs. Q3
2867  const Int_t n3pCumulantTerms = 7;
2868  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++)
2869  const Int_t n2pCumulantTerms = 3;
2870  TString s2pCumulantTerms[n2pCumulantTerms] = {"#LTX_{1}#GT","#LTX_{2}#GT","#LTX_{1}X_{2}#GT"};
2871  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"};
2872  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
2873  TString sXYZ[3] = {"x","y","z"};
2874  TString sQs[nQs] = {"Q_{2}","Q_{3}"};
2875  for(Int_t t=0;t<nTests;t++)
2876  {
2877  if(fFillBackgroundTEST[t])
2878  {
2879 
2880  // Background correlations vs. Q2:
2881  for(Int_t ct=0;ct<n2pCumulantTerms;ct++)
2882  {
2883  for(Int_t xyz=0;xyz<3;xyz++)
2884  {
2885  fBackgroundTEST[t][0][ct][xyz] = new TProfile(Form("fBackgroundTEST[%d][0][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ2bins,fnQ2min,fnQ2max);
2886  fBackgroundTEST[t][0][ct][xyz]->SetStats(kFALSE);
2887  fBackgroundTEST[t][0][ct][xyz]->SetMarkerStyle(kFullSquare);
2888  fBackgroundTEST[t][0][ct][xyz]->SetMarkerColor(kBlue);
2889  fBackgroundTEST[t][0][ct][xyz]->SetLineColor(kBlue);
2890  fBackgroundTEST[t][0][ct][xyz]->GetXaxis()->SetTitle(sQs[0].Data());
2891  fBackgroundTEST[t][0][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",s2pCumulantTerms[ct].Data(),sXYZ[xyz].Data()));
2892  fBackgroundTESTSublist[t]->Add(fBackgroundTEST[t][0][ct][xyz]);
2893  } // for(Int_t xyz=0;xyz<3;xyz++)
2894  } // for(Int_t ct=0;ct<n2pCumulantTerms;ct++)
2895 
2896  // Background cumulants vs. Q2:
2897  for(Int_t xyz=0;xyz<3;xyz++)
2898  {
2899  fBackgroundCumulantsTEST[t][0][0][xyz] = new TProfile(Form("fBackgroundCumulantsTEST[%d][0][%d][%d]",t,0,xyz),Form("TEST: %d",t),fnQ2bins,fnQ2min,fnQ2max);
2900  fBackgroundCumulantsTEST[t][0][0][xyz]->SetStats(kFALSE);
2901  fBackgroundCumulantsTEST[t][0][0][xyz]->SetMarkerStyle(kFullSquare);
2902  fBackgroundCumulantsTEST[t][0][0][xyz]->SetMarkerColor(kGreen+2);
2903  fBackgroundCumulantsTEST[t][0][0][xyz]->SetLineColor(kGreen+2);
2904  fBackgroundCumulantsTEST[t][0][0][xyz]->GetXaxis()->SetTitle(sQs[0].Data());
2905  fBackgroundCumulantsTEST[t][0][0][xyz]->GetYaxis()->SetTitle(Form("#LTX_{1}X_{2}#GT_{c,%s}",sXYZ[xyz].Data()));
2906  fBackgroundTESTSublist[t]->Add(fBackgroundCumulantsTEST[t][0][0][xyz]);
2907  } // for(Int_t xyz=0;xyz<3;xyz++)
2908 
2909  // Background correlations vs. Q3: TBI this clearly can be unified with the above.
2910  for(Int_t ct=0;ct<n3pCumulantTerms;ct++)
2911  {
2912  for(Int_t xyz=0;xyz<3;xyz++)
2913  {
2914  fBackgroundTEST[t][1][ct][xyz] = new TProfile(Form("fBackgroundTEST[%d][1][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ3bins,fnQ3min,fnQ3max);
2915  fBackgroundTEST[t][1][ct][xyz]->SetStats(kFALSE);
2916  fBackgroundTEST[t][1][ct][xyz]->SetMarkerStyle(kFullSquare);
2917  fBackgroundTEST[t][1][ct][xyz]->SetMarkerColor(kBlue);
2918  fBackgroundTEST[t][1][ct][xyz]->SetLineColor(kBlue);
2919  fBackgroundTEST[t][1][ct][xyz]->GetXaxis()->SetTitle(sQs[1].Data());
2920  fBackgroundTEST[t][1][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",s3pCumulantTerms[ct].Data(),sXYZ[xyz].Data()));
2921  fBackgroundTESTSublist[t]->Add(fBackgroundTEST[t][1][ct][xyz]);
2922  } // for(Int_t xyz=0;xyz<3;xyz++)
2923  } // for(Int_t ct=0;ct<n3pCumulantTerms;ct++)
2924 
2925  // Background cumulants vs. Q3:
2926  for(Int_t ct=0;ct<nCumulants;ct++) // TBI this applies only to vs. Q3
2927  {
2928  for(Int_t xyz=0;xyz<3;xyz++)
2929  {
2930  fBackgroundCumulantsTEST[t][1][ct][xyz] = new TProfile(Form("fBackgroundCumulantsTEST[%d][1][%d][%d]",t,ct,xyz),Form("TEST: %d",t),fnQ3bins,fnQ3min,fnQ3max);
2931  fBackgroundCumulantsTEST[t][1][ct][xyz]->SetStats(kFALSE);
2932  fBackgroundCumulantsTEST[t][1][ct][xyz]->SetMarkerStyle(kFullSquare);
2933  fBackgroundCumulantsTEST[t][1][ct][xyz]->SetMarkerColor(kGreen+2);
2934  fBackgroundCumulantsTEST[t][1][ct][xyz]->SetLineColor(kGreen+2);
2935  fBackgroundCumulantsTEST[t][1][ct][xyz]->GetXaxis()->SetTitle(sQs[1].Data());
2936  fBackgroundCumulantsTEST[t][1][ct][xyz]->GetYaxis()->SetTitle(Form("%s_{%s}",sCumulants[ct].Data(),sXYZ[xyz].Data()));
2937  fBackgroundTESTSublist[t]->Add(fBackgroundCumulantsTEST[t][1][ct][xyz]);
2938  } // for(Int_t xyz=0;xyz<3;xyz++)
2939  } // for(Int_t ct=0;ct<nCumulants;ct++)
2940  } // if(fFillCorrelationFunctionsTEST[t])
2941 
2942  } // for(Int_t t=0;t<nTests;t++)
2943 
2944 
2945  // c) Book fMixedEventsTEST[3] and fGlobalTracksAODTEST[3]:
2946  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
2947  {
2948  fMixedEventsTEST[me] = new TClonesArray("AliAODTrack",10000);
2949  fGlobalTracksAODTEST[me] = new TExMap();
2950  }
2951 
2952  // d) Store temporary the yield as well:
2953  // vs. Q2:
2954  fBackgroundYieldTEST[0] = new TH1F("fBackgroundYieldTEST[0]","dN(X_{1}X_{2})/Q_{2} for mixed-events (TEST_0)",fnQ2bins,fnQ2min,fnQ2max);
2955  fBackgroundYieldTEST[0]->SetStats(kFALSE);
2956  fBackgroundYieldTEST[0]->SetLineColor(kRed);
2957  fBackgroundYieldTEST[0]->SetFillColor(kRed-10);
2958  fBackgroundYieldTEST[0]->GetXaxis()->SetTitle("Q_{2}");
2959  fBackgroundYieldTEST[0]->GetYaxis()->SetTitle("counts");
2961  // vs. Q3:
2962  fBackgroundYieldTEST[1] = new TH1F("fBackgroundYieldTEST[1]","dN(X_{1}X_{2}X_{3})/Q_{3} for mixed-events (TEST_1)",fnQ3bins,fnQ3min,fnQ3max);
2963  fBackgroundYieldTEST[1]->SetStats(kFALSE);
2964  fBackgroundYieldTEST[1]->SetLineColor(kRed);
2965  fBackgroundYieldTEST[1]->SetFillColor(kRed-10);
2966  fBackgroundYieldTEST[1]->GetXaxis()->SetTitle("Q_{3}");
2967  fBackgroundYieldTEST[1]->GetYaxis()->SetTitle("counts");
2969 
2970 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForBackgroundTEST()
2971 
2972 //=======================================================================================================================
2973 
2975 {
2976  // Book everything for 'hybrid approach'.
2977 
2978  // a) Book the profile holding all the flags;
2979  // b) Book fMixedEventsHA[3] and fGlobalTracksAODHA[3];
2980  // c) ...
2981 
2982  // a) Book the profile holding all the flags:
2983  fHybridApproachFlagsPro = new TProfile("fHybridApproachFlagsPro","Flags and settings for 'hybrid approach'",1,0,1);
2984  fHybridApproachFlagsPro->SetTickLength(-0.01,"Y");
2985  fHybridApproachFlagsPro->SetMarkerStyle(25);
2986  fHybridApproachFlagsPro->SetLabelSize(0.04);
2987  fHybridApproachFlagsPro->SetLabelOffset(0.02,"Y");
2988  fHybridApproachFlagsPro->SetStats(kFALSE);
2989  fHybridApproachFlagsPro->SetFillColor(kGray);
2990  fHybridApproachFlagsPro->SetLineColor(kBlack);
2991  fHybridApproachFlagsPro->GetXaxis()->SetBinLabel(1,"fDoHybridApproach"); fHybridApproachFlagsPro->Fill(0.5,fDoHybridApproach);
2993 
2994  if(!fDoHybridApproach){return;}
2995 
2996  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})"};
2997  for(Int_t t=0;t<5;t++) // five distinct terms in the numerator of Eq. (18)
2998  {
2999  fDistributionsVsQ3[t] = new TH1F(Form("fDistributionsVsQ3[%d]",t),sTerm[t].Data(),5*fnQ3bins,fnQ3min,fnQ3max); // TBI hardwired refined binning, remove eventually
3000  //fDistributionsVsQ3[t]->SetStats(kFALSE);
3001  fDistributionsVsQ3[t]->SetMarkerStyle(kFullSquare);
3002  fDistributionsVsQ3[t]->SetMarkerColor(kBlue);
3003  fDistributionsVsQ3[t]->SetLineColor(kBlue);
3004  fDistributionsVsQ3[t]->SetFillColor(kBlue-10);
3005  fDistributionsVsQ3[t]->GetXaxis()->SetTitle("Q_{3}");
3006  fDistributionsVsQ3[t]->GetYaxis()->SetTitle("counts");
3008  }
3009 
3010  // b) Book fMixedEventsHA[3] and fGlobalTracksAODHA[3]:
3011  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
3012  {
3013  fMixedEventsHA[me] = new TClonesArray("AliAODTrack",10000);
3014  fGlobalTracksAODHA[me] = new TExMap();
3015  }
3016 
3017  // c) ...
3018 
3019 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForHybridApproach()
3020 
3021 //=======================================================================================================================
3022 
3024 {
3025  // Book and nest all lists nested in the base list fHistList.
3026 
3027  // a) Book and nest lists for control histograms;
3028  // b) Book and nest lists for eveny-by-event histograms;
3029  // c) Book and nest lists for correlation functions;
3030  // d) Book and nest lists for background;
3031  // e) Book and nest lists for buffers;
3032  // f) Book and nest lists for QA;
3033  // g) Book and nest lists for common global track cuts;
3034  // h) Book and nest lists for test correlation functions;
3035  // i) Book and nest lists for test background;
3036  // j) Book and nest lists for 'hybrid approach'.
3037 
3038  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::BookAndNestAllLists()";
3039  if(!fHistList){Fatal(sMethodName.Data(),"fHistList is NULL");}
3040 
3041  // a) Book and nest lists for control histograms:
3042  fControlHistogramsList = new TList();
3043  fControlHistogramsList->SetName("Control_histograms");
3044  fControlHistogramsList->SetOwner(kTRUE);
3047  fControlHistogramsEventList->SetName("Event");
3048  fControlHistogramsEventList->SetOwner(kTRUE);
3051  fControlHistogramsNonIdentifiedParticlesList->SetName("Non-identified_particles");
3055  fControlHistogramsNonIdentifiedParticlesFTSFList->SetName("Non-identified_particles (f.t.s.f.)"); // for the specified filterbit
3059  fControlHistogramsIdentifiedParticlesList->SetName("Identified_particles");
3063  fControlHistogramsV0sList->SetName("V0s");
3064  fControlHistogramsV0sList->SetOwner(kTRUE);
3066 
3067  // b) Book and nest lists for eveny-by-event histograms:
3068  fEBEHistogramsList = new TList();
3069  fEBEHistogramsList->SetName("Event-by-event_histograms");
3070  fEBEHistogramsList->SetOwner(kTRUE);
3072 
3073  // c) Book and nest lists for correlation functions:
3075  fCorrelationFunctionsList->SetName("CorrelationFunctions");
3076  fCorrelationFunctionsList->SetOwner(kTRUE);
3078  // Correlation functions sublists:
3079  for(Int_t cfs=0;cfs<3;cfs++)
3080  {
3081  fCorrelationFunctionsSublist[cfs] = new TList();
3082  fCorrelationFunctionsSublist[cfs]->SetName(Form("f%dpCorrelationFunctions",cfs+2));
3083  fCorrelationFunctionsSublist[cfs]->SetOwner(kTRUE);
3085  } // for(Int_t cfs=0;cfs<3;cfs++)
3086 
3087  // d) Book and nest lists for background:
3088  fBackgroundList = new TList();
3089  fBackgroundList->SetName("Background");
3090  fBackgroundList->SetOwner(kTRUE);
3091  fHistList->Add(fBackgroundList);
3092  // Background sublists:
3093  for(Int_t bs=0;bs<3;bs++)
3094  {
3095  fBackgroundSublist[bs] = new TList();
3096  fBackgroundSublist[bs]->SetName(Form("f%dpBackground",bs+2));
3097  fBackgroundSublist[bs]->SetOwner(kTRUE);
3099  } // for(Int_t bs=0;bs<3;bs++)
3100 
3101  // e) Book and nest lists for buffers:
3102  fBuffersList = new TList();
3103  fBuffersList->SetName("Buffers");
3104  fBuffersList->SetOwner(kTRUE);
3105  fHistList->Add(fBuffersList);
3106 
3107  // f) Book and nest lists for QA:
3108  fQAList = new TList();
3109  fQAList->SetName("QA");
3110  fQAList->SetOwner(kTRUE);
3111  fHistList->Add(fQAList);
3112  if(fFillQAEvents)
3113  {
3114  fQAEventsList = new TList();
3115  fQAEventsList->SetName("Events");
3116  fQAEventsList->SetOwner(kTRUE);
3117  fQAList->Add(fQAEventsList);
3118  }
3119  if(fFillQAParticles)
3120  {
3121  fQAParticlesList = new TList();
3122  fQAParticlesList->SetName("Particles");
3123  fQAParticlesList->SetOwner(kTRUE);
3124  fQAList->Add(fQAParticlesList);
3125  }
3126 
3127  // g) Book and nest lists for common global track cuts:
3128  fGlobalTrackCutsList = new TList();
3129  fGlobalTrackCutsList->SetName("Global_track_cuts");
3130  fGlobalTrackCutsList->SetOwner(kTRUE);
3132 
3133  // h) Book and nest lists for test correlation functions:
3135  fCorrelationFunctionsTESTList->SetName("CorrelationFunctionsTEST");
3136  fCorrelationFunctionsTESTList->SetOwner(kTRUE);
3138  const Int_t nTests = 10;
3139  for(Int_t t=0;t<nTests;t++)
3140  {
3141  if(!fFillCorrelationFunctionsTEST[t]){continue;}
3143  fCorrelationFunctionsTESTSublist[t]->SetName(Form("TEST_%d",t));
3144  fCorrelationFunctionsTESTSublist[t]->SetOwner(kTRUE);
3146  } // for(Int_t t=0;t<nTests;t++)
3147 
3148  // i) Book and nest lists for test background:
3149  fBackgroundTESTList = new TList();
3150  fBackgroundTESTList->SetName("BackgroundTEST");
3151  fBackgroundTESTList->SetOwner(kTRUE);
3153  for(Int_t t=0;t<nTests;t++)
3154  {
3155  if(!fFillBackgroundTEST[t]){continue;}
3156  fBackgroundTESTSublist[t] = new TList();
3157  fBackgroundTESTSublist[t]->SetName(Form("TEST_%d",t));
3158  fBackgroundTESTSublist[t]->SetOwner(kTRUE);
3160  } // for(Int_t t=0;t<nTests;t++)
3161 
3162  // j) Book and nest lists for 'hybrid approach':
3163  fHybridApproachList = new TList();
3164  fHybridApproachList->SetName("HybridApproach");
3165  fHybridApproachList->SetOwner(kTRUE);
3167 
3168 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookAndNestAllLists()
3169 
3170 //=======================================================================================================================
3171 
3173 {
3174  // Book all unclassified objects temporary here. TBI
3175 
3176  fAnalysisType = new TString();
3177 
3178  fPIDResponse = new AliPIDResponse();
3179 
3180  for(Int_t index=0;index<fMaxNoGlobalTracksAOD;index++)
3181  {
3182  fGlobalTracksAOD[index] = new TExMap();
3183  }
3184 
3185 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverything()
3186 
3187 //=======================================================================================================================
3188 
3190 {
3191  // Book all the stuff for control histograms.
3192 
3193  // a) Book the profile holding all the flags for control histograms;
3194  // b) Common variables;
3195  // c) Book all control histograms...
3196  // c0) Event;
3197  // c1) Non-identified particles (for AOD these are "normal global" tracks);
3198  // c2) Non-identified particles for the specified filterbit (f.t.s.f.) (by default TPC-only);
3199  // c3) Identified particles;
3200  // c4) V0s.
3201 
3202  // a) Book the profile holding all the flags for control histograms: TBI stil incomplete
3203  fControlHistogramsFlagsPro = new TProfile("fControlHistogramsFlagsPro","Flags and settings for control histograms",1,0,1);
3204  fControlHistogramsFlagsPro->SetTickLength(-0.01,"Y");
3205  fControlHistogramsFlagsPro->SetMarkerStyle(25);
3206  fControlHistogramsFlagsPro->SetLabelSize(0.04);
3207  fControlHistogramsFlagsPro->SetLabelOffset(0.02,"Y");
3208  fControlHistogramsFlagsPro->SetStats(kFALSE);
3209  fControlHistogramsFlagsPro->SetFillColor(kGray);
3210  fControlHistogramsFlagsPro->SetLineColor(kBlack);
3211  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistograms"); fControlHistogramsFlagsPro->Fill(0.5,fFillControlHistograms);
3213 
3214  if(!fFillControlHistograms){return;}
3215 
3216  // b) Common variables:
3217  TString sParticleLabel[5] = {"e","#mu","#pi","K","p"};
3218  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()};
3219  // ...
3220 
3221  // c0) Event:
3222  // Book the profile holding all the flags for control histograms for global event observables:
3223  fControlHistogramsEventFlagsPro = new TProfile("fControlHistogramsEventFlagsPro","Flags and settings for TBI",1,0,1);
3224  fControlHistogramsEventFlagsPro->SetTickLength(-0.01,"Y");
3225  fControlHistogramsEventFlagsPro->SetMarkerStyle(25);
3226  fControlHistogramsEventFlagsPro->SetLabelSize(0.04);
3227  fControlHistogramsEventFlagsPro->SetLabelOffset(0.02,"Y");
3228  fControlHistogramsEventFlagsPro->SetStats(kFALSE);
3229  fControlHistogramsEventFlagsPro->SetFillColor(kGray);
3230  fControlHistogramsEventFlagsPro->SetLineColor(kBlack);
3231  fControlHistogramsEventFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsEvent"); fControlHistogramsEventFlagsPro->Fill(0.5,fFillControlHistogramsEvent);
3234  {
3235  fGetNumberOfTracksHist = new TH1I("fGetNumberOfTracksHist","aAOD->GetNumberOfTracks() (Remark: Not all of tracks are unique.)",10000,0,10000);
3236  //fGetNumberOfTracksHist->SetStats(kFALSE);
3237  fGetNumberOfTracksHist->SetFillColor(kBlue-10);
3239  fGetNumberOfGlobalTracksHist = new TH1I("fGetNumberOfGlobalTracksHist","fGlobalTracksAOD[0]->GetSize()",10000,0,10000);
3240  //fGetNumberOfGlobalTracksHist->SetStats(kFALSE);
3241  fGetNumberOfGlobalTracksHist->SetFillColor(kBlue-10);
3243  fGetNumberOfV0sHist = new TH1I("fGetNumberOfV0sHist","aAOD->GetNumberOfV0s() (Remark: Some V0s share the daughter.)",10000,0,10000);
3244  fGetNumberOfV0sHist->SetStats(kFALSE);
3245  fGetNumberOfV0sHist->SetFillColor(kBlue-10);
3247  fGetNumberOfCascadesHist = new TH1I("fGetNumberOfCascadesHist","aAOD->GetNumberOfCascades() (TBI: Not validated.)",10000,0,10000);
3248  fGetNumberOfCascadesHist->SetStats(kFALSE);
3249  fGetNumberOfCascadesHist->SetFillColor(kBlue-10);
3251  fGetMagneticFieldHist = new TH1D("fGetMagneticFieldHist","aAOD->GetMagneticField()",200,-10.,10.);
3252  fGetMagneticFieldHist->SetFillColor(kBlue-10);
3254  fGetEventTypeHist = new TH1I("fGetEventTypeHist","aAOD->GetEventType()",1000,0,10);
3255  fGetEventTypeHist->SetFillColor(kBlue-10);
3257  fGetCentralityHist = new TH1D("fGetCentralityHist","aAOD->GetCentrality()",100,0.,100.);
3258  fGetCentralityHist->SetFillColor(kBlue-10);
3260  TString sxyz[3] = {"X","Y","Z"};
3261  for(Int_t xyz=0;xyz<3;xyz++)
3262  {
3263  fVertexXYZ[xyz] = new TH1F(Form("fVertex%s",sxyz[xyz].Data()),Form("avtz->Get%s()",sxyz[xyz].Data()),100000,-50.,50);
3264  fVertexXYZ[xyz]->SetStats(kFALSE);
3266  }
3267  fGetNContributorsHist = new TH1I("fGetNContributorsHist","avtx->GetNContributors()",10000,0,10000);
3268  fGetNContributorsHist->SetStats(kFALSE);
3269  fGetNContributorsHist->SetFillColor(kBlue-10);
3271  fGetChi2perNDFHist = new TH1F("fGetChi2perNDFHist","avtx->GetChi2perNDF()",5000,0.,50.);
3272  fGetChi2perNDFHist->SetStats(kFALSE);
3273  fGetChi2perNDFHist->SetFillColor(kBlue-10);
3275  fGetNDaughtersHist = new TH1I("GetNDaughtersHist","avtx->GetNDaughters()",10000,0,10000);
3276  fGetNDaughtersHist->SetStats(kFALSE);
3277  fGetNDaughtersHist->SetFillColor(kBlue-10);
3279  // ...
3280  } // if(fFillControlHistogramsEvent)
3281 
3282  // c1) Non-identified particles:
3283  // Book the profile holding all the flags for control histograms for non-identified particles:
3284  fControlHistogramsNonIdentifiedParticlesFlagsPro = new TProfile("fControlHistogramsNonIdentifiedParticlesFlagsPro","Flags and settings for TBI",1,0,1);
3285  fControlHistogramsNonIdentifiedParticlesFlagsPro->SetTickLength(-0.01,"Y");
3288  fControlHistogramsNonIdentifiedParticlesFlagsPro->SetLabelOffset(0.02,"Y");
3292  fControlHistogramsNonIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsNonIdentifiedParticles"); fControlHistogramsNonIdentifiedParticlesFlagsPro->Fill(0.5,fFillControlHistogramsNonIdentifiedParticles);
3295  {
3296  fChargeHist = new TH1I("fChargeHist","atrack->Charge()",5,-2,3);
3297  fChargeHist->SetStats(kFALSE);
3298  fChargeHist->SetFillColor(kBlue-10);
3299  fChargeHist->GetXaxis()->SetBinLabel(1,"-2");
3300  fChargeHist->GetXaxis()->SetBinLabel(2,"-1");
3301  fChargeHist->GetXaxis()->SetBinLabel(3,"0");
3302  fChargeHist->GetXaxis()->SetBinLabel(4,"1");
3303  fChargeHist->GetXaxis()->SetBinLabel(5,"2");
3305  fGetTPCNclsHist = new TH1I("fGetTPCNclsHist","atrack->fGetTPCNclsHist()",200,0,200);
3306  fGetTPCNclsHist->SetStats(kFALSE);
3307  fGetTPCNclsHist->SetFillColor(kBlue-10);
3308  fGetTPCNclsHist->GetXaxis()->SetTitle("TPCNcls");
3310  fGetTPCsignalNHist = new TH1I("fGetTPCsignalNHist","atrack->fGetTPCsignalNHist()",200,0,200);
3311  fGetTPCsignalNHist->SetStats(kFALSE);
3312  fGetTPCsignalNHist->SetFillColor(kBlue-10);
3313  fGetTPCsignalNHist->GetXaxis()->SetTitle("TPCsignalN");
3315  fGetITSNclsHist = new TH1I("fGetITSNclsHist","atrack->fGetITSNclsHist()",200,0,200);
3316  fGetITSNclsHist->SetStats(kFALSE);
3317  fGetITSNclsHist->SetFillColor(kBlue-10);
3318  fGetITSNclsHist->GetXaxis()->SetTitle("ITSNcls");
3320  fdEdxVsPtHist = new TH2F("fdEdxVsPtHist","atrack->GetTPCmomentum(),atrack->GetTPCsignal()",1000,0.,20.,1000,-500.,500.);
3321  fdEdxVsPtHist->SetStats(kFALSE);
3323  fPtHist = new TH1F("fPtHist","atrack->Pt()",1000,0.,20.);
3324  fPtHist->SetStats(kFALSE);
3325  fPtHist->SetFillColor(kBlue-10);
3326  fPtHist->SetMinimum(0.);
3328  fEtaHist = new TH1F("fEtaHist","atrack->Eta()",200,-2.,2.);
3329  fEtaHist->SetStats(kFALSE);
3330  fEtaHist->SetFillColor(kBlue-10);
3331  fEtaHist->SetMinimum(0.);
3333  fPhiHist = new TH1F("fPhiHist","atrack->Phi()",360,0.,TMath::TwoPi());
3334  fPhiHist->SetStats(kFALSE);
3335  fPhiHist->SetFillColor(kBlue-10);
3336  fPhiHist->SetMinimum(0.);
3338  // TBI
3339  fMassHist = new TH1F("fMassHist","atrack->M()",10000,0.,10.);
3340  fMassHist->SetStats(kFALSE);
3341  fMassHist->SetFillColor(kBlue-10);
3342  fMassHist->SetMinimum(0.);
3343  for(Int_t nm=0;nm<5;nm++) // nominal masses
3344  {
3345  fMassHist->GetXaxis()->SetBinLabel(fMassHist->FindBin(dNominalMass[nm]),Form("m_{%s}",sParticleLabel[nm].Data()));
3346  }
3348  fGetFilterMap = new TH1I("fGetFilterMap","atrack->fGetFilterMap()",10000,0,10000);
3349  fGetFilterMap->SetStats(kFALSE);
3350  fGetFilterMap->SetFillColor(kBlue-10);
3351  fGetFilterMap->SetMinimum(0.);
3353  fGetPdgCode = new TH1I("fGetPdgCode","atrack->fGetPdgCode()",20000,-10000,10000);
3354  fGetPdgCode->SetStats(kFALSE);
3355  fGetPdgCode->SetFillColor(kBlue-10);
3356  fGetPdgCode->SetMinimum(0.);
3358  } // if(fFillControlHistogramsNonIdentifiedParticles)
3359 
3360  // c2) Non-identified particles for the specified filterbit (f.t.s.f.) (by default TPC-only);
3361  // 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);
3362  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro = new TProfile("fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro","Flags and settings for TBI",2,0,2);
3363  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro->SetTickLength(-0.01,"Y");
3366  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro->SetLabelOffset(0.02,"Y");
3370  fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsNonIdentifiedParticlesFTSF"); fControlHistogramsNonIdentifiedParticlesFTSFFlagsPro->Fill(0.5,fFillControlHistogramsNonIdentifiedParticlesFTSF);
3374  {
3375  fChargeFTSFHist = new TH1I("fChargeFTSFHist",Form("atrack->Charge(), fb = %d",fFilterBitFTSF),5,-2,3);
3376  fChargeFTSFHist->SetStats(kFALSE);
3377  fChargeFTSFHist->SetFillColor(kBlue-10);
3378  fChargeFTSFHist->GetXaxis()->SetBinLabel(1,"-2");
3379  fChargeFTSFHist->GetXaxis()->SetBinLabel(2,"-1");
3380  fChargeFTSFHist->GetXaxis()->SetBinLabel(3,"0");
3381  fChargeFTSFHist->GetXaxis()->SetBinLabel(4,"1");
3382  fChargeFTSFHist->GetXaxis()->SetBinLabel(5,"2");
3384  fGetTPCNclsFTSFHist = new TH1I("fGetTPCNclsFTSFHist",Form("atrack->fGetTPCNclsFTSFHist(), fb = %d",fFilterBitFTSF),200,0,200);
3385  fGetTPCNclsFTSFHist->SetStats(kFALSE);
3386  fGetTPCNclsFTSFHist->SetFillColor(kBlue-10);
3387  fGetTPCNclsFTSFHist->GetXaxis()->SetTitle("TPCNcls");
3389  fGetTPCsignalNFTSFHist = new TH1I("fGetTPCsignalNFTSFHist",Form("atrack->fGetTPCsignalNFTSFHist(), fb = %d",fFilterBitFTSF),200,0,200);
3390  fGetTPCsignalNFTSFHist->SetStats(kFALSE);
3391  fGetTPCsignalNFTSFHist->SetFillColor(kBlue-10);
3392  fGetTPCsignalNFTSFHist->GetXaxis()->SetTitle("TPCsignalN");
3394  fGetITSNclsFTSFHist = new TH1I("fGetITSNclsFTSFHist",Form("atrack->fGetITSNclsFTSFHist(), fb = %d",fFilterBitFTSF),200,0,200);
3395  fGetITSNclsFTSFHist->SetStats(kFALSE);
3396  fGetITSNclsFTSFHist->SetFillColor(kBlue-10);
3397  fGetITSNclsFTSFHist->GetXaxis()->SetTitle("ITSNcls");
3399  fdEdxVsPtFTSFHist = new TH2F("fdEdxVsPtFTSFHist",Form("atrack->GetTPCmomentum(),atrack->GetTPCsignal(), fb = %d",fFilterBitFTSF),1000,0.,20.,1000,-500.,500.);
3400  fdEdxVsPtFTSFHist->SetStats(kFALSE);
3402  fPtFTSFHist = new TH1F("fPtFTSFHist",Form("atrack->Pt(), fb = %d",fFilterBitFTSF),1000,0.,20.);
3403  fPtFTSFHist->SetStats(kFALSE);
3404  fPtFTSFHist->SetFillColor(kBlue-10);
3405  fPtFTSFHist->SetMinimum(0.);
3407  fEtaFTSFHist = new TH1F("fEtaFTSFHist",Form("atrack->Eta(), fb = %d",fFilterBitFTSF),200,-2.,2.);
3408  fEtaFTSFHist->SetStats(kFALSE);
3409  fEtaFTSFHist->SetFillColor(kBlue-10);
3410  fEtaFTSFHist->SetMinimum(0.);
3412  fPhiFTSFHist = new TH1F("fPhiFTSFHist",Form("atrack->Phi(), fb = %d",fFilterBitFTSF),360,0.,TMath::TwoPi());
3413  fPhiFTSFHist->SetStats(kFALSE);
3414  fPhiFTSFHist->SetFillColor(kBlue-10);
3415  fPhiFTSFHist->SetMinimum(0.);
3417  // TBI
3418  fMassFTSFHist = new TH1F("fMassFTSFHist",Form("atrack->M(), fb = %d",fFilterBitFTSF),10000,0.,10.);
3419  fMassFTSFHist->SetStats(kFALSE);
3420  fMassFTSFHist->SetFillColor(kBlue-10);
3421  fMassFTSFHist->SetMinimum(0.);
3422  for(Int_t nm=0;nm<5;nm++) // nominal masses
3423  {
3424  fMassFTSFHist->GetXaxis()->SetBinLabel(fMassFTSFHist->FindBin(dNominalMass[nm]),Form("m_{%s}",sParticleLabel[nm].Data()));
3425  }
3427  fGetFilterMap = new TH1I("fGetFilterMap",Form("atrack->fGetFilterMap(), fb = %d",fFilterBitFTSF),10000,0,10000);
3428  fGetFilterMap->SetStats(kFALSE);
3429  fGetFilterMap->SetFillColor(kBlue-10);
3430  fGetFilterMap->SetMinimum(0.);
3432  fGetPdgCode = new TH1I("fGetPdgCode",Form("atrack->fGetPdgCode(), fb = %d",fFilterBitFTSF),20000,-10000,10000);
3433  fGetPdgCode->SetStats(kFALSE);
3434  fGetPdgCode->SetFillColor(kBlue-10);
3435  fGetPdgCode->SetMinimum(0.);
3437  } // if(fFillControlFTSFHistogramsNonIdentifiedParticlesFTSF)
3438 
3439  // c3) Identified particles:
3440  fControlHistogramsIdentifiedParticlesFlagsPro = new TProfile("fControlHistogramsIdentifiedParticlesFlagsPro","Flags and settings for identified particles",4,0,4);
3441  fControlHistogramsIdentifiedParticlesFlagsPro->SetTickLength(-0.01,"Y");
3444  fControlHistogramsIdentifiedParticlesFlagsPro->SetLabelOffset(0.02,"Y");
3446  fControlHistogramsIdentifiedParticlesFlagsPro->SetFillColor(kGray);
3447  fControlHistogramsIdentifiedParticlesFlagsPro->SetLineColor(kBlack);
3448  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsIdentifiedParticles"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(0.5,fFillControlHistogramsIdentifiedParticles);
3449  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(2,"fUseDefaultInclusiveSigmaCuts"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(1.5,(Int_t)fUseDefaultInclusiveSigmaCuts);
3450  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(3,"fUseDefaultExclusiveSigmaCuts"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(2.5,(Int_t)fUseDefaultExclusiveSigmaCuts);
3451  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(4,"fFillControlHistogramsWithGlobalTrackInfo"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(3.5,(Int_t)fFillControlHistogramsWithGlobalTrackInfo);
3453 
3455  {
3456  // Default inclusive sigma cuts:
3457  if(fUseDefaultInclusiveSigmaCuts)
3458  {
3459  fInclusiveSigmaCuts[2] = 3.; // i.e. in function Pion(...) the inclusive cut for pions is 2 sigma
3460  fInclusiveSigmaCuts[3] = 3.; // i.e. in function Kaon(...) the inclusive cut for kaons is 2 sigma
3461  fInclusiveSigmaCuts[4] = 3.; // i.e. in function Proton(...) the inclusive cut for protons is 2 sigma
3462  }
3463  const Int_t nPidFunctions = 5; //
3464  TString sPidFunctions[nPidFunctions] = {"Electron(...)","Muon(...)","Pion(...)","Kaon(...)","Proton(...)"};
3465  const Int_t nParticleSpecies = 5;
3466  TString sParticleSpecies[nParticleSpecies] = {"e","#mu","#pi","K","p"};
3467  fInclusiveSigmaCutsPro = new TProfile("fInclusiveSigmaCutsPro","Inclusive sigma cuts",nPidFunctions,0.,nPidFunctions);
3468  fInclusiveSigmaCutsPro->GetYaxis()->SetTitle("n#sigma");
3469  for(Int_t pidFunction=0;pidFunction<5;pidFunction++)
3470  {
3471  fInclusiveSigmaCutsPro->SetTickLength(-0.01,"Y");
3472  fInclusiveSigmaCutsPro->SetMarkerStyle(25);
3473  fInclusiveSigmaCutsPro->SetLabelSize(0.04);
3474  fInclusiveSigmaCutsPro->SetLabelOffset(0.02,"Y");
3475  fInclusiveSigmaCutsPro->SetStats(kFALSE);
3476  fInclusiveSigmaCutsPro->SetFillColor(kGray);
3477  fInclusiveSigmaCutsPro->SetLineColor(kBlack);
3478  fInclusiveSigmaCutsPro->GetXaxis()->SetBinLabel(pidFunction+1,sPidFunctions[pidFunction].Data());
3479  fInclusiveSigmaCutsPro->Fill(pidFunction+0.5,fInclusiveSigmaCuts[pidFunction]);
3480  } // for(Int_t pidFunction=0;pidFunction<5;pidFunction++)
3482 
3483  // Default exclusive sigma cuts:
3484  if(fUseDefaultExclusiveSigmaCuts)
3485  {
3486  // Pion(...)
3487  fExclusiveSigmaCuts[2][3] = 4.; // i.e. in function Pion(...) the exclusive cut for kaons is 4 sigma
3488  fExclusiveSigmaCuts[2][4] = 4.; // i.e. in function Pion(...) the exclusive cut for protons is 4 sigma
3489  // Kaon(...)
3490  fExclusiveSigmaCuts[3][2] = 4.; // i.e. in function Kaon(...) the exclusive cut for pions is 4 sigma
3491  fExclusiveSigmaCuts[3][4] = 4.; // i.e. in function Kaon(...) the exclusive cut for protons is 4 sigma
3492  // Proton(...)
3493  fExclusiveSigmaCuts[4][2] = 4.; // i.e. in function Proton(...) the exclusive cut for pions is 4 sigma
3494  fExclusiveSigmaCuts[4][3] = 4.; // i.e. in function Proton(...) the exclusive cut for kaons is 4 sigma
3495  } // if(fUseDefaultExclusiveSigmaCuts)
3496  fExclusiveSigmaCutsPro = new TProfile2D("fExclusiveSigmaCutsPro","Exclusive sigma cuts",nPidFunctions,0.,nPidFunctions,nParticleSpecies,0.,nParticleSpecies);
3497  fExclusiveSigmaCutsPro->SetTickLength(-0.01,"Y");
3498  fExclusiveSigmaCutsPro->SetMarkerStyle(25);
3499  fExclusiveSigmaCutsPro->SetLabelSize(0.04);
3500  //fExclusiveSigmaCutsPro->SetLabelOffset(0.01,"X");
3501  //fExclusiveSigmaCutsPro->SetLabelOffset(0.01,"Y");
3502  fExclusiveSigmaCutsPro->SetTitleOffset(0.9,"Z");
3503  fExclusiveSigmaCutsPro->SetStats(kFALSE);
3504  fExclusiveSigmaCutsPro->SetFillColor(kGray);
3505  fExclusiveSigmaCutsPro->SetLineColor(kBlack);
3506  fExclusiveSigmaCutsPro->GetZaxis()->SetTitle("n#sigma");
3507  for(Int_t pidFunction=0;pidFunction<nPidFunctions;pidFunction++)
3508  {
3509  fExclusiveSigmaCutsPro->GetXaxis()->SetBinLabel(pidFunction+1,sPidFunctions[pidFunction].Data());
3510  for(Int_t particleSpecies=0;particleSpecies<nParticleSpecies;particleSpecies++)
3511  {
3512  if(0==pidFunction){fExclusiveSigmaCutsPro->GetYaxis()->SetBinLabel(particleSpecies+1,sParticleSpecies[particleSpecies].Data());}
3513  fExclusiveSigmaCutsPro->Fill(pidFunction+0.5,particleSpecies+0.5,fExclusiveSigmaCuts[pidFunction][particleSpecies]);
3514  } // for(Int_t pidFunction=0;pidFunction<nPidFunctions;pidFunction++)
3515  } // for(Int_t pidFunction=0;pidFunction<nPidFunctions;pidFunction++)
3517 
3518  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
3519  {
3520  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
3521  {
3522  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
3523  {
3524  if(fFillControlHistogramsWithGlobalTrackInfo)
3525  {
3526  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.);
3527  }
3528  else
3529  {
3530  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.);
3531  }
3532  fMassPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("m [GeV/c^{2}]");
3533  for(Int_t nm=0;nm<5;nm++) // nominal masses
3534  {
3535  fMassPIDHist[pid][pa][ps]->GetXaxis()->SetBinLabel(fMassPIDHist[pid][pa][ps]->FindBin(dNominalMass[nm]),Form("m_{%s}",sParticleLabel[nm].Data()));
3536  }
3537  fMassPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
3539  if(fFillControlHistogramsWithGlobalTrackInfo)
3540  {
3541  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.);
3542  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.);
3543  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.);
3544  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.);
3545  }
3546  else
3547  {
3548  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.);
3549  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.);
3550  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.);
3551  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.);
3552  }
3553  fPtPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
3554  fPtPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
3556  fPPIDHist[pid][pa][ps][0]->GetXaxis()->SetTitle("p_{x} [GeV/c]");
3557  fPPIDHist[pid][pa][ps][0]->SetFillColor(kBlue-10);
3559  fPPIDHist[pid][pa][ps][1]->GetXaxis()->SetTitle("p_{y} [GeV/c]");
3560  fPPIDHist[pid][pa][ps][1]->SetFillColor(kBlue-10);
3562  fPPIDHist[pid][pa][ps][2]->GetXaxis()->SetTitle("p_{z} [GeV/c]");
3563  fPPIDHist[pid][pa][ps][2]->SetFillColor(kBlue-10);
3565 
3566  if(fFillControlHistogramsWithGlobalTrackInfo)
3567  {
3568  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.);
3569  }
3570  else
3571  {
3572  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.);
3573  }
3574  fEtaPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
3576  if(fFillControlHistogramsWithGlobalTrackInfo)
3577  {
3578  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());
3579  }
3580  else
3581  {
3582  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());
3583  }
3584  fPhiPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("#phi");
3585  fPhiPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
3587  } // for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
3588  } // for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
3589  } // for(Int_t pid=0;pid<4;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
3590  } // if(fFillControlHistogramsIdentifiedParticles)
3591 
3592  // c4) V0s:
3593  // Book the profile holding all the flags for V0s:
3594  fControlHistogramsV0sFlagsPro = new TProfile("fControlHistogramsV0sFlagsPro","Flags and settings for V0s",1,0,1);
3595  fControlHistogramsV0sFlagsPro->SetTickLength(-0.01,"Y");
3596  fControlHistogramsV0sFlagsPro->SetMarkerStyle(25);
3597  fControlHistogramsV0sFlagsPro->SetLabelSize(0.04);
3598  fControlHistogramsV0sFlagsPro->SetLabelOffset(0.02,"Y");
3599  fControlHistogramsV0sFlagsPro->SetStats(kFALSE);
3600  fControlHistogramsV0sFlagsPro->SetFillColor(kGray);
3601  fControlHistogramsV0sFlagsPro->SetLineColor(kBlack);
3602  fControlHistogramsV0sFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsV0s"); fControlHistogramsV0sFlagsPro->Fill(0.5,fFillControlHistogramsV0s);
3605  {
3606  fGetNProngsHist = new TH1I("fGetNProngsHist","aAODv0->GetNProngs()",10,0,10);
3607  fGetNProngsHist->SetStats(kFALSE);
3608  fGetNProngsHist->SetFillColor(kBlue-10);
3610  // TBI
3611  fMassK0ShortHist = new TH1F("fMassK0ShortHist","aAODv0->MassK0Short()",1000000,0.,100.);
3612  //fMassK0ShortHist->SetStats(kFALSE);
3613  fMassK0ShortHist->SetFillColor(kBlue-10);
3614  Double_t dMassK0Short = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // nominal mass
3615  //fMassK0ShortHist->GetXaxis()->SetBinLabel(fMassK0ShortHist->FindBin(dMassK0Short),Form("m_{K_{S}^{0}} = %f",dMassK0Short));
3616  fMassK0ShortHist->SetBinContent(fMassK0ShortHist->FindBin(dMassK0Short),1e6);
3618  // TBI
3619  fMassLambdaHist = new TH1F("fMassLambdaHist","aAODv0->MassLambda()",1000000,0.,100.);
3620  //fMassLambdaHist->SetStats(kFALSE);
3621  fMassLambdaHist->SetFillColor(kBlue-10);
3622  Double_t dMassLambda = TDatabasePDG::Instance()->GetParticle(3122)->Mass(); // nominal mass
3623  //fMassLambdaHist->GetXaxis()->SetBinLabel(fMassLambdaHist->FindBin(dMassLambda),Form("m_{#Lambda^{0}} = %f",dMassLambda));
3624  fMassLambdaHist->SetBinContent(fMassLambdaHist->FindBin(dMassLambda),1e6);
3626  // TBI
3627  fMassAntiLambdaHist = new TH1F("fMassAntiLambdaHist","aAODv0->MassAntiLambda()",1000000,0.,100.);
3628  //fMassAntiLambdaHist->SetStats(kFALSE);
3629  fMassAntiLambdaHist->SetFillColor(kBlue-10);
3630  Double_t dMassAntiLambda = TDatabasePDG::Instance()->GetParticle(-3122)->Mass(); // nominal mass
3631  //fMassAntiLambdaHist->GetXaxis()->SetBinLabel(fMassAntiLambdaHist->FindBin(dMassAntiLambda),Form("m_{#bar{Lambda}^{0}} = %f",dMassAntiLambda));
3632  fMassAntiLambdaHist->SetBinContent(fMassAntiLambdaHist->FindBin(dMassAntiLambda),1e6);
3634  // TBI
3635  fOpenAngleV0Hist = new TH1F("fOpenAngleV0Hist","aAODv0->fOpenAngleV0()",10000,-0.044,TMath::Pi()+0.044);
3636  fOpenAngleV0Hist->SetStats(kFALSE);
3637  fOpenAngleV0Hist->SetFillColor(kBlue-10);
3639  // TBI
3640  fRadiusV0Hist = new TH1F("fRadiusV0Hist","aAODv0->fRadiusV0()",10000,0.,1000.);
3641  fRadiusV0Hist->SetStats(kFALSE);
3642  fRadiusV0Hist->SetFillColor(kBlue-10);
3644  // TBI
3645  fDcaV0ToPrimVertexHist = new TH1F("fDcaV0ToPrimVertexHist","aAODv0->fDcaV0ToPrimVertex()",10000,0.,1000.);
3646  fDcaV0ToPrimVertexHist->SetStats(kFALSE);
3647  fDcaV0ToPrimVertexHist->SetFillColor(kBlue-10);
3649  // TBI
3650  fMomV0XHist = new TH1F("fMomV0XHist","aAODv0->fMomV0X() = px(+) + px(-)",10000,-1000.,1000.);
3651  fMomV0XHist->SetStats(kFALSE);
3652  fMomV0XHist->SetFillColor(kBlue-10);
3654  // TBI
3655  fMomV0YHist = new TH1F("fMomV0YHist","aAODv0->fMomV0Y() = py(+) + py(-)",10000,-1000.,1000.);
3656  fMomV0YHist->SetStats(kFALSE);
3657  fMomV0YHist->SetFillColor(kBlue-10);
3659  // TBI
3660  fMomV0ZHist = new TH1F("fMomV0ZHist","aAODv0->fMomV0Z() = pz(+) + pz(-)",10000,-1000.,1000.);
3661  fMomV0ZHist->SetStats(kFALSE);
3662  fMomV0ZHist->SetFillColor(kBlue-10);
3664  // TBI
3665  fPtV0Hist = new TH1F("fPtV0Hist","pow(aAODv0->fPt2V0(),0.5)",10000,0.,100.);
3666  fPtV0Hist->SetStats(kFALSE);
3667  fPtV0Hist->SetFillColor(kBlue-10);
3669  // TBI
3670  fPseudoRapV0Hist = new TH1F("fPseudoRapV0Hist","aAODv0->PseudoRapV0()",1000,-10.,10.);
3671  fPseudoRapV0Hist->SetStats(kFALSE);
3672  fPseudoRapV0Hist->SetFillColor(kBlue-10);
3674  // TBI
3675  fPAHist = new TH2F("fPAHist","TBI",100,-2.,2.,100,0.,1.);
3676  fPAHist->SetStats(kFALSE);
3677  fPAHist->GetXaxis()->SetTitle("#alpha");
3678  fPAHist->GetYaxis()->SetTitle("p_{T}");
3680  } // if(fFillControlHistogramsV0s)
3681 
3682 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForControlHistograms()
3683 
3684 //=======================================================================================================================
3685 
3687 {
3688  // Book all the stuff for event-by-event objects.
3689 
3690  // a) Book the profile holding all the flags for event-by-event objects;
3691  // b) Book all event-by-event objects.
3692 
3693  // a) Book the profile holding all the flags for EBE objects:
3694  fEBEObjectsFlagsPro = new TProfile("fEBEObjectsFlagsPro","Flags and settings for event-by-event histograms",1,0,1);
3695  fEBEObjectsFlagsPro->SetTickLength(-0.01,"Y");
3696  fEBEObjectsFlagsPro->SetMarkerStyle(25);
3697  fEBEObjectsFlagsPro->SetLabelSize(0.04);
3698  fEBEObjectsFlagsPro->SetLabelOffset(0.02,"Y");
3699  fEBEObjectsFlagsPro->SetStats(kFALSE);
3700  fEBEObjectsFlagsPro->SetFillColor(kGray);
3701  fEBEObjectsFlagsPro->SetLineColor(kBlack);
3702  //fEBEObjectsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillEBEHistograms"); fEBEObjectsFlagsPro->Fill(0.5,fFillEBEHistograms);
3704 
3705  //if(!fFillEBEHistograms){return;} // TBI rethink
3706 
3707  // TBI
3708  fUniqueIDHistEBE = new TH1I("fUniqueIDHistEBE","TBI",40000,-20000,20000);
3709  fUniqueIDHistEBE->SetStats(kFALSE);
3710  fUniqueIDHistEBE->SetFillColor(kBlue-10);
3711  // TBI I do not want to store this histogram, right?
3712 
3713  // TBI add comment
3714  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
3715  {
3716  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
3717  {
3718  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
3719  {
3720  fPIDCA[pid][pa][ps] = new TClonesArray("AliAODTrack",10000);
3721  }
3722  }
3723  }
3724 
3725  // TBI add comment
3726  for(Int_t pid=0;pid<1;pid++) // [0=Lambda,1=...]
3727  {
3728  fPIDV0sCA[pid] = new TClonesArray("AliAODv0",10000);
3729  }
3730 
3731 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForEBEObjects()
3732 
3733 //=======================================================================================================================
3734 
3736 {
3737  // Book all the stuff for correlation functions objects.
3738 
3739  // a) Book the profile holding all the flags for correlation functions objects;
3740  // b) Book all 2p correlation functions;
3741  // c) Book TExMap *fCorrelationFunctionsIndices;
3742  // d) Book all 3p correlation functions;
3743  // e) Book all 4p correlation functions.
3744 
3745  // a) Book the profile holding all the flags for correlation functions objects:
3746  fCorrelationFunctionsFlagsPro = new TProfile("fCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",8,0,8);
3747  fCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
3748  fCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
3749  fCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
3750  fCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
3751  fCorrelationFunctionsFlagsPro->SetStats(kFALSE);
3752  fCorrelationFunctionsFlagsPro->SetFillColor(kGray);
3753  fCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
3754  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(0.5,fFillCorrelationFunctions);
3755  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(2,"fNormalizeCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(1.5,fNormalizeCorrelationFunctions);
3756  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(3,"fFill3pCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(2.5,fFill3pCorrelationFunctions);
3757  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(4,"fFill4pCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(3.5,fFill4pCorrelationFunctions);
3758  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(5,"fNormalizationOption"); fCorrelationFunctionsFlagsPro->Fill(4.5,fNormalizationOption);
3759  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(6,"fNormalizationInterval[0]"); fCorrelationFunctionsFlagsPro->Fill(5.5,fNormalizationInterval[0]);
3760  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(7,"fNormalizationInterval[1]"); fCorrelationFunctionsFlagsPro->Fill(6.5,fNormalizationInterval[1]);
3761  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(8,"fnMergedBins"); fCorrelationFunctionsFlagsPro->Fill(7.5,fnMergedBins);
3763  // 2p:
3764  f2pCorrelationFunctionsFlagsPro = new TProfile("f2pCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",1,0,1);
3765  f2pCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
3766  f2pCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
3767  f2pCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
3768  f2pCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
3769  f2pCorrelationFunctionsFlagsPro->SetStats(kFALSE);
3770  f2pCorrelationFunctionsFlagsPro->SetFillColor(kGray);
3771  f2pCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
3772  f2pCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f2pCorrelationFunctionsFlagsPro->Fill(0.5,0.); // TBI dummy at the moment
3774  // 3p:
3775  f3pCorrelationFunctionsFlagsPro = new TProfile("f3pCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",1,0,1);
3776  f3pCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
3777  f3pCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
3778  f3pCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
3779  f3pCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
3780  f3pCorrelationFunctionsFlagsPro->SetStats(kFALSE);
3781  f3pCorrelationFunctionsFlagsPro->SetFillColor(kGray);
3782  f3pCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
3783  f3pCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f3pCorrelationFunctionsFlagsPro->Fill(0.5,0.); // TBI dummy at the moment
3785  // 4p:
3786  f4pCorrelationFunctionsFlagsPro = new TProfile("f4pCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",1,0,1);
3787  f4pCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
3788  f4pCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
3789  f4pCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
3790  f4pCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
3791  f4pCorrelationFunctionsFlagsPro->SetStats(kFALSE);
3792  f4pCorrelationFunctionsFlagsPro->SetFillColor(kGray);
3793  f4pCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
3794  f4pCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f4pCorrelationFunctionsFlagsPro->Fill(0.5,0.); // TBI dummy at the moment
3796 
3797  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
3798 
3799  const Int_t nParticleSpecies = 5; // Supported at the moment: 0=e,1=mu,2=pi,3=K,4=p
3800  TString sParticles[2*nParticleSpecies] = {"e^{+}","#mu^{+}","#pi^{+}","K^{+}","p^{+}","e^{-}","#mu^{-}","#pi^{-}","K^{-}","p^{-}"};
3801 
3802  // b) Book all 2p correlation functions:
3803  // Remark 0: First particle in the pair is always the one with positive charge.
3804  // Remark 1: Diagonal elements are particle/antiparticle pairs.
3805  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]
3806  {
3807  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]
3808  {
3809  // Correlation functions:
3810  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
3811  fCorrelationFunctions[pid1][pid2]->SetStats(kTRUE);
3812  fCorrelationFunctions[pid1][pid2]->SetFillColor(kBlue-10);
3813  fCorrelationFunctions[pid1][pid2]->SetXTitle("Q_{2}");
3814  fCorrelationFunctions[pid1][pid2]->SetYTitle("C(Q_{2})");
3816  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3817  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3818 
3819  // c) Book TExMap *fCorrelationFunctionsIndices:
3820  // 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]
3821  fCorrelationFunctionsIndices = new TExMap();
3822  fCorrelationFunctionsIndices->Add(11,0);
3823  fCorrelationFunctionsIndices->Add(13,1);
3824  fCorrelationFunctionsIndices->Add(211,2);
3825  fCorrelationFunctionsIndices->Add(321,3);
3826  fCorrelationFunctionsIndices->Add(2212,4);
3827  fCorrelationFunctionsIndices->Add(-11,5);
3828  fCorrelationFunctionsIndices->Add(-13,6);
3829  fCorrelationFunctionsIndices->Add(-211,7);
3830  fCorrelationFunctionsIndices->Add(-321,8);
3831  fCorrelationFunctionsIndices->Add(-2212,9);
3832 
3833  // d) Book all 3p correlation functions:
3835  {
3836  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]
3837  {
3838  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]
3839  {
3840  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]
3841  {
3842  // 3-p correlation functions:
3843  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
3844  f3pCorrelationFunctions[pid1][pid2][pid3]->SetStats(kTRUE);
3845  f3pCorrelationFunctions[pid1][pid2][pid3]->SetFillColor(kBlue-10);
3846  f3pCorrelationFunctions[pid1][pid2][pid3]->SetXTitle("Q_{3}");
3847  f3pCorrelationFunctions[pid1][pid2][pid3]->SetYTitle("C(Q_{3})");
3848  fCorrelationFunctionsSublist[1]->Add(f3pCorrelationFunctions[pid1][pid2][pid3]);
3849  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3850  } // 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]
3851  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3852  } // if(fFill3pCorrelationFunctions)
3853 
3854  // e) Book all 4p correlation functions:
3856  {
3857  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]
3858  {
3859  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]
3860  {
3861  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]
3862  {
3863  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]
3864  {
3865  // 4-p correlation functions:
3866  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
3867  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->SetStats(kTRUE);
3868  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->SetFillColor(kBlue-10);
3869  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->SetXTitle("Q_{4}");
3870  f4pCorrelationFunctions[pid1][pid2][pid3][pid4]->SetYTitle("C(Q_{4})");
3871  fCorrelationFunctionsSublist[2]->Add(f4pCorrelationFunctions[pid1][pid2][pid3][pid4]);
3872  } // 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]
3873  } // 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]
3874  } // 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]
3875  } // 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]
3876  } // if(fFill4pCorrelationFunctions)
3877 
3878 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForCorrelationFunctions()
3879 
3880 //=======================================================================================================================
3881 
3883 {
3884  // Book all the stuff for background objects.
3885 
3886  // a) Book the profile holding all the flags for background objects;
3887  // b) Book all histograms for 2p background;
3888  // c) Book all histograms for 3p background;
3889  // d) Book all histograms for 4p background;
3890  // e) Book buffer objects.
3891 
3892  // a) Book the profile holding all the flags for background objects:
3893  fBackgroundFlagsPro = new TProfile("fBackgroundFlagsPro","Flags and settings for background histograms",4,0,4);
3894  fBackgroundFlagsPro->SetTickLength(-0.01,"Y");
3895  fBackgroundFlagsPro->SetMarkerStyle(25);
3896  fBackgroundFlagsPro->SetLabelSize(0.04);
3897  fBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
3898  fBackgroundFlagsPro->SetStats(kFALSE);
3899  fBackgroundFlagsPro->SetFillColor(kGray);
3900  fBackgroundFlagsPro->SetLineColor(kBlack);
3901  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"fBackgroundOption"); fBackgroundFlagsPro->Fill(0.5,fBackgroundOption);
3902  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(2,"fEstimate2pBackground"); fBackgroundFlagsPro->Fill(1.5,fEstimate2pBackground);
3903  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(3,"fEstimate3pBackground"); fBackgroundFlagsPro->Fill(2.5,fEstimate3pBackground);
3904  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(4,"fEstimate4pBackground"); fBackgroundFlagsPro->Fill(3.5,fEstimate4pBackground);
3906  // 2p:
3907  f2pBackgroundFlagsPro = new TProfile("f2pBackgroundFlagsPro","Flags and settings for 2p background histograms",1,0,1);
3908  f2pBackgroundFlagsPro->SetTickLength(-0.01,"Y");
3909  f2pBackgroundFlagsPro->SetMarkerStyle(25);
3910  f2pBackgroundFlagsPro->SetLabelSize(0.04);
3911  f2pBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
3912  f2pBackgroundFlagsPro->SetStats(kFALSE);
3913  f2pBackgroundFlagsPro->SetFillColor(kGray);
3914  f2pBackgroundFlagsPro->SetLineColor(kBlack);
3915  f2pBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f2pBackgroundFlagsPro->Fill(0.5,0); // TBI dummy at the moment
3917  // 3p:
3918  f3pBackgroundFlagsPro = new TProfile("f3pBackgroundFlagsPro","Flags and settings for 3p background histograms",1,0,1);
3919  f3pBackgroundFlagsPro->SetTickLength(-0.01,"Y");
3920  f3pBackgroundFlagsPro->SetMarkerStyle(25);
3921  f3pBackgroundFlagsPro->SetLabelSize(0.04);
3922  f3pBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
3923  f3pBackgroundFlagsPro->SetStats(kFALSE);
3924  f3pBackgroundFlagsPro->SetFillColor(kGray);
3925  f3pBackgroundFlagsPro->SetLineColor(kBlack);
3926  f3pBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f3pBackgroundFlagsPro->Fill(0.5,0); // TBI dummy at the moment
3928  // 4p:
3929  f4pBackgroundFlagsPro = new TProfile("f4pBackgroundFlagsPro","Flags and settings for 4p background histograms",1,0,1);
3930  f4pBackgroundFlagsPro->SetTickLength(-0.01,"Y");
3931  f4pBackgroundFlagsPro->SetMarkerStyle(25);
3932  f4pBackgroundFlagsPro->SetLabelSize(0.04);
3933  f4pBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
3934  f4pBackgroundFlagsPro->SetStats(kFALSE);
3935  f4pBackgroundFlagsPro->SetFillColor(kGray);
3936  f4pBackgroundFlagsPro->SetLineColor(kBlack);
3937  f4pBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"TBI"); f4pBackgroundFlagsPro->Fill(0.5,0); // TBI dummy at the moment
3939 
3940  //if(!fEstimateBackground){return;} // TBI is this safe?
3941 
3942  const Int_t nParticleSpecies = 5; // Supported at the moment: 0=e,1=mu,2=pi,3=K,4=p
3943  TString sParticles[2*nParticleSpecies] = {"e^{+}","#mu^{+}","#pi^{+}","K^{+}","p^{+}","e^{-}","#mu^{-}","#pi^{-}","K^{-}","p^{-}"};
3944 
3945  // b) Book all histograms for 2p background:
3947  {
3948  // Remark 0: First particle in the pair is always the one with positive charge.
3949  // Remark 1: Diagonal elements are particle/antiparticle pairs.
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=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]
3953  {
3954  // Background:
3955  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
3956  f2pBackground[pid1][pid2]->SetStats(kTRUE);
3957  f2pBackground[pid1][pid2]->SetFillColor(kRed-10);
3958  f2pBackground[pid1][pid2]->SetLineColor(kRed);
3959  f2pBackground[pid1][pid2]->SetXTitle("k");
3960  f2pBackground[pid1][pid2]->SetYTitle("B(k)");
3961  fBackgroundSublist[0]->Add(f2pBackground[pid1][pid2]);
3962  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3963  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3964  } // if(fEstimate2pBackground)
3965 
3966  // c) Book all histograms for 3p background:
3968  {
3969  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]
3970  {
3971  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]
3972  {
3973  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]
3974  {
3975  // Background:
3976  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
3977  f3pBackground[pid1][pid2][pid3]->SetStats(kTRUE);
3978  f3pBackground[pid1][pid2][pid3]->SetFillColor(kRed-10);
3979  f3pBackground[pid1][pid2][pid3]->SetLineColor(kRed);
3980  f3pBackground[pid1][pid2][pid3]->SetXTitle("Q_{3}");
3981  f3pBackground[pid1][pid2][pid3]->SetYTitle("B(Q_{3})");
3982  fBackgroundSublist[1]->Add(f3pBackground[pid1][pid2][pid3]);
3983  } // for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3984  } // for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3985  } // for(Int_t pid=0;pid<2*nParticleSpecies;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3986  } // if(fEstimate2pBackground)
3987 
3988  // d) Book all histograms for 4p background:
3990  {
3991  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]
3992  {
3993  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]
3994  {
3995  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]
3996  {
3997  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]
3998  {
3999  // Background:
4000  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
4001  f4pBackground[pid1][pid2][pid3][pid4]->SetStats(kTRUE);
4002  f4pBackground[pid1][pid2][pid3][pid4]->SetFillColor(kRed-10);
4003  f4pBackground[pid1][pid2][pid3][pid4]->SetLineColor(kRed);
4004  f4pBackground[pid1][pid2][pid3][pid4]->SetXTitle("Q_{4}");
4005  f4pBackground[pid1][pid2][pid3][pid4]->SetYTitle("B(Q_{4})");
4006  fBackgroundSublist[2]->Add(f4pBackground[pid1][pid2][pid3][pid4]);
4007  } // for(Int_t pid4=0;pid4<2*nParticleSpecies;pid4++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
4008  } // for(Int_t pid3=0;pid3<2*nParticleSpecies;pid3++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
4009  } // for(Int_t pid2=0;pid2<2*nParticleSpecies;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
4010  } // for(Int_t pid=0;pid<2*nParticleSpecies;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
4011  } // if(fEstimate2pBackground)
4012 
4013  // e) Book buffer objects:
4014  switch(fBackgroundOption)
4015  {
4016  case 0: // shifting (see documentation in EstimateBackground(...))
4017  for(Int_t me=0;me<3;me++) // [0] is buffer for 1st event; [1] for 2nd, etc.
4018  {
4019  fMixedEvents0[me] = new TClonesArray("AliAODTrack",10000);
4020  }
4021  break; // case 0: // shifting
4022 
4023  case 1: // permutations + binning in vertex z-position
4024  for(Int_t vzr=0;vzr<10;vzr++) // vertez-z range
4025  {
4026  for(Int_t me=0;me<fMaxBufferSize1;me++) // buffer for mixed events
4027  {
4028  fMixedEvents1[vzr][me] = new TClonesArray("AliAODTrack",10000);
4029  fGlobalTracksAOD1[vzr][me] = new TExMap();
4030  }
4031  }
4032  break; // permutations + binning in vertex z-position
4033 
4034  /*
4035  default:
4036  cout<<Form("And the fatal 'fBackgroundOption' value is... %d. Congratulations!!",fBackgroundOption)<<endl;
4037  Fatal(sMethodName.Data(),"switch(fBackgroundOption)");
4038  */
4039  } // switch(fBackgroundOption)
4040 
4041 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForBackground()
4042 
4043 //=======================================================================================================================
4044 
4046 {
4047  // Book all the stuff for buffers.
4048 
4049  // a) Book the profile holding all the flags for buffers;
4050  // b) Book TClonesArray *fChargedParticlesCA[2][10][10000];
4051  // c) Book TExMap *fChargedParticlesEM[10];
4052 
4053  // a) Book the profile holding all the flags for buffers:
4054  fBuffersFlagsPro = new TProfile("fBuffersFlagsPro","Flags and settings for buffers",2,0,2);
4055  fBuffersFlagsPro->SetTickLength(-0.01,"Y");
4056  fBuffersFlagsPro->SetMarkerStyle(25);
4057  fBuffersFlagsPro->SetLabelSize(0.04);
4058  fBuffersFlagsPro->SetLabelOffset(0.02,"Y");
4059  fBuffersFlagsPro->SetStats(kFALSE);
4060  fBuffersFlagsPro->SetFillColor(kGray);
4061  fBuffersFlagsPro->SetLineColor(kBlack);
4062  fBuffersFlagsPro->GetXaxis()->SetBinLabel(1,"fFillBuffers"); fBuffersFlagsPro->Fill(0.5,(Int_t)fFillBuffers);
4063  fBuffersFlagsPro->GetXaxis()->SetBinLabel(2,"fMaxBuffer"); fBuffersFlagsPro->Fill(1.5,fMaxBuffer);
4065 
4066  if(!fFillBuffers){return;} // TBI is this safe? Well, let's hope so for the time being...
4067 
4068  // b) Book TClonesArray *fChargedParticlesCA[2][10][10000]:
4069  for(Int_t am=0;am<1+(Int_t)fProcessBothKineAndReco;am++) // analysis method [0=AOD||ESD,1=MC]
4070  {
4071  for(Int_t e=0;e<fMaxBuffer;e++) // number of events to buffer, max = apparently 10
4072  {
4073  for(Int_t p=0;p<10000;p++) // charged particles, max = apparently 10000
4074  {
4075  fChargedParticlesCA[am][e][p] = NULL; // TBI not sure if I need this really here, this way...
4076  // Example: fChargedParticlesCA[0][4][6] is 6th charged particles in 5th buffered event, in AOD||ESD analysis
4077  } // for(Int_t p=0;p<10000;p++) // charged particles, max = apparently 10000
4078  } // for(Int_t e=0;e<10;e++) // number of events to buffer, max = apparently 10
4079  } // for(Int_t am=0;am<2;am++) // analysis method [0=AOD||ESD,1=MC]
4080 
4081  // c) Book TExMap *fChargedParticlesEM[10]:
4082  for(Int_t e=0;e<fMaxBuffer;e++) // number of events to buffer, max = apparently 10
4083  {
4084  fChargedParticlesEM[e] = new TExMap();
4085  } // for(Int_t e=0;e<fMaxBuffer;e++) // number of events to buffer, max = apparently 10
4086 
4087 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForBuffers()
4088 
4089 //=======================================================================================================================
4090 
4092 {
4093  // Book all the stuff for QA.
4094 
4095  // a) Book the profile holding all the flags for QA;
4096  // b) Book all objects for "QA events";
4097  // c) Book all objects for "QA particles".
4098 
4099  // a) Book the profile holding all the flags for QA:
4100  fQAFlagsPro = new TProfile("fQAFlagsPro","Flags and settings for QA",4,0.,4.);
4101  fQAFlagsPro->SetTickLength(-0.01,"Y");
4102  fQAFlagsPro->SetMarkerStyle(25);
4103  fQAFlagsPro->SetLabelSize(0.04);
4104  fQAFlagsPro->SetLabelOffset(0.02,"Y");
4105  fQAFlagsPro->SetStats(kFALSE);
4106  fQAFlagsPro->SetFillColor(kGray);
4107  fQAFlagsPro->SetLineColor(kBlack);
4108  fQAFlagsPro->GetXaxis()->SetBinLabel(1,"fFillQA"); fQAFlagsPro->Fill(0.5,(Int_t)fFillQA);
4109  fQAFlagsPro->GetXaxis()->SetBinLabel(2,"fBailOutAfterQA"); fQAFlagsPro->Fill(1.5,(Int_t)fBailOutAfterQA);
4110  fQAFlagsPro->GetXaxis()->SetBinLabel(3,"fFillQAEvents"); fQAFlagsPro->Fill(2.5,(Int_t)fFillQAEvents);
4111  fQAFlagsPro->GetXaxis()->SetBinLabel(4,"fFillQAParticles"); fQAFlagsPro->Fill(3.5,(Int_t)fFillQAParticles);
4112  fQAList->Add(fQAFlagsPro);
4113 
4114  if(!fFillQA){return;}
4115 
4116  // b) Book all objects for "QA events":
4117  if(fFillQAEvents)
4118  {
4119  // ...
4120  } // if(fFillQAEvents)
4121 
4122  // c) Book all objects for "QA particles":
4123  if(fFillQAParticles)
4124  {
4125  // c0) Book fQAFilterBitScan:
4126  Int_t nFilterBits = 14;
4127  fQAFilterBitScan = new TH1I("fQAFilterBitScan","fQAFilterBitScan",nFilterBits,0,nFilterBits);
4128  fQAFilterBitScan->SetStats(kFALSE);
4129  fQAFilterBitScan->SetFillColor(kBlue-10);
4130  fQAFilterBitScan->SetXTitle("FilterBit");
4131  fQAFilterBitScan->SetYTitle("# particles");
4132  for(Int_t fb=0;fb<nFilterBits;fb++) // 'fb' is a 'left shifting opetator', i.e. 1<<fb = filterbit
4133  {
4134  fQAFilterBitScan->GetXaxis()->SetBinLabel(fb+1,Form("%d",1<<fb));
4135  }
4137 
4138  // c1) Book fQAIDvsFilterBit:
4139  Int_t nIDsMax = 10000;
4140  fQAIDvsFilterBit = new TH2I("fQAIDvsFilterBit","fQAIDvsFilterBit",nFilterBits,0,nFilterBits,2*nIDsMax,-nIDsMax,nIDsMax);
4141  //fQAIDvsFilterBit->SetStats(kFALSE);
4142  //fQAIDvsFilterBit->SetFillColor(kBlue-10);
4143  fQAIDvsFilterBit->SetXTitle("FilterBit");
4144  fQAIDvsFilterBit->SetYTitle("atrack->GetID()");
4145  for(Int_t fb=0;fb<nFilterBits;fb++) // 'fb' is a 'left shifting opetator', i.e. 1<<fb = filterbit
4146  {
4147  fQAIDvsFilterBit->GetXaxis()->SetBinLabel(fb+1,Form("%d",1<<fb));
4148  }
4150 
4151  // c2) Book fQAParticleHist[2][10][10];
4152  TString sBeforeAfter[2] = {"before","after"};
4153  TString sDistributions[10] = {"p_{T}","#eta","","","","","","","",""}; // Only the ones with non-empty entries are booked and filled, see below if(sDistributions[di].EqualTo("")){continue;}
4154  TString sCuts[10] = {"PID","","","","","","","","",""}; // Only the ones with non-empty entries are booked and filled, see below if(sCuts[ci].EqualTo("")){continue;}
4155  for(Int_t ba=0;ba<2;ba++) // 0="before rain",1="after rain"
4156  {
4157  for(Int_t di=0;di<10;di++) // number of distinct distributions, maximum is 10 (at the moment, as it seems...)
4158  {
4159  if(sDistributions[di].EqualTo("")){continue;}
4160  for(Int_t ci=0;ci<10;ci++) // number of distinct track cuts, maximum is 10 (at the moment, as it seems...)
4161  {
4162  if(sCuts[ci].EqualTo("")){continue;}
4163  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
4164  fQAParticleHist[ba][di][ci]->SetStats(kFALSE);
4165  fQAParticleHist[ba][di][ci]->SetFillColor(kBlue-10);
4166  fQAParticleHist[ba][di][ci]->SetXTitle("TBI");
4167  fQAParticleHist[ba][di][ci]->SetYTitle("TBI");
4168  fQAParticlesList->Add(fQAParticleHist[ba][di][ci]);
4169  } // for(Int_t ci=0;ci<10;ci++) // number of distinct track cuts, maximum is 10 (at the moment)
4170  } // for(Int_t di=0;di<10;di++) // number of distinct distributions, maximum is 10 (at the moment)
4171  } // for(Int_t ba=0;ba<2;ba++) // 0="before rain",1="after rain"
4172  } // if(fFillQAParticles)
4173 
4174 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForQA()
4175 
4176 //=======================================================================================================================
4177 
4179 {
4180  // Book all objects for global track cuts (applied only on "normal" global tracks in AOD).
4181 
4182  // a) Book the profile holding all the flags;
4183  // b) ...
4184 
4185  // a) Book the profile holding all the flags:
4186  fGlobalTrackCutsFlagsPro = new TProfile("fGlobalTrackCutsFlagsPro","Flags and settings for global track cuts",7,0.,7.);
4187  fGlobalTrackCutsFlagsPro->SetTickLength(-0.01,"Y");
4188  fGlobalTrackCutsFlagsPro->SetMarkerStyle(25);
4189  fGlobalTrackCutsFlagsPro->SetLabelSize(0.04);
4190  fGlobalTrackCutsFlagsPro->SetLabelOffset(0.02,"Y");
4191  fGlobalTrackCutsFlagsPro->SetStats(kFALSE);
4192  fGlobalTrackCutsFlagsPro->SetFillColor(kGray);
4193  fGlobalTrackCutsFlagsPro->SetLineColor(kBlack);
4194  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(1,"fApplyGlobalTrackCuts"); fGlobalTrackCutsFlagsPro->Fill(0.5,(Int_t)fApplyGlobalTrackCuts);
4195  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(2,"fPtRange[0]"); fGlobalTrackCutsFlagsPro->Fill(1.5,fPtRange[0]);
4196  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(3,"fPtRange[1]"); fGlobalTrackCutsFlagsPro->Fill(2.5,fPtRange[1]);
4197  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(4,"fEtaRange[0]"); fGlobalTrackCutsFlagsPro->Fill(3.5,fEtaRange[0]);
4198  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(5,"fEtaRange[1]"); fGlobalTrackCutsFlagsPro->Fill(4.5,fEtaRange[1]);
4199  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(6,"fPhiRange[0]"); fGlobalTrackCutsFlagsPro->Fill(5.5,fPhiRange[0]);
4200  fGlobalTrackCutsFlagsPro->GetXaxis()->SetBinLabel(7,"fPhiRange[1]"); fGlobalTrackCutsFlagsPro->Fill(6.5,fPhiRange[1]);
4202 
4203 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForGlobalTrackCuts()
4204 
4205 //=======================================================================================================================
4206 
4208 {
4209  // Filter out unique global tracks in AOD and store them in fGlobalTracksAOD[index].
4210 
4211  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
4212  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
4213  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
4214  // 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
4215  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
4216  // 1: 0: 2.073798 -0.503640 2.935432
4217  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
4218  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
4219  // In fact, this is important, otherwise there is double or even triple counting in some cases.
4220  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
4221  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
4222  // b) All of them apparently have positive ID, i.e. these are global tracks
4223  // c) Clearly, we cannot use TestFilterBit() on them
4224  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
4225  // Remark 3: There is a performance penalty when fGlobalTracksAOD[1] and fGlobalTracksAOD[2] needed for mixed events are calculated.
4226  // Yes, I can get them directly from fGlobalTracksAOD[0], without calling this method for them again. TBI today
4227 
4228  // a) Insanity checks;
4229  // b) Determine the map.
4230 
4231  // a) Insanity checks:
4232  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD, Int_t index)";
4233  if(!fGlobalTracksAOD[index]){Fatal(sMethodName.Data(),"fGlobalTracksAOD[%d]",index);}
4234  if(0 != fGlobalTracksAOD[index]->GetSize()){fGlobalTracksAOD[index]->Delete();} // yes, this method determines mapping from scratch each time
4235 
4236  // b) Determine the map:
4237  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4238  {
4239  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
4240  if(aodTrack)
4241  {
4242  Int_t id = aodTrack->GetID();
4243  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
4244  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
4245  {
4246  fGlobalTracksAOD[index]->Add(id,iTrack); // "key" = id, "value" = iTrack
4247  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
4248  } // if(aodTrack)
4249  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4250 
4251 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD)
4252 
4253 //=======================================================================================================================
4254 
4256 {
4257  // Filter out unique global tracks in AOD and store them in fGlobalTracksAODTEST[index].
4258 
4259  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
4260  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
4261  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
4262  // 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
4263  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
4264  // 1: 0: 2.073798 -0.503640 2.935432
4265  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
4266  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
4267  // In fact, this is important, otherwise there is double or even triple counting in some cases.
4268  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
4269  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
4270  // b) All of them apparently have positive ID, i.e. these are global tracks
4271  // c) Clearly, we cannot use TestFilterBit() on them
4272  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
4273  // Remark 3: There is a performance penalty when fGlobalTracksAODTEST[1] and fGlobalTracksAODTEST[2] needed for mixed events are calculated.
4274  // Yes, I can get them directly from fGlobalTracksAODTEST[0], without calling this method for them again. TBI today
4275 
4276  // a) Insanity checks;
4277  // b) Determine the map.
4278 
4279  // a) Insanity checks:
4280  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAODTEST(AliAODEvent *aAOD, Int_t index)";
4281  if(!fGlobalTracksAODTEST[index]){Fatal(sMethodName.Data(),"fGlobalTracksAODTEST[%d]",index);}
4282  if(0 != fGlobalTracksAODTEST[index]->GetSize()){fGlobalTracksAODTEST[index]->Delete();} // yes, this method determines mapping from scratch each time
4283 
4284  // b) Determine the map:
4285  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4286  {
4287  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
4288  if(aodTrack)
4289  {
4290  Int_t id = aodTrack->GetID();
4291  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
4292  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
4293  {
4294  fGlobalTracksAODTEST[index]->Add(id,iTrack); // "key" = id, "value" = iTrack
4295  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
4296  } // if(aodTrack)
4297  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4298 
4299 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAODTEST(AliAODEvent *aAOD)
4300 
4301 //=======================================================================================================================
4302 
4304 {
4305  // Filter out unique global tracks in AOD and store them in fGlobalTracksAODHA[index].
4306 
4307  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
4308  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
4309  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
4310  // 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
4311  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
4312  // 1: 0: 2.073798 -0.503640 2.935432
4313  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
4314  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
4315  // In fact, this is important, otherwise there is double or even triple counting in some cases.
4316  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
4317  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
4318  // b) All of them apparently have positive ID, i.e. these are global tracks
4319  // c) Clearly, we cannot use TestFilterBit() on them
4320  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
4321  // Remark 3: There is a performance penalty when fGlobalTracksAODHA[1] and fGlobalTracksAODHA[2] needed for mixed events are calculated.
4322  // Yes, I can get them directly from fGlobalTracksAODHA[0], without calling this method for them again. TBI today
4323 
4324  // a) Insanity checks;
4325  // b) Determine the map.
4326 
4327  // a) Insanity checks:
4328  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAODHA(AliAODEvent *aAOD, Int_t index)";
4329  if(!fGlobalTracksAODHA[index]){Fatal(sMethodName.Data(),"fGlobalTracksAODHA[%d]",index);}
4330  if(0 != fGlobalTracksAODHA[index]->GetSize()){fGlobalTracksAODHA[index]->Delete();} // yes, this method determines mapping from scratch each time
4331 
4332  // b) Determine the map:
4333  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4334  {
4335  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
4336  if(aodTrack)
4337  {
4338  Int_t id = aodTrack->GetID();
4339  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
4340  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
4341  {
4342  fGlobalTracksAODHA[index]->Add(id,iTrack); // "key" = id, "value" = iTrack
4343  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
4344  } // if(aodTrack)
4345  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4346 
4347 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAODHA(AliAODEvent *aAOD, Int_t index)
4348 
4349 //=======================================================================================================================
4350 
4352 {
4353  // 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...
4354 
4355  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
4356  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
4357  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
4358  // 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
4359  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
4360  // 1: 0: 2.073798 -0.503640 2.935432
4361  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
4362  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
4363  // In fact, this is important, otherwise there is double or even triple counting in some cases.
4364  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
4365  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
4366  // b) All of them apparently have positive ID, i.e. these are global tracks
4367  // c) Clearly, we cannot use TestFilterBit() on them
4368  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
4369  // Remark 3: Used only for fBackgroundOption = 1, i.e. permutations + binning in vertex z-position
4370 
4371  // a) Insanity checks;
4372  // b) Determine the map.
4373 
4374  // a) Insanity checks:
4375  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD, Int_t indexX, Int_t indexY)";
4376  if(!fGlobalTracksAOD1[indexX][indexY]){Fatal(sMethodName.Data(),"fGlobalTracksAOD[%d][%d]",indexX,indexY);}
4377  if(0 != fGlobalTracksAOD1[indexX][indexY]->GetSize()){fGlobalTracksAOD1[indexX][indexY]->Delete();} // yes, this method determines mapping from scratch each time
4378 
4379  // b) Determine the map:
4380  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4381  {
4382  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
4383  if(aodTrack)
4384  {
4385  Int_t id = aodTrack->GetID();
4386  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
4387  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
4388  {
4389  fGlobalTracksAOD1[indexX][indexY]->Add(id,iTrack); // "key" = id, "value" = iTrack
4390  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
4391  } // if(aodTrack)
4392  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
4393 
4394 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD, Int_t indexX, Int_t indexY)
4395 
4396 //=======================================================================================================================
4397 
4399 {
4400  // 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).
4401 
4402  if(run != fRun) return kFALSE;
4403  else if(bunchCross != fBunchCross) return kFALSE;
4404  else if(orbit != fOrbit) return kFALSE;
4405  else if(period != fPeriod) return kFALSE;
4406 
4407  return kTRUE;
4408 
4409 } // void AliAnalysisTaskMultiparticleFemtoscopy::SpecifiedEvent(UInt_t run, UShort_t bunchCross, UInt_t orbit, UInt_t period)
4410 
4411 //=======================================================================================================================
4412 
4414 {
4415  // Check if the track passes common global track cuts (irrespectively of PID).
4416 
4417  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesGlobalTrackCuts(AliAODTrack *gtrack)";
4418  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
4419 
4420  // To do: add data members and corresponding setters:
4421  // fTPCNclsMin, fTPCNclsMax
4422  // fTPCsignalNMin, fTPCsignalNMax
4423  if(gtrack->Pt()<fPtRange[0]) return kFALSE;
4424  if(gtrack->Pt()>=fPtRange[1]) return kFALSE;
4425  if(gtrack->Eta()<fEtaRange[0]) return kFALSE;
4426  if(gtrack->Eta()>=fEtaRange[1]) return kFALSE;
4427  if(gtrack->Phi()<fPhiRange[0]) return kFALSE;
4428  if(gtrack->Phi()>=fPhiRange[1]) return kFALSE;
4429  //if(gtrack->GetTPCNcls()<70) return kFALSE;
4430  //if(gtrack->GetTPCsignalN()<70) return kFALSE;
4431 
4432  if(fRejectFakeTracks && gtrack->GetLabel()<0) return kFALSE; // TBI is it more precise <=0 ?
4433 
4434  return kTRUE;
4435 
4436 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesGlobalTrackCuts(AliAODTrack *gtrack)
4437 
4438 //=======================================================================================================================
4439 
4441 {
4442  // Check if the track passes common analysis specific track (e.g. TPC-only) cuts.
4443  // Therefore we can cut independetly on global track parameters, and on TPC-only cut parameters.
4444 
4445  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODTrack *atrack)";
4446  if(!atrack){Fatal(sMethodName.Data(),"!atrack");}
4447 
4448  if(!atrack->TestFilterBit(128)) return kFALSE; // TPC-only TBI setter TBI#2 There might be some conflict with FillControlHistogramsNonIdentifiedParticlesFTSF
4449 
4450  // TBI: implement eventually separate cuts for 'atracks' and 'gtracks'
4451  if(atrack->Pt()<fPtRange[0]) return kFALSE;
4452  if(atrack->Pt()>=fPtRange[1]) return kFALSE;
4453  if(atrack->Eta()<fEtaRange[0]) return kFALSE;
4454  if(atrack->Eta()>=fEtaRange[1]) return kFALSE;
4455  if(atrack->Phi()<fPhiRange[0]) return kFALSE;
4456  if(atrack->Phi()>=fPhiRange[1]) return kFALSE;
4457  //if(gtrack->GetTPCNcls()<70) return kFALSE;
4458  //if(gtrack->GetTPCsignalN()<70) return kFALSE;
4459 
4460  if(fRejectFakeTracks && atrack->GetLabel()<0) return kFALSE; // TBI is it more precise <=0 ?
4461 
4462  return kTRUE;
4463 
4464 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODTrack *atrack)
4465 
4466 //=======================================================================================================================
4467 
4469 {
4470  // TBI not validated. this method applies only to MC, make it uniform with AOD.
4471 
4472  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODMCParticle *amcparticle)";
4473  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
4474 
4475  if(TMath::Abs(amcparticle->Charge())<1.e-10) return kFALSE; // skipping neutral particles TBI
4476  if(!amcparticle->IsPhysicalPrimary()) return kFALSE; // skipping secondaries TBI
4477  if(amcparticle->Pt()<0.2) return kFALSE;
4478  if(amcparticle->Pt()>=5.0) return kFALSE;
4479  if(amcparticle->Eta()<-0.8) return kFALSE;
4480  if(amcparticle->Eta()>=0.8) return kFALSE;
4481 
4482  return kTRUE;
4483 
4484 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODMCParticle *amcparticle)
4485 
4486 //=======================================================================================================================
4487 
4489 {
4490  // Check if the event passes common event cuts.
4491 
4492  // TBI: add support for fProcessBothKineAndReco
4493 
4494  // a) Determine Ali{MC,ESD,AOD}Event;
4495 
4496  // a) Determine Ali{MC,ESD,AOD}Event:
4497  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
4498  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
4499  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
4500 
4501  if(aMC && fProcessOnlyKine)
4502  {
4503  // TBI
4504  }
4505  else if(aESD)
4506  {
4507  // TBI
4508  }
4509  else if(aAOD) // shall work both for fProcessBothKineAndReco = kTRUE and fProcessOnlyReco = kTRUE
4510  {
4511  // a) Cuts on AliAODEvent:
4512  if(fRejectEventsWithoutPrimaryVertex && !aAOD->GetPrimaryVertex()) return kFALSE;
4513  if(TMath::Abs(aAOD->GetMagneticField())<fMinMagneticField) return kFALSE;
4515  {
4516  if(aAOD->GetNumberOfTracks() < fMinNumberOfTracks) return kFALSE;
4517  if(aAOD->GetNumberOfTracks() > fMaxNumberOfTracks) return kFALSE;
4518  }
4519  if(fCutOnNumberOfV0s)
4520  {
4521  if(aAOD->GetNumberOfV0s() < fMinNumberOfV0s) return kFALSE;
4522  if(aAOD->GetNumberOfV0s() > fMaxNumberOfV0s) return kFALSE;
4523  }
4525  {
4526  if(aAOD->GetNumberOfCascades() < fMinNumberOfCascades) return kFALSE;
4527  if(aAOD->GetNumberOfCascades() > fMaxNumberOfCascades) return kFALSE;
4528  }
4529 
4530  // b) Cuts on AliAODVertex:
4531  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
4532  if(fCutOnVertexX)
4533  {
4534  if(avtx->GetX() < fMinVertexX) return kFALSE;
4535  if(avtx->GetX() > fMaxVertexX) return kFALSE;
4536  }
4537  if(fCutOnVertexY)
4538  {
4539  if(avtx->GetY() < fMinVertexY) return kFALSE;
4540  if(avtx->GetY() > fMaxVertexY) return kFALSE;
4541  }
4542  if(fCutOnVertexZ)
4543  {
4544  if(avtx->GetZ() < fMinVertexZ) return kFALSE;
4545  if(avtx->GetZ() > fMaxVertexZ) return kFALSE;
4546  }
4548  {
4549  if(avtx->GetNContributors() < fMinNContributors) return kFALSE;
4550  if(avtx->GetNContributors() > fMaxNContributors) return kFALSE;
4551  }
4552  } // else if(aAOD)
4553 
4554  return kTRUE;
4555 
4556 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonEventCuts(AliVEvent *ave)
4557 
4558 //=======================================================================================================================
4559 
4561 {
4562  // Check if the event passes mixed event cuts.
4563 
4564  // Remark: PassesCommonEventCuts() is already checked beforehand in UserExec(), therefore here only some new cuts characteristic only for mixed events shall be implemented.
4565 
4566  // a) Determine Ali{MC,ESD,AOD}Event;
4567 
4568  // a) Determine Ali{MC,ESD,AOD}Event:
4569  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
4570  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
4571  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
4572 
4573  if(aMC)
4574  {
4575  // TBI
4576  }
4577  else if(aESD)
4578  {
4579  // TBI
4580  }
4581  else if(aAOD)
4582  {
4583 
4584  // TBI
4585  /* Check whether the cuts below are already being applied in PassesCommonEventCuts()
4586  if(!aAOD->GetPrimaryVertex()) return kFALSE;
4587  if(TMath::Abs(aAOD->GetMagneticField())<0.001) return kFALSE;
4588  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
4589  //if(TMath::Abs(avtx->GetX())>10.0) return kFALSE;
4590  //if(TMath::Abs(avtx->GetY())>10.0) return kFALSE;
4591  if(TMath::Abs(avtx->GetZ())>10.0) return kFALSE; // TBI setter
4592  if(avtx->GetNContributors()<=2) return kFALSE; // TBI setter
4593  */
4594  }
4595 
4596  return kTRUE;
4597 
4598 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesMixedEventCuts(AliVEvent *ave)
4599 
4600 //=======================================================================================================================
4601 
4603 {
4604  // Calculate correlation functions.
4605 
4606  // a) Insanity checks;
4607  // b) Two nested loops to calculate C(k), just an example; TBI
4608  // c) Three nested loops to calculate C(Q3), just an example; TBI
4609  // d) Four nested loops to calculate C(Q4), just an example. TBI
4610 
4611  // a) Insanity checks:
4612  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::CalculateCorrelationFunctions(AliAODEvent *aAOD)";
4613  if(!aAOD){Fatal(sMethodName.Data(),"!aAOD");}
4614  if(0 == fGlobalTracksAOD[0]->GetSize()){Fatal(sMethodName.Data(),"0 == fGlobalTracksAOD[0]->GetSize()");} // this case shall be already treated in UserExec
4615 
4616  // b) Two nested loops to calculate C(k), just an example:
4617  Int_t nTracks = aAOD->GetNumberOfTracks();
4618  for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
4619  {
4620  AliAODTrack *atrack1 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack1));
4621  // TBI Temporary track insanity checks:
4622  if(!atrack1){Fatal(sMethodName.Data(),"!atrack1");} // TBI keep this for some time, eventually just continue
4623  if(atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4624  if(atrack1->TestFilterBit(128) && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->TestFiletrBit(128) && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4625  if(!PassesCommonTrackCuts(atrack1)){continue;} // TBI re-think
4626  // Corresponding AOD global track:
4627  Int_t id1 = atrack1->GetID();
4628  AliAODTrack *gtrack1 = dynamic_cast<AliAODTrack*>(id1>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id1)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id1+1))));
4629  if(!gtrack1){Fatal(sMethodName.Data(),"!gtrack1");} // TBI keep this for some time, eventually just continue
4630  Int_t gid1 = (id1>=0 ? id1 : -(id1+1)); // ID of corresponding global track
4631  // Common track selection criteria for all "normal" global tracks:
4632  if(!PassesGlobalTrackCuts(gtrack1)){continue;}
4633 
4634  // Loop over the 2nd particle:
4635  for(Int_t iTrack2=iTrack1+1;iTrack2<nTracks;iTrack2++)
4636  {
4637  AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack2));
4638  // TBI Temporary track insanity checks:
4639  if(!atrack2){Fatal(sMethodName.Data(),"!atrack2");} // TBI keep this for some time, eventually just continue
4640  if(atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4641  if(atrack2->TestFilterBit(128) && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->TestFiletrBit(128) && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4642  if(!PassesCommonTrackCuts(atrack2)){continue;} // TBI re-think
4643  // Corresponding AOD global track:
4644  Int_t id2 = atrack2->GetID();
4645  AliAODTrack *gtrack2 = dynamic_cast<AliAODTrack*>(id2>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id2)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id2+1))));
4646  if(!gtrack2){Fatal(sMethodName.Data(),"!gtrack2");} // TBI keep this for some time, eventually just continue
4647  Int_t gid2 = (id2>=0 ? id2 : -(id2+1)); // ID of corresponding global track
4648  if(gid1==gid2){continue;} // Eliminate self-correlations:
4649 
4650  // Common track selection criteria for all "normal" global tracks:
4651  if(!PassesGlobalTrackCuts(gtrack2)){continue;}
4652 
4653  // 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:
4654  AliAODTrack *agtrack1 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4655  AliAODTrack *agtrack2 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4657  {
4658  agtrack1 = gtrack1;
4659  agtrack2 = gtrack2;
4660  }
4661  else
4662  {
4663  agtrack1 = atrack1;
4664  agtrack2 = atrack2;
4665  } // if(fFillControlHistogramsWithGlobalTrackInfo)
4666  if(!agtrack1){Fatal(sMethodName.Data(),"!agtrack1");}
4667  if(!agtrack2){Fatal(sMethodName.Data(),"!agtrack2");}
4668 
4669  // 1.) Same particle species:
4670 
4671  // a) pion-pion:
4672  // a1) pi+pi+ [2][2]:
4673  if(Pion(gtrack1,1,kTRUE) && Pion(gtrack2,1,kTRUE))
4674  {
4675  fCorrelationFunctions[2][2]->Fill(Q2(agtrack1,agtrack2));
4676  }
4677  // a2) pi-pi- [7][7]:
4678  if(Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,-1,kTRUE))
4679  {
4680  fCorrelationFunctions[7][7]->Fill(Q2(agtrack1,agtrack2));
4681  }
4682  // a3) pi+pi- || pi-pi+ [2][7]:
4683  if((Pion(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE)) || (Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,1,kTRUE)))
4684  {
4685  fCorrelationFunctions[2][7]->Fill(Q2(agtrack1,agtrack2));
4686  }
4687 
4688  // b) kaon-kaon:
4689  // b1) K+K+ [3][3]:
4690  if(Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE))
4691  {
4692  fCorrelationFunctions[3][3]->Fill(Q2(agtrack1,agtrack2));
4693  }
4694  // b2) K-K- [8][8]:
4695  if(Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
4696  {
4697  fCorrelationFunctions[8][8]->Fill(Q2(agtrack1,agtrack2));
4698  }
4699  // b3) K+K- || K-K+ [3][8]:
4700  if((Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE)) || (Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,1,kTRUE)))
4701  {
4702  fCorrelationFunctions[3][8]->Fill(Q2(agtrack1,agtrack2));
4703  }
4704 
4705  // c) proton-proton:
4706  // c1) p+p+ [4][4]:
4707  if(Proton(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
4708  {
4709  fCorrelationFunctions[4][4]->Fill(Q2(agtrack1,agtrack2));
4710  }
4711  // c2) p-p- [9][9]:
4712  if(Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4713  {
4714  fCorrelationFunctions[9][9]->Fill(Q2(agtrack1,agtrack2));
4715  }
4716  // c3) p+p- || p-p+ [4][9]:
4717  if((Proton(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE)) || (Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,1,kTRUE)))
4718  {
4719  fCorrelationFunctions[4][9]->Fill(Q2(agtrack1,agtrack2));
4720  }
4721 
4722  // 2.) Mixed particle species:
4723  // a) pion-kaon
4724  // a1) pi+K+ [2][3]:
4725  if(Pion(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE))
4726  {
4727  fCorrelationFunctions[2][3]->Fill(Q2(agtrack1,agtrack2));
4728  }
4729  // a2) pi+K- [2][8]
4730  if(Pion(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
4731  {
4732  fCorrelationFunctions[2][8]->Fill(Q2(agtrack1,agtrack2));
4733  }
4734  // a3) K+pi- [3][7]
4735  if(Kaon(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE))
4736  {
4737  fCorrelationFunctions[3][7]->Fill(Q2(agtrack1,agtrack2));
4738  }
4739  // a4) pi-K- [7][8]
4740  if(Pion(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
4741  {
4742  fCorrelationFunctions[7][8]->Fill(Q2(agtrack1,agtrack2));
4743  }
4744  // b) pion-proton
4745  // b1) pi+p+ [2][4]:
4746  if(Pion(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
4747  {
4748  fCorrelationFunctions[2][4]->Fill(Q2(agtrack1,agtrack2));
4749  }
4750  // b2) pi+p- [2][9]
4751  if(Pion(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4752  {
4753  fCorrelationFunctions[2][9]->Fill(Q2(agtrack1,agtrack2));
4754  }
4755  // b3) p+pi- [4][7]
4756  if(Proton(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE))
4757  {
4758  fCorrelationFunctions[4][7]->Fill(Q2(agtrack1,agtrack2));
4759  }
4760  // b4) pi-p- [7][9]
4761  if(Pion(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4762  {
4763  fCorrelationFunctions[7][9]->Fill(Q2(agtrack1,agtrack2));
4764  }
4765  // c) kaon-proton
4766  // c1) K+p+ [3][4]:
4767  if(Kaon(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
4768  {
4769  fCorrelationFunctions[3][4]->Fill(Q2(agtrack1,agtrack2));
4770  }
4771  // c2) K+p- [3][9]
4772  if(Kaon(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4773  {
4774  fCorrelationFunctions[3][9]->Fill(Q2(agtrack1,agtrack2));
4775  }
4776  // c3) p+K- [4][8]
4777  if(Proton(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
4778  {
4779  fCorrelationFunctions[4][8]->Fill(Q2(agtrack1,agtrack2));
4780  }
4781  // c4) K-p- [8][9]
4782  if(Kaon(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
4783  {
4784  fCorrelationFunctions[8][9]->Fill(Q2(agtrack1,agtrack2));
4785  }
4786  } // for(Int_t iTrack2=0;iTrack2<nTracks;iTrack2++)
4787  } // for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
4788 
4789  // c) Three nested loops to calculate C(Q3), just an example; TBI
4791 
4792  // d) Four nested loops to calculate C(Q4), just an example. TBI
4794 
4795 } // void AliAnalysisTaskMultiparticleFemtoscopy::CalculateCorrelationFunctions(AliAODEvent *aAOD)
4796 
4797 //=======================================================================================================================
4798 
4800 {
4801  // Calculate 3-particle correlation functions.
4802 
4803  // a) Insanity checks;
4804  // b) Three nested loops to calculate C(Q3), just an example.
4805 
4806  // a) Insanity checks:
4807  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::Calculate3pCorrelationFunctions(AliAODEvent *aAOD)";
4808  if(!aAOD){Fatal(sMethodName.Data(),"!aAOD");}
4809  if(0 == fGlobalTracksAOD[0]->GetSize()){Fatal(sMethodName.Data(),"0 == fGlobalTracksAOD[0]->GetSize()");} // this case shall be already treated in UserExec
4810 
4811  // b) Three nested loops to calculate C(Q3), just an example:
4812  Int_t nTracks = aAOD->GetNumberOfTracks();
4813  for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
4814  {
4815  AliAODTrack *atrack1 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack1));
4816  // TBI Temporary track insanity checks:
4817  if(!atrack1){Fatal(sMethodName.Data(),"!atrack1");} // TBI keep this for some time, eventually just continue
4818  if(atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4819  if(atrack1->TestFilterBit(128) && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->TestFiletrBit(128) && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4820  if(!PassesCommonTrackCuts(atrack1)){continue;} // TBI re-think
4821  // Corresponding AOD global track:
4822  Int_t id1 = atrack1->GetID();
4823  AliAODTrack *gtrack1 = dynamic_cast<AliAODTrack*>(id1>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id1)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id1+1))));
4824  if(!gtrack1){Fatal(sMethodName.Data(),"!gtrack1");} // TBI keep this for some time, eventually just continue
4825  Int_t gid1 = (id1>=0 ? id1 : -(id1+1)); // ID of corresponding global track
4826  // Common track selection criteria for all "normal" global tracks:
4827  if(!PassesGlobalTrackCuts(gtrack1)){continue;}
4828 
4829  // Loop over the 2nd particle:
4830  for(Int_t iTrack2=0;iTrack2<nTracks;iTrack2++)
4831  {
4832  if(iTrack2<=iTrack1){continue;} // Eliminate self-evident self-correlations, and permutations as well
4833  AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack2));
4834  // TBI Temporary track insanity checks:
4835  if(!atrack2){Fatal(sMethodName.Data(),"!atrack2");} // TBI keep this for some time, eventually just continue
4836  if(atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4837  if(atrack2->TestFilterBit(128) && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->TestFiletrBit(128) && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4838  if(!PassesCommonTrackCuts(atrack2)){continue;} // TBI re-think
4839  // Corresponding AOD global track:
4840  Int_t id2 = atrack2->GetID();
4841  AliAODTrack *gtrack2 = dynamic_cast<AliAODTrack*>(id2>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id2)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id2+1))));
4842  if(!gtrack2){Fatal(sMethodName.Data(),"!gtrack2");} // TBI keep this for some time, eventually just continue
4843  Int_t gid2 = (id2>=0 ? id2 : -(id2+1)); // ID of corresponding global track
4844  if(gid1==gid2){continue;} // Eliminate not-so-evident self-correlations:
4845  // Common track selection criteria for all "normal" global tracks:
4846  if(!PassesGlobalTrackCuts(gtrack2)){continue;}
4847 
4848  // Loop over the 3rd particle:
4849  for(Int_t iTrack3=0;iTrack3<nTracks;iTrack3++)
4850  {
4851  if(iTrack3<=iTrack2 || iTrack3<=iTrack1){continue;} // Eliminate self-evident self-correlations, and permutations as well
4852  AliAODTrack *atrack3 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack3));
4853  // TBI Temporary track insanity checks:
4854  if(!atrack3){Fatal(sMethodName.Data(),"!atrack3");} // TBI keep this for some time, eventually just continue
4855  if(atrack3->GetID()>=0 && atrack3->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack3->GetID()>=0 && atrack3->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4856  if(atrack3->TestFilterBit(128) && atrack3->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack3->TestFiletrBit(128) && atrack3->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
4857  if(!PassesCommonTrackCuts(atrack3)){continue;} // TBI re-think
4858  // Corresponding AOD global track:
4859  Int_t id3 = atrack3->GetID();
4860  AliAODTrack *gtrack3 = dynamic_cast<AliAODTrack*>(id3>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id3)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id3+1))));
4861  if(!gtrack3){Fatal(sMethodName.Data(),"!gtrack3");} // TBI keep this for some time, eventually just continue
4862  Int_t gid3 = (id3>=0 ? id3 : -(id3+1)); // ID of corresponding global track
4863  if(gid3==gid2 || gid3==gid1){continue;} // Eliminate not-so-evident self-correlations
4864  // Common track selection criteria for all "normal" global tracks:
4865  if(!PassesGlobalTrackCuts(gtrack3)){continue;}
4866 
4867  // 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:
4868  AliAODTrack *agtrack1 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4869  AliAODTrack *agtrack2 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4870  AliAODTrack *agtrack3 = NULL; // either 'atrack' or 'gtrack', depending on the flag fFillControlHistogramsWithGlobalTrackInfo. By default, agtrack = 'atrack'
4872  {
4873  agtrack1 = gtrack1;
4874  agtrack2 = gtrack2;
4875  agtrack3 = gtrack3;
4876  }
4877  else
4878  {
4879  agtrack1 = atrack1;
4880  agtrack2 = atrack2;
4881  agtrack3 = atrack3;
4882  } // if(fFillControlHistogramsWithGlobalTrackInfo)
4883  if(!agtrack1){Fatal(sMethodName.Data(),"!agtrack1");}
4884  if(!agtrack2){Fatal(sMethodName.Data(),"!agtrack2");}
4885  if(!agtrack3){Fatal(sMethodName.Data(),"!agtrack3");}
4886 
4887  // Cases of interest:
4888  // a) Same species and same charge;
4889  // b) Same species but different charge combinations (modulo permutations);
4890  // c) Two pions + something else (modulo permutations);
4891  // d) Two nucleons + something else (modulo permutations).
4892 
4893  // a) Same species:
4894  // pi+pi+pi+
4895  if(Pion(gtrack1,1,kTRUE) && Pion(gtrack2,1,kTRUE) && Pion(gtrack3,1,kTRUE))
4896  {
4897  f3pCorrelationFunctions[2][2][2]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4898  }
4899  // pi-pi-pi-
4900  if(Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,-1,kTRUE) && Pion(gtrack3,-1,kTRUE))
4901  {
4902  f3pCorrelationFunctions[7][7][7]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4903  }
4904  // K+K+K+
4905  if(Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE) && Kaon(gtrack3,1,kTRUE))
4906  {
4907  f3pCorrelationFunctions[3][3][3]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4908  }
4909  // K-K-K-
4910  if(Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE) && Kaon(gtrack3,-1,kTRUE))
4911  {
4912  f3pCorrelationFunctions[8][8][8]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4913  }
4914  // p+p+p+
4915  if(Proton(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE) && Proton(gtrack3,1,kTRUE))
4916  {
4917  f3pCorrelationFunctions[4][4][4]->Fill(Q3(agtrack1,agtrack2,agtrack3));
4918  }
4919  // p-p-p-
4920  if(Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRU