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