HBTAN/AliHBTAnalysis.cxx

   1
   2 #include "AliHBTAnalysis.h"
   3
   4 #include <iostream.h>
   5
   6 #include "AliHBTRun.h"
   7 #include "AliHBTEvent.h"
   8 #include "AliHBTReader.h"
   9 #include "AliHBTParticle.h"
  10 #include "AliHBTParticleCut.h"
  11 #include "AliHBTPair.h"
  12 #include "AliHBTPairCut.h"
  13 #include "AliHBTFunction.h"
  14
  15 #include <TBenchmark.h>
  16 #include <TList.h>
  17
  18
  19
  20 ClassImp(AliHBTAnalysis)
  21
  22 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
  23 const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
  24
  25 AliHBTAnalysis::AliHBTAnalysis()
  26  {
  27    fReader = 0x0;
  28
  29    fTrackFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
  30    fParticleFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
  31    fParticleAndTrackFunctions = new AliHBTTwoPairFctn* [fgkFctnArraySize];
  32
  33    fNTrackFunctions = 0;
  34    fNParticleFunctions = 0;
  35    fNParticleAndTrackFunctions = 0;
  36
  37    fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
  38
  39    fBufferSize = 2;
  40  }
  41 /*************************************************************************************/
  42
  43 AliHBTAnalysis::~AliHBTAnalysis()
  44  {
  45  //destructor
  46  //note that we do not delete functions itself
  47  // they should be deleted by whom where created
  48  //we only store pointers, and we use "new" only for pointers array
  49    delete [] fTrackFunctions;
  50    delete [] fParticleFunctions;
  51    delete [] fParticleAndTrackFunctions;
  52
  53    delete fPairCut; // always have an copy of an object - we create we dstroy
  54  }
  55
  56 /*************************************************************************************/
  57
  58 void AliHBTAnalysis::Process(Option_t* option)
  59 {
  60  //default option  = "TracksAndParticles"
  61  //Main method of the HBT Analysis Package
  62  //It triggers reading with the global cut (default is an empty cut)
  63  //Than it checks options and data which are read
  64  //if everything is OK, then it calls one of the looping methods
  65  //depending on tfReaderhe option
  66  //These methods differs on what thay are looping on
  67  //
  68  //        METHOD                           OPTION
  69  //--------------------------------------------------------------------
  70  //ProcessTracksAndParticles    -     "TracksAndParticles"
  71  //     DEFAULT
  72  //     it loops over both, tracks(reconstructed) and particles(simulated)
  73  //     all function gethered in all 3 lists are called for each (double)pair
  74  //
  75  //ProcessTracks                -          "Tracks"
  76  //     it loops only on tracks(reconstructed),
  77  //     functions ONLY from fTrackFunctions list are called
  78  //
  79  //ProcessParticles             -         "Particles"
  80  //     it loops only on particles(simulated),
  81  //     functions ONLY from fParticleAndTrackFunctions list are called
  82  //
  83  //
  84  if (!fReader)
  85   {
  86    Error("Process","The reader is not set");
  87    return;
  88   }
  89
  90
  91  const char *oT = strstr(option,"Tracks");
  92  const char *oP = strstr(option,"Particles");
  93
  94  Bool_t nonid = IsNonIdentAnalysis();
  95
  96  if(oT && oP)
  97   {
  98     if (fReader->GetNumberOfPartEvents() <1)
  99      {
 100        Error("Process","There is no Particles. Maybe change the option?");
 101        return;
 102      }
 103     if (fReader->GetNumberOfTrackEvents() <1)
 104      {
 105        Error("Process","There is no Tracks. Maybe change the option?");
 106        return;
 107      }
 108
 109     if ( RunCoherencyCheck() )
 110       {
 111         Error("Process",
 112               "Coherency check not passed. Maybe change the option?\n");
 113         return;
 114       }
 115     if (nonid) ProcessTracksAndParticlesNonIdentAnal();
 116     else ProcessTracksAndParticles();
 117     return;
 118   }
 119
 120  if(oT)
 121   {
 122     if (fReader->GetNumberOfTrackEvents() <1)
 123      {
 124        Error("Process","There is no data to analyze.");
 125        return;
 126      }
 127     if (nonid) ProcessTracksNonIdentAnal();
 128     else ProcessTracks();
 129     return;
 130   }
 131
 132  if(oP)
 133   {
 134     if (fReader->GetNumberOfPartEvents() <1)
 135      {
 136        Error("Process","There is no data to analyze.");
 137        return;
 138      }
 139     if (nonid) ProcessParticlesNonIdentAnal();
 140     else ProcessParticles();
 141
 142 //    cout<<"NON ID"<<endl;
 143 //    ProcessParticlesNonIdentAnal();
 144 //    cout<<"\n\n\n NORMAL"<<endl;
 145 //    ProcessParticles();
 146
 147     return;
 148   }
 149
 150 }
 151
 152 /*************************************************************************************/
 153
 154 void AliHBTAnalysis::ProcessTracksAndParticles()
 155 {
 156
 157 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
 158 //the loops are splited
 159
 160
 161   AliHBTParticle * part1, * part2;
 162   AliHBTParticle * track1, * track2;
 163
 164   AliHBTEvent * trackEvent, *partEvent;
 165   AliHBTEvent * trackEvent2,*partEvent2;
 166
 167 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
 168
 169   Int_t Nev = fReader->GetNumberOfTrackEvents();
 170
 171   /***************************************/
 172   /******   Looping same events   ********/
 173   /******   filling numerators    ********/
 174   /***************************************/
 175   AliHBTPair * trackpair = new AliHBTPair();
 176   AliHBTPair * partpair = new AliHBTPair();
 177
 178   AliHBTPair * tmptrackpair;//temprary pointers to pairs
 179   AliHBTPair * tmppartpair;
 180
 181
 182
 183   for (Int_t i = 0;i<Nev;i++)
 184     {
 185       partEvent= fReader->GetParticleEvent(i);
 186       trackEvent = fReader->GetTrackEvent(i);
 187
 188       if (!partEvent) continue;
 189
 190       //N = 0;
 191
 192       for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
 193        {
 194          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
 195
 196          part1= partEvent->GetParticle(j);
 197          track1= trackEvent->GetParticle(j);
 198
 199          for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
 200           {
 201             part2= partEvent->GetParticle(k);
 202             partpair->SetParticles(part1,part2);
 203
 204             track2= trackEvent->GetParticle(k);
 205             trackpair->SetParticles(track1,track2);
 206
 207             if(fPairCut->Pass(partpair) || (fPairCut->Pass(trackpair))) //check pair cut
 208               { //do not meets crietria of the pair cut, try with swaped pairs
 209                 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
 210                   continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
 211                 else
 212                  { //swaped pair meets all the criteria
 213                    tmppartpair = partpair->GetSwapedPair();
 214                    tmptrackpair = trackpair->GetSwapedPair();
 215
 216                  }
 217               }
 218             else
 219              {//meets criteria of the pair cut
 220                tmptrackpair = trackpair;
 221                tmppartpair = partpair;
 222              }
 223             UInt_t ii;
 224             for(ii = 0;ii<fNParticleFunctions;ii++)
 225                    fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
 226
 227             for(ii = 0;ii<fNTrackFunctions;ii++)
 228                    fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
 229
 230             for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
 231                    fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
 232            }
 233        }
 234     }
 235
 236
 237   /***** Filling denominators    *********/
 238   /***************************************/
 239   for (Int_t i = 0;i<Nev-1;i++)   //In each event (but last) ....
 240     {
 241
 242       partEvent= fReader->GetParticleEvent(i);
 243       if (!partEvent) continue;
 244
 245       trackEvent = fReader->GetTrackEvent(i);
 246
 247 //      N=0;
 248
 249       for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
 250        {
 251
 252            part1= partEvent->GetParticle(j);
 253
 254            track1= trackEvent->GetParticle(j);
 255
 256            Int_t NNN;
 257
 258            if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
 259                                                                 //or current event+buffersize is greater
 260                                                                 //than max nuber of events
 261             {
 262              NNN = Nev; //set the max event number
 263             }
 264            else
 265             {
 266              NNN = i+fBufferSize; //set the current event number + fBufferSize
 267             }
 268
 269            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
 270             {
 271
 272               partEvent2= fReader->GetParticleEvent(k);
 273               if (!partEvent2) continue;
 274
 275               trackEvent2 = fReader->GetTrackEvent(k);
 276
 277              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
 278
 279              for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
 280               {
 281
 282                 // if (N>MAXCOMB) break;
 283
 284                 part2= partEvent2->GetParticle(l);
 285                 partpair->SetParticles(part1,part2);
 286
 287                 track2= trackEvent2->GetParticle(l);
 288                 trackpair->SetParticles(track1,track2);
 289
 290                 if( (fPairCut->Pass(partpair)) || (fPairCut->Pass(trackpair)) ) //check pair cut
 291                   { //do not meets crietria of the
 292                     if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
 293            continue;
 294         else
 295          {
 296            tmppartpair = partpair->GetSwapedPair();
 297            tmptrackpair = trackpair->GetSwapedPair();
 298          }
 299                   }
 300                 else
 301                  {//meets criteria of the pair cut
 302                   tmptrackpair = trackpair;
 303                   tmppartpair = partpair;
 304                  }
 305                 UInt_t ii;
 306                 for(ii = 0;ii<fNParticleFunctions;ii++)
 307                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
 308
 309                 for(ii = 0;ii<fNTrackFunctions;ii++)
 310                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
 311
 312                 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
 313                   fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
 314
 315
 316               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
 317             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
 318        }
 319
 320     }
 321
 322   /***************************************/
 323
 324
 325 }
 326 /*************************************************************************************/
 327
 328 void AliHBTAnalysis::ProcessTracks()
 329 {
 330   //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
 331 //the loops are splited
 332   AliHBTParticle * track1, * track2;
 333
 334   AliHBTEvent * trackEvent;
 335   AliHBTEvent * trackEvent2;
 336
 337 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
 338
 339   Int_t Nev = fReader->GetNumberOfTrackEvents();
 340
 341   /***************************************/
 342   /******   Looping same events   ********/
 343   /******   filling numerators    ********/
 344   /***************************************/
 345   AliHBTPair * trackpair = new AliHBTPair();
 346   AliHBTPair * tmptrackpair; //temporary pointer
 347
 348   for (Int_t i = 0;i<Nev;i++)
 349     {
 350       trackEvent = fReader->GetTrackEvent(i);
 351       if (!trackEvent) continue;
 352       //N = 0;
 353
 354       for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
 355        {
 356          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
 357
 358          track1= trackEvent->GetParticle(j);
 359
 360          for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
 361           {
 362             track2= trackEvent->GetParticle(k);
 363             trackpair->SetParticles(track1,track2);
 364             if(fPairCut->Pass(trackpair)) //check pair cut
 365               { //do not meets crietria of the
 366                 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
 367                 else tmptrackpair = trackpair->GetSwapedPair();
 368               }
 369             else
 370               {
 371                 tmptrackpair = trackpair;
 372               }
 373             UInt_t ii;
 374
 375             for(ii = 0;ii<fNTrackFunctions;ii++)
 376                    fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
 377
 378
 379            }
 380        }
 381     }
 382
 383   /***************************************/
 384   /***** Filling diff histogram *********/
 385   /***************************************/
 386   for (Int_t i = 0;i<Nev-1;i++)   //In each event (but last) ....
 387     {
 388       trackEvent = fReader->GetTrackEvent(i);
 389       if (!trackEvent) continue;
 390 //      N=0;
 391
 392       for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
 393        {
 394 //         if (N>MAXCOMB) break;
 395
 396            track1= trackEvent->GetParticle(j);
 397
 398            Int_t NNN;
 399
 400            if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
 401                                                                 //or current event+buffersize is greater
 402                                                                 //than max nuber of events
 403             {
 404              NNN = Nev; //set the max event number
 405             }
 406            else
 407             {
 408              NNN = i+fBufferSize; //set the current event number + fBufferSize
 409             }
 410
 411            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
 412             {
 413
 414              trackEvent2 = fReader->GetTrackEvent(k);
 415              if (!trackEvent2) continue;
 416
 417              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
 418
 419              for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++)   //  ... on all particles
 420               {
 421
 422                 // if (N>MAXCOMB) break;
 423                 track2= trackEvent2->GetParticle(l);
 424                 trackpair->SetParticles(track1,track2);
 425
 426                 if(fPairCut->Pass(trackpair)) //check pair cut
 427                   { //do not meets crietria of the
 428                     if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
 429         else tmptrackpair = trackpair->GetSwapedPair();
 430                   }
 431                 else
 432                   {
 433                     tmptrackpair = trackpair;
 434                   }
 435                 UInt_t ii;
 436                 for(ii = 0;ii<fNTrackFunctions;ii++)
 437                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
 438
 439               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
 440             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
 441        }
 442
 443     }
 444
 445   /***************************************/
 446
 447
 448 }
 449
 450 /*************************************************************************************/
 451
 452 void AliHBTAnalysis::ProcessParticles()
 453 {
 454   //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
 455 //the loops are splited
 456   AliHBTParticle * part1, * part2;
 457
 458   AliHBTEvent * partEvent;
 459   AliHBTEvent * partEvent2;
 460
 461   AliHBTPair * partpair = new AliHBTPair();
 462   AliHBTPair * tmppartpair; //temporary pointer to the pair
 463
 464 //  Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
 465
 466   Int_t Nev = fReader->GetNumberOfPartEvents();
 467
 468  // Nev = 1;
 469   /***************************************/
 470   /******   Looping same events   ********/
 471   /******   filling numerators    ********/
 472   /***************************************/
 473 //  gBenchmark->Start("id");
 474
 475   for (Int_t i = 0;i<Nev;i++)
 476     {
 477       partEvent= fReader->GetParticleEvent(i);
 478       if (!partEvent) continue;
 479       //N = 0;
 480
 481       for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
 482        {
 483          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
 484
 485          part1= partEvent->GetParticle(j);
 486
 487          for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
 488           {
 489             part2= partEvent->GetParticle(k);
 490             partpair->SetParticles(part1,part2);
 491
 492             if( fPairCut->Pass(partpair) ) //check pair cut
 493               { //do not meets crietria of the pair cut, try with swaped pairs
 494                 if(  fPairCut->Pass(partpair->GetSwapedPair() )  )
 495                   continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
 496                 else
 497                  { //swaped pair meets all the criteria
 498                    tmppartpair = partpair->GetSwapedPair();
 499                  }
 500               }
 501             else
 502               {
 503                 tmppartpair = partpair;
 504               }
 505
 506             UInt_t ii;
 507
 508             for(ii = 0;ii<fNParticleFunctions;ii++)
 509                    fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
 510            }
 511        }
 512     }
 513
 514   /***************************************/
 515   /***** Filling diff histogram *********/
 516   /***************************************/
 517   for (Int_t i = 0;i<Nev-1;i++)   //In each event (but last)....
 518     {
 519       partEvent= fReader->GetParticleEvent(i);
 520       if (!partEvent) continue;
 521
 522 //      N=0;
 523
 524       for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
 525        {
 526 //         if (N>MAXCOMB) break;
 527
 528            part1= partEvent->GetParticle(j);
 529
 530            Int_t NNN;
 531
 532            if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
 533                                                                 //or current event+buffersize is greater
 534                                                                 //than max nuber of events
 535             {
 536              NNN = Nev; //set the max event number
 537             }
 538            else
 539             {
 540              NNN = i+fBufferSize; //set the current event number + fBufferSize
 541             }
 542
 543 //           cout<<"NNN = "<<NNN<<endl;
 544            for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
 545             {
 546
 547              partEvent2= fReader->GetParticleEvent(k);
 548              if (!partEvent2) continue;
 549
 550              if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
 551
 552              for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++)   //  ... on all particles
 553               {
 554
 555                 // if (N>MAXCOMB) break;
 556                 part2= partEvent2->GetParticle(l);
 557                 partpair->SetParticles(part1,part2);
 558
 559                 if(fPairCut->Pass(partpair)) //check pair cut
 560                   { //do not meets crietria of the
 561                     if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
 562         else tmppartpair = partpair->GetSwapedPair();
 563                   }
 564                 else
 565                   {
 566                     tmppartpair = partpair;
 567                   }
 568                 UInt_t ii;
 569                 for(ii = 0;ii<fNParticleFunctions;ii++)
 570                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
 571
 572               }//for(Int_t l = 0; l<N2;l++)   //  ... on all particles
 573             }//for (Int_t k = i+1; k<NNN;k++)  // ... Loop over next event
 574        }
 575
 576     }
 577
 578   /***************************************/
 579
 580 //  gBenchmark->Stop("id");
 581 //  gBenchmark->Show("id");
 582
 583 }
 584
 585 /*************************************************************************************/
 586
 587
 588 void AliHBTAnalysis::WriteFunctions()
 589 {
 590   UInt_t ii;
 591   for(ii = 0;ii<fNParticleFunctions;ii++)
 592     fParticleFunctions[ii]->Write();
 593
 594   for(ii = 0;ii<fNTrackFunctions;ii++)
 595     fTrackFunctions[ii]->Write();
 596
 597   for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
 598    fParticleAndTrackFunctions[ii]->Write();
 599 }
 600 /*************************************************************************************/
 601
 602 void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
 603 {
 604   if (cut == 0x0)
 605    {
 606      Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
 607    }
 608   delete fPairCut;
 609   fPairCut = (AliHBTPairCut*)cut->Clone();
 610 }
 611
 612 /*************************************************************************************/
 613
 614 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
 615 {
 616   if (f == 0x0) return;
 617   if (fNTrackFunctions == fgkFctnArraySize)
 618    {
 619     Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
 620    }
 621   fTrackFunctions[fNTrackFunctions] = f;
 622   fNTrackFunctions++;
 623 }
 624 /*************************************************************************************/
 625
 626 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
 627 {
 628   if (f == 0x0) return;
 629
 630   if (fNParticleFunctions == fgkFctnArraySize)
 631    {
 632     Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
 633    }
 634   fParticleFunctions[fNParticleFunctions] = f;
 635   fNParticleFunctions++;
 636
 637
 638 }
 639 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
 640 {
 641   if (f == 0x0) return;
 642   if (fNParticleAndTrackFunctions == fgkFctnArraySize)
 643    {
 644     Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
 645    }
 646   fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
 647   fNParticleAndTrackFunctions++;
 648 }
 649
 650
 651 /*************************************************************************************/
 652
 653
 654 /*************************************************************************************/
 655
 656 Bool_t AliHBTAnalysis::RunCoherencyCheck()
 657 {
 658  //Checks if both HBTRuns are similar
 659  //return true if error found
 660  //if they seem to be OK return false
 661  Int_t i;
 662  cout<<"Checking HBT Runs Coherency"<<endl;
 663
 664  cout<<"Number of events ..."; fflush(stdout);
 665
 666  if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same  number of events
 667   {
 668     cout<<"OK. "<<fReader->GetNumberOfTrackEvents()<<"  found"<<endl;
 669   }
 670  else
 671   { //if not the same - ERROR
 672    Error("AliHBTAnalysis::RunCoherencyCheck()",
 673          "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
 674          fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
 675    return kTRUE;
 676   }
 677
 678  cout<<"Checking number of Particles AND Particles Types in each event ...";fflush(stdout);
 679
 680  AliHBTEvent *partEvent;
 681  AliHBTEvent *trackEvent;
 682  for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
 683   {
 684     partEvent= fReader->GetParticleEvent(i); //gets the "ith" event
 685     trackEvent = fReader->GetTrackEvent(i);
 686
 687     if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
 688     if ( (partEvent == 0x0) || (partEvent == 0x0) )
 689      {
 690        Error("AliHBTAnalysis::RunCoherencyCheck()",
 691              "One event is NULL and the other one not. Event Number %d",i);
 692        return kTRUE;
 693      }
 694
 695     if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
 696      {
 697        Error("AliHBTAnalysis::RunCoherencyCheck()",
 698              "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
 699              i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
 700        return kTRUE;
 701      }
 702     else
 703      for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
 704       {
 705         if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
 706          {
 707            Error("AliHBTAnalysis::RunCoherencyCheck()",
 708                  "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
 709                  i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
 710            return kTRUE;
 711
 712          }
 713       }
 714   }
 715  cout<<"  Done"<<endl;
 716  cout<<"  Everything looks OK"<<endl;
 717  return kFALSE;
 718 }
 719
 720 /*************************************************************************************/
 721
 722 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
 723 {
 724   AliHBTParticle * part1, * part2;
 725   AliHBTParticle * track1, * track2;
 726
 727   AliHBTEvent * trackEvent1,*partEvent1;
 728   AliHBTEvent * trackEvent2,*partEvent2;
 729   AliHBTEvent * trackEvent3,*partEvent3;
 730
 731   AliHBTEvent * rawtrackEvent, *rawpartEvent;
 732 //  AliHBTEvent * rawtrackEvent2,*rawpartEvent2;
 733
 734
 735   Int_t Nev = fReader->GetNumberOfTrackEvents();
 736
 737   AliHBTPair * trackpair = new AliHBTPair();
 738   AliHBTPair * partpair = new AliHBTPair();
 739
 740   TList tbuffer;
 741   TList pbuffer;
 742   Int_t ninbuffer = 0;
 743
 744   trackEvent1 = new AliHBTEvent();
 745   partEvent1 = new AliHBTEvent();
 746   trackEvent1->SetOwner(kFALSE);
 747   partEvent1->SetOwner(kFALSE);;
 748
 749   cout<<"**************************************"<<endl;
 750   cout<<"*****  NON IDENT MODE ****************"<<endl;
 751   cout<<"**************************************"<<endl;
 752   for (Int_t i = 0;i<Nev;i++)
 753     {
 754       rawpartEvent  = fReader->GetParticleEvent(i);
 755       rawtrackEvent = fReader->GetTrackEvent(i);
 756       if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error
 757
 758       /********************************/
 759       /*      Filtering out           */
 760       /********************************/
 761       if (ninbuffer > fBufferSize)
 762        {
 763         partEvent2  = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
 764         trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
 765         ninbuffer--;
 766        }
 767       else
 768        {
 769         partEvent2  = new AliHBTEvent();
 770         trackEvent2 = new AliHBTEvent();
 771         partEvent2->SetOwner(kFALSE);
 772         trackEvent2->SetOwner(kFALSE);
 773        }
 774       FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
 775
 776       for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
 777        {
 778          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
 779
 780          part1= partEvent1->GetParticle(j);
 781          track1= trackEvent1->GetParticle(j);
 782
 783
 784          /***************************************/
 785          /******   filling numerators    ********/
 786          /****** (particles from event2) ********/
 787          /***************************************/
 788          for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
 789           {
 790             part2= partEvent2->GetParticle(k);
 791             partpair->SetParticles(part1,part2);
 792
 793             track2= trackEvent2->GetParticle(k);
 794             trackpair->SetParticles(track1,track2);
 795
 796
 797             if( (fPairCut->PassPairProp(partpair)) || (fPairCut->PassPairProp(trackpair))) //check pair cut
 798               { //do not meets crietria of the pair cut
 799                   continue;
 800               }
 801             else
 802              {//meets criteria of the pair cut
 803               UInt_t ii;
 804               for(ii = 0;ii<fNParticleFunctions;ii++)
 805                      fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
 806
 807               for(ii = 0;ii<fNTrackFunctions;ii++)
 808                      fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
 809
 810               for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
 811                      fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
 812              }
 813            }
 814      /***************************************/
 815      /***** Filling denominators    *********/
 816      /***************************************/
 817      TIter piter(&pbuffer);
 818      TIter titer(&tbuffer);
 819      Int_t nmonitor = 0;
 820
 821      while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
 822       {
 823         trackEvent3 = (AliHBTEvent*)titer.Next();
 824         if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
 825
 826         for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
 827           {
 828             part2= partEvent3->GetParticle(k);
 829             partpair->SetParticles(part1,part2);
 830
 831             track2= trackEvent3->GetParticle(k);
 832             trackpair->SetParticles(track1,track2);
 833
 834
 835             if( (fPairCut->PassPairProp(partpair)) || (fPairCut->PassPairProp(trackpair))) //check pair cut
 836               { //do not meets crietria of the pair cut
 837                   continue;
 838               }
 839             else
 840              {//meets criteria of the pair cut
 841               UInt_t ii;
 842               for(ii = 0;ii<fNParticleFunctions;ii++)
 843                      fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
 844
 845               for(ii = 0;ii<fNTrackFunctions;ii++)
 846                      fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
 847
 848               for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
 849                      fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
 850              }
 851            }// for particles event2
 852          }//while event2
 853        }//for over particles in event1
 854
 855      pbuffer.AddFirst(partEvent2);
 856      tbuffer.AddFirst(trackEvent2);
 857      ninbuffer++;
 858
 859   }//end of loop over events (1)
 860
 861  pbuffer.SetOwner();  //to clean stored events by the way of buffer destruction
 862  tbuffer.SetOwner();
 863 }
 864 /*************************************************************************************/
 865
 866 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
 867 {
 868   AliHBTParticle * track1, * track2;
 869
 870   AliHBTEvent * trackEvent1;
 871   AliHBTEvent * trackEvent2;
 872   AliHBTEvent * trackEvent3;
 873
 874   AliHBTEvent * rawtrackEvent;
 875 //  AliHBTEvent * rawtrackEvent2,*rawpartEvent2;
 876
 877
 878   Int_t Nev = fReader->GetNumberOfTrackEvents();
 879
 880   AliHBTPair * trackpair = new AliHBTPair();
 881
 882   TList tbuffer;
 883   Int_t ninbuffer = 0;
 884
 885   trackEvent1 = new AliHBTEvent();
 886   trackEvent1->SetOwner(kFALSE);
 887
 888   cout<<"**************************************"<<endl;
 889   cout<<"*****  NON IDENT MODE ****************"<<endl;
 890   cout<<"**************************************"<<endl;
 891   for (Int_t i = 0;i<Nev;i++)
 892     {
 893       rawtrackEvent = fReader->GetTrackEvent(i);
 894       if (rawtrackEvent == 0x0)  continue;//in case of any error
 895
 896       /********************************/
 897       /*      Filtering out           */
 898       /********************************/
 899       if (ninbuffer > fBufferSize)
 900        {
 901         trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
 902         ninbuffer--;
 903        }
 904       else
 905        {
 906         trackEvent2 = new AliHBTEvent();
 907         trackEvent2->SetOwner(kFALSE);
 908        }
 909       FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
 910
 911       for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
 912        {
 913          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
 914
 915          track1= trackEvent1->GetParticle(j);
 916
 917
 918          /***************************************/
 919          /******   filling numerators    ********/
 920          /****** (particles from event2) ********/
 921          /***************************************/
 922          for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++)
 923           {
 924             track2= trackEvent2->GetParticle(k);
 925             trackpair->SetParticles(track1,track2);
 926
 927
 928             if( fPairCut->PassPairProp(trackpair)) //check pair cut
 929               { //do not meets crietria of the pair cut
 930                   continue;
 931               }
 932             else
 933              {//meets criteria of the pair cut
 934               UInt_t ii;
 935               for(ii = 0;ii<fNTrackFunctions;ii++)
 936                      fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
 937              }
 938            }
 939      /***************************************/
 940      /***** Filling denominators    *********/
 941      /***************************************/
 942      TIter titer(&tbuffer);
 943      Int_t nmonitor = 0;
 944
 945      while ( (trackEvent3 = (AliHBTEvent*)titer.Next()) != 0x0)
 946       {
 947
 948         if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
 949
 950         for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
 951           {
 952
 953             track2= trackEvent3->GetParticle(k);
 954             trackpair->SetParticles(track1,track2);
 955
 956
 957             if( fPairCut->PassPairProp(trackpair)) //check pair cut
 958               { //do not meets crietria of the pair cut
 959                   continue;
 960               }
 961             else
 962              {//meets criteria of the pair cut
 963               UInt_t ii;
 964               for(ii = 0;ii<fNTrackFunctions;ii++)
 965                      fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
 966
 967              }
 968            }// for particles event2
 969          }//while event2
 970        }//for over particles in event1
 971
 972      tbuffer.AddFirst(trackEvent2);
 973      ninbuffer++;
 974
 975   }//end of loop over events (1)
 976
 977  tbuffer.SetOwner();
 978 }
 979 /*************************************************************************************/
 980
 981 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
 982 {
 983   AliHBTParticle * part1 = 0x0, * part2 = 0x0;
 984
 985   AliHBTEvent * partEvent1 = 0x0;
 986   AliHBTEvent * partEvent2 = 0x0;
 987   AliHBTEvent * partEvent3 = 0x0;
 988
 989   AliHBTEvent * rawpartEvent = 0x0;
 990
 991   Int_t Nev = fReader->GetNumberOfPartEvents();
 992
 993   AliHBTPair * partpair = new AliHBTPair();
 994
 995   TList pbuffer;
 996   Int_t ninbuffer = 0;
 997
 998   partEvent1 = new AliHBTEvent();
 999   partEvent1->SetOwner(kFALSE);;
1000
1001   cout<<"**************************************"<<endl;
1002   cout<<"*****    PART NON IDENT MODE    ******"<<endl;
1003   cout<<"**************************************"<<endl;
1004
1005 //  gBenchmark->Start("non_id");
1006   for (Int_t i = 0;i<Nev;i++)
1007     {
1008       rawpartEvent  = fReader->GetParticleEvent(i);
1009       if ( rawpartEvent == 0x0  ) continue;//in case of any error
1010
1011       /********************************/
1012       /*      Filtering out           */
1013       /********************************/
1014       if (ninbuffer > fBufferSize)
1015        {
1016         partEvent2  = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
1017         ninbuffer--;
1018        }
1019       else
1020        {
1021         partEvent2  = new AliHBTEvent();
1022         partEvent2->SetOwner(kFALSE);
1023        }
1024       FilterOut(partEvent1, partEvent2, rawpartEvent);
1025
1026       for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
1027        {
1028          if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
1029
1030          part1= partEvent1->GetParticle(j);
1031
1032
1033          /***************************************/
1034          /******   filling numerators    ********/
1035          /****** (particles from event2) ********/
1036          /***************************************/
1037          for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
1038           {
1039             part2= partEvent2->GetParticle(k);
1040             partpair->SetParticles(part1,part2);
1041
1042             if(fPairCut->PassPairProp(partpair) ) //check pair cut
1043               { //do not meets crietria of the pair cut
1044                   continue;
1045               }
1046             else
1047              {//meets criteria of the pair cut
1048               UInt_t ii;
1049               for(ii = 0;ii<fNParticleFunctions;ii++)
1050                      fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1051
1052              }
1053            }
1054      /***************************************/
1055      /***** Filling denominators    *********/
1056      /***************************************/
1057      TIter piter(&pbuffer);
1058      Int_t nmonitor = 0;
1059
1060      while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
1061       {
1062         if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
1063
1064         for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1065           {
1066             part2= partEvent3->GetParticle(k);
1067             partpair->SetParticles(part1,part2);
1068
1069
1070             if(fPairCut->PassPairProp(partpair) ) //check pair cut
1071               { //do not meets crietria of the pair cut
1072                   continue;
1073               }
1074             else
1075              {//meets criteria of the pair cut
1076               UInt_t ii;
1077               for(ii = 0;ii<fNParticleFunctions;ii++)
1078                      fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1079
1080              }
1081            }// for particles event2
1082          }//while event2
1083        }//for over particles in event1
1084
1085      pbuffer.AddFirst(partEvent2);
1086      ninbuffer++;
1087
1088   }//end of loop over events (1)
1089
1090 // gBenchmark->Stop("non_id");
1091 // gBenchmark->Show("non_id");
1092  pbuffer.SetOwner();//to delete stered events.
1093 }
1094
1095 /*************************************************************************************/
1096 void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, AliHBTEvent* inpart,
1097                      AliHBTEvent* outtrack1, AliHBTEvent* outtrack2, AliHBTEvent* intrack)
1098 {
1099  //Puts particles accepted as a first particle by global cut in out1
1100  //and as a second particle in out2
1101
1102   AliHBTParticle* part, *track;
1103
1104   outpart1->Reset();
1105   outpart2->Reset();
1106   outtrack1->Reset();
1107   outtrack2->Reset();
1108
1109   AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1110   AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1111
1112   Bool_t in1, in2;
1113
1114   for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1115    {
1116      in1 = in2 = kTRUE;
1117      part = inpart->GetParticle(i);
1118      track = intrack->GetParticle(i);
1119
1120      if ( (cut1->Pass(part)) || ( cut1->Pass(track) ) ) in1 = kFALSE; //if any, part or track, is rejected by cut1, in1 is false
1121      if ( (cut2->Pass(part)) || ( cut2->Pass(track) ) ) in2 = kFALSE; //if any, part or track, is rejected by cut2, in2 is false
1122
1123      if (gDebug)//to be removed in real analysis
1124      if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1125       {
1126       //Particle accpted by both cuts
1127        Error("FilterOut","Particle accepted by both cuts");
1128        continue;
1129       }
1130
1131      if (in1)
1132       {
1133         outpart1->AddParticle(part);
1134         outtrack1->AddParticle(track);
1135         continue;
1136       }
1137
1138      if (in2)
1139       {
1140         outpart2->AddParticle(part);
1141         outtrack2->AddParticle(track);
1142         continue;
1143       }
1144    }
1145
1146 }
1147
1148
1149 /*************************************************************************************/
1150 void AliHBTAnalysis::FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent* in)
1151 {
1152  //Puts particles accepted as a first particle by global cut in out1
1153  //and as a second particle in out2
1154   AliHBTParticle* part;
1155
1156   out1->Reset();
1157   out2->Reset();
1158
1159   AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1160   AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1161
1162   Bool_t in1, in2;
1163
1164   for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1165    {
1166      in1 = in2 = kTRUE;
1167      part = in->GetParticle(i);
1168
1169      if ( cut1->Pass(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1170      if ( cut2->Pass(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1171
1172      if (gDebug)//to be removed in real analysis
1173      if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1174       {
1175       //Particle accpted by both cuts
1176        Error("FilterOut","Particle accepted by both cuts");
1177        continue;
1178       }
1179
1180      if (in1)
1181       {
1182         out1->AddParticle(part);
1183         continue;
1184       }
1185
1186      if (in2)
1187       {
1188         out2->AddParticle(part);
1189         continue;
1190       }
1191    }
1192 }
1193 /*************************************************************************************/
1194
1195 Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1196 {
1197  //checks if it is possible to use special analysis for non identical particles
1198  //it means - in global pair cut first particle id is different than second one
1199  //and both are different from 0
1200  //in the future is possible to perform more sophisticated check
1201  //if cuts have excluding requirements
1202
1203  if (fPairCut->IsEmpty()) return kFALSE;
1204  if (fPairCut->GetFirstPartCut()->IsEmpty()) return kFALSE;
1205  if (fPairCut->GetSecondPartCut()->IsEmpty()) return kFALSE;
1206
1207  Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1208  Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1209
1210  if ( (id1==0) || (id2==0) || (id1==id2) ) return kFALSE;
1211
1212  return kTRUE;
1213 }