HBTAN/AliHBTAnalysis.cxx

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