Bug correction: files closed to early
[u/mrichter/AliRoot.git] / HBTAN / AliHBTAnalysis.cxx
1
2 #include "AliHBTAnalysis.h"
3
4 #include <iostream.h>
5
6 #include "AliHBTRun.h"
7 #include "AliHBTReader.h"
8 #include "AliHBTParticle.h"
9 #include "AliHBTParticleCut.h"
10 #include "AliHBTPair.h"
11 #include "AliHBTPairCut.h"
12 #include "AliHBTFunction.h"
13
14 #include <TList.h>
15
16
17
18 ClassImp(AliHBTAnalysis)
19
20 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
21 const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
22
23 AliHBTAnalysis::AliHBTAnalysis()
24  {
25    fReader = 0x0;
26    
27    fTrackFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
28    fParticleFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
29    fParticleAndTrackFunctions = new AliHBTTwoPairFctn* [fgkFctnArraySize];
30    
31    fNTrackFunctions = 0;
32    fNParticleFunctions = 0;
33    fNParticleAndTrackFunctions = 0;
34    
35    fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
36    
37  }
38 /*************************************************************************************/ 
39
40 AliHBTAnalysis::~AliHBTAnalysis()
41  {
42  //destructor
43  //note that we do not delete functions itself
44  // they should be deleted by whom where created
45  //we only store pointers, and we use "new" only for pointers array
46    delete [] fTrackFunctions;
47    delete [] fParticleFunctions;
48    delete [] fParticleAndTrackFunctions;
49    
50    delete fPairCut; // always have an copy of an object - we create we dstroy
51  }
52
53 /*************************************************************************************/ 
54
55 void AliHBTAnalysis::Process(Option_t* option)
56 {
57  //default option  = "TracksAndParticles"
58  //Main method of the HBT Analysis Package
59  //It triggers reading with the global cut (default is an empty cut)
60  //Than it checks options and data which are read
61  //if everything is OK, then it calls one of the looping methods
62  //depending on tfReaderhe option
63  //These methods differs on what thay are looping on
64  //
65  //        METHOD                           OPTION
66  //--------------------------------------------------------------------
67  //ProcessTracksAndParticles    -     "TracksAndParticles"
68  //     DEFAULT
69  //     it loops over both, tracks(reconstructed) and particles(simulated)
70  //     all function gethered in all 3 lists are called for each (double)pair
71  //                             
72  //ProcessTracks                -          "Tracks" 
73  //     it loops only on tracks(reconstructed),
74  //     functions ONLY from fTrackFunctions list are called
75  //
76  //ProcessParticles             -         "Particles"
77  //     it loops only on particles(simulated),
78  //     functions ONLY from fParticleAndTrackFunctions list are called
79  //
80  //
81  if (!fReader) 
82   {
83    Error("Process","The reader is not set");
84    return;
85   }
86  
87  
88  const char *oT = strstr(option,"Tracks");
89  const char *oP = strstr(option,"Particles");
90  
91  if(oT && oP)
92   { 
93     if (fReader->GetNumberOfPartEvents() <1)
94      {
95        Error("Process","There is no Particles. Maybe change the option?");
96        return;
97      }
98     if (fReader->GetNumberOfTrackEvents() <1)
99      {
100        Error("Process","There is no Tracks. Maybe change the option?");
101        return;
102      }
103     
104     if ( RunCoherencyCheck() )
105       {
106         Error("Process",
107               "Coherency check not passed. Maybe change the option?\n");
108         return;
109       }
110     ProcessTracksAndParticles();
111     return;
112   }
113
114  if(oT)
115   {
116     if (fReader->GetNumberOfTrackEvents() <1)
117      {
118        Error("Process","There is no data to analyze.");
119        return;
120      }
121     ProcessTracks();
122     return;
123   }
124  
125  if(oP)
126   {
127     if (fReader->GetNumberOfPartEvents() <1)
128      {
129        Error("Process","There is no data to analyze.");
130        return;
131      }
132     ProcessParticles();
133     return;
134   }
135  
136 }
137
138 /*************************************************************************************/ 
139
140 void AliHBTAnalysis::ProcessTracksAndParticles()
141 {
142
143 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
144 //the loops are splited
145   
146   
147   AliHBTParticle * part1, * part2;
148   AliHBTParticle * track1, * track2;
149   
150   AliHBTEvent * trackEvent, *partEvent;
151   AliHBTEvent * trackEvent2,*partEvent2;
152   
153 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
154   
155   Int_t Nev = fReader->GetNumberOfTrackEvents();
156   
157   /***************************************/
158   /******   Looping same events   ********/
159   /******   filling numerators    ********/
160   /***************************************/
161   AliHBTPair * trackpair = new AliHBTPair();
162   AliHBTPair * partpair = new AliHBTPair();
163
164   AliHBTPair * tmptrackpair;//temprary pointers to pairs
165   AliHBTPair * tmppartpair;
166   
167   for (Int_t i = 0;i<Nev;i++)
168     {
169       partEvent= fReader->GetParticleEvent(i);
170       trackEvent = fReader->GetTrackEvent(i);
171       
172       if (!partEvent) continue;
173       
174       //N = 0;
175       
176       for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
177        {
178          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
179
180          part1= partEvent->GetParticle(j);
181          track1= trackEvent->GetParticle(j);    
182          
183          for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
184           {
185             part2= partEvent->GetParticle(k);
186             partpair->SetParticles(part1,part2);
187            
188             track2= trackEvent->GetParticle(k); 
189             trackpair->SetParticles(track1,track2);
190
191             if(fPairCut->Pass(partpair) || (fPairCut->Pass(trackpair))) //check pair cut
192               { //do not meets crietria of the pair cut, try with swaped pairs
193                 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) ) 
194                   continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
195                 else 
196                  { //swaped pair meets all the criteria
197                    tmppartpair = partpair->GetSwapedPair();
198                    tmptrackpair = trackpair->GetSwapedPair();
199                    
200                  }
201               }
202             else
203              {//meets criteria of the pair cut
204                tmptrackpair = trackpair;
205                tmppartpair = partpair;
206              }
207             UInt_t ii;
208             for(ii = 0;ii<fNParticleFunctions;ii++)
209                    fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
210                 
211             for(ii = 0;ii<fNTrackFunctions;ii++)
212                    fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
213                  
214             for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
215                    fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
216            }
217        }
218     }
219
220   /***************************************/
221   /***** Filling denominators    *********/
222   /***************************************/
223   for (Int_t i = 0;i<Nev;i++)   //In each event ....
224     {
225       
226       partEvent= fReader->GetParticleEvent(i);
227       if (!partEvent) continue;
228       
229       trackEvent = fReader->GetTrackEvent(i); 
230       
231 //      N=0;
232       
233       for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
234        {
235 //         if (N>MAXCOMB) break;
236            
237            part1= partEvent->GetParticle(j);
238
239            track1= trackEvent->GetParticle(j);
240  
241 //         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
242            Int_t NNN;
243   
244            if ( (i+2) < Nev) NNN = i+2;
245            else NNN = Nev;
246  
247            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
248             {
249              
250               partEvent2= fReader->GetParticleEvent(k);
251               if (!partEvent2) continue;
252               
253               trackEvent2 = fReader->GetTrackEvent(k);
254              
255              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
256              
257              for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
258               {
259                
260                 // if (N>MAXCOMB) break;
261                 
262                 part2= partEvent2->GetParticle(l);
263                 partpair->SetParticles(part1,part2);
264
265                 track2= trackEvent2->GetParticle(l);
266                 trackpair->SetParticles(track1,track2);
267
268                 if( (fPairCut->Pass(partpair)) || (fPairCut->Pass(trackpair)) ) //check pair cut
269                   { //do not meets crietria of the 
270                     if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
271            continue;
272         else 
273          {
274            tmppartpair = partpair->GetSwapedPair();
275            tmptrackpair = trackpair->GetSwapedPair();
276          }
277                   }
278                 else
279                  {//meets criteria of the pair cut
280                   tmptrackpair = trackpair;
281                   tmppartpair = partpair;
282                  }
283                 UInt_t ii;
284                 for(ii = 0;ii<fNParticleFunctions;ii++)
285                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
286                  
287                 for(ii = 0;ii<fNTrackFunctions;ii++)
288                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
289                  
290                 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
291                   fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
292            
293                
294               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
295             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
296        }
297     
298     } 
299
300   /***************************************/
301   
302    
303
304 /*************************************************************************************/
305
306 void AliHBTAnalysis::ProcessTracks()
307 {
308   //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
309 //the loops are splited
310   AliHBTParticle * track1, * track2;
311   
312   AliHBTEvent * trackEvent;
313   AliHBTEvent * trackEvent2;
314   
315 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
316   
317   Int_t Nev = fReader->GetNumberOfTrackEvents();
318   
319   /***************************************/
320   /******   Looping same events   ********/
321   /******   filling numerators    ********/
322   /***************************************/
323   AliHBTPair * trackpair = new AliHBTPair();
324   AliHBTPair * tmptrackpair; //temporary pointer 
325   
326   for (Int_t i = 0;i<Nev;i++)
327     {
328       trackEvent = fReader->GetTrackEvent(i);
329       if (!trackEvent) continue;
330       //N = 0;
331       
332       for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
333        {
334          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
335
336          track1= trackEvent->GetParticle(j);    
337          
338          for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
339           {
340             track2= trackEvent->GetParticle(k); 
341             trackpair->SetParticles(track1,track2);
342             if(fPairCut->Pass(trackpair)) //check pair cut
343               { //do not meets crietria of the 
344                 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
345                 else tmptrackpair = trackpair->GetSwapedPair();
346               }
347             else
348               {
349                 tmptrackpair = trackpair;
350               }
351             UInt_t ii;
352                 
353             for(ii = 0;ii<fNTrackFunctions;ii++)
354                    fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
355                  
356            
357            }
358        }
359     }
360
361   /***************************************/
362   /***** Filling diff histogram *********/
363   /***************************************/
364   for (Int_t i = 0;i<Nev;i++)   //In each event ....
365     {
366       trackEvent = fReader->GetTrackEvent(i);
367       if (!trackEvent) continue;
368 //      N=0;
369       
370       for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
371        {
372 //         if (N>MAXCOMB) break;
373            
374            track1= trackEvent->GetParticle(j);
375  
376 //         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
377            Int_t NNN;
378   
379            if ( (i+2) < Nev) NNN = i+2;
380            else NNN = Nev;
381  
382            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
383             {
384              
385              trackEvent2 = fReader->GetTrackEvent(k);
386              if (!trackEvent2) continue;
387              
388              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
389              
390              for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++)   //  ... on all particles
391               {
392                
393                 // if (N>MAXCOMB) break;
394                 track2= trackEvent2->GetParticle(l);
395                 trackpair->SetParticles(track1,track2);
396                 
397                 if(fPairCut->Pass(trackpair)) //check pair cut
398                   { //do not meets crietria of the 
399                     if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
400         else tmptrackpair = trackpair->GetSwapedPair();
401                   }
402                 else
403                   {
404                     tmptrackpair = trackpair;
405                   }
406                 UInt_t ii;
407                 for(ii = 0;ii<fNTrackFunctions;ii++)
408                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
409                
410               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
411             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
412        }
413     
414     } 
415
416   /***************************************/
417   
418
419 }
420
421 /*************************************************************************************/
422
423 void AliHBTAnalysis::ProcessParticles()
424 {
425   //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
426 //the loops are splited
427   AliHBTParticle * part1, * part2;
428   
429   AliHBTEvent * partEvent;
430   AliHBTEvent * partEvent2;
431
432   AliHBTPair * partpair = new AliHBTPair();
433   AliHBTPair * tmppartpair; //temporary pointer to the pair
434   
435 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
436   
437   Int_t Nev = fReader->GetNumberOfPartEvents();
438   
439   /***************************************/
440   /******   Looping same events   ********/
441   /******   filling numerators    ********/
442   /***************************************/
443   for (Int_t i = 0;i<Nev;i++)
444     {
445       partEvent= fReader->GetParticleEvent(i);
446       if (!partEvent) continue;
447       //N = 0;
448       
449       for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
450        {
451          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
452
453          part1= partEvent->GetParticle(j);
454          
455          for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
456           {
457             part2= partEvent->GetParticle(k);
458             partpair->SetParticles(part1,part2);
459             
460             if( fPairCut->Pass(partpair) ) //check pair cut
461               { //do not meets crietria of the pair cut, try with swaped pairs
462                 if(  fPairCut->Pass(partpair->GetSwapedPair() )  ) 
463                   continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
464                 else 
465                  { //swaped pair meets all the criteria
466                    tmppartpair = partpair->GetSwapedPair();
467                  }
468               }
469             else
470               {
471                 tmppartpair = partpair;
472               }
473
474             UInt_t ii;
475                 
476             for(ii = 0;ii<fNParticleFunctions;ii++)
477                    fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
478            }
479        }
480     }
481
482   /***************************************/
483   /***** Filling diff histogram *********/
484   /***************************************/
485   for (Int_t i = 0;i<Nev;i++)   //In each event ....
486     {
487       partEvent= fReader->GetParticleEvent(i);
488       if (!partEvent) continue;
489
490 //      N=0;
491       
492       for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
493        {
494 //         if (N>MAXCOMB) break;
495            
496            part1= partEvent->GetParticle(j);
497  
498 //         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
499            Int_t NNN;
500   
501            if ( (i+2) < Nev) NNN = i+2; //loop over next event
502            else NNN = Nev;
503  
504            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
505             {
506              
507              partEvent2= fReader->GetParticleEvent(k);
508              if (!partEvent2) continue;
509              
510              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
511              
512              for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
513               {
514                
515                 // if (N>MAXCOMB) break;
516                 part2= partEvent2->GetParticle(l);
517                 partpair->SetParticles(part1,part2);
518                 
519                 if(fPairCut->Pass(partpair)) //check pair cut
520                   { //do not meets crietria of the 
521                     if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
522         else tmppartpair = partpair->GetSwapedPair();
523                   }
524                 else
525                   {
526                     tmppartpair = partpair;
527                   }
528                 UInt_t ii;
529                 for(ii = 0;ii<fNParticleFunctions;ii++)
530                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
531                
532               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
533             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
534        }
535     
536     } 
537
538   /***************************************/
539   
540
541 }
542
543 /*************************************************************************************/
544
545
546 void AliHBTAnalysis::WriteFunctions()
547 {
548   UInt_t ii;
549   for(ii = 0;ii<fNParticleFunctions;ii++)
550     fParticleFunctions[ii]->Write();
551                  
552   for(ii = 0;ii<fNTrackFunctions;ii++)
553     fTrackFunctions[ii]->Write();
554                  
555   for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
556    fParticleAndTrackFunctions[ii]->Write();
557 }
558 /*************************************************************************************/
559
560 void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
561 {
562   if (cut == 0x0)
563    {
564      Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
565    }
566   delete fPairCut;
567   fPairCut = (AliHBTPairCut*)cut->Clone();
568 }
569
570 /*************************************************************************************/
571
572 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
573 {
574   if (f == 0x0) return;
575   if (fNTrackFunctions == fgkFctnArraySize)
576    {
577     Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
578    }
579   fTrackFunctions[fNTrackFunctions] = f;
580   fNTrackFunctions++;
581 }
582 /*************************************************************************************/ 
583
584 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
585 {
586   if (f == 0x0) return;
587   
588   if (fNParticleFunctions == fgkFctnArraySize)
589    {
590     Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
591    }
592   fParticleFunctions[fNParticleFunctions] = f;
593   fNParticleFunctions++;
594   
595   
596 }
597 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
598 {
599   if (f == 0x0) return;
600   if (fNParticleAndTrackFunctions == fgkFctnArraySize)
601    {
602     Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
603    }  
604   fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
605   fNParticleAndTrackFunctions++;
606 }
607
608
609 /*************************************************************************************/ 
610
611
612 /*************************************************************************************/  
613
614 Bool_t AliHBTAnalysis::RunCoherencyCheck()
615 {
616  //Checks if both HBTRuns are similar
617  //return true if error found
618  //if they seem to be OK return false
619  Int_t i;  
620  cout<<"Checking HBT Runs Coherency"<<endl;
621  
622  cout<<"Number of events ..."; fflush(stdout);
623  
624  if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same  number of events
625   {
626     cout<<"OK. "<<fReader->GetNumberOfTrackEvents()<<"  found"<<endl;
627   }
628  else
629   { //if not the same - ERROR
630    Error("AliHBTAnalysis::RunCoherencyCheck()",
631          "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
632          fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
633    return kTRUE;
634   }
635  
636  cout<<"Checking number of Particles AND Particles Types in each event ...";fflush(stdout);
637  
638  AliHBTEvent *partEvent;
639  AliHBTEvent *trackEvent;
640  for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
641   {
642     partEvent= fReader->GetParticleEvent(i); //gets the "ith" event 
643     trackEvent = fReader->GetTrackEvent(i);
644     
645     if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
646     if ( (partEvent == 0x0) || (partEvent == 0x0) )
647      {
648        Error("AliHBTAnalysis::RunCoherencyCheck()",
649              "One event is NULL and the other one not. Event Number %d",i);
650        return kTRUE;    
651      }
652     
653     if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
654      {
655        Error("AliHBTAnalysis::RunCoherencyCheck()",
656              "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
657              i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
658        return kTRUE;
659      }
660     else
661      for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
662       {
663         if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
664          {
665            Error("AliHBTAnalysis::RunCoherencyCheck()",
666                  "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
667                  i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
668            return kTRUE;
669            
670          }
671       }
672   }
673  cout<<"  Done"<<endl;
674  cout<<"  Everything looks OK"<<endl;
675  return kFALSE;
676 }
677
678
679 /*************************************************************************************/  
680  
681
682 /*************************************************************************************/  
683
684