EVGEN/AliGenHBTosl.cxx

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