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