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