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