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