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