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