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