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