Better error reporting
[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 AliHBTTwoPartFctn* [fgkFctnArraySize];
28    fParticleFunctions = new AliHBTTwoPartFctn* [fgkFctnArraySize];
29    fParticleAndTrackFunctions = new AliHBTFourPartFctn* [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     ProcessTracks();
117     return;
118   }
119  
120  if(oP)
121   {
122     ProcessParticles();
123     return;
124   }
125  
126 }
127
128 /*************************************************************************************/ 
129
130 void AliHBTAnalysis::ProcessTracksAndParticles()
131 {
132
133 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
134 //the loops are splited
135   
136   
137   AliHBTParticle * part1, * part2;
138   AliHBTParticle * track1, * track2;
139   
140   AliHBTEvent * trackEvent, *partEvent;
141   AliHBTEvent * trackEvent2,*partEvent2;
142   
143 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
144   
145   Int_t Nev = fReader->GetNumberOfTrackEvents();
146   
147   /***************************************/
148   /******   Looping same events   ********/
149   /******   filling numerators    ********/
150   /***************************************/
151   AliHBTPair * trackpair = new AliHBTPair();
152   AliHBTPair * partpair = new AliHBTPair();
153
154   AliHBTPair * tmptrackpair;//temprary pointers to pairs
155   AliHBTPair * tmppartpair;
156   
157   for (Int_t i = 0;i<Nev;i++)
158     {
159       partEvent= fReader->GetParticleEvent(i);
160       trackEvent = fReader->GetTrackEvent(i);
161       
162       if (!partEvent) continue;
163       
164       //N = 0;
165       
166       for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
167        {
168          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
169
170          part1= partEvent->GetParticle(j);
171          track1= trackEvent->GetParticle(j);    
172          
173          for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
174           {
175             part2= partEvent->GetParticle(k);
176             partpair->SetParticles(part1,part2);
177            
178             track2= trackEvent->GetParticle(k); 
179             trackpair->SetParticles(track1,track2);
180
181             if(fPairCut->Pass(partpair) || (fPairCut->Pass(trackpair))) //check pair cut
182               { //do not meets crietria of the pair cut, try with swaped pairs
183                 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) ) 
184                   continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
185                 else 
186                  { //swaped pair meets all the criteria
187                    tmppartpair = partpair->GetSwapedPair();
188                    tmptrackpair = trackpair->GetSwapedPair();
189                    
190                  }
191               }
192             else
193              {//meets criteria of the pair cut
194                tmptrackpair = trackpair;
195                tmppartpair = partpair;
196              }
197             UInt_t ii;
198             for(ii = 0;ii<fNParticleFunctions;ii++)
199                    fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
200                 
201             for(ii = 0;ii<fNTrackFunctions;ii++)
202                    fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
203                  
204             for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
205                    fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
206            }
207        }
208     }
209
210   /***************************************/
211   /***** Filling denominators    *********/
212   /***************************************/
213   for (Int_t i = 0;i<Nev;i++)   //In each event ....
214     {
215       
216       partEvent= fReader->GetParticleEvent(i);
217       if (!partEvent) continue;
218       
219       trackEvent = fReader->GetTrackEvent(i); 
220       
221 //      N=0;
222       
223       for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
224        {
225 //         if (N>MAXCOMB) break;
226            
227            part1= partEvent->GetParticle(j);
228
229            track1= trackEvent->GetParticle(j);
230  
231 //         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
232            Int_t NNN;
233   
234            if ( (i+2) < Nev) NNN = i+2;
235            else NNN = Nev;
236  
237            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
238             {
239              
240               partEvent2= fReader->GetParticleEvent(k);
241               if (!partEvent2) continue;
242               
243               trackEvent2 = fReader->GetTrackEvent(k);
244              
245              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
246              
247              for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
248               {
249                
250                 // if (N>MAXCOMB) break;
251                 
252                 part2= partEvent2->GetParticle(l);
253                 partpair->SetParticles(part1,part2);
254
255                 track2= trackEvent2->GetParticle(l);
256                 trackpair->SetParticles(track1,track2);
257
258                 if( (fPairCut->Pass(partpair)) || (fPairCut->Pass(trackpair)) ) //check pair cut
259                   { //do not meets crietria of the 
260                     if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
261            continue;
262         else 
263          {
264            tmppartpair = partpair->GetSwapedPair();
265            tmptrackpair = trackpair->GetSwapedPair();
266          }
267                   }
268                 else
269                  {//meets criteria of the pair cut
270                   tmptrackpair = trackpair;
271                   tmppartpair = partpair;
272                  }
273                 UInt_t ii;
274                 for(ii = 0;ii<fNParticleFunctions;ii++)
275                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
276                  
277                 for(ii = 0;ii<fNTrackFunctions;ii++)
278                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
279                  
280                 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
281                   fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
282            
283                
284               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
285             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
286        }
287     
288     } 
289
290   /***************************************/
291   
292    
293
294 /*************************************************************************************/
295
296 void AliHBTAnalysis::ProcessTracks()
297 {
298   //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
299 //the loops are splited
300   AliHBTParticle * track1, * track2;
301   
302   AliHBTEvent * trackEvent;
303   AliHBTEvent * trackEvent2;
304   
305 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
306   
307   Int_t Nev = fReader->GetNumberOfTrackEvents();
308   
309   /***************************************/
310   /******   Looping same events   ********/
311   /******   filling numerators    ********/
312   /***************************************/
313   AliHBTPair * trackpair = new AliHBTPair();
314   AliHBTPair * tmptrackpair; //temporary pointer 
315   
316   for (Int_t i = 0;i<Nev;i++)
317     {
318       trackEvent = fReader->GetTrackEvent(i);
319       if (!trackEvent) continue;
320       //N = 0;
321       
322       for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
323        {
324          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
325
326          track1= trackEvent->GetParticle(j);    
327          
328          for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
329           {
330             track2= trackEvent->GetParticle(k); 
331             trackpair->SetParticles(track1,track2);
332             if(fPairCut->Pass(trackpair)) //check pair cut
333               { //do not meets crietria of the 
334                 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
335                 else tmptrackpair = trackpair->GetSwapedPair();
336               }
337             else
338               {
339                 tmptrackpair = trackpair;
340               }
341             UInt_t ii;
342                 
343             for(ii = 0;ii<fNTrackFunctions;ii++)
344                    fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
345                  
346            
347            }
348        }
349     }
350
351   /***************************************/
352   /***** Filling diff histogram *********/
353   /***************************************/
354   for (Int_t i = 0;i<Nev;i++)   //In each event ....
355     {
356       trackEvent = fReader->GetTrackEvent(i);
357       if (!trackEvent) continue;
358 //      N=0;
359       
360       for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
361        {
362 //         if (N>MAXCOMB) break;
363            
364            track1= trackEvent->GetParticle(j);
365  
366 //         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
367            Int_t NNN;
368   
369            if ( (i+2) < Nev) NNN = i+2;
370            else NNN = Nev;
371  
372            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
373             {
374              
375              trackEvent2 = fReader->GetTrackEvent(k);
376              if (!trackEvent2) continue;
377              
378              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
379              
380              for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++)   //  ... on all particles
381               {
382                
383                 // if (N>MAXCOMB) break;
384                 track2= trackEvent2->GetParticle(l);
385                 trackpair->SetParticles(track1,track2);
386                 
387                 if(fPairCut->Pass(trackpair)) //check pair cut
388                   { //do not meets crietria of the 
389                     if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
390         else tmptrackpair = trackpair->GetSwapedPair();
391                   }
392                 else
393                   {
394                     tmptrackpair = trackpair;
395                   }
396                 UInt_t ii;
397                 for(ii = 0;ii<fNTrackFunctions;ii++)
398                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
399                
400               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
401             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
402        }
403     
404     } 
405
406   /***************************************/
407   
408
409 }
410
411 /*************************************************************************************/
412
413 void AliHBTAnalysis::ProcessParticles()
414 {
415   //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
416 //the loops are splited
417   AliHBTParticle * part1, * part2;
418   
419   AliHBTEvent * partEvent;
420   AliHBTEvent * partEvent2;
421
422   AliHBTPair * partpair = new AliHBTPair();
423   AliHBTPair * tmppartpair; //temporary pointer to the pair
424   
425 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
426   
427   Int_t Nev = fReader->GetNumberOfPartEvents();
428   
429   /***************************************/
430   /******   Looping same events   ********/
431   /******   filling numerators    ********/
432   /***************************************/
433   for (Int_t i = 0;i<Nev;i++)
434     {
435       partEvent= fReader->GetParticleEvent(i);
436       if (!partEvent) continue;
437       //N = 0;
438       
439       for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
440        {
441          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
442
443          part1= partEvent->GetParticle(j);
444          
445          for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
446           {
447             part2= partEvent->GetParticle(k);
448             partpair->SetParticles(part1,part2);
449             
450             if( fPairCut->Pass(partpair) ) //check pair cut
451               { //do not meets crietria of the pair cut, try with swaped pairs
452                 if(  fPairCut->Pass(partpair->GetSwapedPair() )  ) 
453                   continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
454                 else 
455                  { //swaped pair meets all the criteria
456                    tmppartpair = partpair->GetSwapedPair();
457                  }
458               }
459             else
460               {
461                 tmppartpair = partpair;
462               }
463
464             UInt_t ii;
465                 
466             for(ii = 0;ii<fNParticleFunctions;ii++)
467                    fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
468            }
469        }
470     }
471
472   /***************************************/
473   /***** Filling diff histogram *********/
474   /***************************************/
475   for (Int_t i = 0;i<Nev;i++)   //In each event ....
476     {
477       partEvent= fReader->GetParticleEvent(i);
478       if (!partEvent) continue;
479
480 //      N=0;
481       
482       for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
483        {
484 //         if (N>MAXCOMB) break;
485            
486            part1= partEvent->GetParticle(j);
487  
488 //         for (Int_t k = i+1; k<Nev;k++)  //  ... Loop over all proceeding events ...
489            Int_t NNN;
490   
491            if ( (i+2) < Nev) NNN = i+2; //loop over next event
492            else NNN = Nev;
493  
494            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
495             {
496              
497              partEvent2= fReader->GetParticleEvent(k);
498              if (!partEvent2) continue;
499              
500              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
501              
502              for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
503               {
504                
505                 // if (N>MAXCOMB) break;
506                 part2= partEvent2->GetParticle(l);
507                 partpair->SetParticles(part1,part2);
508                 
509                 if(fPairCut->Pass(partpair)) //check pair cut
510                   { //do not meets crietria of the 
511                     if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
512         else tmppartpair = partpair->GetSwapedPair();
513                   }
514                 else
515                   {
516                     tmppartpair = partpair;
517                   }
518                 UInt_t ii;
519                 for(ii = 0;ii<fNParticleFunctions;ii++)
520                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
521                
522               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
523             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
524        }
525     
526     } 
527
528   /***************************************/
529   
530
531 }
532
533 /*************************************************************************************/
534
535
536 void AliHBTAnalysis::WriteFunctions()
537 {
538   UInt_t ii;
539   for(ii = 0;ii<fNParticleFunctions;ii++)
540     fParticleFunctions[ii]->Write();
541                  
542   for(ii = 0;ii<fNTrackFunctions;ii++)
543     fTrackFunctions[ii]->Write();
544                  
545   for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
546    fParticleAndTrackFunctions[ii]->Write();
547 }
548 /*************************************************************************************/
549
550 void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
551 {
552   if (cut == 0x0)
553    {
554      Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
555    }
556   delete fPairCut;
557   fPairCut = (AliHBTPairCut*)cut->Clone();
558 }
559
560 /*************************************************************************************/
561
562 void AliHBTAnalysis::AddTrackFunction(AliHBTTwoPartFctn* f)
563 {
564   if (f == 0x0) return;
565   if (fNTrackFunctions == fgkFctnArraySize)
566    {
567     Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
568    }
569   fTrackFunctions[fNTrackFunctions] = f;
570   fNTrackFunctions++;
571 }
572 /*************************************************************************************/ 
573
574 void AliHBTAnalysis::AddParticleFunction(AliHBTTwoPartFctn* f)
575 {
576   if (f == 0x0) return;
577   
578   if (fNParticleFunctions == fgkFctnArraySize)
579    {
580     Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
581    }
582   fParticleFunctions[fNParticleFunctions] = f;
583   fNParticleFunctions++;
584   
585   
586 }
587 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTFourPartFctn* f)
588 {
589   if (f == 0x0) return;
590   if (fNParticleAndTrackFunctions == fgkFctnArraySize)
591    {
592     Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
593    }  
594   fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
595   fNParticleAndTrackFunctions++;
596 }
597
598
599 /*************************************************************************************/ 
600
601
602 /*************************************************************************************/  
603
604 Bool_t AliHBTAnalysis::RunCoherencyCheck()
605 {
606  //Checks if both HBTRuns are similar
607  //return true if error found
608  //if they seem to be OK return false
609  Int_t i;  
610  cout<<"Checking HBT Runs Coherency"<<endl;
611  
612  cout<<"Number of events ..."; fflush(stdout);
613  
614  if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same  number of events
615   {
616     cout<<"OK. "<<fReader->GetNumberOfTrackEvents()<<"  found"<<endl;
617   }
618  else
619   { //if not the same - ERROR
620    Error("AliHBTAnalysis::RunCoherencyCheck()",
621          "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
622          fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
623    return kTRUE;
624   }
625  
626  cout<<"Checking number of Particles AND Particles Types in each event ...";fflush(stdout);
627  
628  AliHBTEvent *partEvent;
629  AliHBTEvent *trackEvent;
630  for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
631   {
632     partEvent= fReader->GetParticleEvent(i); //gets the "ith" event 
633     trackEvent = fReader->GetTrackEvent(i);
634     
635     if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
636     if ( (partEvent == 0x0) || (partEvent == 0x0) )
637      {
638        Error("AliHBTAnalysis::RunCoherencyCheck()",
639              "One event is NULL and the other one not. Event Number %d",i);
640        return kTRUE;    
641      }
642     
643     if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
644      {
645        Error("AliHBTAnalysis::RunCoherencyCheck()",
646              "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
647              i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
648        return kTRUE;
649      }
650     else
651      for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
652       {
653         if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
654          {
655            Error("AliHBTAnalysis::RunCoherencyCheck()",
656                  "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
657                  i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
658            return kTRUE;
659            
660          }
661       }
662   }
663  cout<<"  Done"<<endl;
664  cout<<"  Everything looks OK"<<endl;
665  return kFALSE;
666 }
667
668
669 /*************************************************************************************/  
670  
671
672 /*************************************************************************************/  
673
674