EVGEN/AliGenHBTosl.cxx

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