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