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