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