AliPhysics  251aa1e (251aa1e)
 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(3),
54  fProcessBothKineAndReco(kFALSE),
55  fProcessOnlyKine(kFALSE),
56  fRejectFakeTracks(kTRUE),
57  fMC(NULL),
58  // 1.) Control histograms:
59  fControlHistogramsList(NULL),
60  fControlHistogramsFlagsPro(NULL),
61  fFillControlHistograms(kTRUE), // re-think
62  fControlHistogramsEventList(NULL),
63  fControlHistogramsEventFlagsPro(NULL),
64  fFillControlHistogramsEvent(kFALSE),
65  fGetNumberOfTracksHist(NULL),
66  fGetNumberOfGlobalTracksHist(NULL),
67  fGetNumberOfV0sHist(NULL),
68  fGetNumberOfCascadesHist(NULL),
69  fGetMagneticFieldHist(NULL),
70  fGetEventTypeHist(NULL),
71  fGetCentralityHist(NULL),
72  fGetNContributorsHist(NULL),
73  fGetChi2perNDFHist(NULL),
74  fGetNDaughtersHist(NULL),
75  fControlHistogramsNonIdentifiedParticlesList(NULL),
76  fControlHistogramsNonIdentifiedParticlesFlagsPro(NULL),
77  fFillControlHistogramsNonIdentifiedParticles(kFALSE),
78  fChargeHist(NULL),
79  fGetTPCNclsHist(NULL),
80  fGetTPCsignalNHist(NULL),
81  fGetITSNclsHist(NULL),
82  fdEdxVsPtHist(NULL),
83  fPtHist(NULL),
84  fEtaHist(NULL),
85  fPhiHist(NULL),
86  fMassHist(NULL),
87  fGetFilterMap(NULL),
88  fGetPdgCode(NULL),
89  fControlHistogramsIdentifiedParticlesList(NULL),
90  fControlHistogramsIdentifiedParticlesFlagsPro(NULL),
91  fFillControlHistogramsIdentifiedParticles(kFALSE),
92  fControlHistogramsV0sList(NULL),
93  fControlHistogramsV0sFlagsPro(NULL),
94  fFillControlHistogramsV0s(kFALSE),
95  fGetNProngsHist(NULL),
96  fMassK0ShortHist(NULL),
97  fMassLambdaHist(NULL),
98  fMassAntiLambdaHist(NULL),
99  fOpenAngleV0Hist(NULL),
100  fRadiusV0Hist(NULL),
101  fDcaV0ToPrimVertexHist(NULL),
102  fMomV0XHist(NULL),
103  fMomV0YHist(NULL),
104  fMomV0ZHist(NULL),
105  fPtV0Hist(NULL),
106  fPseudoRapV0Hist(NULL),
107  fPAHist(NULL),
108  // 2.) Event-by-event histograms:
109  fEBEHistogramsList(NULL),
110  fEBEObjectsFlagsPro(NULL),
111  //fFillEBEHistograms(kTRUE),
112  fUniqueIDHistEBE(NULL),
113  // 3.) Correlation functions:
114  fCorrelationFunctionsList(NULL),
115  fCorrelationFunctionsFlagsPro(NULL),
116  fFillCorrelationFunctions(kFALSE),
117  fNormalizeCorrelationFunctions(kFALSE),
118  fCorrelationFunctionsIndices(NULL),
119  // 4.) Background:
120  fBackgroundList(NULL),
121  fBackgroundFlagsPro(NULL),
122  fEstimateBackground(kFALSE),
123  // *.) Online monitoring:
124  fOnlineMonitoring(kFALSE),
125  fUpdateOutputFile(kFALSE),
126  fUpdateFrequency(-44),
127  fUpdateWhichOutputFile(NULL),
128  fMaxNumberOfEvents(-44),
129  // *.) Debugging:
130  fDoSomeDebugging(kFALSE),
131  fWaitForSpecifiedEvent(kFALSE),
132  fRun(0),
133  fBunchCross(0),
134  fOrbit(0),
135  fPeriod(0)
136  {
137  // Constructor.
138 
139  AliDebug(2,"AliAnalysisTaskMultiparticleFemtoscopy::AliAnalysisTaskMultiparticleFemtoscopy(const char *name, Bool_t useParticleWeights)");
140 
141  // Base list:
142  fHistList = new TList();
143  fHistList->SetName("GMlist"); // TBI
144  fHistList->SetOwner(kTRUE);
145 
146  // Initialize all arrays:
147  this->InitializeArrays();
148  this->InitializeArraysForControlHistograms();
149  this->InitializeArraysForEBEObjects();
150  this->InitializeArraysForCorrelationFunctions();
151  this->InitializeArraysForBackground();
152 
153  // Define input and output slots here
154  // Input slot #0 works with an AliFlowEventSimple
155  //DefineInput(0, AliFlowEventSimple::Class());
156  // Input slot #1 is needed for the weights input file:
157  //if(useParticleWeights)
158  //{
159  // DefineInput(1, TList::Class());
160  //}
161  // Output slot #0 is reserved
162  // Output slot #1 writes into a TList container
163 
164  DefineOutput(1, TList::Class());
165 
166  if(useParticleWeights)
167  {
168  // TBI add support eventually for particle weights
169  }
170 
171 } // AliAnalysisTaskMultiparticleFemtoscopy::AliAnalysisTaskMultiparticleFemtoscopy(const char *name, Bool_t useParticleWeights):
172 
173 //================================================================================================================
174 
177  fHistList(NULL),
178  fAnalysisType(NULL),
179  fPIDResponse(NULL),
180  fMaxNoGlobalTracksAOD(3),
181  fProcessBothKineAndReco(kFALSE),
182  fProcessOnlyKine(kFALSE),
183  fRejectFakeTracks(kTRUE),
184  fMC(NULL),
185  // 1.) Control histograms:
186  fControlHistogramsList(NULL),
187  fControlHistogramsFlagsPro(NULL),
188  fFillControlHistograms(kTRUE), // TBI re-think
189  fControlHistogramsEventList(NULL),
190  fControlHistogramsEventFlagsPro(NULL),
191  fFillControlHistogramsEvent(kFALSE),
192  fGetNumberOfTracksHist(NULL),
193  fGetNumberOfGlobalTracksHist(NULL),
194  fGetNumberOfV0sHist(NULL),
195  fGetNumberOfCascadesHist(NULL),
196  fGetMagneticFieldHist(NULL),
197  fGetEventTypeHist(NULL),
198  fGetCentralityHist(NULL),
199  fGetNContributorsHist(NULL),
200  fGetChi2perNDFHist(NULL),
201  fGetNDaughtersHist(NULL),
202  fControlHistogramsNonIdentifiedParticlesList(NULL),
203  fControlHistogramsNonIdentifiedParticlesFlagsPro(NULL),
204  fFillControlHistogramsNonIdentifiedParticles(kFALSE),
205  fChargeHist(NULL),
206  fGetTPCNclsHist(NULL),
207  fGetTPCsignalNHist(NULL),
208  fGetITSNclsHist(NULL),
209  fdEdxVsPtHist(NULL),
210  fPtHist(NULL),
211  fEtaHist(NULL),
212  fPhiHist(NULL),
213  fMassHist(NULL),
214  fGetFilterMap(NULL),
215  fGetPdgCode(NULL),
216  fControlHistogramsIdentifiedParticlesList(NULL),
217  fControlHistogramsIdentifiedParticlesFlagsPro(NULL),
218  fFillControlHistogramsIdentifiedParticles(kFALSE),
219  fControlHistogramsV0sList(NULL),
220  fControlHistogramsV0sFlagsPro(NULL),
221  fFillControlHistogramsV0s(kFALSE),
222  fGetNProngsHist(NULL),
223  fMassK0ShortHist(NULL),
224  fMassLambdaHist(NULL),
225  fMassAntiLambdaHist(NULL),
226  fOpenAngleV0Hist(NULL),
227  fRadiusV0Hist(NULL),
228  fDcaV0ToPrimVertexHist(NULL),
229  fMomV0XHist(NULL),
230  fMomV0YHist(NULL),
231  fMomV0ZHist(NULL),
232  fPtV0Hist(NULL),
233  fPseudoRapV0Hist(NULL),
234  fPAHist(NULL),
235  // 2.) Event-by-event histograms:
236  fEBEHistogramsList(NULL),
237  fEBEObjectsFlagsPro(NULL),
238  //fFillEBEHistograms(kFALSE),
239  fUniqueIDHistEBE(NULL),
240  // 3.) Correlation functions:
241  fCorrelationFunctionsList(NULL),
242  fCorrelationFunctionsFlagsPro(NULL),
243  fFillCorrelationFunctions(kFALSE),
244  fNormalizeCorrelationFunctions(kFALSE),
245  fCorrelationFunctionsIndices(NULL),
246  // 4.) Background:
247  fBackgroundList(NULL),
248  fBackgroundFlagsPro(NULL),
249  fEstimateBackground(kFALSE),
250  // *.) Online monitoring:
251  fOnlineMonitoring(kFALSE),
252  fUpdateOutputFile(kFALSE),
253  fUpdateFrequency(-44),
254  fUpdateWhichOutputFile(NULL),
255  fMaxNumberOfEvents(-44),
256  // *.) Debugging:
257  fDoSomeDebugging(kFALSE),
258  fWaitForSpecifiedEvent(kFALSE),
259  fRun(0),
260  fBunchCross(0),
261  fOrbit(0),
262  fPeriod(0)
263 {
264  // Dummy constructor.
265 
266  AliDebug(2,"AliAnalysisTaskMultiparticleFemtoscopy::AliAnalysisTaskMultiparticleFemtoscopy()");
267 
268 } // AliAnalysisTaskMultiparticleFemtoscopy::AliAnalysisTaskMultiparticleFemtoscopy():
269 
270 //================================================================================================================
271 
273 {
274  // Destructor.
275 
276  if(fHistList) delete fHistList;
277  if(fPIDResponse) delete fPIDResponse;
278 
279  for(Int_t index=0;index<fMaxNoGlobalTracksAOD;index++)
280  {
281  if(fGlobalTracksAOD[index]) delete fGlobalTracksAOD[index];
282  }
283 
284 } // AliAnalysisTaskMultiparticleFemtoscopy::~AliAnalysisTaskMultiparticleFemtoscopy()
285 
286 //================================================================================================================
287 
289 {
290  // Called at every worker node to initialize.
291 
292  // a) Insanity checks;
293  // b) Trick to avoid name clashes, part 1;
294  // c) Book and nest all lists;
295  // d) Book all objects;
296  // e) Set all flags;
297  // f) Trick to avoid name clashes, part 2.
298 
299  // a) Insanity checks:
301 
302  // b) Trick to avoid name clashes, part 1:
303  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
304  TH1::AddDirectory(kFALSE);
305 
306  // c) Book and nest all lists:
307  this->BookAndNestAllLists();
308 
309  // d) Book all objects:
310  this->BookEverything();
315 
316  // e) Set all flags:
317  // ...
318 
319  // f) Trick to avoid name clashes, part 2:
320  TH1::AddDirectory(oldHistAddStatus);
321 
322  PostData(1,fHistList);
323 
324 } // void AliAnalysisTaskMultiparticleFemtoscopy::UserCreateOutputObjects()
325 
326 //================================================================================================================
327 
329 {
330  // Main loop (called for each event).
331 
332  // a) Determine Ali{MC,ESD,AOD}Event;
333  // b) Insanity checks;
334  // c) Start analysis over MC, ESD or AODs;
335  // d) Reset event-by-event objects;
336  // e) PostData;
337  // f) Online monitoring.
338 
339  // a) Determine Ali{MC,ESD,AOD}Event:
340  AliMCEvent *aMC = MCEvent(); // from TaskSE
341  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
342  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
343 
344  // b) Insanity checks:
345  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::UserExec(Option_t *)";
346  if(aMC && aAOD && fProcessBothKineAndReco){*fAnalysisType="MC_AOD";}
347  else if(aMC && fProcessOnlyKine){*fAnalysisType="MC";} // TBI this is a bit shaky, isnt' it...
348  else if(aESD){*fAnalysisType="ESD";}
349  else if(aAOD){*fAnalysisType="AOD";}
350  else{Fatal(sMethodName.Data(),"Neither MC, nor ESD, nor AOD, nor MC_AOD. Okay...");}
352 
353  // c) Start analysis over MC, ESD or AOD:
354  if(fProcessBothKineAndReco) // TBI: this is a bit shaky, isn't it...
355  {
356  fMC = aMC;
357  AOD(aAOD);
358  }
359  else if(aMC && fProcessOnlyKine){MC(aMC);}
360  else if(aESD){ESD(aESD);}
361  else if(aAOD){AOD(aAOD);}
362  else{Fatal(sMethodName.Data(),"Neither MC, nor ESD, nor AOD. Okay...");} // TBI I have this check already
363 
364  // d) Reset event-by-event objects:
365  this->ResetEBEObjects();
366 
367  // e) PostData:
368  PostData(1,fHistList);
369 
370  // f) Online monitoring:
372 
373 } // void AliAnalysisTaskMultiparticleFemtoscopy::UserExec(Option_t *)
374 
375 //================================================================================================================
376 
378 {
379  // At the end of the day...
380 
381  // a) Get pointer to the grandmother of all lists "fHistList";
382  // b) Get pointers to all objects in "fHistList";
383  // c) Normalize correlation functions;
384  // d) Dump the results.
385 
386  // a) Get pointer to the grandmother of all lists:
387  fHistList = (TList*)GetOutputData(1);
388  if(!fHistList){exit(1);}
389 
390  // b) Get pointers to all objects in "fHistList":
392 
393  // c) Normalize correlation functions:
395 
396  // d) Dump the results:
397  //TDirectoryFile *df = new TDirectoryFile("outputMPFanalysis","");
398  //df->Add(fHistList);
399  TFile *f = new TFile("AnalysisResults.root","RECREATE");
400  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
401 
402  delete f;
403 
404 } // end of void AliAnalysisTaskMultiparticleFemtoscopy::Terminate(Option_t *)
405 
406 //================================================================================================================
407 
409 {
410  // Normalize correlation functions with the background.
411 
412  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]
413  {
414  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]
415  {
416  if(fCorrelationFunctions[pid1][pid2] && fCorrelationFunctions[pid1][pid2]->GetEntries() > 0 &&
417  fBackground[pid1][pid2] && fBackground[pid1][pid2]->GetEntries() > 0)
418  {
419  fCorrelationFunctions[pid1][pid2]->Divide(fBackground[pid1][pid2]);
420  } // if(..)
421  } // 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]
422  } // 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]
423 
424 } // void AliAnalysisTaskMultiparticleFemtoscopy::NormalizeCorrelationFunctions()
425 
426 //================================================================================================================
427 
429 {
430  // Fill control histograms for global event observables.
431 
432  // TBI: add support for fProcessBothKineAndReco
433 
434  // To do:
435  // 1) Add support for MC and ESD.
436 
437  // a) Determine Ali{MC,ESD,AOD}Event;
438  // b) ...
439 
440  // a) Determine Ali{MC,ESD,AOD}Event:
441  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
442  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
443  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
444 
445  if(aMC)
446  {
447  // MC event:
448  fGetNumberOfTracksHist->Fill(aMC->GetNumberOfTracks());
449  } // if(aMC)
450  else if(aESD)
451  {
452  // TBI
453  } // else if(aESD)
454  else if(aAOD)
455  {
456  // AOD event:
457  fGetNumberOfTracksHist->Fill(aAOD->GetNumberOfTracks()); // TBI not all tracks are unique
458  fGetNumberOfGlobalTracksHist->Fill(fGlobalTracksAOD[0]->GetSize()); // this is then my multiplicity...
459  fGetNumberOfV0sHist->Fill(aAOD->GetNumberOfV0s()); // TBI some V0s share the daughter
460  fGetNumberOfCascadesHist->Fill(aAOD->GetNumberOfCascades()); // TBI validate
461  fGetMagneticFieldHist->Fill(aAOD->GetMagneticField());
462  fGetEventTypeHist->Fill(aAOD->GetEventType());
463  if(aAOD->GetCentrality())
464  {
465  fGetCentralityHist->Fill(aAOD->GetCentrality()->GetCentralityPercentile("V0M")); // TBI validate this
466  }
467  // AOD primary vertex:
468  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
469  fVertexXYZ[0]->Fill(avtx->GetX());
470  fVertexXYZ[1]->Fill(avtx->GetY());
471  fVertexXYZ[2]->Fill(avtx->GetZ());
472  fGetNContributorsHist->Fill(avtx->GetNContributors());
473  fGetChi2perNDFHist->Fill(avtx->GetChi2perNDF());
474  fGetNDaughtersHist->Fill(avtx->GetNDaughters());
475  } // else if(aAOD)
476 
477 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsEvent(AliVEvent *ave)
478 
479 //================================================================================================================
480 
482 {
483  // Fill control histograms for particles.
484  // Remark: The idea is to have one loop over the particles and to fill everything which needs to be filled.
485 
486  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsParticle(AliVEvent *ave)";
487 
488  // a) Determine Ali{MC,ESD,AOD}Event:
489  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
490  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
491  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
492 
493  if(aMC)
494  {
495  // MC tracks:
496  for(Int_t iTrack=0;iTrack<aMC->GetNumberOfTracks();iTrack++)
497  {
498  // a) Determine "amctrack";
499  // b) Insanity checks for "amctrack";
500  // c) Fill the control histograms for non-identified particles;
501  // d) Fill the control histograms for identified particles.
502 
503  // a) Determine "amctrack":
504  AliAODMCParticle *amcparticle = (AliAODMCParticle*)aMC->GetTrack(iTrack);
505 
506  // b) Insanity checks for "amctrack":
507  if(!amcparticle){Fatal(sMethodName.Data(),"!amctrack");}
508 
509  // c) Fill the control histograms for non-identified particles:
511 
512  // d) Fill the control histograms for identified particles:
514 
515  } // for(Int_t iTrack=0;iTrack<aMC->GetNumberOfTracks();iTrack++)
516  } // if(aMC)
517 
518  //=================================================================================================================
519 
520  else if(aESD)
521  {
522  // TBI
523  } // else if(aESD)
524 
525  //=================================================================================================================
526 
527  else if(aAOD)
528  {
529  // AOD tracks:
530  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
531  {
532  // a) Determine "atrack";
533  // b) Insanity checks for "atrack";
534  // c) Determine the corresponding "gtrack" (a.k.a. "normal global" track);
535  // d) Fill the control histograms for non-identified particles;
536  // e) Fill the control histograms for identified particles;
537 
538  // a) Determine "atrack":
539  AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
540 
541  // b) Insanity checks for "atrack":
542  if(!atrack){Fatal(sMethodName.Data(),"!atrack");}
543  if(atrack->GetID()>=0 && atrack->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack->GetID()>=0 && atrack->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
544  if(atrack->TestFilterBit(128) && atrack->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack->TestFiletrBit(128) && atrack->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
545  // if(0==atrack->GetFilterMap())
546  // TBI a vast majority of such tracks pass
547  // if(0==atrack->GetFilterMap() && atrack->GetType() == AliAODTrack::kFromDecayVtx)
548  // but there are exceptions which pass
549  // if(0==atrack->GetFilterMap() && atrack->GetType() == AliAODTrack::kPrimary)
550  // So it's a mess...
551  // 1) Clearly, for none of such tracks we can test a filter bit
552  // 2) Apparently all of them have positive ID, meaning they are global (?)
553  // Corresponding AOD global track:
554 
555  // c) Determine the corresponding "gtrack" (a.k.a. "normal global" track):
556  if(0 == fGlobalTracksAOD[0]->GetSize()){GlobalTracksAOD(aAOD,0);}
557  if(0 == fGlobalTracksAOD[0]->GetSize()){return;}
558  Int_t id = atrack->GetID();
559  AliAODTrack *gtrack = dynamic_cast<AliAODTrack*>(id>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id+1))));
560  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
561 
562  // d) Fill the control histograms for non-identified particles:
564 
565  // e) Fill the control histograms for identified particles:
567 
568  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
569 
570  } // else if(aAOD)
571 
572 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsParticle(AliVEvent *ave)
573 
574 //================================================================================================================
575 
577 {
578  // Fill control histograms for non-identified particles.
579  // Uses PassesCommonGlobalTrackCuts(...) TBI
580 
581  // a) Insanity checks;
582  // b) Check cut selection criteria;
583  // c) Fill control histograms.
584 
585  // a) Insanity checks:
586  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticles(AliAODTrack *gtrack)";
587  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");} // TBI keep this for some time, eventually just continue
588 
589  // b) Check cut selection criteria:
590  if(!PassesCommonGlobalTrackCuts(gtrack)){return;} // TBI in the method PassesCommonGlobalTrackCuts track is hardwired to AliAODTrack. Try to generalize
591 
592  // c) Fill control histograms:
593  fChargeHist->Fill(gtrack->Charge()+0.5); // see how this histogram was booked for convention used
594  fGetTPCNclsHist->Fill(gtrack->GetTPCNcls());
595  fGetTPCsignalNHist->Fill(gtrack->GetTPCsignalN());
596  fGetITSNclsHist->Fill(gtrack->GetITSNcls());
597  fdEdxVsPtHist->Fill(gtrack->GetTPCmomentum(),gtrack->GetTPCsignal());
598  fPtHist->Fill(gtrack->Pt());
599  fEtaHist->Fill(gtrack->Eta());
600  fPhiHist->Fill(gtrack->Phi());
601  fMassHist->Fill(gtrack->M());
602 
603 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticles(AliAODTrack *gtrack)
604 
605 //================================================================================================================
606 
608 {
609  // Fill control histograms for all charged MC particles.
610 
611  // a) Insanity checks;
612  // b) Check cut selection criteria;
613  // c) Fill control histograms.
614 
615  // a) Insanity checks:
616  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticles(AliAODMCParticle *amcparticle)";
617  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
618 
619  // b) Check cut selection criteria:
620  //if(!PassesCommonGlobalTrackCuts(gtrack)){return;} // TBI in the method PassesCommonGlobalTrackCuts track is hardwired to AliAODTrack. Try to generalize
621 
622  // c) Fill control histograms:
623  fChargeHist->Fill(amcparticle->Charge()+0.5); // see how this histogram was booked for convention used
624  fPtHist->Fill(amcparticle->Pt());
625  fEtaHist->Fill(amcparticle->Eta());
626  fPhiHist->Fill(amcparticle->Phi());
627  fMassHist->Fill(amcparticle->M());
628  fGetPdgCode->Fill(amcparticle->GetPdgCode());
629 
630 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsNonIdentifiedParticles(AliAODMCParticle *amcparticle)
631 
632 //================================================================================================================
633 
635 {
636  // Fill control histograms for identified particles.
637  // Uses PassesCommonTrackCuts(...) TBI yet still it stored the paremeters of corresponding global track... TBI
638 
639  // To do:
640  // 1) Add support for MC and ESD, now it works only for AliAODTrack.
641 
642  // a) Insanity checks;
643  // b) Check cut selection criteria;
644  // c) Fill control histograms.
645 
646  // a) Insanity checks:
647  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsIdentifiedParticles(AliAODTrack *atrack, AliAODTrack *gtrack)";
648  if(!atrack){Fatal(sMethodName.Data(),"!atrack");} // TBI keep this for some time, eventually just continue
649  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");} // TBI keep this for some time, eventually just continue
650 
651  // b) Check cut selection criteria:
652  if(!PassesCommonTrackCuts(atrack)){return;} // TBI in the method PassesCommonGlobalTrackCuts track is hardwired to AliAODTrack. Try to generalize
653  // TBI do I need some check for atrack?
654 
655  // c) Fill control histograms:
656  Int_t nCharge = 1;
657  Bool_t bPrimary = kTRUE;
658  for(Int_t charge=0;charge<2;charge++) // 0 = +q, 1 = -q
659  {
660  if(1==charge){nCharge = -1;} // TBI this is just ugly
661  for(Int_t ps=0;ps<2;ps++) // 0 = kPrimary, 1 = kFromDecayVtx
662  {
663  if(1==ps){bPrimary = kFALSE;}
664  // c2) Pions:
665  if(Pion(gtrack,nCharge,bPrimary))
666  {
667  fPtPIDHist[2][charge][ps]->Fill(gtrack->Pt());
668  fMassPIDHist[2][charge][ps]->Fill(gtrack->M());
669  fEtaPIDHist[2][charge][ps]->Fill(gtrack->Eta());
670  fPhiPIDHist[2][charge][ps]->Fill(gtrack->Phi());
671  }
672  // c3) Kaons:
673  if(Kaon(gtrack,nCharge,bPrimary))
674  {
675  fPtPIDHist[3][charge][ps]->Fill(gtrack->Pt());
676  fMassPIDHist[3][charge][ps]->Fill(gtrack->M());
677  fEtaPIDHist[3][charge][ps]->Fill(gtrack->Eta());
678  fPhiPIDHist[3][charge][ps]->Fill(gtrack->Phi());
679  }
680  // c4) Protons:
681  if(Proton(gtrack,nCharge,bPrimary))
682  {
683  fPtPIDHist[4][charge][ps]->Fill(gtrack->Pt());
684  fMassPIDHist[4][charge][ps]->Fill(gtrack->M());
685  fEtaPIDHist[4][charge][ps]->Fill(gtrack->Eta());
686  fPhiPIDHist[4][charge][ps]->Fill(gtrack->Phi());
687  }
688  } // for(Int_t ps=0;ps<2;ps++) // 0 = kPrimary, 1 = kFromDecayVtx
689  } // for(Int_t charge=0;charge<2;charge++) // 0 = +q, 1 = -q
690 
691 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsIdentifiedParticles(AliAODTrack *atrack, AliAODTrack *gtrack)
692 
693 
694 //================================================================================================================
695 
697 {
698  // Fill control histograms for identified Monte Carlo particles.
699 
700  // a) Insanity checks;
701  // b) Check cut selection criteria;
702  // c) Fill control histograms.
703 
704  // a) Insanity checks:
705  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsIdentifiedParticles(AliAODMCParticle *amcparticle)";
706  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
707 
708  // b) Check cut selection criteria:
709  if(!PassesCommonTrackCuts(amcparticle)){return;} // TBI have a look at implementation
710 
711  // c) Fill control histograms:
712  Int_t index = -44;
713  // c2) Pions:
714  if(211==TMath::Abs(amcparticle->GetPdgCode()))
715  {
716  index = 2;
717  } // if(211==TMath::Abs(amcparticle->GetPdgCode()))
718  // c3) Kaons:
719  else if(321==TMath::Abs(amcparticle->GetPdgCode()))
720  {
721  index = 3;
722  } // if(321==TMath::Abs(amcparticle->GetPdgCode()))
723  // c4) Protons:
724  else if(2212==TMath::Abs(amcparticle->GetPdgCode()))
725  {
726  index = 4;
727  } // if(2212==TMath::Abs(amcparticle->GetPdgCode()))
728  if(-44 != index)
729  {
730  Int_t charge = ((amcparticle->GetPdgCode()>0. ? 0 : 1)); // Okay... TBI
731  Int_t isPhysicalPrimary = ((amcparticle->IsPhysicalPrimary() ? 0 : 1)); // Okay... TBI
732  fPtPIDHist[index][charge][isPhysicalPrimary]->Fill(amcparticle->Pt());
733  fMassPIDHist[index][charge][isPhysicalPrimary]->Fill(amcparticle->M()); // in this context, clearly this is just a control histogram
734  fEtaPIDHist[index][charge][isPhysicalPrimary]->Fill(amcparticle->Eta());
735  fPhiPIDHist[index][charge][isPhysicalPrimary]->Fill(amcparticle->Phi());
736  } // if(-44 != index)
737 
738 } // void AliAnalysisTaskMultiparticleFemtoscopy::FillControlHistogramsIdentifiedParticles(AliAODMCParticle *amcparticle)
739 
740 //================================================================================================================
741 
743 {
744  // Estimate background.
745 
746  // To do:
747  // 1) Add support for MC and ESD.
748 
749  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::EstimateBackground(AliVEvent *ave)";
750 
751  // a) Determine Ali{MC,ESD,AOD}Event;
752  // b) ...
753 
754  // a) Determine Ali{MC,ESD,AOD}Event:
755  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
756  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
757  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
758 
759  if(aMC)
760  {
761  TString pattern = ".11.13.211.321.2212."; // TBI this can be done programatically TBI only for these particles I need to estimate background
762  if(!PassesMixedEventCuts(aMC)){return;} // TBI this is empty at the moment for MC
763  Int_t nTracks = aMC->GetNumberOfTracks();
764  if(!fMixedEvents[0])
765  {
766  TClonesArray ca0("AliAODMCParticle"); // temporary holder, its entries will be copied to fMixedEvents[0]
767  Int_t ca0Counter = 0;
768  for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
769  {
770  AliAODMCParticle *amcparticle = dynamic_cast<AliAODMCParticle*>(aMC->GetTrack(iTrack));
771  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
772  if(!PassesCommonTrackCuts(amcparticle)){continue;} // TBI re-think, see implemntation of this method
773  if(!(pattern.Contains(Form(".%d.",TMath::Abs(amcparticle->GetPdgCode()))))){continue;}
774  // Fill fMixedEvents[0] with tracks which survived all checks:
775  ca0[ca0Counter++] = amcparticle;
776  } // for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
777  fMixedEvents[0] = (TClonesArray*)ca0.Clone();
778  if(!fMixedEvents[0]){Fatal(sMethodName.Data(),"!fMixedEvents[0]");}
779  }
780  else if(!fMixedEvents[1])
781  {
782  TClonesArray ca1("AliAODMCParticle"); // temporary holder, its entries will be copied to fMixedEvents[1]
783  Int_t ca1Counter = 0;
784  for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
785  {
786  AliAODMCParticle *amcparticle = dynamic_cast<AliAODMCParticle*>(aMC->GetTrack(iTrack));
787  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
788  if(!PassesCommonTrackCuts(amcparticle)){continue;} // TBI re-think, see implemntation of this method
789  if(!(pattern.Contains(Form(".%d.",TMath::Abs(amcparticle->GetPdgCode()))))){continue;}
790  // Fill fMixedEvents[0] with tracks which survived all checks:
791  ca1[ca1Counter++] = amcparticle;
792  } // for(Int_t iTrack=0;iTrack<nTracks;iTrack++)
793  fMixedEvents[1] = (TClonesArray*)ca1.Clone();
794  if(!fMixedEvents[1]){Fatal(sMethodName.Data(),"!fMixedEvents[1]");}
795  }
796  // Shall we do something?
797  if(fMixedEvents[0] && fMixedEvents[1])
798  {
800  }
801  } // if(aMC)
802  else if(aESD)
803  {
804  // TBI
805  } // else if(aESD)
806  else if(aAOD)
807  {
808  if(!PassesMixedEventCuts(aAOD)){return;}
809  if(!fMixedEvents[0])
810  {
811  if(fGlobalTracksAOD[1]) fGlobalTracksAOD[1]->Delete();
812  if(0 != fGlobalTracksAOD[1]->GetSize()){Fatal(sMethodName.Data(),"0 != fGlobalTracksAOD[1]->GetSize()");}
813  GlobalTracksAOD(aAOD,1);
814  if(0 == fGlobalTracksAOD[1]->GetSize()){return;}
815  fMixedEvents[0] = (TClonesArray*)aAOD->GetTracks()->Clone();
816  if(!fMixedEvents[0]){Fatal(sMethodName.Data(),"!fMixedEvents[0]");}
817  }
818  else if(!fMixedEvents[1])
819  {
820  if(fGlobalTracksAOD[2]) fGlobalTracksAOD[2]->Delete();
821  if(0 != fGlobalTracksAOD[2]->GetSize()){Fatal(sMethodName.Data(),"0 != fGlobalTracksAOD[2]->GetSize()");}
822  GlobalTracksAOD(aAOD,2);
823  if(0 == fGlobalTracksAOD[2]->GetSize()){return;}
824  fMixedEvents[1] = (TClonesArray*)aAOD->GetTracks()->Clone();
825  if(!fMixedEvents[1]){Fatal(sMethodName.Data(),"!fMixedEvents[1]");}
826  }
827  // Shall we do something?
828  if(fMixedEvents[0] && fMixedEvents[1])
829  {
831  }
832  } // else if(aAOD)
833 
834 } // void AliAnalysisTaskMultiparticleFemtoscopy::EstimateBackground(AliVEvent *ave)
835 
836 //================================================================================================================
837 
839 {
840  // Do all debugging in this function.
841 
842  // TBI: add support for fProcessBothKineAndReco
843 
844  // a) Determine Ali{MC,ESD,AOD}Event;
845  // b) Wait for specified event.
846 
847  // a) Determine Ali{MC,ESD,AOD}Event:
848  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
849  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
850  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
851 
852  // b) Wait for specified event:
854  {
855  if(aMC)
856  {
857  // TBI
858  } // if(aMC)
859  else if(aESD)
860  {
861  // TBI
862  } // else if(aESD)
863  else if(aAOD)
864  {
865  cout<<Form("aAOD->GetRunNumber() = %d",aAOD->GetRunNumber())<<endl;
866  cout<<Form("aAOD->GetBunchCrossNumber() = %d",aAOD->GetBunchCrossNumber())<<endl;
867  cout<<Form("aAOD->GetOrbitNumber() = %d",aAOD->GetOrbitNumber())<<endl;
868  cout<<Form("aAOD->GetPeriodNumber() = %d",aAOD->GetPeriodNumber())<<endl;
869  if(!SpecifiedEvent(aAOD->GetRunNumber(),aAOD->GetBunchCrossNumber(),aAOD->GetOrbitNumber(),aAOD->GetPeriodNumber())){return;}
870  } // else if(aAOD)
871  } // if(fWaitForSpecifiedEvent)
872 
873 } // void AliAnalysisTaskMultiparticleFemtoscopy::DoSomeDebugging(AliVEvent *ave)
874 
875 //================================================================================================================
876 
878 {
879  // Determine the current event number.
880 
881  Int_t currentEventNumber = -44;
882 
883  if(fAnalysisType->EqualTo("MC"))
884  {
885  currentEventNumber = (Int_t)fGetNumberOfTracksHist->GetEntries(); // TBI there is an issue clearly if fFillControlHistogramsEvent is not enabled. Yes, this is really shaky...
886  }
887 
888  else if(fAnalysisType->EqualTo("ESD"))
889  {
890  // TBI
891  }
892 
893  else if(fAnalysisType->Contains("AOD")) // TBI: Okay...?
894  {
895  currentEventNumber = (Int_t)fGetNumberOfTracksHist->GetEntries(); // TBI there is an issue clearly if fFillControlHistogramsEvent is not enabled. Yes, this is really shaky...
896  }
897 
898  return currentEventNumber;
899 
900 } // Int_t AliAnalysisTaskMultiparticleFemtoscopy::CurrentEventNumber()
901 
902 //================================================================================================================
903 
905 {
906  // Per request, do some online monitoring.
907 
908  // a) Update regularly the output file;
909  // b) Bail out after specified number of events.
910 
911  // a) Update regularly the output file:
913  {
914  Int_t currentEventNo = this->CurrentEventNumber(); // TBI not supported yet for MC and ESD
915  if(0 == currentEventNo % fUpdateFrequency)
916  {
917  cout<<Form("nEvts: %d",currentEventNo)<<endl;
918  TFile *f = new TFile(fUpdateWhichOutputFile->Data(),"recreate");
919  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
920  f->Close();
921  }
922  } // if(fUpdateOutputFile)
923 
924  // b) Bail out after specified number of events:
925  if(fMaxNumberOfEvents > 0)
926  {
927  Int_t currentEventNo = this->CurrentEventNumber(); // TBI not supported yet for MC and ESD
928  if(fMaxNumberOfEvents == currentEventNo)
929  {
930  cout<<Form("\nPer request, bailing out after %d events in the file %s .\n",fMaxNumberOfEvents,fUpdateWhichOutputFile->Data())<<endl;
931  TFile *f = new TFile(fUpdateWhichOutputFile->Data(),"recreate");
932  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
933  f->Close();
934  exit(0);
935  }
936  } // if(fMaxNumberOfEvents > 0)
937 
938 } // void AliAnalysisTaskMultiparticleFemtoscopy::OnlineMonitoring()
939 
940 //================================================================================================================
941 
943 {
944  // Insanity checks for UserCreateOutputObjects() method.
945 
946  // TBI
947 
948  //
949 
950 } // void AliAnalysisTaskMultiparticleFemtoscopy::InsanityChecksUserCreateOutputObjects()
951 
952 //================================================================================================================
953 
955 {
956  // Insanity checks for UserExec() method.
957 
958  // a) Insanity checks specific for all analyses types;
959  // b) Insanity checks specific only for MC analyses;
960  // c) Insanity checks specific only for ESD analyses;
961  // d) Insanity checks specific only for AOD analyses.
962 
963  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::InsanityChecksUserExec()";
964 
965  // a) Insanity checks specific for all analyses types: // TBI is this really specific for all analyses types
966  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
967  {
968  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
969  {
970  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
971  {
972  if(!fPIDCA[pid][pa][ps]) {Fatal(sMethodName.Data(),"!fPIDCA[pid][pa][ps]");}
973  if(0 != fPIDCA[pid][pa][ps]->GetEntriesFast()){Fatal(sMethodName.Data(),"0 != fPIDCA[pid][pa][ps]->GetEntriesFast()"); }
974  }
975  }
976  }
977  for(Int_t pid=0;pid<1;pid++) // [0=Lambda,1=...] // TBI is this really specific for all analyses types
978  {
979  if(!fPIDV0sCA[pid]){Fatal(sMethodName.Data(),"!fPIDV0sCA[pid]");}
980  if(0 != fPIDV0sCA[pid]->GetEntriesFast()){Fatal(sMethodName.Data(),"0 != fPIDV0sCA[pid]->GetEntriesFast()");}
981  }
982 
983  // b) Insanity checks specific only for MC analyses:
984  if(fAnalysisType->EqualTo("MC"))
985  {
986  // TBI
987  }
988 
989  // c) Insanity checks specific only for ESD analyses:
990  if(fAnalysisType->EqualTo("ESD")) // TBI 'else if' ?
991  {
992  // TBI
993  }
994 
995  // d) Insanity checks specific only for AOD analyses:
996  if(fAnalysisType->Contains("AOD")) // TBI 'else if' ?
997  {
998  if(!fGlobalTracksAOD[0]){Fatal(sMethodName.Data(),"!fGlobalTracksAOD[0]");}
999  if(0 != fGlobalTracksAOD[0]->GetSize()){Fatal(sMethodName.Data(),"0 != fGlobalTracksAOD[0]->GetSize()");}
1000  if(fProcessBothKineAndReco && fEstimateBackground){Fatal(sMethodName.Data(),"fProcessBothKineAndReco && fEstimateBackground");}
1001  // Remark: Note that the similar check is not needed for instance for fGlobalTracksAOD[1] and fGlobalTracksAOD[2],
1002  // which are used to estimate background, since some events can be stored in a buffer, until a suitable co-event is found,
1003  // to calculate background.
1004  } // if(fAnalysisType->Contains("AOD")
1005 
1006 } // void AliAnalysisTaskMultiparticleFemtoscopy::InsanityChecksUserExec()
1007 
1008 //=======================================================================================================================
1009 
1011 {
1012  // Monte Carlo analysis is performed in this method.
1013 
1014  // a) Debugging;
1015  // b) Common event selection criteria;
1016  // c) Fill control histogram for global event observables;
1017  // d) Fill control histogram for particles;
1018  // e) Calculate correlation functions;
1019  // f) Calculate background;
1020  // g) V0s.
1021 
1022  // a) Debugging:
1023  if(fDoSomeDebugging){this->DoSomeDebugging(aMC);}
1024 
1025  // b) Common event selection criteria:
1026  if(!this->PassesCommonEventCuts(aMC)){return;}
1027 
1028  // c) Fill control histogram for global event observables:
1030 
1031  // d) Fill control histograms for particles:
1033 
1034  // e) Calculate correlation functions:
1036 
1037  // f) Calculate background:
1038  if(fEstimateBackground){this->EstimateBackground(aMC);}
1039 
1040  return;
1041 
1042  // g) V0s:
1043  // TBI
1044 
1045 } // void AliAnalysisTaskMultiparticleFemtoscopy::MC(AliMCEvent *aMC)
1046 
1047 //=======================================================================================================================
1048 
1050 {
1051  // ESD analysis is performed in this method.
1052 
1053  // ...
1054  if(aESD)
1055  {
1056  aESD = NULL; // TBI eliminating warnings temporarily
1057  }
1058 
1059 } // void AliAnalysisTaskMultiparticleFemtoscopy::ESD(AliESDEvent *aESD)
1060 
1061 //=======================================================================================================================
1062 
1064 {
1065  // AOD analysis is performed in this method.
1066 
1067  // a) Debugging;
1068  // b) Common event selection criteria;
1069  // c) Filter our "normal global" tracks;
1070  // d) Fill control histogram for global event observables;
1071  // e) Fill control histogram for particles;
1072  // f) Calculate correlation functions;
1073  // g) Calculate background;
1074  // h) V0s.
1075 
1076  // a) Debugging:
1077  if(fDoSomeDebugging){this->DoSomeDebugging(aAOD);}
1078 
1079  // b) Common event selection criteria:
1080  if(!this->PassesCommonEventCuts(aAOD)){return;}
1081 
1082  // c) Filter out "normal global" tracks:
1083  this->GlobalTracksAOD(aAOD,0); // [0] stands for default analysis
1084  //if(0 == fGlobalTracksAOD[0]->GetSize()) return;
1085 
1086  // d) Fill control histogram for global event observables:
1088 
1089  // e) Fill control histograms for particles:
1091 
1092  // f) Calculate correlation functions:
1094 
1095  // g) Calculate background:
1096  if(fEstimateBackground){this->EstimateBackground(aAOD);}
1097 
1098 
1099  return;
1100 
1101 
1102  // h) V0s:
1103  V0s(aAOD); // TBI implement flag to enable/disable this call
1104 
1105 } // void AliAnalysisTaskMultiparticleFemtoscopy::AOD(AliAODEvent *aAOD)
1106 
1107 //=======================================================================================================================
1108 
1110 {
1111  // Analysis with V0s.
1112 
1113  // a) Determine Ali{MC,ESD,AOD}Event;
1114 
1115  // a) Determine Ali{MC,ESD,AOD}Event:
1116  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
1117  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
1118  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
1119 
1120  if(aMC)
1121  {
1122  // TBI
1123  }
1124  else if(aESD)
1125  {
1126  // TBI
1127  }
1128  else if(aAOD)
1129  {
1130  // AOD V0s:
1131  TClonesArray *caV0s = aAOD->GetV0s();
1132  Int_t index = 0;
1133  AliAODv0 *aAODv0 = NULL;
1134  Int_t nProngs = -44;
1135  while(caV0s->At(index))
1136  {
1137  aAODv0 = (AliAODv0*) caV0s->At(index++);
1138  if(!aAODv0) break;
1139 
1140 
1141  //AliAODv0 *temp = (AliAODv0*)fPIDV0sCA[0]->ConstructedAt(index-1);
1142  //temp = (AliAODv0*)aAODv0->Clone();
1143 
1144  //cout<<Form("indddex = %d, fPIDV0sCA[0]->GetEntries() = %d",index,fPIDV0sCA[0]->GetEntries())<<endl;
1145 
1146  // TBI...
1147  continue;
1148 
1149  nProngs = aAODv0->GetNProngs();
1150  fGetNProngsHist->Fill(nProngs);
1151  fMassK0ShortHist->Fill(aAODv0->MassK0Short());
1152  fMassLambdaHist->Fill(aAODv0->MassLambda());
1153  fMassAntiLambdaHist->Fill(aAODv0->MassAntiLambda());
1154  fOpenAngleV0Hist->Fill(aAODv0->OpenAngleV0());
1155  fRadiusV0Hist->Fill(aAODv0->RadiusV0());
1156  fDcaV0ToPrimVertexHist->Fill(aAODv0->DcaV0ToPrimVertex());
1157  fMomV0XHist->Fill(aAODv0->MomV0X());
1158  fMomV0YHist->Fill(aAODv0->MomV0Y());
1159  fMomV0ZHist->Fill(aAODv0->MomV0Z());
1160  fPtV0Hist->Fill(pow(aAODv0->Pt2V0(),0.5));
1161  fPseudoRapV0Hist->Fill(aAODv0->PseudoRapV0());
1162  fPAHist->Fill(aAODv0->Alpha(),aAODv0->PtArmV0());
1163  // Check sharing:
1164 
1165  fUniqueIDHistEBE->Fill(aAODv0->GetPosID());
1166  fUniqueIDHistEBE->Fill(aAODv0->GetNegID());
1167 
1168  /*
1169  //cout<<Form("V0PosID: %d , V0NegID: %d",aAODv0->GetPosID(),aAODv0->GetNegID())<<endl;
1170  Int_t trackPos = aAODv0->GetPosID()>=0 ? fGlobalTracksAOD[0]->GetValue(aAODv0->GetPosID()) : fGlobalTracksAOD[0]->GetValue(-(aAODv0->GetPosID()+1));
1171  Int_t trackNeg = aAODv0->GetNegID()>=0 ? fGlobalTracksAOD[0]->GetValue(aAODv0->GetNegID()) : fGlobalTracksAOD[0]->GetValue(-(aAODv0->GetNegID()+1));
1172  //cout<<Form("global : %d, global : %d", trackPos, trackNeg)<<endl;
1173  */
1174 
1175  if(-1 != fUniqueIDHistEBE->FindFirstBinAbove(1.44,1)) // TBI
1176  {
1177  cout<<Form("fUniqueIDHistEBE->FindFirstBinAbove(1.44) %d:",(Int_t)fUniqueIDHistEBE->FindFirstBinAbove(1.44,1))<<endl; //
1178  }
1179  } // while(caV0s->At(index))
1180 
1181 
1182  /*
1183  Int_t index2 = 0;
1184  while(caV0s->At(index2))
1185  {
1186  cout<<((AliAODv0*) caV0s->At(index2))->Pt()<<endl;
1187  cout<<((AliAODv0*) fPIDV0sCA[0]->At(index2))->Pt()<<endl;
1188  index2++;
1189  cout<<endl;
1190  }
1191  */
1192 
1193  // AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,globaltrack);
1194  // AliPIDResponse::EDetPidStatus statusTOF = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,globaltrack);
1195 
1196  }
1197 
1198 } // void AliAnalysisTaskMultiparticleFemtoscopy::V0s(AliVEvent *ave)
1199 
1200 //=======================================================================================================================
1201 
1203 {
1204  // Reset all event-by-event objects.
1205 
1206  // a) Reset event-by-event objects specific for all analyses types;
1207  // b) Reset event-by-event objects specific only for MC analyses;
1208  // c) Reset event-by-event objects specific only for ESD analyses;
1209  // d) Reset event-by-event objects specific only for AOD analyses.
1210 
1211  // a) Reset event-by-event objects specific for all analyses types:
1212  if(fUniqueIDHistEBE) fUniqueIDHistEBE->Reset();
1213  // TBI add comment
1214  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1215  {
1216  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
1217  {
1218  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
1219  {
1220  if(fPIDCA[pid][pa][ps]) fPIDCA[pid][pa][ps]->Delete();
1221  }
1222  }
1223  }
1224  // TBI add comment
1225  for(Int_t pid=0;pid<1;pid++) // [0=Lambda,1=...]
1226  {
1227  if(fPIDV0sCA[pid]) fPIDV0sCA[pid]->Delete();
1228  }
1229 
1230  // b) Reset event-by-event objects specific only for MC analyses:
1231  if(fAnalysisType->EqualTo("MC"))
1232  {
1233  // TBI
1234  }
1235 
1236  // c) Reset event-by-event objects specific only for ESD analyses:
1237  if(fAnalysisType->EqualTo("ESD")) // TBI 'else if' ?
1238  {
1239  // TBI
1240  }
1241 
1242  // d) Reset event-by-event objects specific only for AOD analyses:
1243  if(fAnalysisType->Contains("AOD")) // TBI 'else if' ?
1244  {
1245  if(fGlobalTracksAOD[0]) fGlobalTracksAOD[0]->Delete();
1246  // Remark: Note that fGlobalTracksAOD[1] and fGlobalTracksAOD[2] used for event mixing do not have
1247  // to be reset e-b-e. Where then I reset them? TBI
1248  }
1249 
1250 } // void AliAnalysisTaskMultiparticleFemtoscopy::ResetEBEObjects()
1251 
1252 //=======================================================================================================================
1253 
1255 {
1256  // Is this primary/secondary (anti-)pion?
1257 
1258  // a) Insanity checks;
1259  // b) Trivial checks; // TBI
1260  // c) Track quality cuts;
1261  // d) PID;
1262  // e) PID MC (if available).
1263 
1264  // a) Insanity checks
1265  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Pion(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)";
1266  if(!(1 == charge || -1 == charge)){Fatal(sMethodName.Data(),"!(1 == charge || -1 == charge)");}
1267  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
1268 
1269  // b) Trivial checks:
1270  if(charge != (Int_t)gtrack->Charge()) return kFALSE;
1271  if(bPrimary && gtrack->GetType() != AliAODTrack::kPrimary)
1272  {
1273  return kFALSE;
1274  }
1275  else if(!bPrimary && gtrack->GetType() != AliAODTrack::kFromDecayVtx)
1276  {
1277  return kFALSE;
1278  }
1279 
1280  // c) Track quality cuts:
1281 /*
1282  if(bPrimary)
1283  {
1284  // TBI
1285  }
1286 */
1287 
1288  // d) PID:
1289  // For pT < 0.75 use only TPC
1290  if(gtrack->GetTPCmomentum() < 0.75) // TBI hardwired 0.75
1291  {
1292  AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,gtrack);
1293  if(!statusTPC) return kFALSE;
1294  // sigmas:
1295  Double_t dSigmaTPCPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kPion)));
1296  Double_t dSigmaTPCKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kKaon)));
1297  Double_t dSigmaTPCProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kProton)));
1298  // inclusive cuts and exclusive cuts:
1299  if(!(dSigmaTPCPion < fInclusiveSigmaCuts[2] && dSigmaTPCProton > fExclusiveSigmaCuts[2][4] && dSigmaTPCKaon > fExclusiveSigmaCuts[2][3])) return kFALSE;
1300  } // if(gtrack->GetTPCmomentum() <= 0.75) // TBI hardwired 0.75
1301  else if(gtrack->GetTPCmomentum() >= 0.75 )
1302  {
1303  // TBI use combined TPC and TOf
1304  return kFALSE; // TBI remove eventually
1305  }
1306 
1307  // e) PID MC (if available):
1309  {
1310  if(!fMC){Fatal(sMethodName.Data(),"!fMC");}
1311  AliAODMCParticle *mcParticle = dynamic_cast<AliAODMCParticle*>(fMC->GetTrack(TMath::Abs(gtrack->GetLabel())));
1312  if(charge < 0 && mcParticle->GetPdgCode() == -211) return kTRUE;
1313  else if(charge > 0 && mcParticle->GetPdgCode() == 211) return kTRUE;
1314  else return kFALSE;
1315  } // if(fProcessBothKineAndReco)
1316 
1317  return kTRUE;
1318 
1319 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Pion(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)
1320 
1321 
1322 //=======================================================================================================================
1323 
1325 {
1326  // Is this primary/secondary (anti-)kaon?
1327 
1328  // a) Insanity checks;
1329  // b) Trivial checks;
1330  // c) Track quality cuts;
1331  // d) PID;
1332  // e) PID MC (if available).
1333 
1334  // a) Insanity checks:
1335  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Kaon(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)";
1336  if(!(1 == charge || -1 == charge)){Fatal(sMethodName.Data(),"!(1 == charge || -1 == charge)");}
1337  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
1338 
1339  // b) Trivial checks:
1340  if(charge != (Int_t)gtrack->Charge()) return kFALSE;
1341  if(bPrimary && gtrack->GetType() != AliAODTrack::kPrimary)
1342  {
1343  return kFALSE;
1344  }
1345  else if(!bPrimary && gtrack->GetType() != AliAODTrack::kFromDecayVtx)
1346  {
1347  return kFALSE;
1348  }
1349 
1350  // c) Track quality cuts:
1351  /*
1352  if(bPrimary)
1353  {
1354  // TBI
1355  }
1356  */
1357 
1358  // d) PID:
1359  // For pT < 0.75 use only TPC
1360  if(gtrack->GetTPCmomentum() < 0.75) // TBI hardwired 0.75
1361  {
1362  AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,gtrack);
1363  if(!statusTPC) return kFALSE;
1364  // sigmas:
1365  Double_t dSigmaTPCKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kKaon)));
1366  Double_t dSigmaTPCProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kProton)));
1367  Double_t dSigmaTPCPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kPion)));
1368  // inclusive and exclusive cuts:
1369  if(!(dSigmaTPCKaon < fInclusiveSigmaCuts[3] && dSigmaTPCProton > fExclusiveSigmaCuts[3][4] && dSigmaTPCPion > fExclusiveSigmaCuts[3][2])) return kFALSE;
1370  } // if(gtrack->GetTPCmomentum() <= 0.75) // TBI hardwired 0.75
1371  else if(gtrack->GetTPCmomentum() >= 0.75 )
1372  {
1373  // TBI use combined TPC and TOf
1374  return kFALSE; // TBI remove eventually
1375  }
1376 
1377  // e) PID MC (if available):
1379  {
1380  if(!fMC){Fatal(sMethodName.Data(),"!fMC");}
1381  AliAODMCParticle *mcParticle = dynamic_cast<AliAODMCParticle*>(fMC->GetTrack(TMath::Abs(gtrack->GetLabel())));
1382  if(charge < 0 && mcParticle->GetPdgCode() == -321) return kTRUE;
1383  else if(charge > 0 && mcParticle->GetPdgCode() == 321) return kTRUE;
1384  else return kFALSE;
1385  } // if(fProcessBothKineAndReco)
1386 
1387  return kTRUE;
1388 
1389 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Kaon(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)
1390 
1391 //=======================================================================================================================
1392 
1394 {
1395  // Is this primary/secondary (anti-)proton?
1396 
1397  // a) Insanity checks;
1398  // b) Trivial checks;
1399  // c) Track quality cuts;
1400  // d) PID;
1401  // e) PID MC (if available).
1402 
1403  // a) Insanity checks:
1404  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Proton(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)";
1405  if(!(1 == charge || -1 == charge)){Fatal(sMethodName.Data(),"!(1 == charge || -1 == charge)");}
1406  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
1407 
1408  // b) Trivial checks:
1409  if(charge != (Int_t)gtrack->Charge()) return kFALSE;
1410  if(bPrimary && gtrack->GetType() != AliAODTrack::kPrimary)
1411  {
1412  return kFALSE;
1413  }
1414  else if(!bPrimary && gtrack->GetType() != AliAODTrack::kFromDecayVtx)
1415  {
1416  return kFALSE;
1417  }
1418 
1419  // c) Track quality cuts:
1420 /*
1421  if(bPrimary)
1422  {
1423  // TBI
1424  }
1425 */
1426 
1427  // d) PID:
1428  // For pT < 0.75 use only TPC
1429  if(gtrack->GetTPCmomentum() < 0.75) // TBI hardwired 0.75
1430  {
1431  AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,gtrack);
1432  if(!statusTPC) return kFALSE;
1433  // sigmas:
1434  Double_t dSigmaTPCProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kProton)));
1435  Double_t dSigmaTPCPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kPion)));
1436  Double_t dSigmaTPCKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(gtrack,(AliPID::kKaon)));
1437  if(!(dSigmaTPCProton < fInclusiveSigmaCuts[4] && dSigmaTPCPion > fExclusiveSigmaCuts[4][2] && dSigmaTPCKaon > fExclusiveSigmaCuts[4][3])) return kFALSE;
1438  } // if(gtrack->GetTPCmomentum() <= 0.75) // TBI hardwired 0.75
1439  else if(gtrack->GetTPCmomentum() >= 0.75 )
1440  {
1441  // TBI use combined TPC and TOf
1442  return kFALSE; // TBI remove eventually
1443  }
1444 
1445  // e) PID MC (if available):
1447  {
1448  if(!fMC){Fatal(sMethodName.Data(),"!fMC");}
1449  AliAODMCParticle *mcParticle = dynamic_cast<AliAODMCParticle*>(fMC->GetTrack(TMath::Abs(gtrack->GetLabel())));
1450  if(charge < 0 && mcParticle->GetPdgCode() == -2212) return kTRUE;
1451  else if(charge > 0 && mcParticle->GetPdgCode() == 2212) return kTRUE;
1452  else return kFALSE;
1453  } // if(fProcessBothKineAndReco)
1454 
1455  return kTRUE;
1456 
1457 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::Proton(AliAODTrack *gtrack, Int_t charge, Bool_t bPrimary)
1458 
1459 //=======================================================================================================================
1460 
1462 {
1463  // Initialize arrays for all objects not classified yet.
1464 
1465  for(Int_t index=0;index<10;index++) // [0] is used in the default analysis, [1] and [2] for event mixing, etc.
1466  {
1467  fGlobalTracksAOD[index] = NULL;
1468  }
1469 
1470 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArrays()
1471 
1472 //=======================================================================================================================
1473 
1475 {
1476  // Initialize all arrays for control histograms.
1477 
1478  for(Int_t xyz=0;xyz<3;xyz++)
1479  {
1480  fVertexXYZ[xyz] = NULL;
1481  }
1482 
1483  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1484  {
1485  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
1486  {
1487  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
1488  {
1489  fMassPIDHist[pid][pa][ps] = NULL;
1490  fPtPIDHist[pid][pa][ps] = NULL;
1491  fEtaPIDHist[pid][pa][ps] = NULL;
1492  fPhiPIDHist[pid][pa][ps] = NULL;
1493  }
1494  }
1495  } // for(Int_t pid=0;pid<4;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1496 
1497  // Inclusive sigma cuts:
1498  for(Int_t pf=0;pf<5;pf++) // PID function [0=Electron(...),1=Muon(...),2=Pion(...),3=Kaon(...),4=Proton(...)]
1499  {
1500  fInclusiveSigmaCuts[pf] = 0.;
1501  }
1502 
1503  // Exclusive sigma cuts:
1504  for(Int_t pf=0;pf<5;pf++) // PID function [0=Electron(...),1=Muon(...),2=Pion(...),3=Kaon(...),4=Proton(...)]
1505  {
1506  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1507  {
1508  fExclusiveSigmaCuts[pf][pid] = 0.;
1509  }
1510  }
1511 
1512 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForControlHistograms()
1513 
1514 //=======================================================================================================================
1515 
1517 {
1518  // Initialize all arrays for e-b-e objects.
1519 
1520  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1521  {
1522  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
1523  {
1524  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
1525  {
1526  fPIDCA[pid][pa][ps] = NULL;
1527  }
1528  }
1529  } // for(Int_t pid=0;pid<4;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1530 
1531 
1532  for(Int_t pid=0;pid<1;pid++)
1533  {
1534  fPIDV0sCA[pid] = NULL;
1535  }
1536 
1537 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForEBEObjects()
1538 
1539 //=======================================================================================================================
1540 
1542 {
1543  // Initialize all arrays for correlation functions.
1544 
1545  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]
1546  {
1547  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]
1548  {
1549  fCorrelationFunctions[pid1][pid2] = NULL;
1550  }
1551  }
1552 
1553 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForCorrelationFunctions()
1554 
1555 //=======================================================================================================================
1556 
1558 {
1559  // Initialize all arrays for background.
1560 
1561  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]
1562  {
1563  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]
1564  {
1565  fBackground[pid1][pid2] = NULL;
1566  }
1567  }
1568 
1569  for(Int_t me=0;me<2;me++) // [0] is buffer for 1st event; [1] for 2nd
1570  {
1571  fMixedEvents[me] = NULL;
1572  }
1573 
1574 } // void AliAnalysisTaskMultiparticleFemtoscopy::InitializeArraysForBackground()
1575 
1576 //=======================================================================================================================
1577 
1579 {
1580  // Book and nest all lists nested in the base list fHistList.
1581 
1582  // a) Book and nest lists for control histograms;
1583  // b) Book and nest lists for eveny-by-event histograms;
1584  // c) Book and nest lists for correlation functions;
1585  // d) Book and nest lists for background.
1586 
1587  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::BookAndNestAllLists()";
1588  if(!fHistList){Fatal(sMethodName.Data(),"fHistList is NULL");}
1589 
1590  // a) Book and nest lists for control histograms:
1591  fControlHistogramsList = new TList();
1592  fControlHistogramsList->SetName("Control_histograms");
1593  fControlHistogramsList->SetOwner(kTRUE);
1596  fControlHistogramsEventList->SetName("Event");
1597  fControlHistogramsEventList->SetOwner(kTRUE);
1600  fControlHistogramsNonIdentifiedParticlesList->SetName("Non-identified_particles");
1604  fControlHistogramsIdentifiedParticlesList->SetName("Identified_particles");
1608  fControlHistogramsV0sList->SetName("V0s");
1609  fControlHistogramsV0sList->SetOwner(kTRUE);
1611 
1612  // b) Book and nest lists for eveny-by-event histograms:
1613  fEBEHistogramsList = new TList();
1614  fEBEHistogramsList->SetName("Event-by-event_histograms");
1615  fEBEHistogramsList->SetOwner(kTRUE);
1617 
1618  // c) Book and nest lists for correlation functions:
1620  fCorrelationFunctionsList->SetName("Correlation_Functions");
1621  fCorrelationFunctionsList->SetOwner(kTRUE);
1623 
1624  // d) Book and nest lists for background:
1625  fBackgroundList = new TList();
1626  fBackgroundList->SetName("Background");
1627  fBackgroundList->SetOwner(kTRUE);
1628  fHistList->Add(fBackgroundList);
1629 
1630 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookAndNestAllLists()
1631 
1632 //=======================================================================================================================
1633 
1635 {
1636  // Book all unclassified objects temporary here. TBI
1637 
1638  fAnalysisType = new TString();
1639 
1640  fPIDResponse = new AliPIDResponse();
1641 
1642  for(Int_t index=0;index<fMaxNoGlobalTracksAOD;index++)
1643  {
1644  fGlobalTracksAOD[index] = new TExMap();
1645  }
1646 
1647  // Inclusive sigma cuts:
1648  fInclusiveSigmaCuts[2] = 2.; // i.e. in function Pion(...) the inclusive cut for pions is 2 sigma
1649  fInclusiveSigmaCuts[3] = 2.; // i.e. in function Kaon(...) the inclusive cut for kaons is 2 sigma
1650  fInclusiveSigmaCuts[4] = 2.; // i.e. in function Proton(...) the inclusive cut for protons is 2 sigma
1651 
1652  // Exclusive sigma cuts:
1653  // Pion(...)
1654  fExclusiveSigmaCuts[2][3] = 4.; // i.e. in function Pion(...) the exclusive cut for kaons is 4 sigma
1655  fExclusiveSigmaCuts[2][4] = 4.; // i.e. in function Pion(...) the exclusive cut for protons is 4 sigma
1656  // Kaon(...)
1657  fExclusiveSigmaCuts[3][2] = 4.; // i.e. in function Kaon(...) the exclusive cut for pions is 4 sigma
1658  fExclusiveSigmaCuts[3][4] = 4.; // i.e. in function Kaon(...) the exclusive cut for protons is 4 sigma
1659  // Proton(...)
1660  fExclusiveSigmaCuts[4][2] = 4.; // i.e. in function Proton(...) the exclusive cut for pions is 4 sigma
1661  fExclusiveSigmaCuts[4][3] = 4.; // i.e. in function Proton(...) the exclusive cut for kaons is 4 sigma
1662 
1663 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverything()
1664 
1665 //=======================================================================================================================
1666 
1668 {
1669  // Book all the stuff for control histograms.
1670 
1671  // a) Book the profile holding all the flags for control histograms;
1672  // b) Common vaiables;
1673  // c) Book all control histograms...
1674  // c0) Event;
1675  // c1) Non-identified particles;
1676  // c2) Identified particles;
1677  // c3) V0s.
1678 
1679  // a) Book the profile holding all the flags for control histograms: TBI stil incomplete
1680  fControlHistogramsFlagsPro = new TProfile("fControlHistogramsFlagsPro","Flags and settings for control histograms",1,0,1);
1681  fControlHistogramsFlagsPro->SetTickLength(-0.01,"Y");
1682  fControlHistogramsFlagsPro->SetMarkerStyle(25);
1683  fControlHistogramsFlagsPro->SetLabelSize(0.04);
1684  fControlHistogramsFlagsPro->SetLabelOffset(0.02,"Y");
1685  fControlHistogramsFlagsPro->SetStats(kFALSE);
1686  fControlHistogramsFlagsPro->SetFillColor(kGray);
1687  fControlHistogramsFlagsPro->SetLineColor(kBlack);
1688  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistograms"); fControlHistogramsFlagsPro->Fill(0.5,fFillControlHistograms);
1690 
1691  if(!fFillControlHistograms){return;} // TBI is this safe?
1692 
1693  // b) Common vaiables:
1694  TString sParticleLabel[5] = {"e","#mu","#pi","K","p"};
1695  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()};
1696  // ...
1697 
1698  // c0) Event:
1699  // Book the profile holding all the flags for control histograms for global event observables:
1700  fControlHistogramsEventFlagsPro = new TProfile("fControlHistogramsEventFlagsPro","Flags and settings for TBI",1,0,1);
1701  fControlHistogramsEventFlagsPro->SetTickLength(-0.01,"Y");
1702  fControlHistogramsEventFlagsPro->SetMarkerStyle(25);
1703  fControlHistogramsEventFlagsPro->SetLabelSize(0.04);
1704  fControlHistogramsEventFlagsPro->SetLabelOffset(0.02,"Y");
1705  fControlHistogramsEventFlagsPro->SetStats(kFALSE);
1706  fControlHistogramsEventFlagsPro->SetFillColor(kGray);
1707  fControlHistogramsEventFlagsPro->SetLineColor(kBlack);
1708  fControlHistogramsEventFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsEvent"); fControlHistogramsEventFlagsPro->Fill(0.5,fFillControlHistogramsEvent);
1711  {
1712  fGetNumberOfTracksHist = new TH1I("fGetNumberOfTracksHist","aAOD->GetNumberOfTracks() (Remark: Not all of tracks are unique.)",10000,0,10000);
1713  //fGetNumberOfTracksHist->SetStats(kFALSE);
1714  fGetNumberOfTracksHist->SetFillColor(kBlue-10);
1716  fGetNumberOfGlobalTracksHist = new TH1I("fGetNumberOfGlobalTracksHist","fGlobalTracksAOD[0]->GetSize()",10000,0,10000);
1717  //fGetNumberOfGlobalTracksHist->SetStats(kFALSE);
1718  fGetNumberOfGlobalTracksHist->SetFillColor(kBlue-10);
1720  fGetNumberOfV0sHist = new TH1I("fGetNumberOfV0sHist","aAOD->GetNumberOfV0s() (Remark: Some V0s share the daughter.)",10000,0,10000);
1721  fGetNumberOfV0sHist->SetStats(kFALSE);
1722  fGetNumberOfV0sHist->SetFillColor(kBlue-10);
1724  fGetNumberOfCascadesHist = new TH1I("fGetNumberOfCascadesHist","aAOD->GetNumberOfCascades() (TBI: Not validated.)",10000,0,10000);
1725  fGetNumberOfCascadesHist->SetStats(kFALSE);
1726  fGetNumberOfCascadesHist->SetFillColor(kBlue-10);
1728  fGetMagneticFieldHist = new TH1D("fGetMagneticFieldHist","aAOD->GetMagneticField()",20,-10.,10.);
1729  fGetMagneticFieldHist->SetFillColor(kBlue-10);
1731  fGetEventTypeHist = new TH1I("fGetEventTypeHist","aAOD->GetEventType()",1000,0,10);
1732  fGetEventTypeHist->SetFillColor(kBlue-10);
1734  fGetCentralityHist = new TH1D("fGetCentralityHist","aAOD->GetCentrality()",100,0.,100.);
1735  fGetCentralityHist->SetFillColor(kBlue-10);
1737  TString sxyz[3] = {"X","Y","Z"};
1738  for(Int_t xyz=0;xyz<3;xyz++)
1739  {
1740  fVertexXYZ[xyz] = new TH1F(Form("fVertex%s",sxyz[xyz].Data()),Form("avtz->Get%s()",sxyz[xyz].Data()),100000,-50.,50);
1741  fVertexXYZ[xyz]->SetStats(kFALSE);
1743  }
1744  fGetNContributorsHist = new TH1I("fGetNContributorsHist","avtx->GetNContributors()",10000,0,10000);
1745  fGetNContributorsHist->SetStats(kFALSE);
1746  fGetNContributorsHist->SetFillColor(kBlue-10);
1748  fGetChi2perNDFHist = new TH1F("fGetChi2perNDFHist","avtx->GetChi2perNDF()",5000,0.,50.);
1749  fGetChi2perNDFHist->SetStats(kFALSE);
1750  fGetChi2perNDFHist->SetFillColor(kBlue-10);
1752  fGetNDaughtersHist = new TH1I("GetNDaughtersHist","avtx->GetNDaughters()",10000,0,10000);
1753  fGetNDaughtersHist->SetStats(kFALSE);
1754  fGetNDaughtersHist->SetFillColor(kBlue-10);
1756  // ...
1757  } // if(fFillControlHistogramsEvent)
1758 
1759  // c1) Non-identified particles:
1760  // Book the profile holding all the flags for control histograms for non-identified particles:
1761  fControlHistogramsNonIdentifiedParticlesFlagsPro = new TProfile("fControlHistogramsNonIdentifiedParticlesFlagsPro","Flags and settings for TBI",1,0,1);
1762  fControlHistogramsNonIdentifiedParticlesFlagsPro->SetTickLength(-0.01,"Y");
1765  fControlHistogramsNonIdentifiedParticlesFlagsPro->SetLabelOffset(0.02,"Y");
1769  fControlHistogramsNonIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsNonIdentifiedParticles"); fControlHistogramsNonIdentifiedParticlesFlagsPro->Fill(0.5,fFillControlHistogramsNonIdentifiedParticles);
1772  {
1773  fChargeHist = new TH1I("fChargeHist","atrack->Charge()",5,-2,3);
1774  fChargeHist->SetStats(kFALSE);
1775  fChargeHist->SetFillColor(kBlue-10);
1776  fChargeHist->GetXaxis()->SetBinLabel(1,"-2");
1777  fChargeHist->GetXaxis()->SetBinLabel(2,"-1");
1778  fChargeHist->GetXaxis()->SetBinLabel(3,"0");
1779  fChargeHist->GetXaxis()->SetBinLabel(4,"1");
1780  fChargeHist->GetXaxis()->SetBinLabel(5,"2");
1782  fGetTPCNclsHist = new TH1I("fGetTPCNclsHist","atrack->fGetTPCNclsHist()",200,0,200);
1783  fGetTPCNclsHist->SetStats(kFALSE);
1784  fGetTPCNclsHist->SetFillColor(kBlue-10);
1785  fGetTPCNclsHist->GetXaxis()->SetTitle("TPCNcls");
1787  fGetTPCsignalNHist = new TH1I("fGetTPCsignalNHist","atrack->fGetTPCsignalNHist()",200,0,200);
1788  fGetTPCsignalNHist->SetStats(kFALSE);
1789  fGetTPCsignalNHist->SetFillColor(kBlue-10);
1790  fGetTPCsignalNHist->GetXaxis()->SetTitle("TPCsignalN");
1792  fGetITSNclsHist = new TH1I("fGetITSNclsHist","atrack->fGetITSNclsHist()",200,0,200);
1793  fGetITSNclsHist->SetStats(kFALSE);
1794  fGetITSNclsHist->SetFillColor(kBlue-10);
1795  fGetITSNclsHist->GetXaxis()->SetTitle("ITSNcls");
1797  fdEdxVsPtHist = new TH2F("fdEdxVsPtHist","atrack->GetTPCmomentum(),atrack->GetTPCsignal()",1000,0.,20.,1000,-500.,500.);
1798  fdEdxVsPtHist->SetStats(kFALSE);
1800  fPtHist = new TH1F("fPtHist","atrack->Pt()",1000,0.,20.);
1801  fPtHist->SetStats(kFALSE);
1802  fPtHist->SetFillColor(kBlue-10);
1803  fPtHist->SetMinimum(0.);
1805  fEtaHist = new TH1F("fEtaHist","atrack->Eta()",200,-2.,2.);
1806  fEtaHist->SetStats(kFALSE);
1807  fEtaHist->SetFillColor(kBlue-10);
1808  fEtaHist->SetMinimum(0.);
1810  fPhiHist = new TH1F("fPhiHist","atrack->Phi()",360,0.,TMath::TwoPi());
1811  fPhiHist->SetStats(kFALSE);
1812  fPhiHist->SetFillColor(kBlue-10);
1813  fPhiHist->SetMinimum(0.);
1815  // TBI
1816  fMassHist = new TH1F("fMassHist","atrack->M()",10000,0.,10.);
1817  fMassHist->SetStats(kFALSE);
1818  fMassHist->SetFillColor(kBlue-10);
1819  fMassHist->SetMinimum(0.);
1820  for(Int_t nm=0;nm<5;nm++) // nominal masses
1821  {
1822  fMassHist->GetXaxis()->SetBinLabel(fMassHist->FindBin(dNominalMass[nm]),Form("m_{%s}",sParticleLabel[nm].Data()));
1823  }
1825  fGetFilterMap = new TH1I("fGetFilterMap","atrack->fGetFilterMap()",10000,0,10000);
1826  fGetFilterMap->SetStats(kFALSE);
1827  fGetFilterMap->SetFillColor(kBlue-10);
1828  fGetFilterMap->SetMinimum(0.);
1830  fGetPdgCode = new TH1I("fGetPdgCode","atrack->fGetPdgCode()",20000,-10000,10000);
1831  fGetPdgCode->SetStats(kFALSE);
1832  fGetPdgCode->SetFillColor(kBlue-10);
1833  fGetPdgCode->SetMinimum(0.);
1835  } // if(fFillControlHistogramsNonIdentifiedParticles)
1836 
1837  // c2) Identified particles:
1838  // Book the profile holding all the flags for TBI:
1839  fControlHistogramsIdentifiedParticlesFlagsPro = new TProfile("fControlHistogramsIdentifiedParticlesFlagsPro","Flags and settings for TBI",1,0,1);
1840  fControlHistogramsIdentifiedParticlesFlagsPro->SetTickLength(-0.01,"Y");
1843  fControlHistogramsIdentifiedParticlesFlagsPro->SetLabelOffset(0.02,"Y");
1845  fControlHistogramsIdentifiedParticlesFlagsPro->SetFillColor(kGray);
1846  fControlHistogramsIdentifiedParticlesFlagsPro->SetLineColor(kBlack);
1847  fControlHistogramsIdentifiedParticlesFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsIdentifiedParticles"); fControlHistogramsIdentifiedParticlesFlagsPro->Fill(0.5,fFillControlHistogramsIdentifiedParticles);
1850  {
1851  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1852  {
1853  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
1854  {
1855  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
1856  {
1857  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.);
1858  fMassPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("m [GeV/c^{2}]");
1859  for(Int_t nm=0;nm<5;nm++) // nominal masses
1860  {
1861  fMassPIDHist[pid][pa][ps]->GetXaxis()->SetBinLabel(fMassPIDHist[pid][pa][ps]->FindBin(dNominalMass[nm]),Form("m_{%s}",sParticleLabel[nm].Data()));
1862  }
1863  fMassPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
1865  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.);
1866  fPtPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("p_{T} [TBI units]");
1867  fPtPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
1869  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.);
1870  fEtaPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
1872  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());
1873  fPhiPIDHist[pid][pa][ps]->GetXaxis()->SetTitle("#phi");
1874  fPhiPIDHist[pid][pa][ps]->SetFillColor(kBlue-10);
1876  } // for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
1877  } // for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
1878  } // for(Int_t pid=0;pid<4;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
1879  } // if(fFillControlHistogramsIdentifiedParticles)
1880 
1881  // c3) V0s:
1882  // Book the profile holding all the flags for V0s:
1883  fControlHistogramsV0sFlagsPro = new TProfile("fControlHistogramsV0sFlagsPro","Flags and settings for V0s",1,0,1);
1884  fControlHistogramsV0sFlagsPro->SetTickLength(-0.01,"Y");
1885  fControlHistogramsV0sFlagsPro->SetMarkerStyle(25);
1886  fControlHistogramsV0sFlagsPro->SetLabelSize(0.04);
1887  fControlHistogramsV0sFlagsPro->SetLabelOffset(0.02,"Y");
1888  fControlHistogramsV0sFlagsPro->SetStats(kFALSE);
1889  fControlHistogramsV0sFlagsPro->SetFillColor(kGray);
1890  fControlHistogramsV0sFlagsPro->SetLineColor(kBlack);
1891  fControlHistogramsV0sFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistogramsV0s"); fControlHistogramsV0sFlagsPro->Fill(0.5,fFillControlHistogramsV0s);
1894  {
1895  fGetNProngsHist = new TH1I("fGetNProngsHist","aAODv0->GetNProngs()",10,0,10);
1896  fGetNProngsHist->SetStats(kFALSE);
1897  fGetNProngsHist->SetFillColor(kBlue-10);
1899  // TBI
1900  fMassK0ShortHist = new TH1F("fMassK0ShortHist","aAODv0->MassK0Short()",1000000,0.,100.);
1901  //fMassK0ShortHist->SetStats(kFALSE);
1902  fMassK0ShortHist->SetFillColor(kBlue-10);
1903  Double_t dMassK0Short = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // nominal mass
1904  //fMassK0ShortHist->GetXaxis()->SetBinLabel(fMassK0ShortHist->FindBin(dMassK0Short),Form("m_{K_{S}^{0}} = %f",dMassK0Short));
1905  fMassK0ShortHist->SetBinContent(fMassK0ShortHist->FindBin(dMassK0Short),1e6);
1907  // TBI
1908  fMassLambdaHist = new TH1F("fMassLambdaHist","aAODv0->MassLambda()",1000000,0.,100.);
1909  //fMassLambdaHist->SetStats(kFALSE);
1910  fMassLambdaHist->SetFillColor(kBlue-10);
1911  Double_t dMassLambda = TDatabasePDG::Instance()->GetParticle(3122)->Mass(); // nominal mass
1912  //fMassLambdaHist->GetXaxis()->SetBinLabel(fMassLambdaHist->FindBin(dMassLambda),Form("m_{#Lambda^{0}} = %f",dMassLambda));
1913  fMassLambdaHist->SetBinContent(fMassLambdaHist->FindBin(dMassLambda),1e6);
1915  // TBI
1916  fMassAntiLambdaHist = new TH1F("fMassAntiLambdaHist","aAODv0->MassAntiLambda()",1000000,0.,100.);
1917  //fMassAntiLambdaHist->SetStats(kFALSE);
1918  fMassAntiLambdaHist->SetFillColor(kBlue-10);
1919  Double_t dMassAntiLambda = TDatabasePDG::Instance()->GetParticle(-3122)->Mass(); // nominal mass
1920  //fMassAntiLambdaHist->GetXaxis()->SetBinLabel(fMassAntiLambdaHist->FindBin(dMassAntiLambda),Form("m_{#bar{Lambda}^{0}} = %f",dMassAntiLambda));
1921  fMassAntiLambdaHist->SetBinContent(fMassAntiLambdaHist->FindBin(dMassAntiLambda),1e6);
1923  // TBI
1924  fOpenAngleV0Hist = new TH1F("fOpenAngleV0Hist","aAODv0->fOpenAngleV0()",10000,-0.044,TMath::Pi()+0.044);
1925  fOpenAngleV0Hist->SetStats(kFALSE);
1926  fOpenAngleV0Hist->SetFillColor(kBlue-10);
1928  // TBI
1929  fRadiusV0Hist = new TH1F("fRadiusV0Hist","aAODv0->fRadiusV0()",10000,0.,1000.);
1930  fRadiusV0Hist->SetStats(kFALSE);
1931  fRadiusV0Hist->SetFillColor(kBlue-10);
1933  // TBI
1934  fDcaV0ToPrimVertexHist = new TH1F("fDcaV0ToPrimVertexHist","aAODv0->fDcaV0ToPrimVertex()",10000,0.,1000.);
1935  fDcaV0ToPrimVertexHist->SetStats(kFALSE);
1936  fDcaV0ToPrimVertexHist->SetFillColor(kBlue-10);
1938  // TBI
1939  fMomV0XHist = new TH1F("fMomV0XHist","aAODv0->fMomV0X() = px(+) + px(-)",10000,-1000.,1000.);
1940  fMomV0XHist->SetStats(kFALSE);
1941  fMomV0XHist->SetFillColor(kBlue-10);
1943  // TBI
1944  fMomV0YHist = new TH1F("fMomV0YHist","aAODv0->fMomV0Y() = py(+) + py(-)",10000,-1000.,1000.);
1945  fMomV0YHist->SetStats(kFALSE);
1946  fMomV0YHist->SetFillColor(kBlue-10);
1948  // TBI
1949  fMomV0ZHist = new TH1F("fMomV0ZHist","aAODv0->fMomV0Z() = pz(+) + pz(-)",10000,-1000.,1000.);
1950  fMomV0ZHist->SetStats(kFALSE);
1951  fMomV0ZHist->SetFillColor(kBlue-10);
1953  // TBI
1954  fPtV0Hist = new TH1F("fPtV0Hist","pow(aAODv0->fPt2V0(),0.5)",10000,0.,100.);
1955  fPtV0Hist->SetStats(kFALSE);
1956  fPtV0Hist->SetFillColor(kBlue-10);
1958  // TBI
1959  fPseudoRapV0Hist = new TH1F("fPseudoRapV0Hist","aAODv0->PseudoRapV0()",1000,-10.,10.);
1960  fPseudoRapV0Hist->SetStats(kFALSE);
1961  fPseudoRapV0Hist->SetFillColor(kBlue-10);
1963  // TBI
1964  fPAHist = new TH2F("fPAHist","TBI",100,-2.,2.,100,0.,1.);
1965  fPAHist->SetStats(kFALSE);
1966  fPAHist->GetXaxis()->SetTitle("#alpha");
1967  fPAHist->GetYaxis()->SetTitle("p_{T}");
1969  } // if(fFillControlHistogramsV0s)
1970 
1971 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForControlHistograms()
1972 
1973 //=======================================================================================================================
1974 
1976 {
1977  // Book all the stuff for event-by-event objects.
1978 
1979  // a) Book the profile holding all the flags for event-by-event objects;
1980  // b) Book all event-by-event objects.
1981 
1982  // a) Book the profile holding all the flags for EBE objects:
1983  fEBEObjectsFlagsPro = new TProfile("fEBEObjectsFlagsPro","Flags and settings for event-by-event histograms",1,0,1);
1984  fEBEObjectsFlagsPro->SetTickLength(-0.01,"Y");
1985  fEBEObjectsFlagsPro->SetMarkerStyle(25);
1986  fEBEObjectsFlagsPro->SetLabelSize(0.04);
1987  fEBEObjectsFlagsPro->SetLabelOffset(0.02,"Y");
1988  fEBEObjectsFlagsPro->SetStats(kFALSE);
1989  fEBEObjectsFlagsPro->SetFillColor(kGray);
1990  fEBEObjectsFlagsPro->SetLineColor(kBlack);
1991  //fEBEObjectsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillEBEHistograms"); fEBEObjectsFlagsPro->Fill(0.5,fFillEBEHistograms);
1993 
1994  //if(!fFillEBEHistograms){return;} // TBI rethink
1995 
1996  // TBI
1997  fUniqueIDHistEBE = new TH1I("fUniqueIDHistEBE","TBI",40000,-20000,20000);
1998  fUniqueIDHistEBE->SetStats(kFALSE);
1999  fUniqueIDHistEBE->SetFillColor(kBlue-10);
2000  // TBI I do not want to store this histogram, right?
2001 
2002  // TBI add comment
2003  for(Int_t pid=0;pid<5;pid++) // [0=e,1=mu,2=pi,3=K,4=p]
2004  {
2005  for(Int_t pa=0;pa<2;pa++) // particle(+q)/antiparticle(-q)
2006  {
2007  for(Int_t ps=0;ps<2;ps++) // kPrimary/kFromDecayVtx
2008  {
2009  fPIDCA[pid][pa][ps] = new TClonesArray("AliAODTrack",10000);
2010  }
2011  }
2012  }
2013 
2014  // TBI add comment
2015  for(Int_t pid=0;pid<1;pid++) // [0=Lambda,1=...]
2016  {
2017  fPIDV0sCA[pid] = new TClonesArray("AliAODv0",10000);
2018  }
2019 
2020 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForEBEObjects()
2021 
2022 //=======================================================================================================================
2023 
2025 {
2026  // Book all the stuff for correlation functions objects.
2027 
2028  // a) Book the profile holding all the flags for correlation functions objects;
2029  // b) Book all correlation functions;
2030  // c) Book TExMap *fCorrelationFunctionsIndices.
2031 
2032  // a) Book the profile holding all the flags for correlation functions objects:
2033  fCorrelationFunctionsFlagsPro = new TProfile("fCorrelationFunctionsFlagsPro","Flags and settings for correlation functions histograms",2,0,2);
2034  fCorrelationFunctionsFlagsPro->SetTickLength(-0.01,"Y");
2035  fCorrelationFunctionsFlagsPro->SetMarkerStyle(25);
2036  fCorrelationFunctionsFlagsPro->SetLabelSize(0.04);
2037  fCorrelationFunctionsFlagsPro->SetLabelOffset(0.02,"Y");
2038  fCorrelationFunctionsFlagsPro->SetStats(kFALSE);
2039  fCorrelationFunctionsFlagsPro->SetFillColor(kGray);
2040  fCorrelationFunctionsFlagsPro->SetLineColor(kBlack);
2041  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(0.5,fFillCorrelationFunctions);
2042  fCorrelationFunctionsFlagsPro->GetXaxis()->SetBinLabel(2,"fNormalizeCorrelationFunctions"); fCorrelationFunctionsFlagsPro->Fill(1.5,fNormalizeCorrelationFunctions);
2044 
2045  if(!fFillCorrelationFunctions){return;} // TBI is this safe?
2046 
2047  const Int_t nParticleSpecies = 5; // Supported at the moment: 0=e,1=mu,2=pi,3=K,4=p
2048  TString sParticles[2*nParticleSpecies] = {"e^{+}","#mu^{+}","#pi^{+}","K^{+}","p^{+}","e^{-}","#mu^{-}","#pi^{-}","K^{-}","p^{-}"};
2049 
2050  // b) Book all correlation functions:
2051  // Remark 0: First particle in the pair is always the one with positive charge.
2052  // Remark 1: Diagonal elements are particle/antiparticle pairs.
2053  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]
2054  {
2055  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]
2056  {
2057  // Correlation functions:
2058  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
2059  fCorrelationFunctions[pid1][pid2]->SetStats(kFALSE);
2060  fCorrelationFunctions[pid1][pid2]->SetFillColor(kBlue-10);
2061  fCorrelationFunctions[pid1][pid2]->SetXTitle("k");
2062  fCorrelationFunctions[pid1][pid2]->SetYTitle("C(k)");
2064  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
2065  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
2066 
2067  // c) Book TExMap *fCorrelationFunctionsIndices:
2068  // 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]
2069  fCorrelationFunctionsIndices = new TExMap();
2070  fCorrelationFunctionsIndices->Add(11,0);
2071  fCorrelationFunctionsIndices->Add(13,1);
2072  fCorrelationFunctionsIndices->Add(211,2);
2073  fCorrelationFunctionsIndices->Add(321,3);
2074  fCorrelationFunctionsIndices->Add(2212,4);
2075  fCorrelationFunctionsIndices->Add(-11,5);
2076  fCorrelationFunctionsIndices->Add(-13,6);
2077  fCorrelationFunctionsIndices->Add(-211,7);
2078  fCorrelationFunctionsIndices->Add(-321,8);
2079  fCorrelationFunctionsIndices->Add(-2212,9);
2080 
2081 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForCorrelationFunctions()
2082 
2083 //=======================================================================================================================
2084 
2086 {
2087  // Book all the stuff for background objects.
2088 
2089  // a) Book the profile holding all the flags for background objects;
2090  // b) Book all histograms for background;
2091  // c) Book buffer objects.
2092 
2093  // a) Book the profile holding all the flags for correlation functions objects:
2094  fBackgroundFlagsPro = new TProfile("fBackgroundFlagsPro","Flags and settings for background histograms",1,0,1);
2095  fBackgroundFlagsPro->SetTickLength(-0.01,"Y");
2096  fBackgroundFlagsPro->SetMarkerStyle(25);
2097  fBackgroundFlagsPro->SetLabelSize(0.04);
2098  fBackgroundFlagsPro->SetLabelOffset(0.02,"Y");
2099  fBackgroundFlagsPro->SetStats(kFALSE);
2100  fBackgroundFlagsPro->SetFillColor(kGray);
2101  fBackgroundFlagsPro->SetLineColor(kBlack);
2102  fBackgroundFlagsPro->GetXaxis()->SetBinLabel(1,"fEstimateBackground"); fBackgroundFlagsPro->Fill(0.5,fEstimateBackground);
2104 
2105  if(!fEstimateBackground){return;} // TBI is this safe?
2106 
2107  const Int_t nParticleSpecies = 5; // Supported at the moment: 0=e,1=mu,2=pi,3=K,4=p
2108  TString sParticles[2*nParticleSpecies] = {"e^{+}","#mu^{+}","#pi^{+}","K^{+}","p^{+}","e^{-}","#mu^{-}","#pi^{-}","K^{-}","p^{-}"};
2109 
2110  // b) Book all histograms for background:
2111  // Remark 0: First particle in the pair is always the one with positive charge.
2112  // Remark 1: Diagonal elements are particle/antiparticle pairs.
2113  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]
2114  {
2115  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]
2116  {
2117  // Background:
2118  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
2119  fBackground[pid1][pid2]->SetStats(kFALSE);
2120  fBackground[pid1][pid2]->SetFillColor(kBlue-10);
2121  fBackground[pid1][pid2]->SetXTitle("k");
2122  fBackground[pid1][pid2]->SetYTitle("B(k)");
2123  fBackgroundList->Add(fBackground[pid1][pid2]);
2124  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
2125  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
2126 
2127  // c) Book buffer objects:
2128  // TBI not sure why I need this here again. Re-think...
2129  for(Int_t me=0;me<2;me++) // [0] is buffer for 1st event; [1] for 2nd
2130  {
2131  fMixedEvents[me] = NULL;
2132  }
2133 
2134 } // void AliAnalysisTaskMultiparticleFemtoscopy::BookEverythingForBackground()
2135 
2136 //=======================================================================================================================
2137 
2139 {
2140  // Filter out unique global tracks in AOD and store them in fGlobalTracksAOD[index].
2141 
2142  // Remark 0: All global tracks have positive ID, the duplicated TPC-only tracks have -(ID+1);
2143  // Remark 1: The issue here is that there are apparently two sets of global tracks: a) "normal" and b) constrained to primary vertex.
2144  // However, only the "normal" global tracks come with positive ID, additionally they can be discriminated simply via: aodTrack->IsGlobalConstrained()
2145  // 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
2146  // iTrack: atrack->GetID(): atrack->Pt() atrack->Eta() atrack->Phi()
2147  // 1: 0: 2.073798 -0.503640 2.935432
2148  // 19: -1: 2.075537 -0.495988 2.935377 => this is TPC-only
2149  // 35: -1: 2.073740 -0.493576 2.935515 => this is IsGlobalConstrained()
2150  // In fact, this is important, otherwise there is double or even triple counting in some cases.
2151  // Remark 2: There are tracks for which: 0 == aodTrack->GetFilterMap()
2152  // a) Basically all of them pass: atrack->GetType() == AliAODTrack::kFromDecayVtx , but few exceptions also pass atrack->GetType() == AliAODTrack::kPrimary
2153  // b) All of them apparently have positive ID, i.e. these are global tracks
2154  // c) Clearly, we cannot use TestFilterBit() on them
2155  // d) None of them apparently satisfies: atrack->IsGlobalConstrained()
2156  // Remark 3: There is a performance penalty when fGlobalTracksAOD[1] and fGlobalTracksAOD[2] needed for mixed events are calculated.
2157  // Yes, I can get them directly from fGlobalTracksAOD[0], without calling this method for them again. TBI today
2158 
2159  for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
2160  {
2161  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack));
2162  if(aodTrack)
2163  {
2164  Int_t id = aodTrack->GetID();
2165  //if(id>=0 && aodTrack->GetFilterMap()>0 && !aodTrack->IsGlobalConstrained()) // TBI rethink this
2166  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
2167  {
2168  fGlobalTracksAOD[index]->Add(id,iTrack);
2169  } // if(id>=0 && !aodTrack->IsGlobalConstrained())
2170  } // if(aodTrack)
2171  } // for(Int_t iTrack=0;iTrack<aAOD->GetNumberOfTracks();iTrack++)
2172 
2173 } // void AliAnalysisTaskMultiparticleFemtoscopy::GlobalTracksAOD(AliAODEvent *aAOD)
2174 
2175 //=======================================================================================================================
2176 
2178 {
2179  // 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).
2180 
2181  if(run != fRun) return kFALSE;
2182  else if(bunchCross != fBunchCross) return kFALSE;
2183  else if(orbit != fOrbit) return kFALSE;
2184  else if(period != fPeriod) return kFALSE;
2185 
2186  return kTRUE;
2187 
2188 } // void AliAnalysisTaskMultiparticleFemtoscopy::SpecifiedEvent(UInt_t run, UShort_t bunchCross, UInt_t orbit, UInt_t period)
2189 
2190 //=======================================================================================================================
2192 {
2193  // Check if the track passes common global track cuts (irrespectively of PID).
2194 
2195  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonGlobalTrackCuts(AliAODTrack *gtrack)";
2196  if(!gtrack){Fatal(sMethodName.Data(),"!gtrack");}
2197 
2198  // To do: add data members and corresponding setters:
2199  // fPtMin, fPtMax
2200  // fEtaMin, fEtaMax
2201  // fPhiMin, fPhiMax
2202  // fTPCNclsMin, fTPCNclsMax
2203  // fTPCsignalNMin, fTPCsignalNMax
2204 
2205  if(gtrack->Pt()<0.2) return kFALSE;
2206  if(gtrack->Pt()>=5.0) return kFALSE;
2207  if(gtrack->Eta()<-0.8) return kFALSE;
2208  if(gtrack->Eta()>=0.8) return kFALSE;
2209  //if(gtrack->Phi()<-0.6) return kFALSE;
2210  //if(gtrack->Phi()>=0.6) return kFALSE;
2211  //if(gtrack->GetTPCNcls()<70) return kFALSE;
2212  //if(gtrack->GetTPCsignalN()<70) return kFALSE;
2213 
2214  if(fRejectFakeTracks && gtrack->GetLabel()<0) return kFALSE; // TBI is it more precise <=0 ?
2215 
2216  return kTRUE;
2217 
2218 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonGlobalTrackCuts(AliAODTrack *gtrack)
2219 
2220 //=======================================================================================================================
2221 
2223 {
2224  // Check if the track passes common analysis specific track (e.g. TPC-only) cuts.
2225  // Thereforem we can cut independetly on global track parameters, and on TPC-only cut parameters.
2226 
2227  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODTrack *atrack)";
2228  if(!atrack){Fatal(sMethodName.Data(),"!atrack");}
2229 
2230  if(!atrack->TestFilterBit(128)) return kFALSE; // TPC-only TBI setter
2231 
2232  if(fRejectFakeTracks && atrack->GetLabel()<0) return kFALSE; // TBI is it more precise <=0 ?
2233 
2234  return kTRUE;
2235 
2236 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODTrack *atrack)
2237 
2238 //=======================================================================================================================
2239 
2241 {
2242  // TBI this method applies only to MC, make it uniform with AOD.
2243 
2244  TString sMethodName = "Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODMCParticle *amcparticle)";
2245  if(!amcparticle){Fatal(sMethodName.Data(),"!amcparticle");}
2246 
2247  if(TMath::Abs(amcparticle->Charge())<1.e-10) return kFALSE; // skipping neutral particles TBI
2248  if(!amcparticle->IsPhysicalPrimary()) return kFALSE; // skipping secondaries TBI
2249  if(amcparticle->Pt()<0.2) return kFALSE;
2250  if(amcparticle->Pt()>=5.0) return kFALSE;
2251  if(amcparticle->Eta()<-0.8) return kFALSE;
2252  if(amcparticle->Eta()>=0.8) return kFALSE;
2253 
2254  return kTRUE;
2255 
2256 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonTrackCuts(AliAODMCParticle *amcparticle)
2257 
2258 //=======================================================================================================================
2259 
2261 {
2262  // Check if the event passes common event cuts.
2263 
2264  // TBI: add support for fProcessBothKineAndReco
2265 
2266  // a) Determine Ali{MC,ESD,AOD}Event;
2267 
2268  // a) Determine Ali{MC,ESD,AOD}Event:
2269  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
2270  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
2271  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
2272 
2273  if(aMC)
2274  {
2275  // TBI
2276  }
2277  else if(aESD)
2278  {
2279  // TBI
2280  }
2281  else if(aAOD)
2282  {
2283  if(!aAOD->GetPrimaryVertex()) return kFALSE;
2284  if(TMath::Abs(aAOD->GetMagneticField())<0.001) return kFALSE;
2285  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
2286  //if(TMath::Abs(avtx->GetX())>10.0) return kFALSE;
2287  //if(TMath::Abs(avtx->GetY())>10.0) return kFALSE;
2288  if(TMath::Abs(avtx->GetZ())>10.0) return kFALSE; // TBI setter
2289  if(avtx->GetNContributors()<=2) return kFALSE; // TBI setter
2290  } // else if(aAOD)
2291 
2292  return kTRUE;
2293 
2294 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesCommonEventCuts(AliAODEvent *aAOD)
2295 
2296 //=======================================================================================================================
2297 
2299 {
2300  // Check if the event passes mixed event cuts.
2301 
2302  // a) Determine Ali{MC,ESD,AOD}Event;
2303 
2304  // a) Determine Ali{MC,ESD,AOD}Event:
2305  AliMCEvent *aMC = dynamic_cast<AliMCEvent*>(ave);
2306  AliESDEvent *aESD = dynamic_cast<AliESDEvent*>(ave);
2307  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(ave);
2308 
2309  if(aMC)
2310  {
2311  // TBI
2312  }
2313  else if(aESD)
2314  {
2315  // TBI
2316  }
2317  else if(aAOD)
2318  {
2319  if(!aAOD->GetPrimaryVertex()) return kFALSE;
2320  if(TMath::Abs(aAOD->GetMagneticField())<0.001) return kFALSE;
2321  AliAODVertex *avtx = (AliAODVertex*)aAOD->GetPrimaryVertex();
2322  //if(TMath::Abs(avtx->GetX())>10.0) return kFALSE;
2323  //if(TMath::Abs(avtx->GetY())>10.0) return kFALSE;
2324  if(TMath::Abs(avtx->GetZ())>10.0) return kFALSE; // TBI setter
2325  if(avtx->GetNContributors()<=2) return kFALSE; // TBI setter
2326  }
2327 
2328  return kTRUE;
2329 
2330 } // Bool_t AliAnalysisTaskMultiparticleFemtoscopy::PassesMixedEventCuts(AliVEvent *ave)
2331 
2332 //=======================================================================================================================
2333 
2335 {
2336  // Calculate correlation functions.
2337 
2338  // a) Insanity checks;
2339  // b) Two nested loops to calculate C(k), just an example.
2340 
2341  // a) Insanity checks:
2342  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::CalculateCorrelationFunctions(AliAODEvent *aAOD)";
2343  if(!aAOD){Fatal(sMethodName.Data(),"!aAOD");}
2344  if(0 == fGlobalTracksAOD[0]->GetSize()){GlobalTracksAOD(aAOD,0);}
2345  if(0 == fGlobalTracksAOD[0]->GetSize()){return;}
2346 
2347  // b) Two nested loops to calculate C(k), just an example:
2348  Int_t nTracks = aAOD->GetNumberOfTracks();
2349  for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
2350  {
2351  AliAODTrack *atrack1 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack1));
2352  // TBI Temporary track insanity checks:
2353  if(!atrack1){Fatal(sMethodName.Data(),"!atrack1");} // TBI keep this for some time, eventually just continue
2354  if(atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
2355  if(atrack1->TestFilterBit(128) && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->TestFiletrBit(128) && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
2356  if(!PassesCommonTrackCuts(atrack1)){continue;} // TBI re-think
2357  // Corresponding AOD global track:
2358  Int_t id1 = atrack1->GetID();
2359  AliAODTrack *gtrack1 = dynamic_cast<AliAODTrack*>(id1>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id1)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id1+1))));
2360  if(!gtrack1){Fatal(sMethodName.Data(),"!gtrack1");} // TBI keep this for some time, eventually just continue
2361  Int_t gid1 = (id1>=0 ? id1 : -(id1+1)); // ID of corresponding global track
2362  // Common track selection criteria for all "normal" global tracks:
2363  if(!PassesCommonGlobalTrackCuts(gtrack1)){continue;}
2364 
2365  // Loop over the 2nd particle:
2366  for(Int_t iTrack2=iTrack1+1;iTrack2<nTracks;iTrack2++)
2367  {
2368  AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack2));
2369  // TBI Temporary track insanity checks:
2370  if(!atrack2){Fatal(sMethodName.Data(),"!atrack2");} // TBI keep this for some time, eventually just continue
2371  if(atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
2372  if(atrack2->TestFilterBit(128) && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->TestFiletrBit(128) && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
2373  if(!PassesCommonTrackCuts(atrack2)){continue;} // TBI re-think
2374  // Corresponding AOD global track:
2375  Int_t id2 = atrack2->GetID();
2376  AliAODTrack *gtrack2 = dynamic_cast<AliAODTrack*>(id2>=0 ? aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(id2)) : aAOD->GetTrack(fGlobalTracksAOD[0]->GetValue(-(id2+1))));
2377  if(!gtrack2){Fatal(sMethodName.Data(),"!gtrack2");} // TBI keep this for some time, eventually just continue
2378  Int_t gid2 = (id2>=0 ? id2 : -(id2+1)); // ID of corresponding global track
2379  if(gid1==gid2){continue;} // Eliminate self-correlations:
2380 
2381  // Common track selection criteria for all "normal" global tracks:
2382  if(!PassesCommonGlobalTrackCuts(gtrack2)){continue;}
2383 
2384  // Okay, so we have two tracks, let's check PID, and fill the correlation functions:
2385 
2386  // 1.) Same particle species:
2387 
2388  // a) pion-pion:
2389  // a1) pi+pi+ [2][2]:
2390  if(Pion(gtrack1,1,kTRUE) && Pion(gtrack2,1,kTRUE))
2391  {
2392  fCorrelationFunctions[2][2]->Fill(RelativeMomenta(gtrack1,gtrack2));
2393  }
2394  // a2) pi-pi- [7][7]:
2395  if(Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,-1,kTRUE))
2396  {
2397  fCorrelationFunctions[7][7]->Fill(RelativeMomenta(gtrack1,gtrack2));
2398  }
2399  // a3) pi+pi- || pi-pi+ [2][7]:
2400  if((Pion(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE)) || (Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,1,kTRUE)))
2401  {
2402  fCorrelationFunctions[2][7]->Fill(RelativeMomenta(gtrack1,gtrack2));
2403  }
2404 
2405  // b) kaon-kaon:
2406  // b1) K+K+ [3][3]:
2407  if(Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE))
2408  {
2409  fCorrelationFunctions[3][3]->Fill(RelativeMomenta(gtrack1,gtrack2));
2410  }
2411  // b2) K-K- [8][8]:
2412  if(Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
2413  {
2414  fCorrelationFunctions[8][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2415  }
2416  // b3) K+K- || K-K+ [3][8]:
2417  if((Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE)) || (Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,1,kTRUE)))
2418  {
2419  fCorrelationFunctions[3][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2420  }
2421 
2422  // c) proton-proton:
2423  // c1) p+p+ [4][4]:
2424  if(Proton(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
2425  {
2426  fCorrelationFunctions[4][4]->Fill(RelativeMomenta(gtrack1,gtrack2));
2427  }
2428  // c2) p-p- [9][9]:
2429  if(Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2430  {
2431  fCorrelationFunctions[9][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2432  }
2433  // c3) p+p- || p-p+ [4][9]:
2434  if((Proton(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE)) || (Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,1,kTRUE)))
2435  {
2436  fCorrelationFunctions[4][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2437  }
2438 
2439  // 2.) Mixed particle species:
2440  // a) pion-kaon
2441  // a1) pi+K+ [2][3]:
2442  if(Pion(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE))
2443  {
2444  fCorrelationFunctions[2][3]->Fill(RelativeMomenta(gtrack1,gtrack2));
2445  }
2446  // a2) pi+K- [2][8]
2447  if(Pion(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
2448  {
2449  fCorrelationFunctions[2][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2450  }
2451  // a3) K+pi- [3][7]
2452  if(Kaon(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE))
2453  {
2454  fCorrelationFunctions[3][7]->Fill(RelativeMomenta(gtrack1,gtrack2));
2455  }
2456  // a4) pi-K- [7][8]
2457  if(Pion(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
2458  {
2459  fCorrelationFunctions[7][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2460  }
2461  // b) pion-proton
2462  // b1) pi+p+ [2][4]:
2463  if(Pion(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
2464  {
2465  fCorrelationFunctions[2][4]->Fill(RelativeMomenta(gtrack1,gtrack2));
2466  }
2467  // b2) pi+p- [2][9]
2468  if(Pion(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2469  {
2470  fCorrelationFunctions[2][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2471  }
2472  // b3) p+pi- [4][7]
2473  if(Proton(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE))
2474  {
2475  fCorrelationFunctions[4][7]->Fill(RelativeMomenta(gtrack1,gtrack2));
2476  }
2477  // b4) pi-p- [7][9]
2478  if(Pion(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2479  {
2480  fCorrelationFunctions[7][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2481  }
2482  // c) kaon-proton
2483  // c1) K+p+ [3][4]:
2484  if(Kaon(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
2485  {
2486  fCorrelationFunctions[3][4]->Fill(RelativeMomenta(gtrack1,gtrack2));
2487  }
2488  // c2) K+p- [3][9]
2489  if(Kaon(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2490  {
2491  fCorrelationFunctions[3][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2492  }
2493  // c3) p+K- [4][8]
2494  if(Proton(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
2495  {
2496  fCorrelationFunctions[4][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2497  }
2498  // c4) K-p- [8][9]
2499  if(Kaon(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2500  {
2501  fCorrelationFunctions[8][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2502  }
2503  } // for(Int_t iTrack2=0;iTrack2<nTracks;iTrack2++)
2504  } // for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
2505 
2506 } // void AliAnalysisTaskMultiparticleFemtoscopy::CalculateCorrelationFunctions()
2507 
2508 //=======================================================================================================================
2509 
2511 {
2512  // Calculate correlation functions for Monte Carlo.
2513 
2514  // a) Insanity checks;
2515  // b) Pattern from supported PDG codes;
2516  // c) Two nested loops to calculate C(k), just an example.
2517 
2518  // b) Pattern from supported PDG codes:
2519  TString pattern = ".11.13.211.321.2212."; // TBI this can be done programatically
2520 
2521  // a) Insanity checks:
2522  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::CalculateCorrelationFunctions(AliMCEvent *aMC)";
2523  if(0 == aMC->GetNumberOfTracks()){return;} // TBI re-think
2524 
2525  // c) Two nested loops to calculate C(k), just an example:
2526  Int_t nTracks = aMC->GetNumberOfTracks();
2527  for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
2528  {
2529  AliAODMCParticle *amcparticle1 = dynamic_cast<AliAODMCParticle*>(aMC->GetTrack(iTrack1));
2530  if(!amcparticle1){Fatal(sMethodName.Data(),"!amcparticle1");}
2531 
2532  if(!PassesCommonTrackCuts(amcparticle1)){continue;} // TBI re-think, see implemntation of this method
2533 
2534  // Loop over the 2nd particle:
2535  for(Int_t iTrack2=iTrack1+1;iTrack2<nTracks;iTrack2++)
2536  {
2537  AliAODMCParticle *amcparticle2 = dynamic_cast<AliAODMCParticle*>(aMC->GetTrack(iTrack2));
2538  if(!amcparticle2){Fatal(sMethodName.Data(),"!amcparticle2");}
2539 
2540  if(!PassesCommonTrackCuts(amcparticle2)){continue;} // TBI re-think, see implemntation of this method
2541 
2542  // Okay, so we have two tracks, let's check PID, and fill the correlation functions:
2543  // Check if this PID is supported in current implementation:
2544  if(!(pattern.Contains(Form(".%d.",TMath::Abs(amcparticle1->GetPdgCode()))) && pattern.Contains(Form(".%d.",TMath::Abs(amcparticle2->GetPdgCode()))))){continue;}
2545 
2546  // Determine the indices of correlation function to be filled:
2547  Int_t index1 = -44;
2548  Int_t index2 = -44;
2549  if(fCorrelationFunctionsIndices->GetValue(amcparticle1->GetPdgCode())<=fCorrelationFunctionsIndices->GetValue(amcparticle2->GetPdgCode()))
2550  {
2551  index1 = fCorrelationFunctionsIndices->GetValue(amcparticle1->GetPdgCode());
2552  index2 = fCorrelationFunctionsIndices->GetValue(amcparticle2->GetPdgCode());
2553  }
2554  else
2555  {
2556  index1 = fCorrelationFunctionsIndices->GetValue(amcparticle2->GetPdgCode());
2557  index2 = fCorrelationFunctionsIndices->GetValue(amcparticle1->GetPdgCode());
2558  }
2559  if(-44 == index1){Fatal(sMethodName.Data(),"-44 == index1");}
2560  if(-44 == index2){Fatal(sMethodName.Data(),"-44 == index2");}
2561 
2562  // Fill or die:
2563  fCorrelationFunctions[index1][index2]->Fill(RelativeMomenta(amcparticle1,amcparticle2)); // for the relative momenta ordering doesn't matter
2564 
2565  } // for(Int_t iTrack2=iTrack1+1;iTrack2<nTracks;iTrack2++)
2566  } // for(Int_t iTrack1=0;iTrack1<nTracks;iTrack1++)
2567 
2568 } // void AliAnalysisTaskMultiparticleFemtoscopy::CalculateCorrelationFunctions(AliMCEvent *aMC)
2569 
2570 //=======================================================================================================================
2571 
2572 Double_t AliAnalysisTaskMultiparticleFemtoscopy::RelativeMomenta(AliAODTrack *agtrack1, AliAODTrack *agtrack2)
2573 {
2574  // Comment the weather here. TBI
2575 
2576  Double_t k = -44.; // relative momenta k = \frac{1}{2}|\vec{p_1}-\vec{p_2}|
2577 
2578  // \vec{p_1}:
2579  Double_t p1x = agtrack1->Px();
2580  Double_t p1y = agtrack1->Py();
2581  Double_t p1z = agtrack1->Pz();
2582  // \vec{p_2}:
2583  Double_t p2x = agtrack2->Px();
2584  Double_t p2y = agtrack2->Py();
2585  Double_t p2z = agtrack2->Pz();
2586  // k:
2587  k = 0.5*pow(pow(p1x-p2x,2.)+pow(p1y-p2y,2.)+pow(p1z-p2z,2.),0.5);
2588 
2589  if(k<1.e-14) // TBI rething this constraint
2590  {
2591  cout<<Form("\nSelf-correlation !!!!")<<endl;
2592  cout<<Form("p1: %f %f %f",p1x,p1y,p1z)<<endl;
2593  cout<<Form("p2: %f %f %f\n",p2x,p2y,p2z)<<endl;
2594  exit(0);
2595  }
2596 
2597  return k;
2598 
2599 } // Double_t AliAnalysisTaskMultiparticleFemtoscopy::RelativeMomenta(AliAODTrack *agtrack1, AliAODTrack *agtrack2)
2600 
2601 //=======================================================================================================================
2602 
2603 Double_t AliAnalysisTaskMultiparticleFemtoscopy::RelativeMomenta(AliAODMCParticle *amcparticle1, AliAODMCParticle *amcparticle2)
2604 {
2605  // Comment the weather here. TBI
2606 
2607  Double_t k = -44.; // relative momenta k = \frac{1}{2}|\vec{p_1}-\vec{p_2}|
2608 
2609  // \vec{p_1}:
2610  Double_t p1x = amcparticle1->Px();
2611  Double_t p1y = amcparticle1->Py();
2612  Double_t p1z = amcparticle1->Pz();
2613  // \vec{p_2}:
2614  Double_t p2x = amcparticle2->Px();
2615  Double_t p2y = amcparticle2->Py();
2616  Double_t p2z = amcparticle2->Pz();
2617  // k:
2618  k = 0.5*pow(pow(p1x-p2x,2.)+pow(p1y-p2y,2.)+pow(p1z-p2z,2.),0.5);
2619 
2620  if(k<1.e-14) // TBI rething this constraint
2621  {
2622  cout<<Form("\nSelf-correlation !!!!")<<endl;
2623  cout<<Form("p1: %f %f %f",p1x,p1y,p1z)<<endl;
2624  cout<<Form("p2: %f %f %f\n",p2x,p2y,p2z)<<endl;
2625  exit(0);
2626  }
2627 
2628  return k;
2629 
2630 } // Double_t AliAnalysisTaskMultiparticleFemtoscopy::RelativeMomenta(AliAODMCParticle *amcparticle1, AliAODMCParticle *amcparticle2)
2631 
2632 //=======================================================================================================================
2633 
2634 void AliAnalysisTaskMultiparticleFemtoscopy::CalculateBackground(TClonesArray *ca1, TClonesArray *ca2)
2635 {
2636  // Calculate background.
2637 
2638  // TBI this method is not really validated
2639 
2640  // a) Insanity checks;
2641  // b) Temporary primitive algorithm: correlate particles from 'previous event' to particles from 'current event';
2642  // c) Two nested loops to calculate B(k);
2643  // d) Shift [1] -> [0];
2644  // e) Clean [1].
2645 
2646  // a) Insanity checks:
2647  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::CalculateBackground(TClonesArray *ca1, TClonesArray *ca2)";
2648  if(!ca1){Fatal(sMethodName.Data(),"!ca1");}
2649  if(!ca2){Fatal(sMethodName.Data(),"!ca2");}
2650 
2651  // b) Temporary primitive algorithm: correlate particles from 'previous event' to particles from 'current event'
2652  // ...
2653 
2654  // c) Two nested loops to calculate B(k):
2655  Int_t nTracks1 = ca1->GetEntries();
2656  Int_t nTracks2 = ca2->GetEntries();
2657 
2658  // Start loop over tracks in the 1st event:
2659  for(Int_t iTrack1=0;iTrack1<nTracks1;iTrack1++)
2660  {
2661  AliAODTrack *atrack1 = dynamic_cast<AliAODTrack*>(ca1->UncheckedAt(iTrack1)); // TBI cross-check UncheckedAt
2662  // TBI Temporary track insanity checks:
2663  if(!atrack1){Fatal(sMethodName.Data(),"!atrack1");} // TBI keep this for some time, eventually just continue
2664  if(atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->GetID()>=0 && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
2665  if(atrack1->TestFilterBit(128) && atrack1->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack1->TestFiletrBit(128) && atrack1->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
2666  if(!PassesCommonTrackCuts(atrack1)){continue;} // TBI re-think
2667  // Corresponding AOD global track:
2668  Int_t id1 = atrack1->GetID();
2669  AliAODTrack *gtrack1 = dynamic_cast<AliAODTrack*>(id1>=0 ? ca1->UncheckedAt(fGlobalTracksAOD[1]->GetValue(id1)) : ca1->UncheckedAt(fGlobalTracksAOD[1]->GetValue(-(id1+1))));
2670  if(!gtrack1){Fatal(sMethodName.Data(),"!gtrack1");} // TBI keep this for some time, eventually just continue
2671  // Common track selection criteria for all "normal" global tracks:
2672  if(!PassesCommonGlobalTrackCuts(gtrack1)){continue;}
2673 
2674  // Start loop over tracks in the 2nd event:
2675  for(Int_t iTrack2=0;iTrack2<nTracks2;iTrack2++)
2676  {
2677  AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(ca2->UncheckedAt(iTrack2));
2678  // TBI Temporary track insanity checks:
2679  if(!atrack2){Fatal(sMethodName.Data(),"!atrack2");} // TBI keep this for some time, eventually just continue
2680  if(atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->GetID()>=0 && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
2681  if(atrack2->TestFilterBit(128) && atrack2->IsGlobalConstrained()){Fatal(sMethodName.Data(),"atrack2->TestFiletrBit(128) && atrack2->IsGlobalConstrained()");} // TBI keep this for some time, eventually just continue
2682  if(!PassesCommonTrackCuts(atrack2)){continue;} // TBI re-think
2683  // Corresponding AOD global track:
2684  Int_t id2 = atrack2->GetID();
2685  AliAODTrack *gtrack2 = dynamic_cast<AliAODTrack*>(id2>=0 ? ca2->UncheckedAt(fGlobalTracksAOD[2]->GetValue(id2)) : ca2->UncheckedAt(fGlobalTracksAOD[2]->GetValue(-(id2+1))));
2686  if(!gtrack2){Fatal(sMethodName.Data(),"!gtrack2");} // TBI keep this for some time, eventually just continue
2687  // Common track selection criteria for all "normal" global tracks:
2688  if(!PassesCommonGlobalTrackCuts(gtrack2)){continue;}
2689 
2690  // Okay... So we have two tracks from two different events ready. Let's check PID, and calculate the background:
2691 
2692  // 1.) Same particle species:
2693 
2694  // a) pion-pion:
2695  // a1) pi+pi+ [2][2]:
2696  if(Pion(gtrack1,1,kTRUE) && Pion(gtrack2,1,kTRUE))
2697  {
2698  fBackground[2][2]->Fill(RelativeMomenta(gtrack1,gtrack2));
2699  }
2700  // a2) pi-pi- [7][7]:
2701  if(Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,-1,kTRUE))
2702  {
2703  fBackground[7][7]->Fill(RelativeMomenta(gtrack1,gtrack2));
2704  }
2705  // a3) pi+pi- || pi-pi+ [2][7]:
2706  if((Pion(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE)) || (Pion(gtrack1,-1,kTRUE) && Pion(gtrack2,1,kTRUE)))
2707  {
2708  fBackground[2][7]->Fill(RelativeMomenta(gtrack1,gtrack2));
2709  }
2710 
2711  // b) kaon-kaon:
2712  // b1) K+K+ [3][3]:
2713  if(Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE))
2714  {
2715  fBackground[3][3]->Fill(RelativeMomenta(gtrack1,gtrack2));
2716  }
2717  // b2) K-K- [8][8]:
2718  if(Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
2719  {
2720  fBackground[8][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2721  }
2722  // b3) K+K- || K-K+ [3][8]:
2723  if((Kaon(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE)) || (Kaon(gtrack1,-1,kTRUE) && Kaon(gtrack2,1,kTRUE)))
2724  {
2725  fBackground[3][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2726  }
2727 
2728  // c) proton-proton:
2729  // c1) p+p+ [4][4]:
2730  if(Proton(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
2731  {
2732  fBackground[4][4]->Fill(RelativeMomenta(gtrack1,gtrack2));
2733  }
2734  // c2) p-p- [9][9]:
2735  if(Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2736  {
2737  fBackground[9][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2738  }
2739  // c3) p+p- || p-p+ [4][9]:
2740  if((Proton(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE)) || (Proton(gtrack1,-1,kTRUE) && Proton(gtrack2,1,kTRUE)))
2741  {
2742  fBackground[4][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2743  }
2744 
2745  // 2.) Mixed particle species:
2746  // a) pion-kaon
2747  // a1) pi+K+ [2][3]:
2748  if(Pion(gtrack1,1,kTRUE) && Kaon(gtrack2,1,kTRUE))
2749  {
2750  fBackground[2][3]->Fill(RelativeMomenta(gtrack1,gtrack2));
2751  }
2752  // a2) pi+K- [2][8]
2753  if(Pion(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
2754  {
2755  fBackground[2][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2756  }
2757  // a3) K+pi- [3][7]
2758  if(Kaon(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE))
2759  {
2760  fBackground[3][7]->Fill(RelativeMomenta(gtrack1,gtrack2));
2761  }
2762  // a4) pi-K- [7][8]
2763  if(Pion(gtrack1,-1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
2764  {
2765  fBackground[7][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2766  }
2767  // b) pion-proton
2768  // b1) pi+p+ [2][4]:
2769  if(Pion(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
2770  {
2771  fBackground[2][4]->Fill(RelativeMomenta(gtrack1,gtrack2));
2772  }
2773  // b2) pi+p- [2][9]
2774  if(Pion(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2775  {
2776  fBackground[2][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2777  }
2778  // b3) p+pi- [4][7]
2779  if(Proton(gtrack1,1,kTRUE) && Pion(gtrack2,-1,kTRUE))
2780  {
2781  fBackground[4][7]->Fill(RelativeMomenta(gtrack1,gtrack2));
2782  }
2783  // b4) pi-p- [7][9]
2784  if(Pion(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2785  {
2786  fBackground[7][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2787  }
2788  // c) kaon-proton
2789  // c1) K+p+ [3][4]:
2790  if(Kaon(gtrack1,1,kTRUE) && Proton(gtrack2,1,kTRUE))
2791  {
2792  fBackground[3][4]->Fill(RelativeMomenta(gtrack1,gtrack2));
2793  }
2794  // c2) K+p- [3][9]
2795  if(Kaon(gtrack1,1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2796  {
2797  fBackground[3][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2798  }
2799  // c3) p+K- [4][8]
2800  if(Proton(gtrack1,1,kTRUE) && Kaon(gtrack2,-1,kTRUE))
2801  {
2802  fBackground[4][8]->Fill(RelativeMomenta(gtrack1,gtrack2));
2803  }
2804  // c4) K-p- [8][9]
2805  if(Kaon(gtrack1,-1,kTRUE) && Proton(gtrack2,-1,kTRUE))
2806  {
2807  fBackground[8][9]->Fill(RelativeMomenta(gtrack1,gtrack2));
2808  }
2809 
2810  } // for(Int_t iTrack2=0;iTrack2<nTracks2;iTrack2++)
2811  } // for(Int_t iTrack1=0;iTrack1<nTracks1;iTrack1++)
2812 
2813  // d) Shift [1] -> [0]:
2814  // TBI re-think the lines below, there should be a better way...
2815  fMixedEvents[0] = (TClonesArray*)fMixedEvents[1]->Clone();
2816  fGlobalTracksAOD[1] = (TExMap*)fGlobalTracksAOD[2]->Clone();
2817 
2818  // e) Clean [1]:
2819  fMixedEvents[1] = NULL;
2820  fGlobalTracksAOD[2]->Delete(); // TBI or = NULL ?
2821 
2822 } // void AliAnalysisTaskMultiparticleFemtoscopy::CalculateBackground(TClonesArray *ca1, TClonesArray *ca2)
2823 
2824 //=======================================================================================================================
2825 
2826 void AliAnalysisTaskMultiparticleFemtoscopy::CalculateBackground(TClonesArray *ca1, TClonesArray *ca2, Bool_t bMC)
2827 {
2828  // Calculate background for MC events.
2829 
2830  // TBI this method is not really validated
2831  // TBI unify with previous method
2832 
2833  // a) Insanity checks;
2834  // b) Temporary primitive algorithm: correlate particles from 'previous event' to particles from 'current event';
2835  // c) Two nested loops to calculate B(k);
2836  // d) Shift [1] -> [0];
2837  // e) Clean [1].
2838 
2839  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
2840 
2841  // a) Insanity checks:
2842  if(!bMC) return; // TBI this is not really needed...
2843  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::CalculateBackground(TClonesArray *ca1, TClonesArray *ca2, Bool_t bMC)";
2844  if(!ca1){Fatal(sMethodName.Data(),"!ca1");}
2845  if(!ca2){Fatal(sMethodName.Data(),"!ca2");}
2846 
2847  // b) Temporary primitive algorithm: correlate particles from 'previous event' to particles from 'current event'
2848  // ...
2849 
2850  // c) Two nested loops to calculate B(k):
2851  Int_t nTracks1 = ca1->GetEntries();
2852  Int_t nTracks2 = ca2->GetEntries();
2853 
2854  // Start loop over tracks in the 1st event:
2855  for(Int_t iTrack1=0;iTrack1<nTracks1;iTrack1++)
2856  {
2857  AliAODMCParticle *amcparticle1 = dynamic_cast<AliAODMCParticle*>(ca1->UncheckedAt(iTrack1)); // TBI cross-check UncheckedAt
2858  // TBI Temporary track insanity checks:
2859  if(!amcparticle1){Fatal(sMethodName.Data(),"!amcparticle1");} // TBI keep this for some time, eventually just continue
2860  if(!PassesCommonTrackCuts(amcparticle1)){continue;} // TBI re-think
2861 
2862  // Start loop over tracks in the 2nd event:
2863  for(Int_t iTrack2=0;iTrack2<nTracks2;iTrack2++)
2864  {
2865  AliAODMCParticle *amcparticle2 = dynamic_cast<AliAODMCParticle*>(ca2->UncheckedAt(iTrack2)); // TBI cross-check UncheckedAt
2866  // TBI Temporary track insanity checks:
2867  if(!amcparticle2){Fatal(sMethodName.Data(),"!amcparticle2");} // TBI keep this for some time, eventually just continue
2868  if(!PassesCommonTrackCuts(amcparticle2)){continue;} // TBI re-think
2869 
2870  // Okay... So we have two tracks from two different events ready. Let's check PID, and calculate the background:
2871  if(!(pattern.Contains(Form(".%d.",TMath::Abs(amcparticle1->GetPdgCode()))) && pattern.Contains(Form(".%d.",TMath::Abs(amcparticle2->GetPdgCode()))))){continue;}
2872 
2873  // Determine the indices of correlation function to be filled:
2874  Int_t index1 = -44;
2875  Int_t index2 = -44;
2876  if(fCorrelationFunctionsIndices->GetValue(amcparticle1->GetPdgCode())<=fCorrelationFunctionsIndices->GetValue(amcparticle2->GetPdgCode()))
2877  {
2878  index1 = fCorrelationFunctionsIndices->GetValue(amcparticle1->GetPdgCode());
2879  index2 = fCorrelationFunctionsIndices->GetValue(amcparticle2->GetPdgCode());
2880  }
2881  else
2882  {
2883  index1 = fCorrelationFunctionsIndices->GetValue(amcparticle2->GetPdgCode());
2884  index2 = fCorrelationFunctionsIndices->GetValue(amcparticle1->GetPdgCode());
2885  }
2886  if(-44 == index1){Fatal(sMethodName.Data(),"-44 == index1");}
2887  if(-44 == index2){Fatal(sMethodName.Data(),"-44 == index2");}
2888 
2889  // Fill or die:
2890  fBackground[index1][index2]->Fill(RelativeMomenta(amcparticle1,amcparticle2)); // for the relative momenta ordering doesn't matter
2891 
2892  } // for(Int_t iTrack2=0;iTrack2<nTracks2;iTrack2++)
2893  } // for(Int_t iTrack1=0;iTrack1<nTracks1;iTrack1++)
2894 
2895  // d) Shift [1] -> [0]:
2896  // TBI re-think the lines below, there should be a better way...
2897  fMixedEvents[0] = (TClonesArray*)fMixedEvents[1]->Clone();
2898 
2899  // e) Clean [1]:
2900  fMixedEvents[1] = NULL;
2901 
2902 } // void AliAnalysisTaskMultiparticleFemtoscopy::CalculateBackground(TClonesArray *ca1, TClonesArray *ca2, Bool_t bMC)
2903 
2904 //=======================================================================================================================
2905 
2907 {
2908  // Get pointers for everything in the base list "fHistList".
2909 
2910  // a) Get pointer for base list fHistList;
2911  // b) Get pointer for profile holding internal flags and set again all flags TBI this profile is not implemented yet;
2912  // *) TBI ...
2913  // *) Get pointers for correlation functions;
2914  // *) Get pointers for background.
2915 
2916  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GetOutputHistograms(TList *histList)";
2917 
2918  // a) Get pointer for base list fHistList and profile holding internal flags;
2919  fHistList = histList;
2920  if(!fHistList){Fatal(sMethodName.Data(),"fHistList is not around today...");}
2921 
2922  // b) Get pointer for profile holding internal flags and set again all flags TBI this profile is not implemented yet
2923  //fInternalFlagsPro = dynamic_cast<TProfile*>(fHistList->FindObject("fInternalFlagsPro")); TBI this was example from MPC
2924  //if(!fInternalFlagsPro){Fatal(sMethodName.Data(),"fInternalFlagsPro");} TBI this was example from MPC
2925 
2926  // *) TBI
2927  //this->GetPointersFor...;
2928 
2929  // *) Get pointers for correlation functions:
2931 
2932  // *) Get pointers for background:
2933  this->GetPointersForBackground();
2934 
2935 } // void AliAnalysisTaskMultiparticleFemtoscopy::GetOutputHistograms(TList *histList)
2936 
2937 //=======================================================================================================================
2938 
2940 {
2941  // Get pointers for correlation functions.
2942 
2943  // a) Get pointer for fCorrelationFunctionsList;
2944  // b) Get pointer for fCorrelationFunctionsFlagsPro;
2945  // c) Set again all flags;
2946  // d) Get pointers for fCorrelationFunctions[10][10].
2947 
2948  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GetPointersForCorrelationFunctions()";
2949 
2950  // a) Get pointer for fCorrelationFunctionsList:
2951  fCorrelationFunctionsList = dynamic_cast<TList*>(fHistList->FindObject("Correlation_Functions"));
2952  if(!fCorrelationFunctionsList){Fatal(sMethodName.Data(),"fCorrelationFunctionsList");}
2953 
2954  // b) Get pointer for fCorrelationFunctionsFlagsPro:
2955  fCorrelationFunctionsFlagsPro = dynamic_cast<TProfile*>(fCorrelationFunctionsList->FindObject("fCorrelationFunctionsFlagsPro"));
2956  if(!fCorrelationFunctionsFlagsPro){Fatal(sMethodName.Data(),"fCorrelationFunctionsFlagsPro");}
2957 
2958  // c) Set again all flags:
2961 
2962  if(!fFillCorrelationFunctions){return;} // TBI is this safe enough
2963 
2964  // d) Get pointers for fCorrelationFunctions[10][10]:
2965  const Int_t nParticleSpecies = 5; // Supported at the moment: 0=e,1=mu,2=pi,3=K,4=p
2966  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]
2967  {
2968  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]
2969  {
2970  fCorrelationFunctions[pid1][pid2] = dynamic_cast<TH1F*>(fCorrelationFunctionsList->FindObject(Form("fCorrelationFunctions[%d][%d]",pid1,pid2)));
2971  if(!fCorrelationFunctions[pid1][pid2]){Fatal(sMethodName.Data(),"fCorrelationFunctions[%d][%d]",pid1,pid2);}
2972  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
2973  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
2974 
2975 } // void AliAnalysisTaskMultiparticleFemtoscopy::GetPointersForCorrelationFunctions()
2976 
2977 //=======================================================================================================================
2978 
2980 {
2981  // Get pointers for background.
2982 
2983  // a) Get pointer for fBackgroundList;
2984  // b) Get pointer for fBackgroundFlagsPro;
2985  // c) Set again all flags;
2986  // d) Get pointers for fBackground[10][10].
2987 
2988  TString sMethodName = "void AliAnalysisTaskMultiparticleFemtoscopy::GetPointersForBackground()";
2989 
2990  // a) Get pointer for fBackgroundList:
2991  fBackgroundList = dynamic_cast<TList*>(fHistList->FindObject("Background"));
2992  if(!fBackgroundList){Fatal(sMethodName.Data(),"fBackgroundList");}
2993 
2994  // b) Get pointer for fBackgroundFlagsPro:
2995  fBackgroundFlagsPro = dynamic_cast<TProfile*>(fBackgroundList->FindObject("fBackgroundFlagsPro"));
2996  if(!fBackgroundFlagsPro){Fatal(sMethodName.Data(),"fBackgroundFlagsPro");}
2997 
2998  // c) Set again all flags:
2999  fEstimateBackground = (Bool_t) fBackgroundFlagsPro->GetBinContent(1);
3000 
3001  if(!fEstimateBackground){return;} // TBI is this safe enough
3002 
3003  // d) Get pointers for fBackground[10][10]:
3004  const Int_t nParticleSpecies = 5; // Supported at the moment: 0=e,1=mu,2=pi,3=K,4=p
3005  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]
3006  {
3007  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]
3008  {
3009  fBackground[pid1][pid2] = dynamic_cast<TH1F*>(fBackgroundList->FindObject(Form("fBackground[%d][%d]",pid1,pid2)));
3010  if(!fBackground[pid1][pid2]){Fatal(sMethodName.Data(),"fBackground[%d][%d]",pid1,pid2);}
3011  } // for(Int_t pid2=0;pid2<5;pid2++) // anti-particle [0=e,1=mu,2=pi,3=K,4=p]
3012  } // for(Int_t pid=0;pid<5;pid++) // particle [0=e,1=mu,2=pi,3=K,4=p]
3013 
3014 } // void AliAnalysisTaskMultiparticleFemtoscopy::GetPointersForBackground()
3015 
3016 //=======================================================================================================================
3017 
3018 
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
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]...
Definition: External.C:212
Bool_t Pion(AliAODTrack *atrack, Int_t charge=1, Bool_t bPrimary=kTRUE)
UShort_t fBunchCross
do something only for the specified event
Bool_t Kaon(AliAODTrack *atrack, Int_t charge=1, Bool_t bPrimary=kTRUE)
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)
virtual void FillControlHistogramsNonIdentifiedParticles(AliAODTrack *atrack)
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)
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
TClonesArray * fMixedEvents[2]
[particle(+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