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