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