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