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