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