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