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