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 <TSystem.h>
  24 #include <TFile.h>
  25
  26 #include "AliAOD.h"
  27 #include "AliAODParticle.h"
  28 #include "AliAODPairCut.h"
  29 #include "AliEventCut.h"
  30
  31 #include "AliEventBuffer.h"
  32
  33 #include "AliReader.h"
  34 #include "AliHBTPair.h"
  35 #include "AliHBTFunction.h"
  36 #include "AliHBTMonitorFunction.h"
  37
  38
  39 ClassImp(AliHBTAnalysis)
  40
  41 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
  42 const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
  43 const Int_t  AliHBTAnalysis::fgkDefaultBufferSize = 5;
  44
  45 AliHBTAnalysis::AliHBTAnalysis():
  46   fProcEvent(0x0),
  47   fReader(0x0),
  48   fNTrackFunctions(0),
  49   fNParticleFunctions(0),
  50   fNParticleAndTrackFunctions(0),
  51   fNTrackMonitorFunctions(0),
  52   fNParticleMonitorFunctions(0),
  53   fNParticleAndTrackMonitorFunctions(0),
  54   fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
  55   fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
  56   fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
  57   fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
  58   fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
  59   fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),
  60   fBkgEventCut(0x0),
  61   fBufferSize(2),
  62   fDisplayMixingInfo(fgkDefaultMixingInfo),
  63   fIsOwner(kFALSE),
  64   fPartBuffer(0x0),
  65   fTrackBuffer(0x0),
  66   fProcessOption(kSimulatedAndReconstructed),
  67   fNoCorrfctns(kFALSE),
  68   fOutputFileName(0x0)
  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   fBkgEventCut(0x0),
  92   fBufferSize(fgkDefaultBufferSize),
  93   fDisplayMixingInfo(fgkDefaultMixingInfo),
  94   fIsOwner(kFALSE),
  95   fPartBuffer(0x0),
  96   fTrackBuffer(0x0),
  97   fProcessOption(kSimulatedAndReconstructed),
  98   fNoCorrfctns(kFALSE),
  99   fOutputFileName(0x0)
 100  {
 101 //copy constructor
 102    Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
 103  }
 104 /*************************************************************************************/
 105 AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
 106  {
 107 //operator =
 108    Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
 109    return *this;
 110  }
 111 /*************************************************************************************/
 112 AliHBTAnalysis::~AliHBTAnalysis()
 113  {
 114  //destructor
 115  //note that we do not delete functions itself
 116  // they should be deleted by whom where created
 117  //we only store pointers, and we use "new" only for pointers array
 118
 119    if (fIsOwner)
 120     {
 121       if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Is Owner: Attempting to delete functions");
 122       DeleteFunctions();
 123       if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
 124     }
 125    delete [] fTrackFunctions;
 126    delete [] fParticleFunctions;
 127    delete [] fParticleAndTrackFunctions;
 128
 129    delete [] fParticleMonitorFunctions;
 130    delete [] fTrackMonitorFunctions;
 131    delete [] fParticleAndTrackMonitorFunctions;
 132
 133    delete fBkgEventCut;
 134    delete fOutputFileName;
 135  }
 136
 137 /*************************************************************************************/
 138
 139 Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
 140 {
 141   //Processes one event
 142   if (fProcEvent == 0x0)
 143    {
 144      Error("ProcessEvent","Analysis option not specified");
 145      return 1;
 146    }
 147   if ( Pass(aodrec,aodsim) ) return 0;
 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   TFile* oututfile = 0x0;
1261   if (fOutputFileName)
1262    {
1263      oututfile = TFile::Open(*fOutputFileName,"update");
1264    }
1265   UInt_t ii;
1266   for(ii = 0;ii<fNParticleFunctions;ii++)
1267    {
1268      if (AliVAODParticle::GetDebug()>5)
1269       {
1270         Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
1271         Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
1272       }
1273      fParticleFunctions[ii]->Write();
1274    }
1275
1276   for(ii = 0;ii<fNTrackFunctions;ii++)
1277    {
1278      if (AliVAODParticle::GetDebug()>5)
1279       {
1280         Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
1281         Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
1282       }
1283      fTrackFunctions[ii]->Write();
1284    }
1285
1286   for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1287    {
1288      if (AliVAODParticle::GetDebug()>5)
1289       {
1290         Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
1291         Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
1292       }
1293      fParticleAndTrackFunctions[ii]->Write();
1294    }
1295
1296   for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
1297    {
1298      if (AliVAODParticle::GetDebug()>5)
1299       {
1300         Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
1301         Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
1302       }
1303      fParticleMonitorFunctions[ii]->Write();
1304    }
1305
1306   for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
1307    {
1308      if (AliVAODParticle::GetDebug()>5)
1309       {
1310         Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
1311         Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
1312       }
1313      fTrackMonitorFunctions[ii]->Write();
1314    }
1315
1316   for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
1317    {
1318      if (AliVAODParticle::GetDebug()>5)
1319       {
1320        Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
1321        Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
1322       }
1323      fParticleAndTrackMonitorFunctions[ii]->Write();
1324    }
1325   delete oututfile;
1326 }
1327 /*************************************************************************************/
1328
1329 void AliHBTAnalysis::SetOutputFileName(const char* fname)
1330 {
1331   //Sets fiele name where to dump results,
1332   //if not specified reults are written to gDirectory
1333   if (fname == 0x0)
1334    {
1335      delete fOutputFileName;
1336      fOutputFileName = 0x0;
1337      return;
1338    }
1339   if ( strcmp(fname,"") == 0 )
1340    {
1341      delete fOutputFileName;
1342      fOutputFileName = 0x0;
1343      return;
1344    }
1345   if (fOutputFileName == 0x0) fOutputFileName = new TString(fname);
1346   else *fOutputFileName = fname;
1347 }
1348 /*************************************************************************************/
1349
1350 void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
1351 {
1352 //Sets the global cut
1353   if (cut == 0x0)
1354    {
1355      Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
1356    }
1357   delete fPairCut;
1358   fPairCut = (AliAODPairCut*)cut->Clone();
1359 }
1360
1361 /*************************************************************************************/
1362
1363 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1364 {
1365 //Adds track function
1366   if (f == 0x0) return;
1367   if (fNTrackFunctions == fgkFctnArraySize)
1368    {
1369     Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
1370    }
1371   fTrackFunctions[fNTrackFunctions] = f;
1372   fNTrackFunctions++;
1373 }
1374 /*************************************************************************************/
1375
1376 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1377 {
1378 //adds particle function
1379   if (f == 0x0) return;
1380
1381   if (fNParticleFunctions == fgkFctnArraySize)
1382    {
1383     Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
1384    }
1385   fParticleFunctions[fNParticleFunctions] = f;
1386   fNParticleFunctions++;
1387 }
1388 /*************************************************************************************/
1389
1390 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1391 {
1392 //add resolution function
1393   if (f == 0x0) return;
1394   if (fNParticleAndTrackFunctions == fgkFctnArraySize)
1395    {
1396     Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
1397    }
1398   fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
1399   fNParticleAndTrackFunctions++;
1400 }
1401 /*************************************************************************************/
1402
1403 void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
1404 {
1405 //add particle monitoring function
1406   if (f == 0x0) return;
1407
1408   if (fNParticleMonitorFunctions == fgkFctnArraySize)
1409     {
1410       Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
1411    }
1412   fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
1413   fNParticleMonitorFunctions++;
1414 }
1415 /*************************************************************************************/
1416
1417 void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
1418 {
1419 //add track monitoring function
1420   if (f == 0x0) return;
1421
1422   if (fNTrackMonitorFunctions == fgkFctnArraySize)
1423    {
1424     Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
1425    }
1426   fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
1427   fNTrackMonitorFunctions++;
1428 }
1429 /*************************************************************************************/
1430
1431 void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
1432 {
1433 //add resolution monitoring function
1434   if (f == 0x0) return;
1435   if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
1436     {
1437       Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
1438     }
1439   fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
1440   fNParticleAndTrackMonitorFunctions++;
1441 }
1442
1443
1444 /*************************************************************************************/
1445 /*************************************************************************************/
1446
1447 Bool_t AliHBTAnalysis::RunCoherencyCheck()
1448 {
1449  //Checks if both HBTRuns are similar
1450  //return true if error found
1451  //if they seem to be OK return false
1452 /*
1453
1454  Int_t i;
1455  Info("RunCoherencyCheck","Checking HBT Runs Coherency");
1456
1457 //When we use non-buffering reader this is a big waste of time -> We need to read all data to check it
1458 //and reader is implemented safe in this case anyway
1459 // Info("RunCoherencyCheck","Number of events ...");
1460 // if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same  number of events
1461 //  {
1462 //    Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
1463 //  }
1464 // else
1465 //  { //if not the same - ERROR
1466 //   Error("RunCoherencyCheck",
1467 //         "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
1468 //         fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
1469 //   return kTRUE;
1470 //  }
1471
1472  Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
1473
1474  AliAOD *partEvent;
1475  AliAOD *trackEvent;
1476  for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
1477   {
1478     partEvent= fReader->GetEventSim(i); //gets the "ith" event
1479     trackEvent = fReader->GetEventRec(i);
1480
1481     if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
1482     if ( (partEvent == 0x0) || (partEvent == 0x0) )
1483      {
1484        Error("RunCoherencyCheck",
1485              "One event is NULL and the other one not. Event Number %d",i);
1486        return kTRUE;
1487      }
1488
1489     if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
1490      {
1491        Error("RunCoherencyCheck",
1492              "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
1493              i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
1494        return kTRUE;
1495      }
1496     else
1497      for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
1498       {
1499         if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
1500          {
1501            Error("RunCoherencyCheck",
1502                  "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
1503                  i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
1504            return kTRUE;
1505
1506          }
1507       }
1508   }
1509  Info("RunCoherencyCheck","  Done");
1510  Info("RunCoherencyCheck","  Everything looks OK");
1511 */
1512  return kFALSE;
1513 }
1514
1515 /*************************************************************************************/
1516
1517 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
1518 {
1519 //Performs analysis for both, tracks and particles
1520   AliAOD* rawtrackEvent, * rawpartEvent;
1521   fReader->Rewind();
1522
1523   Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1524   Info("ProcessTracksAndParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
1525   Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1526
1527   for (Int_t i = 0;;i++)//infinite loop
1528     {
1529       if (fReader->Next()) break; //end when no more events available
1530
1531       rawpartEvent  = fReader->GetEventSim();
1532       rawtrackEvent = fReader->GetEventRec();
1533
1534       ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
1535   }//end of loop over events (1)
1536 }
1537 /*************************************************************************************/
1538
1539 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1540 {
1541 //Process Tracks only with non identical mode
1542   AliAOD * rawtrackEvent;
1543   fReader->Rewind();
1544
1545   Info("ProcessTracksNonIdentAnal","**************************************");
1546   Info("ProcessTracksNonIdentAnal","*****  NON IDENT MODE ****************");
1547   Info("ProcessTracksNonIdentAnal","**************************************");
1548
1549   for (Int_t i = 0;;i++)//infinite loop
1550     {
1551       if (fReader->Next()) break; //end when no more events available
1552       rawtrackEvent = fReader->GetEventRec();
1553       ProcessRecNonId(rawtrackEvent,0x0);
1554   }//end of loop over events (1)
1555 }
1556 /*************************************************************************************/
1557
1558 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1559 {
1560 //process paricles only with non identical mode
1561   AliAOD * rawpartEvent = 0x0;
1562   fReader->Rewind();
1563
1564   Info("ProcessParticlesNonIdentAnal","**************************************");
1565   Info("ProcessParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
1566   Info("ProcessParticlesNonIdentAnal","**************************************");
1567
1568   for (Int_t i = 0;;i++)//infinite loop
1569     {
1570       if (fReader->Next()) break; //end when no more events available
1571
1572       rawpartEvent  = fReader->GetEventSim();
1573       ProcessSimNonId(0x0,rawpartEvent);
1574   }//end of loop over events (1)
1575 }
1576
1577 /*************************************************************************************/
1578 void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
1579                      AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
1580 {
1581  //Puts particles accepted as a first particle by global cut in out1
1582  //and as a second particle in out2
1583
1584   AliVAODParticle* part, *track;
1585
1586   outpart1->Reset();
1587   outpart2->Reset();
1588   outtrack1->Reset();
1589   outtrack2->Reset();
1590
1591   Bool_t in1, in2;
1592
1593   for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1594    {
1595      in1 = in2 = kTRUE;
1596      part = inpart->GetParticle(i);
1597      track = intrack->GetParticle(i);
1598
1599      if ( ((this->*fkPass1)(part,track))  ) in1 = kFALSE; //if part  is rejected by cut1, in1 is false
1600      if ( ((this->*fkPass2)(part,track))  ) in2 = kFALSE; //if part  is rejected by cut2, in2 is false
1601
1602      if (gDebug)//to be removed in real analysis
1603      if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1604       {
1605       //Particle accpted by both cuts
1606        Error("FilterOut","Particle accepted by both cuts");
1607        continue;
1608       }
1609
1610      if (in1)
1611       {
1612         outpart1->AddParticle(part);
1613         outtrack1->AddParticle(track);
1614         continue;
1615       }
1616
1617      if (in2)
1618       {
1619         outpart2->AddParticle(new AliAODParticle(*part));
1620         outtrack2->AddParticle(new AliAODParticle(*track));
1621         continue;
1622       }
1623    }
1624 }
1625 /*************************************************************************************/
1626
1627 void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
1628 {
1629  //Puts particles accepted as a first particle by global cut in out1
1630  //and as a second particle in out2
1631   AliVAODParticle* part;
1632
1633   out1->Reset();
1634   out2->Reset();
1635
1636   AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
1637   AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
1638
1639   Bool_t in1, in2;
1640
1641   for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1642    {
1643      in1 = in2 = kTRUE;
1644      part = in->GetParticle(i);
1645
1646      if ( cut1->Pass(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1647      if ( cut2->Pass(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1648
1649      if (gDebug)//to be removed in real analysis
1650      if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1651       {
1652       //Particle accpted by both cuts
1653        Error("FilterOut","Particle accepted by both cuts");
1654        continue;
1655       }
1656
1657      if (in1)
1658       {
1659         out1->AddParticle(part);
1660         continue;
1661       }
1662
1663      if (in2)
1664       {
1665         out2->AddParticle(part);
1666         continue;
1667       }
1668    }
1669 }
1670 /*************************************************************************************/
1671
1672 Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1673 {
1674  //checks if it is possible to use special analysis for non identical particles
1675  //it means - in global pair cut first particle id is different than second one
1676  //and both are different from 0
1677  //in the future is possible to perform more sophisticated check
1678  //if cuts have excluding requirements
1679
1680  if (fPairCut->IsEmpty())
1681    return kFALSE;
1682
1683  if (fPairCut->GetFirstPartCut()->IsEmpty())
1684    return kFALSE;
1685
1686  if (fPairCut->GetSecondPartCut()->IsEmpty())
1687    return kFALSE;
1688
1689  Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1690  Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1691
1692  if ( (id1==0) || (id2==0) || (id1==id2) )
1693    return kFALSE;
1694
1695  return kTRUE;
1696 }
1697 /*************************************************************************************/
1698
1699 void AliHBTAnalysis::PressAnyKey()
1700 {
1701  //small utility function that helps to make comfortable macros
1702   char c;
1703   int nread = -1;
1704   fcntl(0,  F_SETFL, O_NONBLOCK);
1705   ::Info("","Press Any Key to continue ...");
1706   while (nread<1)
1707    {
1708      nread = read(0, &c, 1);
1709      gSystem->ProcessEvents();
1710    }
1711 }
1712
1713 /*************************************************************************************/
1714