HBTAN/AliHBTAnalysis.cxx

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