HBTAN/AliHBTAnalysis.cxx

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