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