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