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