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