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