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