HBTAN/AliHBTAnalysis.cxx

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