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