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