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