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