EVGEN/AliGenHBTosl.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 #include "AliGenHBTosl.h"
  19 #include "AliLog.h"
  20
  21 //__________________________________________________________
  22 /////////////////////////////////////////////////////////////
  23 //                                                         //
  24 //  class AliGenHBTosl                                     //
  25 //                                                         //
  26 //  Genarator simulating particle correlations             //
  27 //                                                         //
  28 //  The main idea of the generator is to produce particles //
  29 //  according to some distribution of two particle         //
  30 //  property. In HBT they are qout,qsie and qlong.         //
  31 //  In order to be able to generate signal that produces   //
  32 //  given two particle correlation background must be      //
  33 //  known before in order to produce the shape of signal   //
  34 //  to randomize given distribution from.                  //
  35 //                                                         //
  36 //  The generator works as follows:                        //
  37 //  1. Coarse Background (fQCoarseBackground) is generated //
  38 //     ade  from the particles                             //
  39 //     given by the external generator (variable           //
  40 //     fGenerator) by the mixing technique.                //
  41 //  2. Coarse signal is prduced by multiplying Coarse      //
  42 //     background by a required function                   //
  43 //     See method FillCoarseSignal                         //
  44 //  3. Signal is randomized out of the coarse signal       //
  45 //     histogram (two particle property). First particle   //
  46 //     is taken from the external generator, and the       //
  47 //     second one is CALCULATED on the basis of the first  //
  48 //     one and the two particle property (qout,qside,qlong)//
  49 //     Background is made by the mixing out of the         //
  50 //     genereted signal events.                            //
  51 //     This step is cotinued up to the moment signal       //
  52 //     histogram has enough statistics (data member        //
  53 //     fMinFill)                                           //
  54 //     See method StartSignalPass1()                       //
  55 //  4. chi is calculated for each bin (chiarray variqable) //
  56 //     (not the chi2 because sign is important)            //
  57 //     Two particle prioperty                              //
  58 //     (qout,qside,qlong) is chosen at the points that     //
  59 //     chi is the smallest. First particle is taken from   //
  60 //     the the external generator (fGenerator) and second's /
  61 //     momenta are caclulated out of their momenta and     //
  62 //     (qout,qside,qlong). Background is updated           //
  63 //     continuesely for all the events. This step is       //
  64 //     continued until stability conditions are fullfiled  //
  65 //     or maximum number of iteration is reached.          //
  66 //  5. The same as step 4 but events are stored.           //
  67 //                                                         //
  68 ////////////////////////////////////////////////////////////
  69
  70 #include <TCanvas.h>
  71
  72
  73 #include <TH3D.h>
  74 #include <TList.h>
  75 #include <TPDGCode.h>
  76 #include <TParticle.h>
  77 #include <AliStack.h>
  78 #include <TMath.h>
  79 #include <TVector3.h>
  80 #include <TStopwatch.h>
  81 #include <TFile.h>
  82
  83 #include "AliGenCocktailAfterBurner.h"
  84 #include "AliGeVSimParticle.h"
  85 #include "AliGenGeVSim.h"
  86 #include "AliGenHIJINGpara.h"
  87
  88
  89 /***********************************************************/
  90 ClassImp(AliGenHBTosl)
  91
  92 AliGenHBTosl::AliGenHBTosl():
  93  AliGenerator(),
  94  fQCoarseBackground(0x0),
  95  fQCoarseSignal(0x0),
  96  fQSignal(0x0),
  97  fQBackground(0x0),
  98  fQSecondSignal(0x0),
  99  fQSecondBackground(0x0),
 100  fQRange(0.06),
 101  fQNBins(60),
 102  fGenerator(0x0),
 103  fStackBuffer(0),
 104  fBufferSize(5),
 105  fNBinsToScale(Int_t(fQNBins*0.1)),
 106  fDebug(0),
 107  fSignalShapeCreated(0),
 108  fMaxIterations(1000),
 109  fMaxChiSquereChange(0.01),
 110  fMaxChiSquerePerNDF(1.5),
 111  fQRadius(8.0),
 112  fPID(kPiPlus),
 113  fSamplePhiMin(-0.01),
 114  fSamplePhiMax(TMath::TwoPi()+0.01),
 115  fSignalRegion(0.0),
 116  fMinFill(1000),
 117  fSwapped(0),
 118  fLogFile(0x0)
 119 {
 120 //default constructor
 121 }
 122 /***********************************************************/
 123
 124 AliGenHBTosl::AliGenHBTosl(Int_t n, Int_t pid):
 125  AliGenerator(n),
 126  fQCoarseBackground(0x0),
 127  fQCoarseSignal(0x0),
 128  fQSignal(0x0),
 129  fQBackground(0x0),
 130  fQSecondSignal(0x0),
 131  fQSecondBackground(0x0),
 132  fQRange(0.06),
 133  fQNBins(60),
 134  fGenerator(0x0),
 135  fStackBuffer(0),
 136  fBufferSize(5),
 137  fNBinsToScale(Int_t(fQNBins*0.1)),
 138  fDebug(0),
 139  fSignalShapeCreated(kFALSE),
 140  fMaxIterations(1000),
 141  fMaxChiSquereChange(0.01),
 142  fMaxChiSquerePerNDF(1.5),
 143  fQRadius(8.0),
 144  fPID(pid),
 145  fSamplePhiMin(-0.01),
 146  fSamplePhiMax(TMath::TwoPi()+0.01),
 147  fSignalRegion(0.0),
 148  fMinFill(1000),
 149  fSwapped(0),
 150  fLogFile(0x0)
 151 {
 152 //default constructor
 153 }
 154
 155 AliGenHBTosl::AliGenHBTosl(const AliGenHBTosl & hbt):
 156  AliGenerator(-1),
 157  fQCoarseBackground(0x0),
 158  fQCoarseSignal(0x0),
 159  fQSignal(0x0),
 160  fQBackground(0x0),
 161  fQSecondSignal(0x0),
 162  fQSecondBackground(0x0),
 163  fQRange(0.06),
 164  fQNBins(60),
 165  fGenerator(0x0),
 166  fStackBuffer(0),
 167  fBufferSize(5),
 168  fNBinsToScale(Int_t(fQNBins*0.1)),
 169  fDebug(0),
 170  fSignalShapeCreated(kFALSE),
 171  fMaxIterations(1000),
 172  fMaxChiSquereChange(0.01),
 173  fMaxChiSquerePerNDF(1.5),
 174  fQRadius(8.0),
 175  fPID(kPiPlus),
 176  fSamplePhiMin(-0.01),
 177  fSamplePhiMax(TMath::TwoPi()+0.01),
 178  fSignalRegion(0.0),
 179  fMinFill(1000),
 180  fSwapped(0),
 181  fLogFile(0x0)
 182 {
 183 // Copy constructor
 184     hbt.Copy(*this);
 185 }
 186 /***********************************************************/
 187
 188 AliGenHBTosl::~AliGenHBTosl()
 189 {
 190 //destructor
 191  delete fQCoarseSignal;
 192  delete fQCoarseBackground;
 193  delete fQSignal;
 194  delete fQBackground;
 195  delete fGenerator;
 196  delete fQSecondSignal;
 197  delete fQSecondBackground;
 198  delete fLogFile;
 199 }
 200 /***********************************************************/
 201
 202 void AliGenHBTosl::Init()
 203 {
 204   //Initializes generator
 205   if (fGenerator == 0x0)
 206    {
 207
 208      AliGenHIJINGpara* bkggen = new AliGenHIJINGpara(fNpart*4);
 209      fGenerator = bkggen;
 210
 211 /*
 212      AliGenGeVSim * gevsim = new AliGenGeVSim(0.0);
 213      AliGeVSimParticle* kplus = new AliGeVSimParticle(fPID,1,fNpart, 0.17, 0.9);
 214      gevsim->AddParticleType(kplus);
 215
 216      fGenerator = gevsim;
 217 */
 218
 219 /*
 220
 221     AliMevSimConfig *c = new AliMevSimConfig(1);
 222     c->SetRectPlane(1);                                 // reaction plane control, model 4
 223     c->SetGrid(80,80);
 224
 225     AliGenMevSim *mevsim = new AliGenMevSim(c);
 226     mevsim->SetPtRange(0.001, 3);
 227     mevsim->SetMomentumRange(0.1, 3);
 228     mevsim->SetTrackingFlag(0);
 229     mevsim->SetOrigin(0.0, 0.0, 0.0);
 230     mevsim->SetSigma(0.0, 0.0, 0.0);
 231     AliMevSimParticle *kplus = new AliMevSimParticle(kKPlus, fNpart, 0, 0.25, 0.0, 2, 0.15, 0.0, 0.0 );
 232     mevsim->AddParticleType(kplus);
 233     fGenerator = mevsim;
 234 */
 235
 236     fGenerator->SetOrigin(fOrigin[0],fOrigin[1],fOrigin[2]);
 237     static const Double_t kDegToRadCF = 180./TMath::Pi();
 238     fGenerator->SetMomentumRange(fPtMin,fPtMax);
 239     fGenerator->SetPhiRange(kDegToRadCF*fPhiMin,kDegToRadCF*fPhiMax);
 240     fGenerator->SetYRange(fYMin,fYMax);
 241     fGenerator->SetThetaRange(kDegToRadCF*fThetaMin,kDegToRadCF*fThetaMax);
 242     fGenerator->Init();
 243
 244    }
 245
 246 //  fQCoarseBackground =  new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 247 //  fQCoarseSignal =  new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 248 //  fQCoarseBackground->Sumw2();
 249 //  fQCoarseSignal->Sumw2();
 250
 251   fQSignal =  new TH3D("fQSignal1","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 252   fQBackground =  new TH3D("fQBackground1","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 253
 254   fQSecondSignal =  new TH3D("fQSignal2","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 255   fQSecondBackground =  new TH3D("fQBackground2","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 256
 257   fQSignal->Sumw2();
 258   fQBackground->Sumw2();
 259   fQSecondSignal->Sumw2();
 260   fQSecondBackground->Sumw2();
 261
 262   fLogFile = new ofstream("BadEvent",ios::out);
 263
 264 }
 265 /***********************************************************/
 266
 267 void AliGenHBTosl::Generate()
 268 {
 269  //the main method
 270
 271  ofstream& logfile = *fLogFile;
 272  logfile<<"Generate"<<"Attempts to generate "<<fNpart<<" particles.";
 273
 274
 275  if (fStackBuffer == 0x0) fStackBuffer = new TList();
 276  //Here is initialization level
 277  if (fSignalShapeCreated == kFALSE)
 278   {
 279    TH3D *hs = 0x0, *hb = 0x0;
 280    TFile* file;
 281
 282    file = TFile::Open("QTSignal.root");
 283    if (file)
 284     {
 285      hs = (TH3D*)file->Get("fQSignal1");
 286      if (hs) hs->SetDirectory(0x0);
 287     }
 288    delete file;
 289
 290    file = TFile::Open("QTBackground.root");
 291    if (file)
 292     {
 293      hb = (TH3D*)file->Get("fQBackground1");
 294      if (hb) hb->SetDirectory(0x0);
 295     }
 296    delete file;
 297
 298    if (hs && hb)
 299     {
 300       Info("Generate","**********************************");
 301       Info("Generate","Found starting histograms in files");
 302       Info("Generate","**********************************");
 303       delete fQSignal;
 304       delete fQBackground;
 305       fQSignal = hs;
 306       fQBackground = hb;
 307     }
 308    else
 309     {
 310       TH3D *cs = 0x0, *cb = 0x0;
 311       file = TFile::Open("QTCoarseBackground.root");
 312       if (file)
 313        {
 314         cb = (TH3D*)file->Get("fQCoarseBackground");
 315         if (cb) cb->SetDirectory(0x0);
 316        }
 317       delete file;
 318
 319       file = TFile::Open("QTCoarseSignal.root");
 320       if (file)
 321        {
 322         cs = (TH3D*)file->Get("fQCoarseSignal");
 323         if (cs) cs->SetDirectory(0x0);
 324        }
 325       delete file;
 326
 327       if (cs && cb)
 328        {
 329
 330          Info("Generate","Got Coarse signal and bkg from files");
 331          delete fQCoarseBackground;
 332          delete fQCoarseSignal;
 333          fQCoarseSignal = cs;
 334          fQCoarseBackground = cb;
 335        }
 336       else
 337        {
 338          if (cb)
 339            {
 340              Info("Generate","Got Coarse bkg from file");
 341              delete fQCoarseBackground;
 342              fQCoarseBackground = cb;
 343            }
 344          else
 345            {
 346              fQCoarseBackground =  new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 347              fQCoarseBackground->Sumw2();
 348
 349              FillCoarse();      //create coarse background - just to know the spectrum
 350            }
 351
 352          fQCoarseSignal =  new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 353          fQCoarseSignal->Sumw2();
 354          FillCoarseSignal();//create first coarse signal by brutal multplication coarse background and required function shape
 355        }
 356
 357       StartSignal();     //Initilizes the stack that is used for generation
 358     }
 359    fSignalShapeCreated = kTRUE;
 360   }
 361
 362  AliStack* stack = RotateStack();
 363
 364  AliStack* genstack = fGenerator->GetStack();
 365  if (genstack == 0x0)
 366   {
 367    genstack = new AliStack(fNpart);
 368    fGenerator->SetStack(genstack);
 369   }
 370  else
 371   {
 372    genstack->Reset();
 373   }
 374
 375  fGenerator->Generate();
 376  Int_t ntr = 0;
 377  if ( genstack->GetNtrack() < fNpart/2)
 378   {
 379     Warning("Generate","************************************************************");
 380     Warning("Generate","Generator generated (%d) less particles then expected (%d).",
 381              stack->GetNtrack(),fNpart/2);
 382     Warning("Generate","************************************************************");
 383   }
 384
 385  TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
 386  work->SetDirectory(0x0);
 387  work->Sumw2();
 388
 389  Double_t*** chiarray = new Double_t** [fQNBins+1];
 390  Double_t*** sigarray = new Double_t** [fQNBins+1];
 391
 392  for (Int_t i = 1; i<=fQNBins; i++)
 393    {
 394       chiarray[i] =  new Double_t* [fQNBins+1];
 395       sigarray[i] =  new Double_t* [fQNBins+1];
 396
 397       for (Int_t k = 1; k<=fQNBins; k++)
 398        {
 399          chiarray[i][k] =  new Double_t [fQNBins+1];
 400          sigarray[i][k] =  new Double_t [fQNBins+1];
 401        }
 402    }
 403
 404
 405  Double_t scale = Scale(fQSignal,fQBackground);
 406  work->Divide(fQSignal,fQBackground,scale);
 407
 408  Double_t binwdh = work->GetBinWidth(1)/2.;
 409
 410  for (Int_t k = 1; k<=fQNBins; k++)
 411    {
 412     Double_t z = work->GetZaxis()->GetBinCenter(k);
 413     for (Int_t j = 1; j<=fQNBins; j++)
 414       {
 415        Double_t y = work->GetYaxis()->GetBinCenter(j);
 416        for (Int_t i = 1; i<=fQNBins; i++)
 417          {
 418               sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
 419               Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qinv)
 420               Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
 421               Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
 422               chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
 423          }
 424        }
 425     }
 426
 427  char msg[1000];
 428  logfile<<endl;
 429  sprintf(msg,"\n");
 430  Int_t middlebin = fQNBins/2;
 431
 432  for (Int_t k = middlebin-5; k < middlebin+5; k++)
 433    {
 434      Double_t tx = work->GetXaxis()->GetBinCenter(30);
 435      Double_t ty = work->GetYaxis()->GetBinCenter(30);
 436      Double_t tz = work->GetZaxis()->GetBinCenter(k);
 437      sprintf(msg,"% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
 438      logfile<<msg;
 439    }
 440  logfile<<endl;
 441
 442  for (Int_t k = middlebin-5; k < middlebin+5; k++)
 443   {
 444     sprintf(msg,"% 6.5f ",work->GetBinContent(30,30,k));
 445     logfile<<msg;
 446   }
 447  logfile<<endl;
 448
 449  for (Int_t k = middlebin-5; k < middlebin+5; k++)
 450   {
 451     sprintf(msg,"% 6.5f ",chiarray[30][30][k]);
 452     logfile<<msg;
 453   }
 454  logfile<<endl;
 455
 456  TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
 457  TParticle* second = &particle;
 458
 459  Bool_t shortloop = kTRUE;
 460  Int_t sc = 0;//security check against infinite loop
 461
 462  while ( (ntr+1) < fNpart)
 463    {
 464      Int_t xmax = 1;
 465      Int_t ymax = 1;
 466      Int_t zmax = 1;
 467      Double_t qout;
 468      Double_t qside;
 469      Double_t qlong;
 470
 471
 472      Int_t loopmin;
 473      Int_t loopmax;
 474
 475      if (shortloop)
 476       {
 477         shortloop = kFALSE;
 478         loopmax = fQNBins;
 479         loopmin = 1;
 480       }
 481      else
 482       {
 483         shortloop = kTRUE;
 484         loopmax = fQNBins/2+fQNBins/4;
 485         loopmin = fQNBins/2-fQNBins/4;
 486       }
 487
 488
 489      for (Int_t k = loopmin; k <=loopmax; k++ )
 490       {
 491        qlong = work->GetZaxis()->GetBinCenter(k);
 492        for (Int_t j = loopmin; j<=loopmax; j++)
 493          {
 494           qside = work->GetYaxis()->GetBinCenter(j);
 495           for (Int_t i = loopmin; i<=loopmax; i++)
 496             {
 497               qout  = work->GetXaxis()->GetBinCenter(i);
 498               if (chiarray[xmax][ymax][zmax] < chiarray[i][j][k])
 499                {
 500                  xmax = i;
 501                  ymax = j;
 502                  zmax = k;
 503                }
 504
 505 //              Double_t qdist = TMath::Sqrt(qout*qout + qside*qside + qlong*qlong);
 506
 507 //              Double_t fact = chiarray[i][j][k];//chiarray is chi2
 508 //              if (fact > work->GetBinError(i,j,k))//if differece between what we want and
 509 //               {                                  //what we have is bigger than stat. error
 510 //                 xmax = i;                        //we force to fill that bin
 511 //                 ymax = j;
 512 //                 zmax = k;
 513 //                 break;
 514 //               }
 515             }
 516          }
 517       }
 518      Double_t qlongc = work->GetZaxis()->GetBinCenter(zmax);
 519      Double_t qsidec = work->GetYaxis()->GetBinCenter(ymax);
 520      Double_t qoutc  = work->GetXaxis()->GetBinCenter(xmax);
 521
 522
 523      sprintf(msg,"Generate Fill bin chi2(%d,%d,%d)=%f",xmax,ymax,zmax,chiarray[xmax][ymax][zmax]);
 524      logfile<<msg;
 525
 526      qout  = gRandom->Uniform(qoutc -binwdh, qoutc +binwdh);
 527      qside = gRandom->Uniform(qsidec-binwdh, qsidec+binwdh);
 528      qlong = gRandom->Uniform(qlongc-binwdh, qlongc+binwdh);
 529
 530      TParticle* first = 0;
 531      Int_t jj = 0;
 532
 533      while (jj < genstack->GetNtrack())
 534       {
 535         TParticle* tmpp = genstack->Particle(jj++);
 536         if (tmpp->GetPdgCode() == fPID)
 537          {
 538            if (CheckParticle(tmpp,0x0,stack) == kFALSE)
 539             {
 540               first = tmpp;
 541               break;
 542             }
 543          }
 544       }
 545
 546      if (first == 0x0)
 547       {
 548         if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
 549         break;
 550       }
 551
 552      //Here put the check if this particle do not fall into signal region with other partticle
 553
 554       Int_t retval = GetThreeD(first,second,qout,qside,qlong);
 555       if (retval)
 556        {
 557          //Info("StartSignal","Can not find momenta for this OSL and particle");
 558          continue;
 559        }
 560     //in case this particle is falling into signal area with another
 561     //particle we take a another pair
 562     //it can intruduce artificial correlations
 563      Bool_t checkresult = CheckParticle(second,first,stack);
 564      if ( checkresult  && (sc < 10) )
 565       {
 566         sc++;
 567         continue;
 568       }
 569      sc = 0;
 570
 571      //Put on output stack
 572      SetTrack(first,ntr);
 573      SetTrack(second,ntr);
 574
 575      //Put on internal stack
 576      Int_t etmp;
 577      SetTrack(first,etmp,stack);
 578      SetTrack(second,etmp,stack);
 579
 580      Double_t y = GetQOutQSideQLongCorrTheorValue(qoutc,qsidec,qlongc);
 581
 582      sigarray[xmax][ymax][zmax] ++;
 583      chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
 584      chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
 585
 586     }
 587
 588  Mix(fStackBuffer,fQBackground,fQSecondSignal); //upgrate background
 589  Mix(stack,fQSignal,fQSecondBackground); //upgrate signal
 590
 591  delete work;
 592
 593  for (Int_t i = 1; i<=fQNBins; i++)
 594    {
 595      for (Int_t k = 1; k<=fQNBins; k++)
 596       {
 597         delete [] chiarray[i][k];
 598         delete [] sigarray[i][k];
 599       }
 600      delete [] chiarray[i];
 601      delete [] sigarray[i];
 602    }
 603  delete [] chiarray;
 604  delete [] sigarray;
 605 }
 606 /***********************************************************/
 607
 608 void AliGenHBTosl::GetOneD(TParticle* first, TParticle* second,Double_t qinv)
 609 {
 610 //deprecated method that caclulates momenta of the second particle
 611 // out of qinv and the first particle
 612 //first particle is rotated that only X is non-zero
 613
 614
 615   Double_t m = first->GetMass();
 616   Double_t msqrd = m*m;
 617   Double_t fourmassSquered = 4.*msqrd;
 618
 619   //Condition that R must fullfill to be possible to have qinv less smaller then randomized
 620 //  Double_t rRange = qinv*TMath::Sqrt(qinv*qinv + fourmassSquered)/fourmassSquered;
 621 //  Double_t r = gRandom->Uniform(rRange);
 622
 623   Double_t r = gRandom->Uniform(qinv);
 624   Double_t phi = gRandom->Uniform(TMath::TwoPi());
 625
 626   Double_t firstPx = first->P();//first particle is rotated that only X is non-zero  thus P==Px
 627   Double_t px = 2.*msqrd*firstPx + firstPx*qinv*qinv;
 628   Double_t temp = qinv*qinv*qinv*qinv  + fourmassSquered * (qinv*qinv - r*r );
 629   if (temp < 0.0)
 630    {
 631      Error("GetOneD","temp is less then 0: %f",temp);
 632      return;
 633    }
 634   temp = temp*(msqrd+firstPx*firstPx);
 635
 636   px = (px - TMath::Sqrt(temp))/(2.*msqrd);
 637
 638   Double_t py = r*TMath::Sin(phi);
 639   Double_t pz = r*TMath::Cos(phi);
 640
 641   TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
 642   TVector3 vector(px,py,pz);
 643   Rotate(firstpvector,vector);
 644
 645   Double_t e = TMath::Sqrt(msqrd + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
 646   second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
 647 //  TParticle* f = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, firstPx,0,0,e=TMath::Sqrt(msqrd+firstPx*firstPx),0.0,0.0,0.0,0.0);
 648 //        TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
 649 //                px, py, pz, e, vx, vy, vz, tof);
 650
 651   AliDebug(1,Form("Randomized qinv = %f, obtained = %f",qinv,GetQInv(first,second)));
 652
 653 }
 654 /***********************************************************/
 655
 656 Int_t AliGenHBTosl::GetThreeD(TParticle* first,TParticle* second, Double_t qout, Double_t qside, Double_t qlong)
 657 {
 658 //deprecated method that caclulates momenta of the second particle
 659 //out of  qout qside and qlong and the first particle
 660   Double_t m = first->GetMass();
 661   Double_t m2 = m*m;
 662
 663   Double_t px = first->P();//first particle is rotated that only X is non-zero  thus P==Px
 664   Double_t px2 = px*px;
 665
 666
 667   Double_t qout2 = qout*qout;
 668   Double_t qside2 = qside*qside;
 669   Double_t qlong2 = qlong*qlong;
 670
 671
 672   Double_t util1 = 4.*px2 - qside2;//4*P1x^2 - Y^2
 673   if (util1 < 0)
 674    {
 675      Info("GetThreeD","4.*px2* - qside2 is negative px: %f, qside: %f",px,qside);
 676      return 1;
 677    }
 678   Double_t util2 = TMath::Sqrt(px2*qout2*util1);
 679
 680
 681   Double_t p2x,p2y,p2z;
 682
 683 //  if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))
 684   if (qout >= 0)
 685    {
 686      //p2x
 687      Double_t tmp = px*(2.*px2 - qside2);
 688      tmp -=  util2;
 689      p2x = tmp/(2.*px2);
 690
 691      //p2y
 692      tmp =  qout - TMath::Sqrt(util1);
 693      p2y = - (tmp*qside)/(2.*px);
 694
 695      //p2z
 696      tmp = 4.*m2 + 2.*qout2+qlong2;
 697      tmp *= px;
 698      tmp -= 2.*util2;//!!!
 699      tmp += 4.*px*px2;
 700      tmp *= qlong2;
 701
 702      Double_t m2px2 = m2+px2;
 703      Double_t tmp2 = m2px2*tmp;
 704
 705      tmp  = 4.*(m2px2+qout2) + qlong2;
 706      tmp *= px;
 707      tmp -= 4.*util2;
 708      tmp *= 4.*(m2px2) + qlong2;
 709      tmp *= qlong2*qlong2;
 710      tmp *= m2px2*m2px2;
 711      tmp *= px;
 712      if (tmp < 0)
 713       {
 714         Error("","Argument of sqrt is negative");
 715         return 1;
 716       }
 717
 718      tmp2 += TMath::Sqrt(tmp);
 719
 720      tmp = 8.0*px*m2px2*m2px2;
 721      p2z = -TMath::Sqrt(tmp2/tmp);
 722      if (qlong < 0) p2z = -p2z;
 723    }
 724   else
 725    {
 726      //p2x
 727      Double_t tmp = px*(2.*px2 - qside2);
 728      tmp +=  util2;
 729      p2x = tmp/(2.*px2);
 730
 731      //p2y
 732      tmp =  qout - TMath::Sqrt(util1);
 733      p2y = - (tmp*qside)/(2.*px);
 734
 735      //p2z
 736      tmp = 4.*m2 + 2.*qout2+qlong2;
 737      tmp *= px;
 738      tmp += 2.*util2;//!!!
 739      tmp += 4.*px*px2;
 740      tmp *= qlong2;
 741
 742      Double_t m2px2 = m2+px2;
 743      Double_t tmp2 = m2px2*tmp;
 744
 745      tmp  = 4.*(m2px2+qout2) + qlong2;
 746      tmp *= px;
 747      tmp += 4.*util2;
 748      tmp *= 4.*(m2px2) + qlong2;
 749      tmp *= qlong2*qlong2;
 750      tmp *= m2px2*m2px2;
 751      tmp *= px;
 752      if (tmp < 0)
 753       {
 754         Error("","Argument of sqrt is negative");
 755         return 1;
 756       }
 757
 758      tmp2 += TMath::Sqrt(tmp);
 759
 760      tmp = 8.0*px*m2px2*m2px2;
 761      p2z = -TMath::Sqrt(tmp2/tmp);
 762      if (qlong < 0) p2z = -p2z;
 763     }
 764
 765 //  if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))  p2z = -p2z;
 766
 767   TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
 768   TVector3 vector(p2x,p2y,p2z);
 769   Rotate(firstpvector,vector);
 770
 771   Double_t e = TMath::Sqrt(m2 + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
 772   second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
 773
 774 ////////////
 775   if ( AliDebugLevel() > 3 )
 776    {
 777      e=TMath::Sqrt(m2+px*px);
 778      TParticle* f = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, px , 0.0, 0.0, e,0.0,0.0,0.0,0.0);
 779
 780      e = TMath::Sqrt(m2 + p2x*p2x + p2y*p2y + p2z*p2z);
 781      TParticle* s = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, p2x, p2y, p2z, e, 0.0, 0.0, 0.0, 0.0);
 782
 783      Double_t qo, qs, ql;
 784      GetQOutQSideQLong(f,s,qo, qs, ql);
 785
 786      Info("GetThreeD","TEST");
 787      f->Print();
 788      s->Print();
 789      Info("GetThreeD","Required %f %f %f",qout,qside,qlong);
 790      Info("GetThreeD","Got      %f %f %f",qo,qs,ql);
 791      if ( qout == qo)
 792         if (qside == qs)
 793           if  (qlong == ql)
 794             Info("GetThreeD","!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 795    }
 796 /////////////
 797   return 0;
 798 }
 799 /***********************************************************/
 800
 801 void AliGenHBTosl::StartSignal()
 802 {
 803 //Starts the signal histograms
 804  ofstream& logfile = *fLogFile;
 805  logfile<<"\n\n\n\n";
 806  logfile<<"************************************************"<<endl;
 807  logfile<<"************************************************"<<endl;
 808  logfile<<"               StartSignal                      "<<endl;
 809  logfile<<"************************************************"<<endl;
 810  logfile<<"************************************************"<<endl;
 811
 812  AliStack* stack;
 813
 814  fSwapped = kFALSE;
 815
 816  TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
 817  TParticle* second = &particle;
 818
 819  TIter next(fStackBuffer);
 820  while(( stack=(AliStack*)next() ))
 821   {
 822     stack->Reset();
 823   }
 824
 825  AliStack* genstack = fGenerator->GetStack();
 826  if (genstack == 0x0)
 827   {
 828    genstack = new AliStack(fNpart);
 829    fGenerator->SetStack(genstack);
 830   }
 831  else
 832   {
 833    genstack->Reset();
 834   }
 835
 836  StartSignalPass1();
 837  //We alread have detailed histograms and we do not need Coarse anymore
 838  delete fQCoarseSignal;
 839  delete fQCoarseBackground;
 840  fQCoarseSignal = 0x0;
 841  fQCoarseBackground = 0x0;
 842
 843
 844  const Double_t kNDF = fQNBins*fQNBins*fQNBins;
 845
 846  TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
 847  work->Sumw2();
 848  work->SetDirectory(0x0);
 849
 850
 851  Double_t binwdh = work->GetBinWidth(1)/2.;
 852
 853  Double_t*** chiarray = new Double_t** [fQNBins+1];
 854  Double_t*** sigarray = new Double_t** [fQNBins+1];
 855
 856  for (Int_t i = 1; i<=fQNBins; i++)
 857    {
 858       chiarray[i] =  new Double_t* [fQNBins+1];
 859       sigarray[i] =  new Double_t* [fQNBins+1];
 860
 861       for (Int_t k = 1; k<=fQNBins; k++)
 862        {
 863          chiarray[i][k] =  new Double_t [fQNBins+1];
 864          sigarray[i][k] =  new Double_t [fQNBins+1];
 865        }
 866    }
 867
 868
 869  Float_t chisqrchange = fMaxChiSquereChange + 1.;
 870  Float_t chisqrPerDF = fMaxChiSquerePerNDF;
 871  Float_t chisqrold = 0.0;
 872
 873  Int_t counter = 1;
 874  Int_t niterations = 1;
 875  Int_t rotaxisorder = 1;//defines order of looping over 3D histo (X,Y,Z or Y,Z,X or Z,X,Y)
 876
 877  Bool_t flag = kTRUE;
 878  Bool_t shortloop = kTRUE;
 879  TCanvas* c1 = new TCanvas();
 880
 881
 882  printf("\n");
 883  Info("StartSignal","\n\n\n\nSecond Pass\n\n\n\n");
 884
 885  while ( ( (chisqrPerDF > fMaxChiSquereChange) || flag) && (niterations++ < fMaxIterations)  )
 886   {
 887
 888    logfile<<"StartSignal\n";
 889    logfile<<" Row 1 Theory, 2 current value, 3 Chi2 \n";
 890
 891    Double_t chisqrnew = 0.0;
 892    flag = kFALSE;
 893
 894    Double_t scale = Scale(fQSignal,fQBackground);
 895    work->Divide(fQSignal,fQBackground,scale);
 896
 897    if ( (counter%100) == 0)
 898     {
 899       c1->cd();
 900       char buff[50];
 901       sprintf(buff,"QTWorkPass2.%3d.root",counter);
 902       TFile* file = TFile::Open(buff,"update");
 903       work->Write();
 904       work->SetDirectory(0x0);
 905       file->Close();
 906       delete file;
 907
 908       sprintf(buff,"QTBackgroundPass2.%3d.root",counter);
 909       file = TFile::Open(buff,"update");
 910       fQBackground->Write();
 911       fQBackground->SetDirectory(0x0);
 912       file->Close();
 913       delete file;
 914
 915       sprintf(buff,"QTSignalPass2.%3d.root",counter);
 916       file = TFile::Open(buff,"update");
 917       fQSignal->Write();
 918       fQSignal->SetDirectory(0x0);
 919       file->Close();
 920       delete file;
 921     }
 922
 923    counter++;
 924    Int_t novertresh = 0;
 925    for (Int_t k = 1; k<=fQNBins; k++)
 926     {
 927       Double_t z = work->GetZaxis()->GetBinCenter(k);
 928       for (Int_t j = 1; j<=fQNBins; j++)
 929         {
 930           Double_t y = work->GetYaxis()->GetBinCenter(j);
 931           for (Int_t i = 1; i<=fQNBins; i++)
 932             {
 933               Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qout)
 934               sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
 935               Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
 936               Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
 937               chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
 938               Double_t be = work->GetBinError(i,j,k);
 939               chisqrnew += diff*diff/(be*be);//add up chisq
 940
 941               //even if algorithm is stable (chi sqr change less then threshold)
 942               //and any bin differs more then 5% from expected value we continue
 943               Double_t fact = diff;
 944               if (TMath::Abs(fact) > 0.1)
 945                {
 946                  flag = kTRUE;
 947                  novertresh++;
 948                }
 949             }
 950          }
 951      }
 952
 953
 954    char msg[1000];
 955
 956    printf("\n");
 957
 958    for (Int_t k = 25; k < 36; k++)
 959     {
 960       Double_t tx = work->GetXaxis()->GetBinCenter(30);
 961       Double_t ty = work->GetYaxis()->GetBinCenter(30);
 962       Double_t tz = work->GetZaxis()->GetBinCenter(k);
 963       sprintf(msg,"% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
 964       logfile<<msg;
 965     }
 966    logfile<<endl;
 967
 968    for (Int_t k = 25; k < 36; k++)
 969     {
 970       sprintf(msg,"%6.5f ",work->GetBinContent(30,30,k));
 971       logfile<<msg;
 972     }
 973    logfile<<endl;
 974
 975    for (Int_t k = 25; k < 36; k++)
 976     {
 977       sprintf(msg,"% 6.5f ",chiarray[30][30][k]);
 978       logfile<<msg;
 979     }
 980    logfile<<endl;
 981
 982    chisqrchange = TMath::Abs(chisqrnew - chisqrold)/chisqrnew;
 983    chisqrold = chisqrnew;
 984
 985    chisqrPerDF = chisqrnew/kNDF;
 986
 987    Info("StartSignal","Iteration %d Chi-sq change %f%%",niterations,chisqrchange*100.0);
 988    Info("StartSignal","ChiSq = %f, NDF = %f, ChiSq/NDF = %f",chisqrnew, kNDF, chisqrPerDF );
 989    Info("StartSignal","novertresh = %d",novertresh);
 990
 991
 992    stack = RotateStack();
 993    genstack->Reset();
 994    fGenerator->Generate();
 995    Int_t ninputparticle = 0, ntr = 0;
 996    if ( genstack->GetNtrack() < fNpart/2)
 997     {
 998       Warning("StartSignal","**********************************");
 999       Warning("StartSignal","Generator generated (%d) less particles then expected (%d).",
1000                genstack->GetNtrack(),fNpart/2);
1001       Warning("StartSignal","**********************************");
1002     }
1003
1004    Int_t sc = 0; //security check against infinite loop
1005    while ( (ntr+1) < fNpart)//ntr is number of track generated up to now
1006     {
1007      Int_t xmax = 1;
1008      Int_t ymax = 1;
1009      Int_t zmax = 1;
1010      Double_t qout;
1011      Double_t qside;
1012      Double_t qlong;
1013
1014      Int_t i,j,k;
1015
1016      Int_t* cx = 0x0;
1017      Int_t* cy = 0x0;
1018      Int_t* cz = 0x0;
1019
1020      switch (rotaxisorder)
1021       {
1022         case 1:
1023           cx = &i;
1024           cy = &j;
1025           cz = &k;
1026           break;
1027         case 2:
1028           cx = &j;
1029           cy = &k;
1030           cz = &i;
1031           break;
1032         case 3:
1033           cx = &k;
1034           cy = &i;
1035           cz = &j;
1036           break;
1037       }
1038      rotaxisorder++;
1039      if (rotaxisorder > 3) rotaxisorder = 1;
1040      Int_t nrange;
1041
1042      if (shortloop)
1043       {
1044         shortloop = kFALSE;
1045         nrange = fQNBins;
1046       }
1047      else
1048       {
1049         shortloop = kTRUE;
1050         nrange = fQNBins/4;
1051       }
1052
1053 //     Bool_t force = kFALSE;
1054      for ( k = 1; k <=nrange;k++ )
1055       {
1056        for ( j = 1; j<=nrange; j++)
1057          {
1058           for ( i = 1; i<=nrange; i++)
1059             {
1060               if ( (chiarray[*cx][*cy][*cz]) > (chiarray[xmax][ymax][zmax]) )
1061                {
1062                  xmax = *cx;
1063                  ymax = *cy;
1064                  zmax = *cz;
1065                }
1066
1067 //              Double_t fact = chiarray[*cx][*cy][*cz];//chiarray is chi2*qinv
1068 //              if (fact > work->GetBinError(*cx,*cy,*cz))//if differece between what we want and
1069 //               {                                  //what we have is bigger than stat. error
1070 //                                                  //we force to fill that bin
1071 //                 qout  = work->GetXaxis()->GetBinCenter(*cx);
1072 //                 qside = work->GetYaxis()->GetBinCenter(*cy);
1073 //                 qlong = work->GetZaxis()->GetBinCenter(*cz);
1074
1075 //                 Info("StartSignal"," bin: (%d,%d,%d) loop status (%d,%d,%d) \nUsing Force: chiarray: %f \nq(o,s,l): (%f,%f,%f)  signal: %d background: %d binerror: %f",
1076 //        *cx,*cy,*cz,i,j,k,fact,qout,qside,qlong,
1077 //        (Int_t)sigarray[*cx][*cy][*cz],(Int_t)fQBackground->GetBinContent(*cx,*cy,*cz),work->GetBinError(*cx,*cy,*cz));
1078 //                 force = kTRUE;
1079 //                 break;
1080 //               }
1081
1082             }
1083 //           if (force) break;
1084           }
1085 //         if (force) break;
1086         }
1087
1088
1089       qout  = work->GetXaxis()->GetBinCenter(xmax);
1090       qside = work->GetYaxis()->GetBinCenter(ymax);
1091       qlong = work->GetZaxis()->GetBinCenter(zmax);
1092
1093 //      Info("StartSignal"," bin: (%d,%d,%d) chiarray: %f \nq(o,s,l): (%f,%f,%f)  signal: %d background: %d binerror: %f",
1094 //            xmax,ymax,zmax,chiarray[xmax][ymax][zmax],qout,qside,qlong,
1095 //            (Int_t)sigarray[xmax][ymax][zmax],
1096 //            (Int_t)fQBackground->GetBinContent(xmax,ymax,zmax),
1097 //            work->GetBinError(xmax,ymax,zmax));
1098
1099       qout  = gRandom->Uniform(qout-binwdh,qout+binwdh);
1100       qside = gRandom->Uniform(qside-binwdh,qside+binwdh);
1101       qlong = gRandom->Uniform(qlong-binwdh,qlong+binwdh);
1102
1103       TParticle* first = 0;
1104       while (ninputparticle < genstack->GetNtrack())
1105        {
1106          TParticle* tmpp = genstack->Particle(ninputparticle++);
1107          if (tmpp->GetPdgCode() == fPID)
1108           {
1109            if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1110             {
1111               first = tmpp;
1112               break;
1113             }
1114           }
1115        }
1116
1117       if (first == 0x0)
1118        {
1119          if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
1120          break;
1121        }
1122
1123       Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1124       if (retval)
1125        {
1126          Info("StartSignal","Can not find momenta for this OSL and particle OSL = %f %f %f",qout,qside,qlong);
1127          first->Print();
1128          second->Print();
1129
1130          continue;
1131        }
1132      //in case this particle is falling into signal area with another
1133      //particle we take a another pair
1134      //it can intruduce artificial correlations
1135       Bool_t checkresult = CheckParticle(second,first,stack);
1136       if ( checkresult  && (sc < 10) )
1137        {
1138          sc++;
1139          continue;
1140        }
1141       sc = 0;
1142
1143       //Put on output stack
1144       SetTrack(first,ntr,stack);
1145       SetTrack(second,ntr,stack);
1146
1147       Double_t y = GetQOutQSideQLongCorrTheorValue(qout,qside,qlong);
1148
1149       sigarray[xmax][ymax][zmax] ++;
1150       chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
1151       chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
1152
1153     }
1154    Info("StartSignal","Mixing background...");
1155    Mix(fStackBuffer,fQBackground,fQSecondBackground); //upgrate background
1156    Info("StartSignal","Mixing signal...");
1157    Mix(stack,fQSignal,fQSecondSignal); //upgrate background
1158
1159
1160 //   if ( (chisqrPerDF < 2.0) && (fSwapped == kFALSE) )
1161 //     {
1162 //         SwapGeneratingHistograms();
1163 //     }
1164
1165  }
1166  TFile* filef = TFile::Open("QTBackground.root","recreate");
1167  fQBackground->Write();
1168  fQBackground->SetDirectory(0x0);
1169  filef->Close();
1170  delete filef;
1171
1172  filef = TFile::Open("QTSignal.root","recreate");
1173  fQSignal->Write();
1174  fQSignal->SetDirectory(0x0);
1175  filef->Close();
1176  delete filef;
1177
1178
1179  delete c1;
1180  delete work;
1181
1182  for (Int_t i = 1; i<=fQNBins; i++)
1183    {
1184      for (Int_t k = 1; k<=fQNBins; k++)
1185       {
1186         delete [] chiarray[i][k];
1187         delete [] sigarray[i][k];
1188       }
1189      delete [] chiarray[i];
1190      delete [] sigarray[i];
1191    }
1192  delete [] chiarray;
1193  delete [] sigarray;
1194
1195 }
1196 /***********************************************************/
1197
1198 void AliGenHBTosl::StartSignalPass1()
1199 {
1200  //This method makes first part of the initialization of working histograms
1201  //It randomizes qout, qside and qlong from the coarse signal histogram
1202
1203  Bool_t flag = kTRUE;
1204  TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
1205  TParticle* second = &particle;
1206  Double_t qout;
1207  Double_t qside;
1208  Double_t qlong;
1209
1210  Info("StartSignalPass1","\n\nFirst Pass\n\n");
1211
1212  while (flag)
1213   {
1214     Info("StartSignalPass1","NextEvent");
1215     AliStack* stack = RotateStack();
1216     AliStack* genstack = fGenerator->GetStack();
1217     genstack->Reset();
1218     fGenerator->Generate();
1219     Int_t j = 0, ntr = 0;
1220     if ( genstack->GetNtrack() < fNpart/2)
1221      {
1222        Warning("StartSignalPass1","**********************************");
1223        Warning("StartSignalPass1","Generator generated (%d) less particles then expected (%d).",
1224                 genstack->GetNtrack(),fNpart/2);
1225        Warning("StartSignalPass1","**********************************");
1226      }
1227
1228     Int_t sc = 0;//security check against infinite loop
1229     while ((ntr+1)<fNpart)
1230      {
1231
1232 //    Info("StartSignal","Number of track on output stack: = %d", ntr);
1233 //    Info("StartSignal","Number of track on input stack: = %d\n", j);
1234
1235       TParticle* first = 0;
1236       while (j < genstack->GetNtrack())
1237        {
1238          TParticle* tmpp = genstack->Particle(j++);
1239          if (tmpp->GetPdgCode() == fPID)
1240           {
1241            if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1242             {
1243               first = tmpp;
1244               break;
1245             }
1246            else
1247             {
1248               Info("StartSignalPass1","Particle did not pass the safety check 1");
1249               tmpp->Print();
1250             }
1251           }
1252        }
1253
1254       if (first == 0x0)
1255        {
1256          if ( fDebug > 2 ) Info("StartSignalPass1","No more particles of that type");
1257
1258          break;
1259        }
1260
1261       SetTrack(first,ntr,stack);
1262
1263       fQCoarseSignal->GetRandom3(qout,qside,qlong);
1264
1265       Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1266       if (retval)
1267        {
1268          //Info("StartSignal","Can not find momenta for this OSL and particle");
1269          continue;
1270        }
1271       //in case this particle is falling into signal area with another
1272       //particle we take a another pair
1273       //it can intruduce artificial correlations
1274        Bool_t checkresult = CheckParticle(second,first,stack);
1275        if ( checkresult  && (sc < 10) )
1276         {
1277           sc++;
1278           Info("StartSignalPass1","Particle did not pass the safety check 2");
1279           second->Print();
1280           continue;
1281         }
1282
1283        sc = 0;
1284
1285       SetTrack(second,ntr,stack);
1286      }
1287
1288     Mix(stack,fQSignal,fQSecondSignal);
1289     Mix(fStackBuffer,fQBackground,fQSecondBackground);
1290
1291     flag = kFALSE;
1292
1293     for (Int_t k = 1; k<=fQNBins; k++)
1294       {
1295        for (j = 1; j<=fQNBins; j++)
1296          {
1297            for (Int_t i = 1; i<=fQNBins; i++)
1298              {
1299                 if ( (fQBackground->GetBinContent(i,j,k) < fMinFill) )
1300                   {
1301                     //(fQSignal->GetBinContent(i,j,k) < fMinFill) ||
1302         Info("StartSignalPass1","bin (%d,%d,%d): signal=%f background=%f",i,j,k,
1303                           fQSignal->GetBinContent(i,j,k),fQBackground->GetBinContent(i,j,k));
1304         flag = kTRUE;//continue while
1305         break;//breakes for not while
1306                   }
1307              }
1308             if (flag == kTRUE) break;
1309           }
1310         if (flag == kTRUE) break;
1311       }
1312
1313   }//while (flag)
1314
1315
1316 }
1317 /***********************************************************/
1318
1319 void AliGenHBTosl::FillCoarseSignal()
1320 {
1321  //Makes coarse signal by multiplying the coarse background and required function
1322  Info("FillCoarseSignal","START");
1323  for (Int_t k = 1; k <=fQNBins ;k++ )
1324    {
1325     Double_t z = fQCoarseBackground->GetZaxis()->GetBinCenter(k);
1326     for (Int_t j = 1; j <=fQNBins; j++)
1327       {
1328        Double_t y = fQCoarseBackground->GetYaxis()->GetBinCenter(j);
1329        for (Int_t i = 1; i <=fQNBins; i++)
1330          {
1331            Double_t x  = fQCoarseBackground->GetXaxis()->GetBinCenter(i);
1332            Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);
1333            Info("FillCoarseSignal","Bin (%d,%d,%d): osl(%f,%f,%f)=%f",i,j,k,x,y,z,v);
1334            fQCoarseSignal->SetBinContent(i,j,k,v*fQCoarseBackground->GetBinContent(i,j,k));
1335          }
1336        }
1337     }
1338
1339   //if (AliDebugLevel())
1340   TestCoarseSignal();
1341
1342   Info("FillCoarseSignal","DONE");
1343 }
1344 /***********************************************************/
1345
1346 void AliGenHBTosl::FillCoarse()
1347 {
1348   //creates the statistical background histogram on the base of input from
1349   //fGenerator
1350   Info("FillCoarse","START");
1351
1352   AliStack* stack;
1353   Int_t niter = 0;
1354
1355   Bool_t cont;
1356   TH3D tmph("tmph","tmph",2,0,1,2,0,1,2,0,1);
1357   printf("\n");
1358
1359   do
1360    {
1361 //     if (niter > 20) break;
1362
1363      cout<<niter++<<"  bincont "<<fQCoarseBackground->GetBinContent(30,30,28)
1364                   <<" "<<fQCoarseBackground->GetBinContent(30,30,29)
1365                   <<" "<<fQCoarseBackground->GetBinContent(30,30,30)
1366                   <<" "<<fQCoarseBackground->GetBinContent(30,30,31)
1367                   <<" "<<fQCoarseBackground->GetBinContent(30,30,32)
1368                   <<"\n";
1369      fflush(0);
1370
1371      stack = RotateStack();
1372      fGenerator->SetStack(stack);
1373      fGenerator->Init();
1374      fGenerator->Generate();
1375
1376      Mix(fStackBuffer,fQCoarseBackground,&tmph);
1377
1378      cont = kFALSE;
1379
1380      Info("FillCoarse","fMinFill = %d",fMinFill);
1381      for (Int_t k = 1; k<=fQNBins; k++)
1382        {
1383         for (Int_t j = 1; j<=fQNBins; j++)
1384           {
1385             for (Int_t i = 1; i<=fQNBins; i++)
1386               {
1387                 if ( fQCoarseBackground->GetBinContent(i,j,k) < fMinFill )
1388                 {
1389                   cont = kTRUE;
1390                   Info("FillCoarse","bin (%d,%d,%d)=%f",i,j,k,fQCoarseBackground->GetBinContent(i,j,k));
1391                   break;
1392                 }
1393
1394               }
1395             if (cont) break;
1396           }
1397          if (cont) break;
1398         }
1399    }while(cont);
1400
1401   printf("\n");
1402   fGenerator->SetStack(0x0);
1403   Info("FillCoarse","DONE");
1404
1405 }
1406 /***********************************************************/
1407
1408 void AliGenHBTosl::Mix(TList* eventbuffer,TH3D* denominator,TH3D* denominator2)
1409 {
1410   //Fills denominators
1411   //Mixes events stored in the eventbuffer and fills the background histograms
1412   static TStopwatch stoper;
1413
1414   if (eventbuffer == 0x0)
1415    {
1416     Error("Mix","Buffer List is null.");
1417     return;
1418    }
1419
1420   if (denominator == 0x0)
1421    {
1422     Error("Mix","Denominator histogram is null.");
1423     return;
1424    }
1425
1426   if (denominator2 == 0x0)
1427    {
1428     Error("Mix","Denominator2 histogram is null.");
1429     return;
1430    }
1431
1432   Info("Mix","%s",denominator->GetName());
1433   stoper.Start();
1434
1435   TIter next(eventbuffer);
1436   AliStack* firstevent;
1437   AliStack* secondevent = 0x0;
1438
1439   while(( firstevent=(AliStack*)next() ))
1440    {
1441     if (secondevent == 0x0)
1442      {
1443        secondevent = firstevent;
1444        continue;
1445      }
1446 //    Info("Mix","Mixing %#x with %#x",firstevent,secondevent);
1447     for(Int_t j = 0; j < firstevent->GetNtrack(); j++ )
1448      {
1449        TParticle* firstpart = firstevent->Particle(j);
1450
1451        Float_t phi = firstpart->Phi();
1452        if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1453
1454 //       Info("Mix","Mixing %d phi %f min %f max %f",j,phi,fSamplePhiMin,fSamplePhiMax);
1455
1456        for(Int_t i = 0; i < secondevent->GetNtrack(); i++ )
1457          {
1458            TParticle* secondpart = secondevent->Particle(i);
1459            phi = secondpart->Phi();
1460            if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1461
1462            Double_t qout;
1463            Double_t qside;
1464            Double_t qlong;
1465            GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1466            denominator->Fill(qout,qside,qlong);
1467            denominator2->Fill(qout,qside,qlong);
1468          }
1469      }
1470
1471     secondevent = firstevent;
1472    }
1473   stoper.Stop();
1474   stoper.Print();
1475
1476 }
1477 /***********************************************************/
1478
1479 void AliGenHBTosl::Mix(AliStack* stack, TH3D* numerator, TH3D* numerator2)
1480 {
1481 //fils numerator with particles from stack
1482   static TStopwatch stoper;
1483   if (stack == 0x0)
1484    {
1485     Error("Mix","Stack is null.");
1486     return;
1487    }
1488
1489   if ( (numerator == 0x0) || (numerator2 == 0x0) )
1490    {
1491     Error("Mix","Numerator histogram is null.");
1492     return;
1493    }
1494
1495   Info("Mix","%s",numerator->GetName());
1496   stoper.Start();
1497
1498   for(Int_t j = 0; j < stack->GetNtrack(); j++ )
1499     {
1500       TParticle* firstpart = stack->Particle(j);
1501       Float_t phi = firstpart->Phi();
1502       if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1503
1504       for(Int_t i = j+1; i < stack->GetNtrack(); i++ )
1505        {
1506          TParticle* secondpart = stack->Particle(i);
1507          phi = secondpart->Phi();
1508          if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1509          Double_t qout;
1510          Double_t qside;
1511          Double_t qlong;
1512          GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1513          numerator->Fill(qout,qside,qlong);
1514          numerator2->Fill(qout,qside,qlong);
1515        }
1516     }
1517   stoper.Stop();
1518   stoper.Print();
1519
1520 }
1521 /***********************************************************/
1522
1523 Double_t AliGenHBTosl::GetQInv(TParticle* f, TParticle* s)
1524 {
1525 //calculates qinv
1526 // cout<<f->Px()<<"   "<<s->Px()<<endl;
1527  Double_t pxdiff = f->Px() - s->Px();
1528  Double_t pydiff = f->Py() - s->Py();
1529  Double_t pzdiff = f->Pz() - s->Pz();
1530  Double_t ediff  = f->Energy() - s->Energy();
1531
1532  Double_t qinvl = ediff*ediff - ( pxdiff*pxdiff + pydiff*pydiff + pzdiff*pzdiff );
1533  Double_t qinv = TMath::Sqrt(TMath::Abs(qinvl));
1534  return qinv;
1535 }
1536 /***********************************************************/
1537
1538 void  AliGenHBTosl::GetQOutQSideQLong(TParticle* f, TParticle* s,Double_t& out, Double_t& side, Double_t& lon)
1539 {
1540  //returns qout,qside and qlong of the pair of particles
1541  out = side = lon = 10e5;
1542
1543  Double_t pxsum = f->Px() + s->Px();
1544  Double_t pysum = f->Py() + s->Py();
1545  Double_t pzsum = f->Pz() + s->Pz();
1546  Double_t esum  = f->Energy() + s->Energy();
1547  Double_t pxdiff = f->Px() - s->Px();
1548  Double_t pydiff = f->Py() - s->Py();
1549  Double_t pzdiff = f->Pz() - s->Pz();
1550  Double_t ediff  = f->Energy() - s->Energy();
1551  Double_t kt =  0.5*TMath::Hypot(pxsum,pysum);
1552
1553  Double_t k2 = pxsum*pxdiff+pysum*pydiff;
1554
1555  if (kt == 0.0)
1556   {
1557      f->Print();
1558      s->Print();
1559      kt = 10e5;
1560   }
1561  else
1562   {
1563     out = 0.5*k2/kt;
1564     side = (f->Px()*s->Py()-s->Px()*f->Py())/kt;
1565   }
1566
1567  Double_t beta = pzsum/esum;
1568  Double_t gamma = 1.0/TMath::Sqrt((1.-beta)*(1.+beta));
1569
1570  lon = gamma * ( pzdiff - beta*ediff );
1571
1572 // out = TMath::Abs(out);
1573 // side = TMath::Abs(side);
1574 // lon = TMath::Abs(lon);
1575 }
1576
1577 /***********************************************************/
1578
1579 Double_t AliGenHBTosl::Scale(TH3D* num, TH3D* den)
1580 {
1581  //Calculates the factor that should be used to scale
1582  //quatience of num and den to 1 at tail
1583
1584   AliDebug(1,"Entered");
1585   if(!num)
1586    {
1587      AliError("No numerator");
1588      return 0.0;
1589    }
1590   if(!den)
1591    {
1592      AliError("No denominator");
1593      return 0.0;
1594    }
1595
1596   if(fNBinsToScale < 1)
1597    {
1598     AliError("Number of bins for scaling is smaller than 1");
1599     return 0.0;
1600    }
1601   Int_t fNBinsToScaleX = fNBinsToScale;
1602   Int_t fNBinsToScaleY = fNBinsToScale;
1603   Int_t fNBinsToScaleZ = fNBinsToScale;
1604
1605   Int_t nbinsX = num->GetNbinsX();
1606   if (fNBinsToScaleX > nbinsX)
1607    {
1608     AliError("Number of X bins for scaling is bigger thnan number of bins in histograms");
1609     return 0.0;
1610    }
1611
1612   Int_t nbinsY = num->GetNbinsX();
1613   if (fNBinsToScaleY > nbinsY)
1614    {
1615     AliError("Number of Y bins for scaling is bigger thnan number of bins in histograms");
1616     return 0.0;
1617    }
1618
1619   Int_t nbinsZ = num->GetNbinsZ();
1620   if (fNBinsToScaleZ > nbinsZ)
1621    {
1622     AliError("Number of Z bins for scaling is bigger thnan number of bins in histograms");
1623     return 0.0;
1624    }
1625
1626   AliDebug(1,"No errors detected");
1627
1628   Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
1629   Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
1630   Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
1631
1632   Double_t densum = 0.0;
1633   Double_t numsum = 0.0;
1634
1635   for (Int_t k = offsetZ; k<nbinsZ; k++)
1636     for (Int_t j = offsetY; j<nbinsY; j++)
1637       for (Int_t i = offsetX; i<nbinsX; i++)
1638        {
1639         if ( num->GetBinContent(i,j,k) > 0.0 )
1640          {
1641
1642            densum += den->GetBinContent(i,j,k);
1643            numsum += num->GetBinContent(i,j,k);
1644          }
1645        }
1646
1647   AliDebug(1,Form("numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
1648           numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ));
1649
1650   if (numsum == 0) return 0.0;
1651   Double_t ret = densum/numsum;
1652
1653   AliDebug(1,Form("returning %f",ret));
1654   return ret;
1655
1656 }
1657 /***********************************************************/
1658
1659 void AliGenHBTosl::TestCoarseSignal()
1660 {
1661 //Tests how works filling from generated histogram shape
1662   TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
1663
1664 //  for (Int_t i = 0; i < fQCoarseBackground->GetEntries() ;i++)
1665 //   {
1666 //     Double_t x,y,z;
1667 //     fQCoarseSignal->GetRandom3(x,y,z);
1668 //     work->Fill(x,y,z);
1669 //   }
1670
1671   TCanvas* c1 = new TCanvas();
1672   c1->cd();
1673   work->Draw();
1674   c1->SaveAs("QTwork.root");
1675   TFile* file = TFile::Open("QTwork.root","update");
1676 //  work->Write();
1677   work->SetDirectory(0x0);
1678   file->Close();
1679
1680   fQCoarseSignal->Draw();
1681   c1->SaveAs("QTCoarseSignal.root");
1682   file = TFile::Open("QTCoarseSignal.root","update");
1683   fQCoarseSignal->Write();
1684   fQCoarseSignal->SetDirectory(0x0);
1685   file->Close();
1686
1687   fQCoarseBackground->Draw();
1688   c1->SaveAs("QTCoarseBackground.root");
1689   file = TFile::Open("QTCoarseBackground.root","update");
1690   fQCoarseBackground->Write();
1691   fQCoarseBackground->SetDirectory(0x0);
1692   file->Close();
1693
1694   TH1 *result = (TH1*)fQCoarseBackground->Clone("ratio");
1695   result->SetTitle("ratio");
1696   Float_t normfactor = Scale(work,fQCoarseBackground);
1697   result->Divide(work,fQCoarseBackground,normfactor);//work
1698
1699
1700   c1->cd();
1701   result->Draw();
1702   c1->SaveAs("QTresult.root");
1703   file = TFile::Open("QTresult.root","update");
1704   result->Write();
1705   result->SetDirectory(0x0);
1706   file->Close();
1707
1708   delete work;
1709   delete c1;
1710 }
1711 /***********************************************************/
1712
1713 void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr)
1714 {
1715 //Shortcut to PushTrack(bla,bla,bla,bla.............)
1716    if (p->P() == 0.0)
1717     {
1718       Error("SetTrack(TParticle*,Int_t&)","Particle has zero momentum");
1719       return;
1720     }
1721
1722
1723    Int_t pdg = p->GetPdgCode();
1724    Double_t px = p->Px();
1725    Double_t py = p->Py();
1726    Double_t pz = p->Pz();
1727    Double_t e  = p->Energy();
1728    Double_t vx = p->Vx();
1729    Double_t vy = p->Vy();
1730    Double_t vz = p->Vz();
1731    Double_t tof = p->T();
1732
1733    TVector3 pol;
1734    p->GetPolarisation(pol);
1735
1736    Double_t polx = pol.X();
1737    Double_t poly = pol.Y();
1738    Double_t polz = pol.Z();
1739    TMCProcess mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
1740    Float_t weight = p->GetWeight();
1741
1742    AliGenerator::PushTrack(fTrackIt, -1, pdg, px, py, pz, e, vx, vy, vz, tof,polx, poly, polz, mech, ntr, weight);
1743 }
1744 /***********************************************************/
1745
1746 void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr, AliStack* stack) const
1747 {
1748 //Shortcut to SetTrack(bla,bla,bla,bla.............)
1749    if (p->P() == 0.0)
1750     {
1751       Error("SetTrack(TParticle*,Int_t&,AliStack*)","Particle has zero momentum");
1752       return;
1753     }
1754
1755    Int_t pdg = p->GetPdgCode();
1756    Double_t px = p->Px();
1757    Double_t py = p->Py();
1758    Double_t pz = p->Pz();
1759    Double_t e  = p->Energy();
1760
1761    stack->PushTrack(fTrackIt, -1, pdg, px, py, pz, e, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, kPPrimary, ntr,1,0);
1762 }
1763 /***********************************************************/
1764
1765 void AliGenHBTosl::Rotate(TVector3& relvector, TVector3& vector)
1766 {
1767 //This method rotates vector about the angeles that are needed to rotate
1768 //relvector from postion (firstPx,0,0) to its actual positon
1769 //In other words: To make equations easier
1770
1771   static TVector3 first;
1772   if (AliDebugLevel()>=1)
1773    {
1774      first.SetXYZ(relvector.x(),relvector.y(),relvector.z());
1775    }
1776
1777   Double_t firstPx = TMath::Sqrt( relvector.x()*relvector.x() +
1778                                   relvector.y()*relvector.y() +
1779                                   relvector.z()*relvector.z()   );
1780
1781   Double_t rotAngleZ = -TMath::ATan2(relvector.y(),relvector.x());//calculating rot angles
1782   relvector.RotateZ(rotAngleZ);
1783   rotAngleZ = -rotAngleZ;
1784   Double_t rotAngleY = -TMath::ATan2(relvector.z(),relvector.x());
1785
1786   vector.RotateY(rotAngleY);
1787   vector.RotateZ(rotAngleZ);
1788
1789   if (AliDebugLevel()>5)
1790    {
1791     TVector3 test(firstPx,0.0,0.0);
1792     test.RotateY(rotAngleY);
1793     test.RotateZ(rotAngleZ);
1794     AliInfo(Form("Rotation test: px %f %f",first.x(),test.x()));
1795     AliInfo(Form("Rotation test: py %f %f",first.y(),test.y()));
1796     AliInfo(Form("Rotation test: pz %f %f",first.z(),test.z()));
1797    }
1798 }
1799 /***********************************************************/
1800
1801 Double_t AliGenHBTosl::Rotate(Double_t x,Double_t y,Double_t z)
1802 {
1803 //Rotates vector to base where only x - coordinate is no-zero, and returns that
1804
1805   Double_t xylength = TMath::Hypot(x,y);
1806   Double_t sinphi = -y/xylength;
1807   Double_t cosphi = x/xylength;
1808
1809   Double_t xprime = cosphi*x - sinphi*y;
1810   Double_t yprime = sinphi*x + cosphi*y;
1811
1812   TVector3 v(x,y,z);
1813   Double_t a1 = -TMath::ATan2(v.Y(),v.X());
1814
1815   if (AliDebugLevel()>5)
1816    {
1817      AliInfo(Form("Xpr = %f  Ypr = %f",xprime,yprime));
1818      AliInfo(Form("Calc sin = %f, and %f",sinphi,TMath::Sin(a1)));
1819      AliInfo(Form("Calc cos = %f, and %f",cosphi,TMath::Cos(a1)));
1820    }
1821
1822   Double_t xprimezlength = TMath::Hypot(xprime,z);
1823
1824   Double_t sintheta = z/xprimezlength;
1825   Double_t costheta = xprime/xprimezlength;
1826
1827
1828   Double_t xbis = sintheta*z + costheta*(cosphi*x - sinphi*y);
1829
1830   AliInfo(Form("Calculated rot %f, modulus %f",xbis,TMath::Sqrt(x*x+y*y+z*z)));
1831   return xbis;
1832 }
1833 /***********************************************************/
1834
1835 AliStack* AliGenHBTosl::RotateStack()
1836 {
1837 //swaps to next stack last goes to first and is reseted
1838
1839  AliStack* stack;
1840  if ( fStackBuffer->GetSize() >= fBufferSize )
1841   {
1842     stack = (AliStack*)fStackBuffer->Remove(fStackBuffer->Last());
1843   }
1844  else
1845   {
1846     stack = new AliStack(fNpart);
1847   }
1848
1849  fStackBuffer->AddFirst(stack);
1850  stack->Reset();
1851  return stack;
1852 }
1853 /***********************************************************/
1854
1855 Double_t AliGenHBTosl::GetQInvCorrTheorValue(Double_t qinv) const
1856 {
1857 //Function (deprecated)
1858  static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1859
1860  return 1.0 + 0.5*TMath::Exp(-qinv*qinv*fQRadius*fQRadius/kFactorsqrd);
1861 }
1862 /***********************************************************/
1863
1864 Double_t AliGenHBTosl::GetQOutQSideQLongCorrTheorValue(Double_t& out, Double_t& side, Double_t& lon) const
1865 {
1866  //Theoretical function. Wa want to get correlation of the shape of this function
1867  static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1868  return 1.0 + 0.7*TMath::Exp(-fQRadius*fQRadius*(out*out+side*side+lon*lon)/kFactorsqrd);
1869 }
1870 /***********************************************************/
1871
1872 Bool_t AliGenHBTosl::CheckParticle(TParticle* p, TParticle* aupair ,AliStack* stack)
1873 {
1874  //Checks if a given particle is falling into signal region with any other particle
1875  //already existing on stack
1876   //PH return kFALSE;
1877
1878  if (fSignalRegion <=0) return kFALSE;
1879
1880  for (Int_t i = 0; i < stack->GetNtrack(); i++)
1881   {
1882     TParticle* part = stack->Particle(i);
1883     if (part == aupair) continue;
1884     Double_t qout = 10e5;
1885     Double_t qside= 10e5;
1886     Double_t qlong= 10e5;
1887     GetQOutQSideQLong(p,part,qout,qside,qlong);
1888
1889     if (TMath::Abs(qout)  < fSignalRegion)
1890       if (TMath::Abs(qside) < fSignalRegion)
1891        if (TMath::Abs(qlong) < fSignalRegion)
1892          return kTRUE;
1893   }
1894   return kFALSE;
1895 }
1896 /***********************************************************/
1897
1898 void AliGenHBTosl::SwapGeneratingHistograms()
1899 {
1900   //Checks if it is time to swap signal and background histograms
1901   //if yes it swaps them
1902   Int_t threshold = fMinFill;
1903   for (Int_t k = 1; k<=fQNBins; k++)
1904    {
1905      for (Int_t j = 1; j<=fQNBins; j++)
1906        {
1907          for (Int_t i = 1; i<=fQNBins; i++)
1908            {
1909              if ( fQSecondBackground->GetBinContent(i,j,k) < threshold) return;
1910            }
1911        }
1912
1913    }
1914
1915
1916   Info("SwapGeneratingHistograms","*******************************************");
1917   Info("SwapGeneratingHistograms","*******************************************");
1918   Info("SwapGeneratingHistograms","*******************************************");
1919   Info("SwapGeneratingHistograms","****   SWAPPING HISTOGRAMS             ****");
1920   Info("SwapGeneratingHistograms","*******************************************");
1921   Info("SwapGeneratingHistograms","*******************************************");
1922   Info("SwapGeneratingHistograms","*******************************************");
1923
1924
1925   TH3D* h = fQSignal;
1926   fQSignal = fQSecondSignal;
1927   fQSecondSignal = h;
1928   fQSecondSignal->Reset();
1929   fQSecondSignal->SetDirectory(0x0);
1930
1931   h = fQBackground;
1932   fQBackground = fQSecondBackground;
1933   fQSecondBackground = h;
1934   fQSecondBackground->Reset();
1935   fQSecondBackground->SetDirectory(0x0);
1936
1937   fSwapped = kTRUE;
1938
1939 }
1940
1941 AliGenHBTosl& AliGenHBTosl::operator=(const  AliGenHBTosl& rhs)
1942 {
1943 // Assignment operator
1944     rhs.Copy(*this);
1945     return *this;
1946 }
1947
1948 void AliGenHBTosl::Copy(TObject&) const
1949 {
1950     //
1951     // Copy
1952     //
1953     Fatal("Copy","Not implemented!\n");
1954 }
1955
1956