HBTAN/AliHBTAnalysis.cxx

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