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