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