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