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