Bug Correction
[u/mrichter/AliRoot.git] / HBTAN / AliHBTAnalysis.cxx
1 #include "AliHBTAnalysis.h"
2 //_________________________________________________________
3 ////////////////////////////////////////////////////////////////////////////
4 //
5 // class AliHBTAnalysis
6 //
7 // Central Object Of HBTAnalyser: 
8 // This class performs main looping within HBT Analysis
9 // User must plug a reader of Type AliReader
10 // User plugs in coorelation and monitor functions
11 // as well as monitor functions
12 //
13 // HBT Analysis Tool, which is integral part of AliRoot,
14 // ALICE Off-Line framework:
15 //
16 // Piotr.Skowronski@cern.ch
17 // more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
18 //
19 ////////////////////////////////////////////////////////////////////////////
20 //_________________________________________________________
21
22
23 #include <TSystem.h>
24 #include <TFile.h>
25
26 #include "AliAOD.h"
27 #include "AliAODParticle.h"
28 #include "AliAODPairCut.h"
29 #include "AliEventCut.h"
30
31 #include "AliEventBuffer.h"
32
33 #include "AliReader.h"
34 #include "AliHBTPair.h"
35 #include "AliHBTFunction.h"
36 #include "AliHBTMonitorFunction.h"
37  
38
39 ClassImp(AliHBTAnalysis)
40
41 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
42 const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
43 const Int_t  AliHBTAnalysis::fgkDefaultBufferSize = 5;
44
45 AliHBTAnalysis::AliHBTAnalysis():
46   fProcEvent(0x0),
47   fReader(0x0),
48   fNTrackFunctions(0),
49   fNParticleFunctions(0),
50   fNParticleAndTrackFunctions(0),
51   fNTrackMonitorFunctions(0),
52   fNParticleMonitorFunctions(0), 
53   fNParticleAndTrackMonitorFunctions(0),
54   fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
55   fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
56   fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
57   fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),    
58   fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),    
59   fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),    
60   fBkgEventCut(0x0),
61   fPartBuffer(0x0),
62   fTrackBuffer(0x0),
63   fBufferSize(2),
64   fDisplayMixingInfo(fgkDefaultMixingInfo),
65   fIsOwner(kFALSE),
66   fProcessOption(kSimulatedAndReconstructed),
67   fNoCorrfctns(kFALSE),
68   fOutputFileName(0x0),
69   fVertexX(0.0),
70   fVertexY(0.0),
71   fVertexZ(0.0)
72  {
73    //default constructor
74    
75  }
76 /*************************************************************************************/ 
77
78 AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
79   AliAnalysis(in),
80   fProcEvent(0x0),
81   fReader(0x0),
82   fNTrackFunctions(0),
83   fNParticleFunctions(0),
84   fNParticleAndTrackFunctions(0),
85   fNTrackMonitorFunctions(0),
86   fNParticleMonitorFunctions(0), 
87   fNParticleAndTrackMonitorFunctions(0),
88   fTrackFunctions(0x0),
89   fParticleFunctions(0x0),
90   fParticleAndTrackFunctions(0x0),
91   fParticleMonitorFunctions(0x0),
92   fTrackMonitorFunctions(0x0),
93   fParticleAndTrackMonitorFunctions(0x0),
94   fBkgEventCut(0x0),
95   fPartBuffer(0x0),
96   fTrackBuffer(0x0),
97   fBufferSize(fgkDefaultBufferSize),
98   fDisplayMixingInfo(fgkDefaultMixingInfo),
99   fIsOwner(kFALSE),
100   fProcessOption(kSimulatedAndReconstructed),
101   fNoCorrfctns(kFALSE),
102   fOutputFileName(0x0),
103   fVertexX(0.0),
104   fVertexY(0.0),
105   fVertexZ(0.0)
106  {
107 //copy constructor
108    Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
109  }
110 /*************************************************************************************/ 
111 AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
112  {
113 //operator =
114    Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
115    return *this;
116  }
117 /*************************************************************************************/ 
118 AliHBTAnalysis::~AliHBTAnalysis()
119  {
120  //destructor
121  //note that we do not delete functions itself
122  // they should be deleted by whom where created
123  //we only store pointers, and we use "new" only for pointers array
124     
125    if (fIsOwner)
126     {
127       if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Is Owner: Attempting to delete functions");
128       DeleteFunctions();
129       if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
130     }
131    delete [] fTrackFunctions;
132    delete [] fParticleFunctions;
133    delete [] fParticleAndTrackFunctions;
134    
135    delete [] fParticleMonitorFunctions; 
136    delete [] fTrackMonitorFunctions; 
137    delete [] fParticleAndTrackMonitorFunctions;
138
139    delete fBkgEventCut;
140    delete fOutputFileName;
141  }
142
143 /*************************************************************************************/ 
144
145 Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
146 {
147   //Processes one event
148   if (fProcEvent == 0x0)
149    {
150      Error("ProcessEvent","Analysis <<%s>> option not specified.",GetName());
151      return 1;
152    }
153   if ( Rejected(aodrec,aodsim) ) 
154    {
155 //     Double_t x = 0, y = 0, z = 0;
156 //     aodrec->GetPrimaryVertex(x,y,z);
157 //     Info("ProcessEvent","Event has vertex at %f %f %f",x,y,z);
158      Info("ProcessEvent","Nch is %d",aodsim->GetNumberOfCharged());
159      Info("ProcessEvent","%s: Event cut rejected this event",GetName());
160      return 0;
161    } 
162   
163   //Move event to the apparent vertex -> must be after the event cut
164   //It is important for any cut that use any spacial coordiantes, 
165   //f.g. anti-merging cut in order to preserve the same bahavior of variables (f.g. distance between tracks)
166   Double_t dvx = 0, dvy = 0, dvz = 0;
167   if (aodrec)
168    {
169      Double_t pvx,pvy,pvz;
170      aodrec->GetPrimaryVertex(pvx,pvy,pvz);
171      
172      dvx = fVertexX - pvx;
173      dvy = fVertexY - pvy;
174      dvz = fVertexZ - pvz;
175      aodrec->Move(dvx,dvy,dvz);
176    }  
177   
178   Int_t result = (this->*fProcEvent)(aodrec,aodsim);
179
180   if (aodrec) aodrec->Move(-dvx,-dvy,-dvz);//move event back to the same spacial coordinates
181   
182   return  result;
183 }
184 /*************************************************************************************/ 
185
186 Int_t AliHBTAnalysis::Finish()
187 {
188 //Finishes analysis
189   WriteFunctions();
190   return 0;
191 }
192 /*************************************************************************************/ 
193
194 void AliHBTAnalysis::DeleteFunctions()
195 {
196  //Deletes all functions added to analysis
197  
198  UInt_t ii;
199  for(ii = 0;ii<fNParticleFunctions;ii++)
200   { 
201     if (AliVAODParticle::GetDebug()>5)
202      {
203        Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]);
204        Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name());
205      }
206     delete fParticleFunctions[ii];
207   } 
208  fNParticleFunctions = 0;
209                 
210  for(ii = 0;ii<fNTrackFunctions;ii++)
211   { 
212     if (AliVAODParticle::GetDebug()>5)
213      {
214        Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]);
215        Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name());
216      }
217     delete fTrackFunctions[ii];
218   }  
219  fNTrackFunctions = 0;
220  
221  for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
222   { 
223     if (AliVAODParticle::GetDebug()>5)
224      {
225        Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
226        Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
227      } 
228     delete fParticleAndTrackFunctions[ii];
229   }  
230  fNParticleAndTrackFunctions = 0;
231  
232  for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
233   { 
234     if (AliVAODParticle::GetDebug()>5)
235      {
236        Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
237        Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
238      }
239     delete fParticleMonitorFunctions[ii];
240   } 
241  fNParticleMonitorFunctions = 0;
242    
243  for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
244   { 
245     if (AliVAODParticle::GetDebug()>5)
246      {
247        Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
248        Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
249      }
250     delete fTrackMonitorFunctions[ii];
251   } 
252  fNTrackMonitorFunctions = 0;
253    
254  for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
255   { 
256     if (AliVAODParticle::GetDebug()>5)
257      {
258       Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
259       Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
260      }
261     delete fParticleAndTrackMonitorFunctions[ii];
262   } 
263  fNParticleAndTrackMonitorFunctions = 0;
264  
265 }
266 /*************************************************************************************/ 
267
268 Int_t AliHBTAnalysis::Init()
269 {
270 //Initializeation method
271 //calls Init for all functions
272
273  //Depending on option and pair cut assigns proper analysis method
274  Bool_t nonid = IsNonIdentAnalysis();
275  switch(fProcessOption)
276   {
277     case kReconstructed:
278       if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecNonId;
279       else fProcEvent = &AliHBTAnalysis::ProcessRec;
280       SetCutsOnRec();
281       break;
282
283     case kSimulated:
284       if (nonid) fProcEvent = &AliHBTAnalysis::ProcessSimNonId;
285       else fProcEvent = &AliHBTAnalysis::ProcessSim;
286       SetCutsOnSim();
287       break;
288
289     case kSimulatedAndReconstructed:
290       if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecAndSimNonId;
291       else fProcEvent = &AliHBTAnalysis::ProcessRecAndSim;
292       break;
293   }
294   
295  if (fPartBuffer == 0x0)  fPartBuffer = new AliEventBuffer (fBufferSize);
296  else fPartBuffer->Reset();
297
298  if (fTrackBuffer == 0x0)  fTrackBuffer = new AliEventBuffer (fBufferSize);
299  else fTrackBuffer->Reset();
300
301
302  fNoCorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
303  
304  UInt_t ii;
305  for(ii = 0;ii<fNParticleFunctions;ii++)
306    fParticleFunctions[ii]->Init();
307                 
308  for(ii = 0;ii<fNTrackFunctions;ii++)
309    fTrackFunctions[ii]->Init();
310
311  for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
312    fParticleAndTrackFunctions[ii]->Init();
313  
314  for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
315    fParticleMonitorFunctions[ii]->Init();
316    
317  for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
318    fTrackMonitorFunctions[ii]->Init();
319    
320  for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
321    fParticleAndTrackMonitorFunctions[ii]->Init();
322    
323  return 0;
324 }
325 /*************************************************************************************/ 
326
327 void AliHBTAnalysis::ResetFunctions()
328 {
329 //In case fOwner is true, deletes all functions
330 //in other case, just set number of analysis to 0
331  if (fIsOwner) DeleteFunctions();
332  else
333   {
334     fNParticleFunctions = 0;
335     fNTrackFunctions = 0;
336     fNParticleAndTrackFunctions = 0;
337     fNParticleMonitorFunctions = 0;
338     fNTrackMonitorFunctions = 0;
339     fNParticleAndTrackMonitorFunctions = 0;
340   }
341 }
342 /*************************************************************************************/ 
343 Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim)
344 {
345 //Does analysis for both tracks and particles
346 //mainly for resolution study and analysies with weighting algirithms
347   
348 // cut on particles only -- why?
349 // - PID: when we make resolution analysis we want to take only tracks with correct PID
350 // We need cut on tracks because there are data characteristic
351   
352   AliVAODParticle * part1, * part2;
353   AliVAODParticle * track1, * track2;
354   
355   AliAOD * trackEvent = aodrec, *partEvent = aodsim;
356   AliAOD* trackEvent1 = new AliAOD();
357   AliAOD* partEvent1 = new AliAOD();
358   partEvent1->SetOwner(kTRUE);
359   trackEvent1->SetOwner(kTRUE);
360
361   AliAOD * trackEvent2,*partEvent2;
362   
363 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
364   
365 //  Int_t nev = fReader->GetNumberOfTrackEvents();
366   static AliHBTPair tpair;
367   static AliHBTPair ppair;
368
369   AliHBTPair* trackpair = &tpair;
370   AliHBTPair* partpair = &ppair;
371  
372   AliHBTPair * tmptrackpair;//temprary pointers to pairs
373   AliHBTPair * tmppartpair;
374   
375   register UInt_t ii;
376   
377   
378
379   if ( !partEvent || !trackEvent ) 
380    {
381      Error("ProcessRecAndSim","<<%s>> Can not get event",GetName());
382      return 1;
383    }
384
385   if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
386    {
387      Error("ProcessRecAndSim",
388            "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d). Skipping Event.",
389             partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
390      return 2;
391    }
392
393
394   for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
395    {
396      /***************************************/
397      /******   Looping same events   ********/
398      /******   filling numerators    ********/
399      /***************************************/
400      if ( (j%fDisplayMixingInfo) == 0)
401         Info("ProcessTracksAndParticles",
402              "Mixing particle %d with particles from the same event",j);
403
404      part1= partEvent->GetParticle(j);
405      track1= trackEvent->GetParticle(j);
406
407      Bool_t firstcut = (this->*fkPass1)(part1,track1);
408      if (fBufferSize != 0) 
409        if ( (firstcut == kFALSE) || ( (this->*fkPass2)(part1,track1) == kFALSE ) )
410         {
411           //accepted by any cut
412           // we have to copy because reader keeps only one event
413
414           partEvent1->AddParticle(new AliAODParticle(*part1));
415           trackEvent1->AddParticle(new AliAODParticle(*track1));
416         }
417
418      if (firstcut) continue;
419
420      for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
421        fParticleMonitorFunctions[ii]->Process(part1);
422      for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
423        fTrackMonitorFunctions[ii]->Process(track1);
424      for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
425        fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
426
427      if (fNoCorrfctns) continue;
428
429      for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
430       {
431         part2= partEvent->GetParticle(k);
432         if (part1->GetUID() == part2->GetUID()) continue;
433         partpair->SetParticles(part1,part2);
434
435         track2= trackEvent->GetParticle(k);
436         trackpair->SetParticles(track1,track2);
437
438         if( (this->*fkPass)(partpair,trackpair) ) //check pair cut 
439           { //do not meets crietria of the pair cut, try with swapped pairs
440             if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
441               continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
442             else 
443              { //swaped pair meets all the criteria
444                tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
445                tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
446              }
447           }
448         else
449          {//meets criteria of the pair cut
450            tmptrackpair = trackpair;
451            tmppartpair = partpair;
452          }
453
454         for(ii = 0;ii<fNParticleFunctions;ii++)
455                fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
456
457         for(ii = 0;ii<fNTrackFunctions;ii++)
458                fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
459
460         for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
461                fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
462        //end of 2nd loop over particles from the same event  
463       }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
464
465      /***************************************/
466      /***** Filling denominators    *********/
467      /***************************************/
468      if (fBufferSize == 0) continue;
469
470      fPartBuffer->ResetIter();
471      fTrackBuffer->ResetIter();
472      Int_t m = 0;
473      while (( partEvent2 = fPartBuffer->Next() ))
474       {
475         trackEvent2 = fTrackBuffer->Next();
476
477         m++;
478         if ( (j%fDisplayMixingInfo) == 0)
479            Info("ProcessTracksAndParticles",
480                 "Mixing particle %d from current event with particles from event %d",j,-m);
481
482         for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
483           {
484             part2= partEvent2->GetParticle(l);
485             partpair->SetParticles(part1,part2);
486
487             track2= trackEvent2->GetParticle(l);
488             trackpair->SetParticles(track1,track2);
489
490             if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
491               { //do not meets crietria of the 
492                 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
493                   continue;
494                 else 
495                  {
496                    tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
497                    tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
498                  }
499               }
500             else
501              {//meets criteria of the pair cut
502               tmptrackpair = trackpair;
503               tmppartpair = partpair;
504              }
505
506             for(ii = 0;ii<fNParticleFunctions;ii++)
507               fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
508
509             for(ii = 0;ii<fNTrackFunctions;ii++)
510               fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
511
512             for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
513               fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
514           }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
515
516       }
517     //end of loop over particles from first event
518    }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
519   delete fPartBuffer->Push(partEvent1);
520   delete fTrackBuffer->Push(trackEvent1);
521  //end of loop over events  
522   return 0;
523 }
524 /*************************************************************************************/ 
525
526 Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim)
527 {
528   //Does analysis of simulated data
529   AliVAODParticle * part1, * part2;
530   
531   AliAOD* partEvent = aodsim;
532   AliAOD* partEvent1 = new AliAOD();
533   partEvent1->SetOwner(kTRUE);
534
535   AliAOD* partEvent2;
536   
537   AliHBTPair ppair;
538
539   AliHBTPair* partpair = &ppair;
540  
541   AliHBTPair * tmppartpair;
542   
543   register UInt_t ii;
544   
545   
546
547   if ( !partEvent )
548    {
549      Error("ProcessRecAndSim","Can not get event");
550      return 1;
551    }
552
553
554   for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
555    {
556      /***************************************/
557      /******   Looping same events   ********/
558      /******   filling numerators    ********/
559      /***************************************/
560      if ( (j%fDisplayMixingInfo) == 0)
561         Info("ProcessTracksAndParticles",
562              "Mixing particle %d with particles from the same event",j);
563
564      part1= partEvent->GetParticle(j);
565
566      Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(part1);
567
568      if (fBufferSize != 0) 
569        if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(part1) == kFALSE ) )
570         {
571           //accepted by any cut
572           // we have to copy because reader keeps only one event
573
574           partEvent1->AddParticle(new AliAODParticle(*part1));
575         }
576
577      if (firstcut) continue;
578
579      for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
580        fParticleMonitorFunctions[ii]->Process(part1);
581
582      if ( fNParticleFunctions == 0 ) continue;
583
584      for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
585       {
586         part2= partEvent->GetParticle(k);
587         if (part1->GetUID() == part2->GetUID()) continue;
588         partpair->SetParticles(part1,part2);
589
590            if(fPairCut->Rejected(partpair)) //check pair cut
591             { //do not meets crietria of the 
592               if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) ) continue;
593               else tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
594             }
595            else
596             {
597               tmppartpair = partpair;
598             }
599
600         for(ii = 0;ii<fNParticleFunctions;ii++)
601                fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
602
603        //end of 2nd loop over particles from the same event  
604       }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
605
606      /***************************************/
607      /***** Filling denominators    *********/
608      /***************************************/
609      if (fBufferSize == 0) continue;
610
611      fPartBuffer->ResetIter();
612          Int_t m = 0;
613          while (( partEvent2 = fPartBuffer->Next() ))
614           {
615             m++;
616             if ( (j%fDisplayMixingInfo) == 0)
617                Info("ProcessParticles",
618                     "Mixing particle %d from current event with particles from event %d",j,-m);
619             for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
620               {
621
622                 part2= partEvent2->GetParticle(l);
623                 partpair->SetParticles(part1,part2);
624
625                 if( fPairCut->Rejected(partpair) ) //check pair cut
626                   { //do not meets crietria of the 
627                     if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) )
628                       continue;
629                     else 
630                      {
631                        tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
632                      }
633                   }
634                 else
635                  {//meets criteria of the pair cut
636                   tmppartpair = partpair;
637                  }
638                  
639                 for(ii = 0;ii<fNParticleFunctions;ii++)
640                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
641                  
642              }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
643           }
644        }
645   delete fPartBuffer->Push(partEvent1);
646  //end of loop over events  
647   return 0;
648 }
649 /*************************************************************************************/ 
650 Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/)
651 {
652   //Does analysis of reconstructed data
653   AliVAODParticle * track1, * track2;
654   
655   AliAOD* trackEvent = aodrec;
656   AliAOD* trackEvent1 = new AliAOD();
657   trackEvent1->SetOwner(kTRUE);
658
659   AliAOD* trackEvent2;
660   
661   AliHBTPair tpair;
662
663   AliHBTPair* trackpair = &tpair;
664  
665   AliHBTPair * tmptrackpair;
666   
667   register UInt_t ii;
668   
669
670   if ( !trackEvent )
671    {
672      Error("ProcessRecAndSim","Can not get event");
673      return 1;
674    }
675
676
677   for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
678    {
679      /***************************************/
680      /******   Looping same events   ********/
681      /******   filling numerators    ********/
682      /***************************************/
683      if ( (j%fDisplayMixingInfo) == 0)
684         Info("ProcessTracksAndParticles",
685              "Mixing Particle %d with Particles from the same event",j);
686
687      track1= trackEvent->GetParticle(j);
688
689      Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(track1);
690
691      if (fBufferSize != 0) 
692        if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(track1) == kFALSE ) )
693         {
694           //accepted by any cut
695           // we have to copy because reader keeps only one event
696
697           trackEvent1->AddParticle(new AliAODParticle(*track1));
698         }
699
700      if (firstcut) continue;
701
702      for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
703        fParticleMonitorFunctions[ii]->Process(track1);
704
705      if ( fNParticleFunctions == 0 ) continue;
706
707      for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
708       {
709         track2= trackEvent->GetParticle(k);
710         if (track1->GetUID() == track2->GetUID()) continue;
711         trackpair->SetParticles(track1,track2);
712
713            if(fPairCut->Rejected(trackpair)) //check pair cut
714             { //do not meets crietria of the 
715               if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) ) continue;
716               else tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
717             }
718            else
719             {
720               tmptrackpair = trackpair;
721             }
722
723         for(ii = 0;ii<fNTrackFunctions;ii++)
724                fParticleFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
725
726        //end of 2nd loop over Particles from the same event  
727       }//for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
728
729      /***************************************/
730      /***** Filling denominators    *********/
731      /***************************************/
732      if (fBufferSize == 0) continue;
733
734      fTrackBuffer->ResetIter();
735          Int_t m = 0;
736          while (( trackEvent2 = fTrackBuffer->Next() ))
737           {
738             m++;
739             if ( (j%fDisplayMixingInfo) == 0)
740                Info("ProcessParticles",
741                     "Mixing Particle %d from current event with Particles from event %d",j,-m);
742             for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++)   //  ... on all Particles
743               {
744
745                 track2= trackEvent2->GetParticle(l);
746                 trackpair->SetParticles(track1,track2);
747
748                 if( fPairCut->Rejected(trackpair) ) //check pair cut
749                   { //do not meets crietria of the 
750                     if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) )
751                       continue;
752                     else 
753                      {
754                        tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
755                      }
756                   }
757                 else
758                  {//meets criteria of the pair cut
759                   tmptrackpair = trackpair;
760                  }
761                  
762                 for(ii = 0;ii<fNTrackFunctions;ii++)
763                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
764                  
765              }//for(Int_t l = 0; l<N2;l++)   //  ... on all Particles
766           }
767        }
768   delete fTrackBuffer->Push(trackEvent1);
769  //end of loop over events  
770   return 0;
771 }
772 /*************************************************************************************/ 
773      
774 Int_t AliHBTAnalysis::ProcessRecAndSimNonId(AliAOD* aodrec, AliAOD* aodsim)
775 {
776 //Analyzes both reconstructed and simulated data
777   if (aodrec == 0x0) 
778    {
779      Error("ProcessTracksAndParticlesNonIdentAnal","Reconstructed event is NULL");
780      return 1;
781    }
782
783   if (aodsim == 0x0) 
784    {
785      Error("ProcessTracksAndParticlesNonIdentAnal","Simulated event is NULL");
786      return 1;
787    }
788
789   if ( aodrec->GetNumberOfParticles() != aodsim->GetNumberOfParticles() )
790    {
791      Error("ProcessTracksAndParticlesNonIdentAnal",
792            "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
793            aodsim->GetNumberOfParticles() , aodrec->GetNumberOfParticles());
794      return 2;
795    }
796
797  
798   AliVAODParticle * part1, * part2;
799   AliVAODParticle * track1, * track2;
800
801   static AliAOD aodrec1;
802   static AliAOD aodsim1;
803
804   AliAOD * trackEvent1=&aodrec1,*partEvent1=&aodsim1;//Particle that passes first particle cut, this event
805   AliAOD * trackEvent2=0x0,*partEvent2=0x0;//Particle that passes second particle cut, this event
806   AliAOD * trackEvent3=0x0,*partEvent3=0x0;//Particle that passes second particle cut, events from buffer
807
808   AliAOD* rawtrackEvent = aodrec;//this we get from Reader
809   AliAOD* rawpartEvent = aodsim;//this we get from Reader
810
811   static AliHBTPair tpair;
812   static AliHBTPair ppair;
813   
814   AliHBTPair* trackpair = &tpair;
815   AliHBTPair* partpair = &ppair;
816
817   register UInt_t ii;
818   
819   
820   trackEvent1 = new AliAOD();
821   partEvent1 = new AliAOD();
822   trackEvent1->SetOwner(kFALSE);
823   partEvent1->SetOwner(kFALSE);;
824   
825   /********************************/
826   /*      Filtering out           */
827   /********************************/
828   if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
829    {
830      partEvent2  = new AliAOD();
831      trackEvent2 = new AliAOD();
832      partEvent2->SetOwner(kTRUE);
833      trackEvent2->SetOwner(kTRUE);
834    }
835
836   FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
837
838   for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
839    {
840      if ( (j%fDisplayMixingInfo) == 0) 
841         Info("ProcessTracksAndParticlesNonIdentAnal",
842              "Mixing particle %d from current event with particles from current event",j);
843
844      part1= partEvent1->GetParticle(j);
845      track1= trackEvent1->GetParticle(j);
846
847
848      for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
849        fParticleMonitorFunctions[ii]->Process(part1);
850      for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
851        fTrackMonitorFunctions[ii]->Process(track1);
852      for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
853        fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
854
855      if (fNoCorrfctns) continue;
856
857      /***************************************/
858      /******   filling numerators    ********/
859      /****** (particles from event2) ********/
860      /***************************************/
861
862      for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups 
863       {
864         part2= partEvent2->GetParticle(k);
865         if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
866         partpair->SetParticles(part1,part2);
867
868         track2= trackEvent2->GetParticle(k);
869         trackpair->SetParticles(track1,track2);
870
871         if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
872          { //do not meets crietria of the pair cut
873           continue; 
874          }
875         else
876          {//meets criteria of the pair cut
877           for(ii = 0;ii<fNParticleFunctions;ii++)
878                  fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
879
880           for(ii = 0;ii<fNTrackFunctions;ii++)
881                  fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
882
883           for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
884                  fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
885          }
886        }
887
888  if ( fBufferSize == 0) continue;//do not mix diff histograms
889  /***************************************/
890  /***** Filling denominators    *********/
891  /***************************************/
892  fPartBuffer->ResetIter();
893  fTrackBuffer->ResetIter();
894
895  Int_t nmonitor = 0;
896
897  while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
898   {
899     trackEvent3 = fTrackBuffer->Next();
900
901     if ( (j%fDisplayMixingInfo) == 0) 
902       Info("ProcessTracksAndParticlesNonIdentAnal",
903            "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
904
905     for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
906       {
907         part2= partEvent3->GetParticle(k);
908         partpair->SetParticles(part1,part2);
909
910         track2= trackEvent3->GetParticle(k);
911         trackpair->SetParticles(track1,track2);
912
913         if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
914          { //do not meets crietria of the pair cut
915           continue; 
916          }
917         else
918          {//meets criteria of the pair cut
919           UInt_t ii;
920           for(ii = 0;ii<fNParticleFunctions;ii++)
921                  fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
922
923           for(ii = 0;ii<fNTrackFunctions;ii++)
924                  fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
925
926           for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
927                  fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
928          }
929        }// for particles event2
930      }//while event2
931    }//for over particles in event1
932
933  delete fPartBuffer->Push(partEvent2);
934  delete fTrackBuffer->Push(trackEvent2);
935  
936  return 0;
937 }
938 /*************************************************************************************/ 
939 Int_t AliHBTAnalysis::ProcessSimNonId(AliAOD* /*aodrec*/, AliAOD* aodsim)
940 {
941 //does analysis of simulated (MC) data in non-identical mode
942 //i.e. when particles selected by first part. cut are a disjunctive set than particles
943 //passed by the second part. cut
944  if (aodsim == 0x0) 
945   {
946     return 1;
947   }
948  
949  
950   AliVAODParticle * part1, * part2;
951
952   static AliAOD aodsim1;
953
954   AliAOD* partEvent1=&aodsim1;//Particle that passes first particle cut, this event
955   AliAOD* partEvent2=0x0;//Particle that passes second particle cut, this event
956   AliAOD* partEvent3=0x0;//Particle that passes second particle cut, events from buffer
957
958   AliAOD* rawpartEvent = aodsim;//this we get from Reader
959
960   static AliHBTPair ppair;
961   
962   AliHBTPair* partpair = &ppair;
963
964   register UInt_t ii;
965   
966   
967   partEvent1 = new AliAOD();
968   partEvent1->SetOwner(kFALSE);;
969   
970
971   /********************************/
972   /*      Filtering out           */
973   /********************************/
974   if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
975    {
976      partEvent2  = new AliAOD();
977      partEvent2->SetOwner(kTRUE);
978    }
979
980   FilterOut(partEvent1, partEvent2, rawpartEvent);
981
982   for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
983    {
984      if ( (j%fDisplayMixingInfo) == 0) 
985         Info("ProcessParticlesNonIdentAnal",
986              "Mixing particle %d from current event with particles from current event",j);
987
988      part1= partEvent1->GetParticle(j);
989
990
991      for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
992        fParticleMonitorFunctions[ii]->Process(part1);
993
994      if (fNParticleFunctions == 0) continue;
995
996      /***************************************/
997      /******   filling numerators    ********/
998      /****** (particles from event2) ********/
999      /***************************************/
1000
1001      for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups 
1002       {
1003         part2= partEvent2->GetParticle(k);
1004         if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
1005         partpair->SetParticles(part1,part2);
1006
1007
1008         if(fPairCut->PassPairProp(partpair) ) //check pair cut
1009           { //do not meets crietria of the pair cut
1010               continue; 
1011           }
1012         else
1013          {//meets criteria of the pair cut
1014           for(ii = 0;ii<fNParticleFunctions;ii++)
1015               fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1016          }
1017        }
1018
1019  if ( fBufferSize == 0) continue;//do not mix diff histograms
1020  /***************************************/
1021  /***** Filling denominators    *********/
1022  /***************************************/
1023  fPartBuffer->ResetIter();
1024
1025  Int_t nmonitor = 0;
1026
1027  while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
1028   {
1029
1030     if ( (j%fDisplayMixingInfo) == 0) 
1031       Info("ProcessParticlesNonIdentAnal",
1032            "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1033
1034     for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1035       {
1036         part2= partEvent3->GetParticle(k);
1037         partpair->SetParticles(part1,part2);
1038
1039
1040         if(fPairCut->PassPairProp(partpair) ) //check pair cut
1041           { //do not meets crietria of the pair cut
1042               continue; 
1043           }
1044         else
1045          {//meets criteria of the pair cut
1046           for(ii = 0;ii<fNParticleFunctions;ii++)
1047            {
1048              fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1049            }
1050          }
1051        }// for particles event2
1052      }//while event2
1053    }//for over particles in event1
1054
1055  delete fPartBuffer->Push(partEvent2);
1056  
1057  return 0;
1058 }
1059 /*************************************************************************************/ 
1060 Int_t AliHBTAnalysis::ProcessRecNonId(AliAOD* aodrec, AliAOD* /*aodsim*/)
1061 {
1062 //Analyzes both reconstructed and simulated data
1063  if (aodrec == 0x0) 
1064   {
1065     return 1;
1066   }
1067  
1068   AliVAODParticle * track1, * track2;
1069
1070   static AliAOD aodrec1;
1071
1072   AliAOD * trackEvent1=&aodrec1;//Particle that passes first particle cut, this event
1073   AliAOD * trackEvent2=0x0;//Particle that passes second particle cut, this event
1074   AliAOD * trackEvent3=0x0;//Particle that passes second particle cut, events from buffer
1075
1076   AliAOD* rawtrackEvent = aodrec;//this we get from Reader
1077
1078   static AliHBTPair tpair;
1079   
1080   AliHBTPair* trackpair = &tpair;
1081
1082   register UInt_t ii;
1083   
1084   
1085   trackEvent1 = new AliAOD();
1086   trackEvent1->SetOwner(kFALSE);
1087   
1088   /********************************/
1089   /*      Filtering out           */
1090   /********************************/
1091   if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
1092    {
1093      trackEvent2 = new AliAOD();
1094      trackEvent2->SetOwner(kTRUE);
1095    }
1096
1097   FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
1098
1099   for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
1100    {
1101      if ( (j%fDisplayMixingInfo) == 0) 
1102         Info("ProcessTracksNonIdentAnal",
1103              "Mixing particle %d from current event with particles from current event",j);
1104
1105      track1= trackEvent1->GetParticle(j);
1106
1107
1108      for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
1109        fTrackMonitorFunctions[ii]->Process(track1);
1110
1111      if (fNTrackFunctions == 0x0) continue;
1112
1113      /***************************************/
1114      /******   filling numerators    ********/
1115      /****** (particles from event2) ********/
1116      /***************************************/
1117
1118      for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups 
1119       {
1120         track2= trackEvent2->GetParticle(k);
1121         if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
1122         trackpair->SetParticles(track1,track2);
1123
1124
1125         if( fPairCut->PassPairProp(trackpair)) //check pair cut
1126          { //do not meets crietria of the pair cut
1127            continue; 
1128          }
1129         else
1130          {//meets criteria of the pair cut
1131            UInt_t ii;
1132            for(ii = 0;ii<fNTrackFunctions;ii++)
1133                   fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1134          }
1135        }
1136
1137  if ( fBufferSize == 0) continue;//do not mix diff histograms
1138  /***************************************/
1139  /***** Filling denominators    *********/
1140  /***************************************/
1141  fTrackBuffer->ResetIter();
1142
1143  Int_t nmonitor = 0;
1144
1145  while ( (trackEvent3 = fTrackBuffer->Next() ) != 0x0)
1146   {
1147     if ( (j%fDisplayMixingInfo) == 0) 
1148       Info("ProcessTracksNonIdentAnal",
1149            "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1150
1151     for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
1152       {
1153         track2= trackEvent3->GetParticle(k);
1154         trackpair->SetParticles(track1,track2);
1155
1156         if( fPairCut->PassPairProp(trackpair)) //check pair cut
1157          { //do not meets crietria of the pair cut
1158            continue; 
1159          }
1160         else
1161          {//meets criteria of the pair cut
1162            for(ii = 0;ii<fNTrackFunctions;ii++)
1163                fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1164          }
1165        }// for particles event2
1166      }//while event2
1167    }//for over particles in event1
1168
1169  delete fTrackBuffer->Push(trackEvent2);
1170  
1171  return 0;
1172 }
1173 /*************************************************************************************/ 
1174
1175 void AliHBTAnalysis::Process(Option_t* option)
1176 {
1177  //default option  = "TracksAndParticles"
1178  //Main method of the HBT Analysis Package
1179  //It triggers reading with the global cut (default is an empty cut)
1180  //Than it checks options and data which are read
1181  //if everything is OK, then it calls one of the looping methods
1182  //depending on tfReaderhe option
1183  //These methods differs on what thay are looping on
1184  //
1185  //        METHOD                           OPTION
1186  //--------------------------------------------------------------------
1187  //ProcessTracksAndParticles    -     "TracksAndParticles"
1188  //     DEFAULT
1189  //     it loops over both, tracks(reconstructed) and particles(simulated)
1190  //     all function gethered in all 3 lists are called for each (double)pair
1191  //                             
1192  //ProcessTracks                -          "Tracks" 
1193  //     it loops only on tracks(reconstructed),
1194  //     functions ONLY from fTrackFunctions list are called
1195  //
1196  //ProcessParticles             -         "Particles"
1197  //     it loops only on particles(simulated),
1198  //     functions ONLY from fParticleAndTrackFunctions list are called
1199  //
1200  //
1201  if (!fReader) 
1202   {
1203    Error("Process","The reader is not set");
1204    return;
1205   }
1206  
1207  const char *oT = strstr(option,"Tracks");
1208  const char *oP = strstr(option,"Particles");
1209  
1210  Bool_t nonid = IsNonIdentAnalysis();
1211  
1212  Init();
1213  
1214  if(oT && oP)
1215   { 
1216     if (nonid) ProcessTracksAndParticlesNonIdentAnal();
1217     else ProcessTracksAndParticles();
1218     return;
1219   }
1220
1221  if(oT)
1222   {
1223     if (nonid) ProcessTracksNonIdentAnal();
1224     else ProcessTracks();
1225     return;
1226   }
1227  
1228  if(oP)
1229   {
1230     if (nonid) ProcessParticlesNonIdentAnal();
1231     else ProcessParticles();
1232     return;
1233   }
1234  
1235 }
1236 /*************************************************************************************/ 
1237
1238 void AliHBTAnalysis::ProcessTracksAndParticles()
1239 {
1240 //Makes analysis for both tracks and particles
1241 //mainly for resolution study and analysies with weighting algirithms
1242 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1243 //the loops are splited
1244   
1245 // cut on particles only -- why?
1246 // - PID: when we make resolution analysis we want to take only tracks with correct PID
1247 // We need cut on tracks because there are data characteristic to 
1248   
1249   AliAOD * trackEvent, *partEvent;
1250
1251   fReader->Rewind();
1252   while (fReader->Next() == kFALSE)
1253     {
1254       partEvent = fReader->GetEventSim();
1255       trackEvent = fReader->GetEventRec();
1256       ProcessRecAndSim(trackEvent,partEvent);
1257
1258     }//while (fReader->Next() == kFALSE)
1259     
1260
1261 /*************************************************************************************/
1262
1263 void AliHBTAnalysis::ProcessTracks()
1264 {
1265 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1266 //the loops are splited
1267   AliAOD * trackEvent;
1268   fReader->Rewind();
1269   while (fReader->Next() == kFALSE)
1270     {
1271       trackEvent = fReader->GetEventRec();
1272       ProcessRec(trackEvent,0x0);
1273     }//while (fReader->Next() == kFALSE)
1274 }
1275
1276 /*************************************************************************************/
1277
1278 void AliHBTAnalysis::ProcessParticles()
1279 {
1280 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1281 //the loops are splited
1282   AliAOD * partEvent;
1283   fReader->Rewind();
1284   while (fReader->Next() == kFALSE)
1285     {
1286       partEvent = fReader->GetEventSim();
1287       ProcessSim(0x0,partEvent);
1288     }//while (fReader->Next() == kFALSE)
1289 }
1290 /*************************************************************************************/
1291
1292 void AliHBTAnalysis::WriteFunctions()
1293 {
1294 //Calls Write for all defined functions in analysis
1295 //== writes all results
1296   TFile* oututfile = 0x0;
1297   if (fOutputFileName)
1298    {
1299      oututfile = TFile::Open(*fOutputFileName,"update");
1300    }
1301   UInt_t ii;
1302   for(ii = 0;ii<fNParticleFunctions;ii++)
1303    {
1304      if (AliVAODParticle::GetDebug()>5)
1305       {
1306         Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
1307         Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
1308       }
1309      fParticleFunctions[ii]->Write();
1310    }
1311                  
1312   for(ii = 0;ii<fNTrackFunctions;ii++)
1313    {
1314      if (AliVAODParticle::GetDebug()>5)
1315       {
1316         Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
1317         Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
1318       }
1319      fTrackFunctions[ii]->Write();
1320    }
1321                  
1322   for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1323    {
1324      if (AliVAODParticle::GetDebug()>5)
1325       {
1326         Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
1327         Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
1328       } 
1329      fParticleAndTrackFunctions[ii]->Write();
1330    }
1331
1332   for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
1333    {
1334      if (AliVAODParticle::GetDebug()>5)
1335       {
1336         Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
1337         Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
1338       }
1339      fParticleMonitorFunctions[ii]->Write();
1340    }
1341
1342   for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
1343    {
1344      if (AliVAODParticle::GetDebug()>5)
1345       {
1346         Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
1347         Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
1348       }
1349      fTrackMonitorFunctions[ii]->Write();
1350    }
1351
1352   for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
1353    {
1354      if (AliVAODParticle::GetDebug()>5)
1355       {
1356        Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
1357        Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
1358       }
1359      fParticleAndTrackMonitorFunctions[ii]->Write();
1360    }
1361   delete oututfile; 
1362 }
1363 /*************************************************************************************/
1364
1365 void AliHBTAnalysis::SetOutputFileName(const char* fname)
1366 {
1367   //Sets fiele name where to dump results, 
1368   //if not specified reults are written to gDirectory
1369   if (fname == 0x0)
1370    {
1371      delete fOutputFileName;
1372      fOutputFileName = 0x0;
1373      return;
1374    }
1375   if ( strcmp(fname,"") == 0 ) 
1376    {
1377      delete fOutputFileName;
1378      fOutputFileName = 0x0;
1379      return;
1380    }
1381   if (fOutputFileName == 0x0) fOutputFileName = new TString(fname);
1382   else *fOutputFileName = fname;
1383 }
1384 /*************************************************************************************/
1385
1386 void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
1387 {
1388 //Sets the global cut
1389   if (cut == 0x0)
1390    {
1391      Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
1392    }
1393   delete fPairCut;
1394   fPairCut = (AliAODPairCut*)cut->Clone();
1395 }
1396
1397 /*************************************************************************************/
1398
1399 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1400 {
1401 //Adds track function
1402   if (f == 0x0) return;
1403   if (fNTrackFunctions == fgkFctnArraySize)
1404    {
1405     Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
1406    }
1407   fTrackFunctions[fNTrackFunctions] = f;
1408   fNTrackFunctions++;
1409 }
1410 /*************************************************************************************/ 
1411
1412 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1413 {
1414 //adds particle function
1415   if (f == 0x0) return;
1416   
1417   if (fNParticleFunctions == fgkFctnArraySize)
1418    {
1419     Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
1420    }
1421   fParticleFunctions[fNParticleFunctions] = f;
1422   fNParticleFunctions++;
1423 }
1424 /*************************************************************************************/ 
1425
1426 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1427 {
1428 //add resolution function
1429   if (f == 0x0) return;
1430   if (fNParticleAndTrackFunctions == fgkFctnArraySize)
1431    {
1432     Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
1433    }  
1434   fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
1435   fNParticleAndTrackFunctions++;
1436 }
1437 /*************************************************************************************/ 
1438
1439 void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
1440 {
1441 //add particle monitoring function
1442   if (f == 0x0) return;
1443
1444   if (fNParticleMonitorFunctions == fgkFctnArraySize)
1445     {
1446       Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
1447    }
1448   fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
1449   fNParticleMonitorFunctions++;
1450 }
1451 /*************************************************************************************/ 
1452
1453 void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
1454 {
1455 //add track monitoring function
1456   if (f == 0x0) return;
1457
1458   if (fNTrackMonitorFunctions == fgkFctnArraySize)
1459    {
1460     Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
1461    }
1462   fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
1463   fNTrackMonitorFunctions++;
1464 }
1465 /*************************************************************************************/ 
1466
1467 void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
1468 {
1469 //add resolution monitoring function
1470   if (f == 0x0) return;
1471   if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
1472     {
1473       Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
1474     }
1475   fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
1476   fNParticleAndTrackMonitorFunctions++;
1477 }
1478
1479
1480 /*************************************************************************************/ 
1481 /*************************************************************************************/  
1482
1483 Bool_t AliHBTAnalysis::RunCoherencyCheck()
1484 {
1485  //Checks if both HBTRuns are similar
1486  //return true if error found
1487  //if they seem to be OK return false
1488 /* 
1489
1490  Int_t i;  
1491  Info("RunCoherencyCheck","Checking HBT Runs Coherency");
1492
1493 //When we use non-buffering reader this is a big waste of time -> We need to read all data to check it
1494 //and reader is implemented safe in this case anyway
1495 // Info("RunCoherencyCheck","Number of events ...");
1496 // if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same  number of events
1497 //  {
1498 //    Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
1499 //  }
1500 // else
1501 //  { //if not the same - ERROR
1502 //   Error("RunCoherencyCheck",
1503 //         "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
1504 //         fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
1505 //   return kTRUE;
1506 //  }
1507  
1508  Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
1509  
1510  AliAOD *partEvent;
1511  AliAOD *trackEvent;
1512  for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
1513   {
1514     partEvent= fReader->GetEventSim(i); //gets the "ith" event 
1515     trackEvent = fReader->GetEventRec(i);
1516     
1517     if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
1518     if ( (partEvent == 0x0) || (partEvent == 0x0) )
1519      {
1520        Error("RunCoherencyCheck",
1521              "One event is NULL and the other one not. Event Number %d",i);
1522        return kTRUE;    
1523      }
1524     
1525     if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
1526      {
1527        Error("RunCoherencyCheck",
1528              "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
1529              i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
1530        return kTRUE;
1531      }
1532     else
1533      for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
1534       {
1535         if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
1536          {
1537            Error("RunCoherencyCheck",
1538                  "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
1539                  i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
1540            return kTRUE;
1541            
1542          }
1543       }
1544   }
1545  Info("RunCoherencyCheck","  Done");
1546  Info("RunCoherencyCheck","  Everything looks OK");
1547 */ 
1548  return kFALSE;
1549 }
1550
1551 /*************************************************************************************/  
1552  
1553 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
1554 {
1555 //Performs analysis for both, tracks and particles
1556   AliAOD* rawtrackEvent, * rawpartEvent;
1557   fReader->Rewind();
1558  
1559   Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1560   Info("ProcessTracksAndParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
1561   Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1562
1563   for (Int_t i = 0;;i++)//infinite loop
1564     {
1565       if (fReader->Next()) break; //end when no more events available
1566       
1567       rawpartEvent  = fReader->GetEventSim();
1568       rawtrackEvent = fReader->GetEventRec();
1569       
1570       ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
1571   }//end of loop over events (1)
1572 }
1573 /*************************************************************************************/  
1574  
1575 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1576 {
1577 //Process Tracks only with non identical mode
1578   AliAOD * rawtrackEvent;
1579   fReader->Rewind();
1580   
1581   Info("ProcessTracksNonIdentAnal","**************************************");
1582   Info("ProcessTracksNonIdentAnal","*****  NON IDENT MODE ****************");
1583   Info("ProcessTracksNonIdentAnal","**************************************");
1584   
1585   for (Int_t i = 0;;i++)//infinite loop
1586     {
1587       if (fReader->Next()) break; //end when no more events available
1588       rawtrackEvent = fReader->GetEventRec();
1589       ProcessRecNonId(rawtrackEvent,0x0);
1590   }//end of loop over events (1)
1591 }
1592 /*************************************************************************************/  
1593
1594 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1595 {
1596 //process paricles only with non identical mode
1597   AliAOD * rawpartEvent = 0x0;
1598   fReader->Rewind();
1599   
1600   Info("ProcessParticlesNonIdentAnal","**************************************");
1601   Info("ProcessParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
1602   Info("ProcessParticlesNonIdentAnal","**************************************");
1603
1604   for (Int_t i = 0;;i++)//infinite loop
1605     {
1606       if (fReader->Next()) break; //end when no more events available
1607       
1608       rawpartEvent  = fReader->GetEventSim();
1609       ProcessSimNonId(0x0,rawpartEvent);
1610   }//end of loop over events (1)
1611 }
1612
1613 /*************************************************************************************/  
1614 void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
1615                      AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
1616 {
1617  //Puts particles accepted as a first particle by global cut in out1
1618  //and as a second particle in out2
1619
1620   AliVAODParticle* part, *track;
1621
1622   outpart1->Reset();
1623   outpart2->Reset();
1624   outtrack1->Reset();
1625   outtrack2->Reset();
1626   
1627   Bool_t in1, in2;
1628   
1629   for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1630    {
1631      in1 = in2 = kTRUE;
1632      part = inpart->GetParticle(i);
1633      track = intrack->GetParticle(i);
1634      
1635      if ( ((this->*fkPass1)(part,track))  ) in1 = kFALSE; //if part  is rejected by cut1, in1 is false
1636      if ( ((this->*fkPass2)(part,track))  ) in2 = kFALSE; //if part  is rejected by cut2, in2 is false
1637      
1638      if (gDebug)//to be removed in real analysis     
1639      if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1640       {
1641       //Particle accpted by both cuts
1642        Error("FilterOut","Particle accepted by both cuts");
1643        continue;
1644       }
1645
1646      if (in1)
1647       {
1648         outpart1->AddParticle(part);
1649         outtrack1->AddParticle(track);
1650         continue;
1651       }
1652      
1653      if (in2)
1654       {
1655         outpart2->AddParticle(new AliAODParticle(*part));
1656         outtrack2->AddParticle(new AliAODParticle(*track));
1657         continue;
1658       }
1659    }
1660 }
1661 /*************************************************************************************/  
1662
1663 void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
1664 {
1665  //Puts particles accepted as a first particle by global cut in out1
1666  //and as a second particle in out2
1667   AliVAODParticle* part;
1668   
1669   out1->Reset();
1670   out2->Reset();
1671   
1672   AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
1673   AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
1674   
1675   Bool_t in1, in2;
1676   
1677   for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1678    {
1679      in1 = in2 = kTRUE;
1680      part = in->GetParticle(i);
1681      
1682      if ( cut1->Rejected(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1683      if ( cut2->Rejected(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1684      
1685      if (gDebug)//to be removed in real analysis     
1686      if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1687       {
1688       //Particle accpted by both cuts
1689        Error("FilterOut","Particle accepted by both cuts");
1690        continue;
1691       }
1692
1693      if (in1)
1694       { 
1695         out1->AddParticle(part);
1696         continue;
1697       }
1698      
1699      if (in2)
1700       {
1701         out2->AddParticle(new AliAODParticle(*part));
1702         continue;
1703       }
1704    }
1705 }
1706 /*************************************************************************************/ 
1707
1708 Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1709 {
1710  //checks if it is possible to use special analysis for non identical particles
1711  //it means - in global pair cut first particle id is different than second one 
1712  //and both are different from 0
1713  //in the future is possible to perform more sophisticated check 
1714  //if cuts have excluding requirements
1715  
1716  if (fPairCut->IsEmpty()) 
1717    return kFALSE;
1718  
1719  if (fPairCut->GetFirstPartCut()->IsEmpty()) 
1720    return kFALSE;
1721
1722  if (fPairCut->GetSecondPartCut()->IsEmpty()) 
1723    return kFALSE;
1724  
1725  Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1726  Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1727
1728  if ( (id1==0) || (id2==0) || (id1==id2) ) 
1729    return kFALSE;
1730  
1731  return kTRUE;
1732 }
1733 /*************************************************************************************/ 
1734
1735 void AliHBTAnalysis::SetApparentVertex(Double_t x, Double_t y, Double_t z)
1736
1737   //Sets apparent vertex
1738   // All events have to be moved to the same vertex position in order to
1739   // to be able to comare any space positions (f.g. anti-merging)
1740   // This method defines this position
1741   
1742   fVertexX = x;
1743   fVertexY = y;
1744   fVertexZ = z;
1745 }
1746 /*************************************************************************************/ 
1747
1748 void AliHBTAnalysis::PressAnyKey()
1749
1750  //small utility function that helps to make comfortable macros
1751   char c;
1752   int nread = -1;
1753   fcntl(0,  F_SETFL, O_NONBLOCK);
1754   ::Info("","Press Any Key to continue ...");
1755   while (nread<1) 
1756    {
1757      nread = read(0, &c, 1);
1758      gSystem->ProcessEvents();
1759    }
1760 }
1761
1762 /*************************************************************************************/ 
1763