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