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