HBTAN/AliHBTAnalysis.cxx

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