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