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