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