Generator for HBT
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 2 Jun 2004 20:05:39 +0000 (20:05 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 2 Jun 2004 20:05:39 +0000 (20:05 +0000)
EVGEN/AliGenHBTosl.cxx [new file with mode: 0644]
EVGEN/AliGenHBTosl.h [new file with mode: 0644]
EVGEN/EVGENLinkDef.h
EVGEN/libEVGEN.pkg

diff --git a/EVGEN/AliGenHBTosl.cxx b/EVGEN/AliGenHBTosl.cxx
new file mode 100644 (file)
index 0000000..77c91d9
--- /dev/null
@@ -0,0 +1,1822 @@
+#include "AliGenHBTosl.h" 
+//__________________________________________________________
+/////////////////////////////////////////////////////////////
+//                                                         //
+//  class AliGenHBTosl                                     //
+//                                                         //
+//  Genarator simulating particle correlations             //
+//                                                         //
+//  The main idea of the generator is to produce particles //
+//  according to some distribution of two particle         //
+//  property. In HBT they are qout,qsie and qlong.         //
+//  In order to be able to generate signal that produces   //
+//  given two particle correlation background must be      //
+//  known before in order to produce the shape of signal   //
+//  to randomize given distribution from.                  //
+//                                                         //
+//  The generator works as follows:                        //
+//  1. Coarse Background (fQCoarseBackground) is generated //
+//     ade  from the particles                             //
+//     given by the external generator (variable           //
+//     fGenerator) by the mixing technique.                //
+//  2. Coarse signal is prduced by multiplying Coarse      //
+//     background by a required function                   //
+//     See method FillCoarseSignal                         //
+//  3. Signal is randomized out of the coarse signal       //
+//     histogram (two particle property). First particle   //
+//     is taken from the external generator, and the       //
+//     second one is CALCULATED on the basis of the first  //
+//     one and the two particle property (qout,qside,qlong)//
+//     Background is made by the mixing out of the         //
+//     genereted signal events.                            //
+//     This step is cotinued up to the moment signal       //
+//     histogram has enough statistics (data member        //
+//     fMinFill)                                           //
+//     See method StartSignalPass1()                       //
+//  4. chi is calculated for each bin (chiarray variqable) // 
+//     (not the chi2 because sign is important)            //
+//     Two particle prioperty                              //
+//     (qout,qside,qlong) is chosen at the points that     //
+//     chi is the smallest. First particle is taken from   //
+//     the the external generator (fGenerator) and second's /
+//     momenta are caclulated out of their momenta and     //
+//     (qout,qside,qlong). Background is updated           //
+//     continuesely for all the events. This step is       //
+//     continued until stability conditions are fullfiled  //
+//     or maximum number of iteration is reached.          //
+//  5. The same as step 4 but events are stored.           //
+//                                                         //
+////////////////////////////////////////////////////////////
+
+#include <Riostream.h>
+#include <TCanvas.h>
+
+
+#include <TH3D.h>
+#include <TList.h>
+#include <TPDGCode.h>
+#include <TParticle.h>
+#include <AliStack.h>
+#include <TMath.h>
+#include <TVector3.h>
+#include <TStopwatch.h>
+#include <TFile.h>
+
+#include "AliGenCocktailAfterBurner.h"
+#include "AliGeVSimParticle.h"
+#include "AliGenGeVSim.h"
+#include "AliGenHIJINGpara.h"
+
+/***********************************************************/
+ClassImp(AliGenHBTosl)
+
+AliGenHBTosl::AliGenHBTosl():
+ AliGenerator(),
+ fQCoarseBackground(0x0),
+ fQCoarseSignal(0x0),
+ fQSignal(0x0),
+ fQBackground(0x0),
+ fQSecondSignal(0x0),
+ fQSecondBackground(0x0),
+ fQRange(0.15),
+ fQNBins(60),
+ fGenerator(0x0),
+ fBufferSize(5),
+ fNBinsToScale(Int_t(fQNBins*0.2)),
+ fDebug(kFALSE),
+ fMaxIterations(300),
+ fMaxChiSquereChange(0.01),
+ fMaxChiSquerePerNDF(1.5), 
+ fQRadius(8.0),
+ fPID(kPiPlus),
+ fSamplePhiMin(-0.01),
+ fSamplePhiMax(TMath::TwoPi()+0.01),
+ fSignalRegion(0.0),
+ fMinFill(300)
+{
+//default constructor
+}
+/***********************************************************/
+
+AliGenHBTosl::AliGenHBTosl(Int_t n, Int_t pid):
+ AliGenerator(n),
+ fQCoarseBackground(0x0),
+ fQCoarseSignal(0x0),
+ fQSignal(0x0),
+ fQBackground(0x0),
+ fQSecondSignal(0x0),
+ fQSecondBackground(0x0),
+ fQRange(0.15),
+ fQNBins(60),
+ fGenerator(0x0),
+ fBufferSize(5),
+ fNBinsToScale(Int_t(fQNBins*0.2)),
+ fDebug(kFALSE),
+ fSignalShapeCreated(kFALSE),
+ fMaxIterations(300),
+ fMaxChiSquereChange(0.01),
+ fMaxChiSquerePerNDF(1.5),
+ fQRadius(8.0),
+ fPID(pid),
+ fSamplePhiMin(-0.01),
+ fSamplePhiMax(TMath::TwoPi()+0.01),
+ fSignalRegion(0.0),
+ fMinFill(300)
+{
+//default constructor
+}
+/***********************************************************/
+
+AliGenHBTosl::~AliGenHBTosl()
+{
+//destructor
+ delete fQCoarseSignal;
+ delete fQCoarseBackground;
+ delete fQSignal;
+ delete fQBackground;
+ delete fGenerator; 
+ delete fQSecondSignal;
+ delete fQSecondBackground;
+
+}
+
+void AliGenHBTosl::Init()
+{
+  //Initializes generator
+  if (fGenerator == 0x0)
+   { 
+   
+     AliGenHIJINGpara* bkggen = new AliGenHIJINGpara(fNpart*4);
+     fGenerator = bkggen;
+
+/*     
+     AliGenGeVSim * gevsim = new AliGenGeVSim(0.0);
+     AliGeVSimParticle* kplus = new AliGeVSimParticle(fPID,1,fNpart, 0.17, 0.9);
+     gevsim->AddParticleType(kplus);
+     
+     fGenerator = gevsim;
+*/
+
+/*   
+   
+    AliMevSimConfig *c = new AliMevSimConfig(1);
+    c->SetRectPlane(1);                                 // reaction plane control, model 4
+    c->SetGrid(80,80);
+
+    AliGenMevSim *mevsim = new AliGenMevSim(c);
+    mevsim->SetPtRange(0.001, 3);
+    mevsim->SetMomentumRange(0.1, 3);
+    mevsim->SetTrackingFlag(0);
+    mevsim->SetOrigin(0.0, 0.0, 0.0);
+    mevsim->SetSigma(0.0, 0.0, 0.0);
+    AliMevSimParticle *kplus = new AliMevSimParticle(kKPlus, fNpart, 0, 0.25, 0.0, 2, 0.15, 0.0, 0.0 );
+    mevsim->AddParticleType(kplus);
+    fGenerator = mevsim;
+*/
+    
+    fGenerator->SetOrigin(fOrigin[0],fOrigin[1],fOrigin[2]);
+    static const Double_t degtoradcf = 180./TMath::Pi();
+    fGenerator->SetMomentumRange(fPtMin,fPtMax);
+    fGenerator->SetPhiRange(degtoradcf*fPhiMin,degtoradcf*fPhiMax);
+    fGenerator->SetYRange(fYMin,fYMax);
+    fGenerator->SetThetaRange(degtoradcf*fThetaMin,degtoradcf*fThetaMax);
+    fGenerator->Init();
+    
+   }
+
+  fQCoarseBackground =  new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
+  fQCoarseSignal =  new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
+  fQSignal =  new TH3D("fQSignal1","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
+  fQBackground =  new TH3D("fQBackground1","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
+
+  fQSecondSignal =  new TH3D("fQSignal2","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
+  fQSecondBackground =  new TH3D("fQBackground2","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
+
+  fQCoarseBackground->Sumw2();
+  fQCoarseSignal->Sumw2();
+  fQSignal->Sumw2();
+  fQBackground->Sumw2();
+  fQSecondSignal->Sumw2();
+  fQSecondBackground->Sumw2();
+}
+/***********************************************************/
+
+void AliGenHBTosl::Generate()
+{
+ //the main method
+ Info("Generate","Attempts to generate %d particles.",fNpart);
+ if (fStackBuffer == 0x0) fStackBuffer = new TList();
+ //Here is initialization level
+ if (fSignalShapeCreated == kFALSE)
+  {
+   TH3D *hs = 0x0, *hb = 0x0;
+   TFile* file;
+
+   file = TFile::Open("QTSignal.root");
+   if (file)
+    {
+     hs = (TH3D*)file->Get("fQSignal1");
+     if (hs) hs->SetDirectory(0x0);
+    }
+   delete file;
+   
+   file = TFile::Open("QTBackground.root");
+   if (file)
+    {
+     hb = (TH3D*)file->Get("fQBackground1");
+     if (hb) hb->SetDirectory(0x0);
+    }
+   delete file;
+   
+   if (hs && hb)
+    {
+      Info("Generate","**********************************");
+      Info("Generate","Found starting histograms in files");
+      Info("Generate","**********************************");
+      delete fQSignal;
+      delete fQBackground;
+      fQSignal = hs;
+      fQBackground = hb;
+    }
+   else
+    { 
+      TH3D *cs = 0x0, *cb = 0x0;
+      file = TFile::Open("QTCoarseBackground.root");
+      if (file)
+       {
+        cb = (TH3D*)file->Get("fQCoarseBackground");
+        if (cb) cb->SetDirectory(0x0);
+       }
+      delete file;
+
+      file = TFile::Open("QTCoarseSignal.root");
+      if (file)
+       {
+        cs = (TH3D*)file->Get("fQCoarseSignal");
+        if (cs) cs->SetDirectory(0x0);
+       }
+      delete file;
+
+      if (cs && cb)
+       {
+
+         Info("Generate","Got Coarse signal and bkg from files");
+         delete fQCoarseBackground;
+         delete fQCoarseSignal;
+         fQCoarseSignal = cs;
+         fQCoarseBackground = cb;
+       }
+      else
+       { 
+         if (cb)
+           {
+             Info("Generate","Got Coarse bkg from file");
+             delete fQCoarseBackground;
+             fQCoarseBackground = cb;
+           }
+         else
+           {  
+             FillCoarse();      //create coarse background - just to know the spectrum
+           }
+         FillCoarseSignal();//create first coarse signal by brutal multplication coarse background and required function shape
+       }
+       
+      StartSignal();     //Initilizes the stack that is used for generation
+    }  
+   fSignalShapeCreated = kTRUE;
+  }
+
+ AliStack* stack = RotateStack();
+
+ AliStack* genstack = fGenerator->GetStack();
+ if (genstack == 0x0)
+  {
+   genstack = new AliStack(fNpart);
+   fGenerator->SetStack(genstack);  
+  }
+ else
+  {
+   genstack->Reset();
+  }
+ fGenerator->Generate();
+ Int_t j = 0, ntr = 0;
+ if ( genstack->GetNtrack() < fNpart/2)
+  {
+    Warning("Generate","************************************************************");
+    Warning("Generate","Generator generated (%d) less particles then expected (%d).",
+             stack->GetNtrack(),fNpart/2);
+    Warning("Generate","************************************************************");
+  }
+
+ TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
+ Double_t*** chiarray = new Double_t** [fQNBins+1];
+ Double_t*** sigarray = new Double_t** [fQNBins+1];
+ for (Int_t i = 1; i<=fQNBins; i++)
+   {
+      chiarray[i] =  new Double_t* [fQNBins+1];
+      sigarray[i] =  new Double_t* [fQNBins+1];
+      
+      for (Int_t k = 1; k<=fQNBins; k++)
+       {
+         chiarray[i][k] =  new Double_t [fQNBins+1];
+         sigarray[i][k] =  new Double_t [fQNBins+1];
+       }
+   }
+
+      
+ Double_t scale = Scale(fQSignal,fQBackground);
+ work->Divide(fQSignal,fQBackground,scale);
+ Double_t binwdh = work->GetBinWidth(1)/2.;
+
+ for (Int_t k = 1; k<=fQNBins; k++)
+   {
+    Double_t z = work->GetZaxis()->GetBinCenter(k);
+    for (Int_t j = 1; j<=fQNBins; j++)
+      {  
+       Double_t y = work->GetYaxis()->GetBinCenter(j);
+       for (Int_t i = 1; i<=fQNBins; i++)
+         {
+              sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
+              Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qinv)
+              Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
+              Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
+              chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
+         }
+       }
+    }
+    
+ printf("\n");
+ Int_t middlebin = fQNBins/2;
+ for (Int_t k = middlebin-5; k < middlebin+5; k++)
+  {
+    printf("% 6.5f ",work->GetBinContent(1,1,k));
+  }
+ printf("\n");
+
+ for (Int_t k = middlebin-5; k < middlebin+5; k++)
+  {
+    printf("% 6.5f ",chiarray[1][1][k]);
+  }
+ printf("\n");
+ TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
+ TParticle* second = &particle;
+
+ Bool_t shortloop = kTRUE;
+ Int_t sc = 0;//security check against infinite loop
+
+ while ( (ntr+1) < fNpart)
+   {
+     Int_t xmax = 1; 
+     Int_t ymax = 1; 
+     Int_t zmax = 1; 
+     Double_t qout;
+     Double_t qside;
+     Double_t qlong;
+     
+     
+     Int_t loopmin;
+     Int_t loopmax;
+     
+     if (shortloop)
+      {
+        shortloop = kFALSE;
+        loopmax = fQNBins;
+        loopmin = 1;
+      }
+     else 
+      {
+        shortloop = kTRUE;
+        loopmax = fQNBins/2+fQNBins/4;
+        loopmin = fQNBins/2-fQNBins/4;
+      }
+     
+     
+     for (Int_t k = loopmin; k <=loopmax; k++ )
+      {
+       qlong = work->GetZaxis()->GetBinCenter(k);
+       for (Int_t j = loopmin; j<=loopmax; j++)
+         {
+          qside = work->GetYaxis()->GetBinCenter(j);
+          for (Int_t i = loopmin; i<=loopmax; i++)
+            {
+              qout  = work->GetXaxis()->GetBinCenter(i);
+              if (chiarray[xmax][ymax][zmax] < chiarray[i][j][k]) 
+               {
+                 xmax = i;
+                 ymax = j;
+                 zmax = k;
+               }  
+              
+//              Double_t qdist = TMath::Sqrt(qout*qout + qside*qside + qlong*qlong);
+              
+//              Double_t fact = chiarray[i][j][k];//chiarray is chi2
+//              if (fact > work->GetBinError(i,j,k))//if differece between what we want and 
+//               {                                  //what we have is bigger than stat. error
+//                 xmax = i;                        //we force to fill that bin
+//                 ymax = j;
+//                 zmax = k;
+//                 break;
+//               }
+            }
+         }
+      }
+     qlong = work->GetZaxis()->GetBinCenter(zmax);
+     qside = work->GetYaxis()->GetBinCenter(ymax);
+     qout  = work->GetXaxis()->GetBinCenter(xmax);
+     
+     Info("Generate","Fill bin (%d,%d,%d)",xmax,ymax,zmax);
+     
+     qout  = gRandom->Uniform(qout-binwdh,qout+binwdh);
+     qside = gRandom->Uniform(qside-binwdh,qside+binwdh);
+     qlong = gRandom->Uniform(qlong-binwdh,qlong+binwdh);
+
+     TParticle* first = 0;
+     while (j < genstack->GetNtrack())
+      {
+        TParticle* tmpp = genstack->Particle(j++);
+        if (tmpp->GetPdgCode() == fPID)
+         {
+           if (CheckParticle(tmpp,0x0,stack) == kFALSE)
+            {
+              first = tmpp;
+              break;
+            } 
+         }
+      } 
+     
+     if (first == 0x0)
+      {
+        if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
+        break;
+      }
+     
+     //Here put the check if this particle do not fall into signal region with other partticle
+     
+      Int_t retval = GetThreeD(first,second,qout,qside,qlong);
+      if (retval)
+       {
+         //Info("StartSignal","Can not find momenta for this OSL and particle");
+         continue;
+       }
+    //in case this particle is falling into signal area with another
+    //particle we take a another pair
+    //it can intruduce artificial correlations 
+     Bool_t checkresult = CheckParticle(second,first,stack);
+     if ( checkresult  && (sc < 10) ) 
+      { 
+        sc++;
+        continue;
+      }  
+     sc = 0;
+     
+     //Put on output stack
+     SetTrack(first,ntr);
+     SetTrack(second,ntr);
+
+     //Put on internal stack   
+     Int_t etmp;
+     SetTrack(first,etmp,stack);
+     SetTrack(second,etmp,stack);
+     
+     Double_t y = GetQOutQSideQLongCorrTheorValue(qout,qside,qlong);
+      
+     sigarray[xmax][ymax][zmax] ++; 
+     chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
+     chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
+      
+    }
+  
+ Mix(fStackBuffer,fQBackground,fQSecondSignal); //upgrate background
+ Mix(stack,fQSignal,fQSecondBackground); //upgrate signal
+  
+ for (Int_t i = 1; i<=fQNBins; i++)
+   {
+     for (Int_t k = 1; k<=fQNBins; k++)
+      {
+        delete [] chiarray[i][k]; 
+        delete [] sigarray[i][k];
+      }
+     delete [] chiarray[i];
+     delete [] sigarray[i];
+   }
+ delete [] chiarray;
+ delete [] sigarray;
+}
+/***********************************************************/
+
+void AliGenHBTosl::GetOneD(TParticle* first, TParticle* second,Double_t qinv)
+{
+//deprecated method that caclulates momenta of the second particle 
+// out of qinv and the first particle    
+//first particle is rotated that only X is non-zero
+
+    
+  Double_t m = first->GetMass();
+  Double_t msqrd = m*m;
+  Double_t fourmassSquered = 4.*msqrd;
+  
+  //Condition that R must fullfill to be possible to have qinv less smaller then randomized
+//  Double_t rRange = qinv*TMath::Sqrt(qinv*qinv + fourmassSquered)/fourmassSquered;
+//  Double_t r = gRandom->Uniform(rRange);
+
+  Double_t r = gRandom->Uniform(qinv);
+  Double_t phi = gRandom->Uniform(TMath::TwoPi());
+  
+  Double_t firstPx = first->P();//first particle is rotated that only X is non-zero  thus P==Px
+  Double_t px = 2.*msqrd*firstPx + firstPx*qinv*qinv;
+  Double_t temp = qinv*qinv*qinv*qinv  + fourmassSquered * (qinv*qinv - r*r );
+  if (temp < 0.0)
+   {
+     Error("GetOneD","temp is less then 0: %f",temp);
+     return;
+   }
+  temp = temp*(msqrd+firstPx*firstPx);
+  
+  px = (px - TMath::Sqrt(temp))/(2.*msqrd);
+  
+  Double_t py = r*TMath::Sin(phi);
+  Double_t pz = r*TMath::Cos(phi);
+  
+  TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
+  TVector3 vector(px,py,pz);
+  Rotate(firstpvector,vector);
+  Double_t e = TMath::Sqrt(msqrd + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
+  second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
+//  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);
+//        TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
+//                px, py, pz, e, vx, vy, vz, tof);
+  if (GetDebug()) Info("GetOneD","Randomized qinv = %f, obtained = %f",qinv,GetQInv(first,second));
+
+}
+/***********************************************************/
+
+Int_t AliGenHBTosl::GetThreeD(TParticle* first,TParticle* second, Double_t qout, Double_t qside, Double_t qlong)
+{
+//deprecated method that caclulates momenta of the second particle 
+//out of  qout qside and qlong and the first particle    
+  Double_t m = first->GetMass();
+  Double_t m2 = m*m;
+  
+  Double_t px = first->P();//first particle is rotated that only X is non-zero  thus P==Px
+  Double_t px2 = px*px;
+
+  
+  Double_t qout2 = qout*qout;
+  Double_t qside2 = qside*qside;
+  Double_t qlong2 = qlong*qlong;
+  
+  
+  Double_t util1 = 4.*px2 - qside2;//4*P1x^2 - Y^2
+  if (util1 < 0) 
+   {
+     Info("GetThreeD","4.*px2* - qside2 is negative px: %f, qside: %f",px,qside);
+     return 1;
+   }  
+  Double_t util2 = TMath::Sqrt(px2*qout2*util1);
+  
+
+  Double_t p2x,p2y,p2z;
+      
+//  if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))
+  if (qout >= 0)
+   {
+     //p2x
+     Double_t tmp = px*(2.*px2 - qside2);
+     tmp -=  util2;
+     p2x = tmp/(2.*px2);
+
+     //p2y
+     tmp =  qout - TMath::Sqrt(util1);
+     p2y = - (tmp*qside)/(2.*px);
+      
+     //p2z
+     tmp = 4.*m2 + 2.*qout2+qlong2;
+     tmp *= px;
+     tmp -= 2.*util2;//!!!
+     tmp += 4.*px*px2;
+     tmp *= qlong2;
+
+     Double_t m2px2 = m2+px2;
+     Double_t tmp2 = m2px2*tmp;
+
+     tmp  = 4.*(m2px2+qout2) + qlong2;
+     tmp *= px;
+     tmp -= 4.*util2;
+     tmp *= 4.*(m2px2) + qlong2;
+     tmp *= qlong2*qlong2;
+     tmp *= m2px2*m2px2;
+     tmp *= px;
+     if (tmp < 0)
+      {
+        Error("","Argument of sqrt is negative");
+        return 1;
+      }
+
+     tmp2 += TMath::Sqrt(tmp); 
+
+     tmp = 8.0*px*m2px2*m2px2;
+     p2z = -TMath::Sqrt(tmp2/tmp);
+     if (qlong < 0) p2z = -p2z;
+   }
+  else
+   {
+     //p2x
+     Double_t tmp = px*(2.*px2 - qside2);
+     tmp +=  util2;
+     p2x = tmp/(2.*px2);
+
+     //p2y
+     tmp =  qout - TMath::Sqrt(util1);
+     p2y = - (tmp*qside)/(2.*px);
+      
+     //p2z
+     tmp = 4.*m2 + 2.*qout2+qlong2;
+     tmp *= px;
+     tmp += 2.*util2;//!!!
+     tmp += 4.*px*px2;
+     tmp *= qlong2;
+
+     Double_t m2px2 = m2+px2;
+     Double_t tmp2 = m2px2*tmp;
+
+     tmp  = 4.*(m2px2+qout2) + qlong2;
+     tmp *= px;
+     tmp += 4.*util2;
+     tmp *= 4.*(m2px2) + qlong2;
+     tmp *= qlong2*qlong2;
+     tmp *= m2px2*m2px2;
+     tmp *= px;
+     if (tmp < 0)
+      {
+        Error("","Argument of sqrt is negative");
+        return 1;
+      }
+
+     tmp2 += TMath::Sqrt(tmp); 
+
+     tmp = 8.0*px*m2px2*m2px2;
+     p2z = -TMath::Sqrt(tmp2/tmp);
+     if (qlong < 0) p2z = -p2z;
+    }
+     
+//  if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))  p2z = -p2z;
+  
+  TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
+  TVector3 vector(p2x,p2y,p2z);
+  Rotate(firstpvector,vector);
+  Double_t e = TMath::Sqrt(m2 + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
+  second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
+  
+////////////  
+  if ( GetDebug() > 3 )
+   {
+     e=TMath::Sqrt(m2+px*px);  
+     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);
+
+     e = TMath::Sqrt(m2 + p2x*p2x + p2y*p2y + p2z*p2z);
+     TParticle* s = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, p2x, p2y, p2z, e, 0.0, 0.0, 0.0, 0.0);
+
+     Double_t qo, qs, ql;
+     GetQOutQSideQLong(f,s,qo, qs, ql);
+
+     Info("GetThreeD","TEST");
+     f->Print();
+     s->Print();
+     Info("GetThreeD","Required %f %f %f",qout,qside,qlong);
+     Info("GetThreeD","Got      %f %f %f",qo,qs,ql);
+     if ( qout == qo)
+        if (qside == qs)
+          if  (qlong == ql)
+            Info("GetThreeD","!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
+   }  
+///////////// 
+  return 0;
+}
+/***********************************************************/
+
+void AliGenHBTosl::StartSignal()
+{ 
+//Starts the signal histograms  
+ AliStack* stack;
+ fSwapped = kFALSE;
+ TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
+ TParticle* second = &particle;
+ TIter next(fStackBuffer);
+ while(( stack=(AliStack*)next() ))
+  {
+    stack->Reset();
+  }
+
+ AliStack* genstack = fGenerator->GetStack();
+ if (genstack == 0x0)
+  {
+   genstack = new AliStack(fNpart);
+   fGenerator->SetStack(genstack);  
+  }
+ else
+  {
+   genstack->Reset();
+  }
+ const Double_t ndf = fQNBins*fQNBins*fQNBins;
+ TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
+ work->Sumw2();
+ work->SetDirectory(0x0);
+ Double_t binwdh = work->GetBinWidth(1)/2.;
+ Double_t*** chiarray = new Double_t** [fQNBins+1];
+ Double_t*** sigarray = new Double_t** [fQNBins+1];
+ for (Int_t i = 1; i<=fQNBins; i++)
+   {
+      chiarray[i] =  new Double_t* [fQNBins+1];
+      sigarray[i] =  new Double_t* [fQNBins+1];
+      
+      for (Int_t k = 1; k<=fQNBins; k++)
+       {
+         chiarray[i][k] =  new Double_t [fQNBins+1];
+         sigarray[i][k] =  new Double_t [fQNBins+1];
+       }
+   }
+
+  
+ Float_t chisqrchange = fMaxChiSquereChange + 1.;
+ Float_t chisqrPerDF = fMaxChiSquerePerNDF;
+ Float_t chisqrold = 0.0;
+ Int_t counter = 1;
+ Int_t niterations = 1;
+ Int_t rotaxisorder = 1;//defines order of looping over 3D histo (X,Y,Z or Y,Z,X or Z,X,Y)
+ Bool_t flag = kTRUE;
+ Bool_t shortloop = kTRUE;
+ TCanvas* c1 = new TCanvas();
+
+ StartSignalPass1();
+
+ printf("\n");
+ Info("StartSignal","\n\n\n\nSecond Pass\n\n\n\n");
+
+ while ( ( (chisqrPerDF > fMaxChiSquereChange) || flag) && (niterations++ < fMaxIterations)  )
+  {
+   
+   printf("StartSignal\n");
+   printf(" Row 1 Theory, 2 current value, 3 Chi2 \n");
+
+   Double_t chisqrnew = 0.0;
+   flag = kFALSE;
+     
+   Double_t scale = Scale(fQSignal,fQBackground);
+   work->Divide(fQSignal,fQBackground,scale);
+   
+   if ( (counter%100) == 0) 
+    {
+      c1->cd();
+      char buff[50];
+      sprintf(buff,"QTWorkPass2.%3d.root",counter);
+      TFile* file = TFile::Open(buff,"update");
+      work->Write();
+      work->SetDirectory(0x0);
+      file->Close();
+      delete file;
+
+      sprintf(buff,"QTBackgroundPass2.%3d.root",counter);
+      file = TFile::Open(buff,"update");
+      fQBackground->Write();
+      fQBackground->SetDirectory(0x0);
+      file->Close();
+      delete file;
+
+      sprintf(buff,"QTSignalPass2.%3d.root",counter);
+      file = TFile::Open(buff,"update");
+      fQSignal->Write();
+      fQSignal->SetDirectory(0x0);
+      file->Close();
+      delete file;
+    }
+
+   counter++;
+   Int_t novertresh = 0;
+   for (Int_t k = 1; k<=fQNBins; k++)
+    {
+      Double_t z = work->GetZaxis()->GetBinCenter(k);
+      for (Int_t j = 1; j<=fQNBins; j++)
+        {  
+          Double_t y = work->GetYaxis()->GetBinCenter(j);
+          for (Int_t i = 1; i<=fQNBins; i++)
+            {
+              Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qout)
+              sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
+              Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
+              Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value 
+              chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
+              Double_t be = work->GetBinError(i,j,k);
+              chisqrnew += diff*diff/(be*be);//add up chisq
+
+              //even if algorithm is stable (chi sqr change less then threshold)
+              //and any bin differs more then 5% from expected value we continue
+              Double_t fact = diff;
+              if (TMath::Abs(fact) > 0.1) 
+               {
+                 flag = kTRUE; 
+                 novertresh++;
+               } 
+            }
+         }   
+     }
+    
+
+   printf("\n");
+
+   for (Int_t k = 1; k < 11; k++)
+    {
+      Double_t tx = work->GetXaxis()->GetBinCenter(1);
+      Double_t ty = work->GetYaxis()->GetBinCenter(1);
+      Double_t tz = work->GetZaxis()->GetBinCenter(k);
+      printf("% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
+    }
+   printf("\n");
+
+   for (Int_t k = 1; k < 11; k++)
+    {
+    printf("% 6.5f ",work->GetBinContent(1,1,k));
+    }
+   printf("\n");
+
+   for (Int_t k = 1; k < 11; k++)
+    {
+      printf("% 6.5f ",chiarray[1][1][k]);
+    }
+   printf("\n");
+    
+   chisqrchange = TMath::Abs(chisqrnew - chisqrold)/chisqrnew;
+   chisqrold = chisqrnew;
+
+   chisqrPerDF = chisqrnew/ndf;
+   
+   Info("StartSignal","Iteration %d Chi-sq change %f%%",niterations,chisqrchange*100.0);
+   Info("StartSignal","ChiSq = %f, NDF = %f, ChiSq/NDF = %f",chisqrnew, ndf, chisqrPerDF );
+   Info("StartSignal","novertresh = %d",novertresh);
+   
+   
+   stack = RotateStack();
+   genstack->Reset();
+   fGenerator->Generate();
+   Int_t ninputparticle = 0, ntr = 0;
+   if ( genstack->GetNtrack() < fNpart/2)
+    {
+      Warning("StartSignal","**********************************");
+      Warning("StartSignal","Generator generated (%d) less particles then expected (%d).",
+               genstack->GetNtrack(),fNpart/2);
+      Warning("StartSignal","**********************************");
+    }
+   
+   Int_t sc = 0; //security check against infinite loop
+   while ( (ntr+1) < fNpart)//ntr is number of track generated up to now
+    {
+     Int_t xmax = 1; 
+     Int_t ymax = 1; 
+     Int_t zmax = 1; 
+     Double_t qout;
+     Double_t qside;
+     Double_t qlong;
+     
+     Int_t i,j,k;
+     
+     Int_t* cx = 0x0;
+     Int_t* cy = 0x0;
+     Int_t* cz = 0x0;
+     
+     switch (rotaxisorder)
+      {
+        case 1:
+          cx = &i;
+          cy = &j;
+          cz = &k;
+          break;
+        case 2:
+          cx = &j;
+          cy = &k;
+          cz = &i;
+          break;
+        case 3:
+          cx = &k;
+          cy = &i;
+          cz = &j;
+          break;
+      }
+     rotaxisorder++;
+     if (rotaxisorder > 3) rotaxisorder = 1;
+     Int_t nrange;
+     
+     if (shortloop)
+      {
+        shortloop = kFALSE;
+        nrange = fQNBins;
+      }
+     else 
+      {
+        shortloop = kTRUE;
+        nrange = fQNBins/4;
+      }
+       
+//     Bool_t force = kFALSE; 
+     for ( k = 1; k <=nrange;k++ )
+      {
+       for ( j = 1; j<=nrange; j++)
+         {
+          for ( i = 1; i<=nrange; i++)
+            {
+              if ( (chiarray[*cx][*cy][*cz]) > (chiarray[xmax][ymax][zmax]) ) 
+               {
+                 xmax = *cx;
+                 ymax = *cy;
+                 zmax = *cz;
+               }  
+        
+//              Double_t fact = chiarray[*cx][*cy][*cz];//chiarray is chi2*qinv
+//              if (fact > work->GetBinError(*cx,*cy,*cz))//if differece between what we want and 
+//               {                                  //what we have is bigger than stat. error
+//                                                  //we force to fill that bin
+//                 qout  = work->GetXaxis()->GetBinCenter(*cx);
+//                 qside = work->GetYaxis()->GetBinCenter(*cy);
+//                 qlong = work->GetZaxis()->GetBinCenter(*cz);
+
+//                 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",
+//       *cx,*cy,*cz,i,j,k,fact,qout,qside,qlong,
+//       (Int_t)sigarray[*cx][*cy][*cz],(Int_t)fQBackground->GetBinContent(*cx,*cy,*cz),work->GetBinError(*cx,*cy,*cz));
+//                 force = kTRUE;
+//                 break;
+//               }
+               
+            }
+//           if (force) break; 
+          }
+//         if (force) break;
+        }
+
+
+      qout  = work->GetXaxis()->GetBinCenter(xmax);
+      qside = work->GetYaxis()->GetBinCenter(ymax);
+      qlong = work->GetZaxis()->GetBinCenter(zmax);
+
+//      Info("StartSignal"," bin: (%d,%d,%d) chiarray: %f \nq(o,s,l): (%f,%f,%f)  signal: %d background: %d binerror: %f",
+//            xmax,ymax,zmax,chiarray[xmax][ymax][zmax],qout,qside,qlong,
+//            (Int_t)sigarray[xmax][ymax][zmax],
+//            (Int_t)fQBackground->GetBinContent(xmax,ymax,zmax),
+//            work->GetBinError(xmax,ymax,zmax));
+          
+      qout  = gRandom->Uniform(qout-binwdh,qout+binwdh);
+      qside = gRandom->Uniform(qside-binwdh,qside+binwdh);
+      qlong = gRandom->Uniform(qlong-binwdh,qlong+binwdh);
+
+      TParticle* first = 0;
+      while (ninputparticle < genstack->GetNtrack())
+       {
+         TParticle* tmpp = genstack->Particle(ninputparticle++);
+         if (tmpp->GetPdgCode() == fPID)
+          {
+           if (CheckParticle(tmpp,0x0,stack) == kFALSE)
+            {
+              first = tmpp;
+              break;
+            } 
+          }
+       } 
+
+      if (first == 0x0)
+       {
+         if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
+         break;
+       }
+      
+      Int_t retval = GetThreeD(first,second,qout,qside,qlong);
+      if (retval)
+       {
+         Info("StartSignal","Can not find momenta for this OSL and particle OSL = %f %f %f",qout,qside,qlong);
+         first->Print();
+         second->Print();
+         
+         continue;
+       }
+     //in case this particle is falling into signal area with another
+     //particle we take a another pair
+     //it can intruduce artificial correlations 
+      Bool_t checkresult = CheckParticle(second,first,stack);
+      if ( checkresult  && (sc < 10) ) 
+       { 
+         sc++;
+         continue;
+       }  
+      sc = 0;
+
+      //Put on output stack
+      SetTrack(first,ntr,stack);
+      SetTrack(second,ntr,stack);
+
+      Double_t y = GetQOutQSideQLongCorrTheorValue(qout,qside,qlong);
+      
+      sigarray[xmax][ymax][zmax] ++; 
+      chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
+      chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
+      
+    }
+   Info("StartSignal","Mixing background...");
+   Mix(fStackBuffer,fQBackground,fQSecondBackground); //upgrate background
+   Info("StartSignal","Mixing signal...");
+   Mix(stack,fQSignal,fQSecondSignal); //upgrate background
+   
+   if ( (chisqrPerDF < 2.0) && (fSwapped == kFALSE) )
+     {
+         SwapGeneratingHistograms();
+     }
+   
+ }
+ TFile* filef = TFile::Open("QTBackground.root","recreate");
+ fQBackground->Write();
+ fQBackground->SetDirectory(0x0);
+ filef->Close();
+ delete filef;
+
+ filef = TFile::Open("QTSignal.root","recreate");
+ fQSignal->Write();
+ fQSignal->SetDirectory(0x0);
+ filef->Close();
+ delete filef;
+
+ delete c1;
+ delete work; 
+
+ for (Int_t i = 1; i<=fQNBins; i++)
+   {
+     for (Int_t k = 1; k<=fQNBins; k++)
+      {
+        delete [] chiarray[i][k]; 
+        delete [] sigarray[i][k];
+      }
+     delete [] chiarray[i];
+     delete [] sigarray[i];
+   }
+ delete [] chiarray;
+ delete [] sigarray;
+}
+/***********************************************************/
+
+void AliGenHBTosl::StartSignalPass1()
+{
+ //This method makes first part of the initialization of working histograms
+ //It randomizes qout, qside and qlong from the coarse signal histogram
+ Bool_t flag = kTRUE;
+ TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
+ TParticle* second = &particle;
+ Double_t qout;
+ Double_t qside;
+ Double_t qlong;
+ Info("StartSignalPass1","\n\nFirst Pass\n\n");
+ while (flag)
+  {
+    Info("StartSignalPass1","NextEvent");
+    AliStack* stack = RotateStack();
+    AliStack* genstack = fGenerator->GetStack();
+    genstack->Reset();
+    fGenerator->Generate();
+    Int_t j = 0, ntr = 0;
+    if ( genstack->GetNtrack() < fNpart/2)
+     {
+       Warning("StartSignalPass1","**********************************");
+       Warning("StartSignalPass1","Generator generated (%d) less particles then expected (%d).",
+                genstack->GetNtrack(),fNpart/2);
+       Warning("StartSignalPass1","**********************************");
+     }
+   
+    Int_t sc = 0;//security check against infinite loop
+    while ((ntr+1)<fNpart)
+     {
+  
+//    Info("StartSignal","Number of track on output stack: = %d", ntr);
+//    Info("StartSignal","Number of track on input stack: = %d\n", j);
+      
+      TParticle* first = 0;
+      while (j < genstack->GetNtrack())
+       {
+         TParticle* tmpp = genstack->Particle(j++);
+         if (tmpp->GetPdgCode() == fPID)
+          {
+           if (CheckParticle(tmpp,0x0,stack) == kFALSE)
+            {
+              first = tmpp;
+              break;
+            }
+           else
+            {
+              Info("StartSignalPass1","Particle did not pass the safety check 1");
+              tmpp->Print();
+            } 
+          }
+       } 
+
+      if (first == 0x0)
+       {
+         if ( fDebug > 2 ) Info("StartSignalPass1","No more particles of that type");
+         
+         break;
+       }
+
+      SetTrack(first,ntr,stack);
+      
+      fQCoarseSignal->GetRandom3(qout,qside,qlong);
+
+      Int_t retval = GetThreeD(first,second,qout,qside,qlong);
+      if (retval)
+       {
+         //Info("StartSignal","Can not find momenta for this OSL and particle");
+         continue;
+       }
+      //in case this particle is falling into signal area with another
+      //particle we take a another pair
+      //it can intruduce artificial correlations 
+       Bool_t checkresult = CheckParticle(second,first,stack);
+       if ( checkresult  && (sc < 10) ) 
+        { 
+          sc++;
+          Info("StartSignalPass1","Particle did not pass the safety check 2");
+          second->Print();
+          continue;
+        }
+       sc = 0;
+      
+      SetTrack(second,ntr,stack);
+     }
+    
+    Mix(stack,fQSignal,fQSecondSignal);
+    Mix(fStackBuffer,fQBackground,fQSecondBackground);
+    
+    flag = kFALSE;
+    
+    for (Int_t k = 1; k<=fQNBins; k++)
+      {
+       for (Int_t j = 1; j<=fQNBins; j++)
+         {  
+           for (Int_t i = 1; i<=fQNBins; i++)
+             {
+                if ( (fQBackground->GetBinContent(i,j,k) < fMinFill) )
+                  {
+                    //(fQSignal->GetBinContent(i,j,k) < fMinFill) || 
+       Info("StartSignalPass1","bin (%d,%d,%d): signal=%f background=%f",i,j,k,
+                          fQSignal->GetBinContent(i,j,k),fQBackground->GetBinContent(i,j,k));
+       flag = kTRUE;//continue while
+       break;//breakes for not while
+                  }
+             }
+            if (flag == kTRUE) break;
+          }
+        if (flag == kTRUE) break;
+      }
+    
+  }//while (flag)
+
+
+}
+/***********************************************************/
+
+void AliGenHBTosl::FillCoarseSignal()
+{
+ //Makes coarse signal by multiplying the coarse background and required function
+ Info("FillCoarseSignal","START");
+ for (Int_t k = 1; k <=fQNBins ;k++ )
+   {
+    Double_t z = fQCoarseBackground->GetZaxis()->GetBinCenter(k);
+    for (Int_t j = 1; j <=fQNBins; j++)
+      {
+       Double_t y = fQCoarseBackground->GetYaxis()->GetBinCenter(j);
+       for (Int_t i = 1; i <=fQNBins; i++)
+         {
+           Double_t x  = fQCoarseBackground->GetXaxis()->GetBinCenter(i);
+           Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);
+           Info("FillCoarseSignal","Bin (%d,%d,%d): osl(%f,%f,%f)=%f",i,j,k,x,y,z,v);
+           fQCoarseSignal->SetBinContent(i,j,k,v*fQCoarseBackground->GetBinContent(i,j,k));
+         }   
+       }  
+    }   
+  //if (GetDebug()) 
+  TestCoarseSignal();
+  Info("FillCoarseSignal","DONE"); 
+}
+/***********************************************************/
+
+void AliGenHBTosl::FillCoarse()
+{
+  //creates the statistical background histogram on the base of input from
+  //fGenerator
+  Info("FillCoarse","START");
+   
+  AliStack* stack;
+  Int_t niter = 0;
+  
+  Bool_t cont;
+  TH3D tmph("tmph","tmph",2,0,1,2,0,1,2,0,1);
+  printf("\n");
+
+  do 
+   { 
+//     if (niter > 20) break;
+     
+     cout<<niter++<<"  bincont "<<fQCoarseBackground->GetBinContent(1,1,1)
+                  <<" "<<fQCoarseBackground->GetBinContent(1,1,2)
+                  <<" "<<fQCoarseBackground->GetBinContent(1,1,3)
+                  <<" "<<fQCoarseBackground->GetBinContent(1,1,4)
+                  <<" "<<fQCoarseBackground->GetBinContent(1,1,5)
+                  <<"\n";
+     fflush(0);
+
+     stack = RotateStack();
+     fGenerator->SetStack(stack);
+     fGenerator->Init();
+     fGenerator->Generate();
+
+     Mix(fStackBuffer,fQCoarseBackground,&tmph);
+     
+     cont = kFALSE;
+     
+     Info("FillCoarse","fMinFill = %d",fMinFill);
+     for (Int_t k = 1; k<=fQNBins; k++)
+       {
+        for (Int_t j = 1; j<=fQNBins; j++)
+          {  
+            for (Int_t i = 1; i<=fQNBins; i++)
+              {
+                if ( fQCoarseBackground->GetBinContent(i,j,k) < fMinFill )
+                {
+                  cont = kTRUE;
+                  Info("FillCoarse","bin (%d,%d,%d)=%f",i,j,k,fQCoarseBackground->GetBinContent(i,j,k));
+                  break;
+                }
+
+              }
+            if (cont) break;
+          }
+         if (cont) break;
+        }
+   }while(cont);
+   
+  printf("\n");
+  fGenerator->SetStack(0x0);
+  Info("FillCoarse","DONE");
+  
+}
+/***********************************************************/
+void AliGenHBTosl::Mix(TList* eventbuffer,TH3D* denominator,TH3D* denominator2)
+{
+  //Fills denominators
+  //Mixes events stored in the eventbuffer and fills the background histograms
+  static TStopwatch stoper;
+  
+  if (eventbuffer == 0x0)
+   {
+    Error("Mix","Buffer List is null.");
+    return;
+   }
+
+  if (denominator == 0x0)
+   {
+    Error("Mix","Denominator histogram is null.");
+    return;
+   }
+
+  if (denominator2 == 0x0)
+   {
+    Error("Mix","Denominator2 histogram is null.");
+    return;
+   }
+
+  Info("Mix","%s",denominator->GetName());
+  stoper.Start();
+  
+  TIter next(eventbuffer);
+  AliStack* firstevent;
+  AliStack* secondevent = 0x0;
+  
+  while(( firstevent=(AliStack*)next() ))
+   {
+    if (secondevent == 0x0) 
+     {
+       secondevent = firstevent;
+       continue;
+     }
+//    Info("Mix","Mixing %#x with %#x",firstevent,secondevent);
+    for(Int_t j = 0; j < firstevent->GetNtrack(); j++ )
+     { 
+       TParticle* firstpart = firstevent->Particle(j);
+       
+       Float_t phi = firstpart->Phi();
+       if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
+       
+//       Info("Mix","Mixing %d phi %f min %f max %f",j,phi,fSamplePhiMin,fSamplePhiMax);
+
+       for(Int_t i = 0; i < secondevent->GetNtrack(); i++ )
+         {
+           TParticle* secondpart = secondevent->Particle(i);
+           phi = secondpart->Phi();
+           if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
+           
+           Double_t qout;
+           Double_t qside;
+           Double_t qlong;
+           GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
+           denominator->Fill(qout,qside,qlong);
+           denominator2->Fill(qout,qside,qlong);
+         }
+     }
+
+    secondevent = firstevent;
+   }
+  stoper.Stop();
+  stoper.Print();
+  
+}
+/***********************************************************/
+
+void AliGenHBTosl::Mix(AliStack* stack, TH3D* numerator, TH3D* numerator2)
+{
+//fils numerator with particles from stack
+  static TStopwatch stoper;
+  if (stack == 0x0)
+   {
+    Error("Mix","Stack is null.");
+    return;
+   }
+
+  if ( (numerator == 0x0) || (numerator2 == 0x0) )
+   {
+    Error("Mix","Numerator histogram is null.");
+    return;
+   }
+
+  Info("Mix","%s",numerator->GetName());
+  stoper.Start();
+
+  for(Int_t j = 0; j < stack->GetNtrack(); j++ )
+    { 
+      TParticle* firstpart = stack->Particle(j);
+      Float_t phi = firstpart->Phi();
+      if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
+       
+      for(Int_t i = j+1; i < stack->GetNtrack(); i++ )
+       {
+         TParticle* secondpart = stack->Particle(i);
+         phi = secondpart->Phi();
+         if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
+         Double_t qout;
+         Double_t qside;
+         Double_t qlong;
+         GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
+         numerator->Fill(qout,qside,qlong);
+         numerator2->Fill(qout,qside,qlong);
+       }
+    }
+  stoper.Stop();
+  stoper.Print();
+  
+}
+/***********************************************************/
+
+Double_t AliGenHBTosl::GetQInv(TParticle* f, TParticle* s)
+{
+//calculates qinv 
+// cout<<f->Px()<<"   "<<s->Px()<<endl;
+ Double_t pxdiff = f->Px() - s->Px();
+ Double_t pydiff = f->Py() - s->Py();
+ Double_t pzdiff = f->Pz() - s->Pz();
+ Double_t ediff  = f->Energy() - s->Energy();
+ Double_t qinvl = ediff*ediff - ( pxdiff*pxdiff + pydiff*pydiff + pzdiff*pzdiff );
+ Double_t qinv = TMath::Sqrt(TMath::Abs(qinvl)); 
+ return qinv;
+}
+/***********************************************************/
+
+void  AliGenHBTosl::GetQOutQSideQLong(TParticle* f, TParticle* s,Double_t& out, Double_t& side, Double_t& lon)
+{
+ //returns qout,qside and qlong of the pair of particles
+ out = side = lon = 10e5;
+ Double_t pxsum = f->Px() + s->Px();
+ Double_t pysum = f->Py() + s->Py();
+ Double_t pzsum = f->Pz() + s->Pz();
+ Double_t esum  = f->Energy() + s->Energy();
+ Double_t pxdiff = f->Px() - s->Px();
+ Double_t pydiff = f->Py() - s->Py();
+ Double_t pzdiff = f->Pz() - s->Pz();
+ Double_t ediff  = f->Energy() - s->Energy();
+ Double_t kt =  0.5*TMath::Hypot(pxsum,pysum);
+
+ Double_t k2 = pxsum*pxdiff+pysum*pydiff;
+ if (kt == 0.0)
+  {
+     f->Print();
+     s->Print();
+     kt = 10e5;
+  }
+ else
+  { 
+    out = 0.5*k2/kt;
+    side = (f->Px()*s->Py()-s->Px()*f->Py())/kt;
+  }
+
+ Double_t beta = pzsum/esum;
+ Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
+
+ lon = gamma * ( pzdiff - beta*ediff );
+
+// out = TMath::Abs(out);
+// side = TMath::Abs(side);
+// lon = TMath::Abs(lon);
+}
+
+/***********************************************************/
+
+Double_t AliGenHBTosl::Scale(TH3D* num, TH3D* den)
+{
+ //Calculates the factor that should be used to scale
+ //quatience of num and den to 1 at tail
+
+  if (GetDebug()) Info("Scale","Entered Scale()");
+  if(!num)
+   {
+     Error("Scale","No numerator");
+     return 0.0;
+   }
+  if(!den)
+   {
+     Error("Scale","No denominator");
+     return 0.0;
+   }
+
+  if(fNBinsToScale < 1)
+   {
+    return 0.0;
+    Error("Scale","Number of bins for scaling is smaller thnan 1");
+   }
+  Int_t fNBinsToScaleX = fNBinsToScale;
+  Int_t fNBinsToScaleY = fNBinsToScale;
+  Int_t fNBinsToScaleZ = fNBinsToScale;
+
+  Int_t nbinsX = num->GetNbinsX();
+  if (fNBinsToScaleX > nbinsX) 
+   {
+    Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
+    return 0.0;
+   }
+   
+  Int_t nbinsY = num->GetNbinsX();
+  if (fNBinsToScaleY > nbinsY) 
+   {
+    Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
+    return 0.0;
+   }
+
+  Int_t nbinsZ = num->GetNbinsZ();
+  if (fNBinsToScaleZ > nbinsZ) 
+   {
+    Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
+    return 0.0;
+   }
+
+  if (GetDebug()>0) Info("Scale","No errors detected");
+
+  Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
+  Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
+  Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
+
+  Double_t densum = 0.0;
+  Double_t numsum = 0.0;
+  
+  for (Int_t k = offsetZ; k<nbinsZ; k++)
+    for (Int_t j = offsetY; j<nbinsY; j++)
+      for (Int_t i = offsetX; i<nbinsX; i++)
+       {
+        if ( num->GetBinContent(i,j,k) > 0.0 )
+         {
+           
+           densum += den->GetBinContent(i,j,k);
+           numsum += num->GetBinContent(i,j,k);
+         }
+       }
+  
+  if(GetDebug() > 0) 
+    Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
+          numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ);
+  
+  if (numsum == 0) return 0.0;
+  Double_t ret = densum/numsum;
+
+  if(GetDebug() > 0) Info("Scale","returning %f",ret);
+  return ret;
+   
+}
+/***********************************************************/
+
+void AliGenHBTosl::TestCoarseSignal()
+{
+//Tests how works filling from generated histogram shape
+  TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
+     
+//  for (Int_t i = 0; i < fQCoarseBackground->GetEntries() ;i++)
+//   {
+//     Double_t x,y,z;
+//     fQCoarseSignal->GetRandom3(x,y,z);
+//     work->Fill(x,y,z);
+//   }
+
+  TCanvas* c1 = new TCanvas();
+  c1->cd();
+  work->Draw();
+  c1->SaveAs("QTwork.root");
+  TFile* file = TFile::Open("QTwork.root","update");
+//  work->Write();
+  work->SetDirectory(0x0);
+  file->Close();
+  
+  fQCoarseSignal->Draw(); 
+  c1->SaveAs("QTCoarseSignal.root");
+  file = TFile::Open("QTCoarseSignal.root","update");
+  fQCoarseSignal->Write();
+  fQCoarseSignal->SetDirectory(0x0);
+  file->Close();
+  
+  fQCoarseBackground->Draw(); 
+  c1->SaveAs("QTCoarseBackground.root");
+  file = TFile::Open("QTCoarseBackground.root","update");
+  fQCoarseBackground->Write();
+  fQCoarseBackground->SetDirectory(0x0);
+  file->Close();
+  
+  TH1 *result = (TH1*)fQCoarseBackground->Clone("ratio");
+  result->SetTitle("ratio");
+  Float_t normfactor = Scale(work,fQCoarseBackground);
+  result->Divide(work,fQCoarseBackground,normfactor);//work
+  
+  
+  c1->cd();
+  result->Draw();
+  c1->SaveAs("QTresult.root");
+  file = TFile::Open("QTresult.root","update");
+  result->Write();
+  result->SetDirectory(0x0);
+  file->Close();
+  
+  delete work;
+  delete c1;
+}
+/***********************************************************/
+
+void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr)
+{
+//Shortcut to PushTrack(bla,bla,bla,bla.............)
+   if (p->P() == 0.0)
+    {
+      Error("SetTrack(TParticle*,Int_t&)","Particle has zero momentum");
+      return;
+    }
+   
+   
+   Int_t pdg = p->GetPdgCode();
+   Double_t px = p->Px();
+   Double_t py = p->Py();
+   Double_t pz = p->Pz();
+   Double_t e  = p->Energy();
+   Double_t vx = p->Vx();
+   Double_t vy = p->Vy();
+   Double_t vz = p->Vz();
+   Double_t tof = p->T();
+
+   TVector3 pol;
+   p->GetPolarisation(pol);
+
+   Double_t polx = pol.X();
+   Double_t poly = pol.Y();
+   Double_t polz = pol.Z();
+   TMCProcess mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
+   Float_t weight = p->GetWeight();
+
+   AliGenerator::PushTrack(fTrackIt, -1, pdg, px, py, pz, e, vx, vy, vz, tof,polx, poly, polz, mech, ntr, weight);
+}
+/***********************************************************/
+
+void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr, AliStack* stack)
+{
+//Shortcut to SetTrack(bla,bla,bla,bla.............)
+   if (p->P() == 0.0)
+    {
+      Error("SetTrack(TParticle*,Int_t&,AliStack*)","Particle has zero momentum");
+      return;
+    }
+   Int_t pdg = p->GetPdgCode();
+   Double_t px = p->Px();
+   Double_t py = p->Py();
+   Double_t pz = p->Pz();
+   Double_t e  = p->Energy();
+
+   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);
+}
+/***********************************************************/
+
+void AliGenHBTosl::Rotate(TVector3& relvector, TVector3& vector)
+{
+//This method rotates vector about the angeles that are needed to rotate 
+//relvector from postion (firstPx,0,0) to its actual positon
+//In other words: To make equations easier
+
+  static TVector3 first;
+  if (GetDebug()) 
+   {
+     first.SetXYZ(relvector.x(),relvector.y(),relvector.z());
+   }
+  
+  Double_t firstPx = TMath::Sqrt( relvector.x()*relvector.x() + 
+                                  relvector.y()*relvector.y() + 
+                                  relvector.z()*relvector.z()   );
+                     
+  Double_t rotAngleZ = -TMath::ATan2(relvector.y(),relvector.x());//calculating rot angles
+  relvector.RotateZ(rotAngleZ);
+  rotAngleZ = -rotAngleZ;
+  Double_t rotAngleY = -TMath::ATan2(relvector.z(),relvector.x());
+  
+  vector.RotateY(rotAngleY);
+  vector.RotateZ(rotAngleZ);
+  
+  if (GetDebug()>5)
+   {
+    TVector3 test(firstPx,0.0,0.0);
+    test.RotateY(rotAngleY);
+    test.RotateZ(rotAngleZ);
+    ::Info("Rotate","Rotation test: px %f %f",first.x(),test.x());
+    ::Info("Rotate","Rotation test: py %f %f",first.y(),test.y());
+    ::Info("Rotate","Rotation test: pz %f %f",first.z(),test.z());
+   }
+}
+/***********************************************************/
+
+Double_t AliGenHBTosl::Rotate(Double_t x,Double_t y,Double_t z)
+{
+//Rotates vector to base where only x - coordinate is no-zero, and returns that 
+
+  Double_t xylength = TMath::Hypot(x,y);
+  Double_t sinphi = -y/xylength;
+  Double_t cosphi = x/xylength;
+  
+  Double_t xprime = cosphi*x - sinphi*y;
+  Double_t yprime = sinphi*x + cosphi*y;
+  
+  TVector3 v(x,y,z);
+  Double_t a1 = -TMath::ATan2(v.Y(),v.X());
+  
+  ::Info("Rotate","Xpr = %f  Ypr = %f",xprime,yprime);
+  ::Info("Rotate","Calc sin = %f, and %f",sinphi,TMath::Sin(a1));
+  ::Info("Rotate","Calc cos = %f, and %f",cosphi,TMath::Cos(a1));
+  
+  
+  Double_t xprimezlength = TMath::Hypot(xprime,z);
+  
+  Double_t sintheta = z/xprimezlength;
+  Double_t costheta = xprime/xprimezlength;
+  
+  
+  Double_t xbis = sintheta*z + costheta*(cosphi*x - sinphi*y);
+  
+  ::Info("Rotate","Calculated rot %f, modulus %f",xbis,TMath::Sqrt(x*x+y*y+z*z));
+  return xbis;
+}
+/***********************************************************/
+
+AliStack* AliGenHBTosl::RotateStack()
+{ 
+//swaps to next stack last goes to first and is reseted
+
+ AliStack* stack;
+ if ( fStackBuffer->GetSize() >= fBufferSize )
+  {
+    stack = (AliStack*)fStackBuffer->Remove(fStackBuffer->Last());
+  }
+ else
+  {
+    stack = new AliStack(fNpart);
+  }
+     
+ fStackBuffer->AddFirst(stack);
+ stack->Reset();
+ return stack;
+}
+/***********************************************************/
+
+Double_t AliGenHBTosl::GetQInvCorrTheorValue(Double_t qinv) const
+{
+//Function (deprecated)
+ static const Double_t factorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
+ return 1.0 + 0.5*TMath::Exp(-qinv*qinv*fQRadius*fQRadius/factorsqrd);
+}
+/***********************************************************/
+
+Double_t AliGenHBTosl::GetQOutQSideQLongCorrTheorValue(Double_t& out, Double_t& side, Double_t& lon) const
+{
+ //Theoretical function. Wa want to get correlation of the shape of this function
+ static const Double_t factorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
+ return 1.0 + 0.7*TMath::Exp(-fQRadius*fQRadius*(out*out+side*side+lon*lon)/factorsqrd);
+}
+/***********************************************************/
+
+Bool_t AliGenHBTosl::CheckParticle(TParticle* p, TParticle* aupair ,AliStack* stack)
+{
+ //Checks if a given particle is falling into signal region with any other particle
+ //already existing on stack
+ return kFALSE;
+  
+ if (fSignalRegion <=0) return kFALSE;
+ for (Int_t i = 0; i < stack->GetNtrack(); i++)
+  {
+    TParticle* part = stack->Particle(i);
+    if (part == aupair) continue;
+    Double_t qout = 10e5;
+    Double_t qside= 10e5;
+    Double_t qlong= 10e5;
+    GetQOutQSideQLong(p,part,qout,qside,qlong);
+    
+    if (TMath::Abs(qout)  < fSignalRegion)
+      if (TMath::Abs(qside) < fSignalRegion)
+       if (TMath::Abs(qlong) < fSignalRegion) 
+         return kTRUE;
+  }
+  return kFALSE; 
+}
+/***********************************************************/
+
+void AliGenHBTosl::SwapGeneratingHistograms()
+{
+  //Checks if it is time to swap signal and background histograms
+  //if yes it swaps them
+  Int_t threshold = fMinFill;
+  for (Int_t k = 1; k<=fQNBins; k++)
+   {
+     for (Int_t j = 1; j<=fQNBins; j++)
+       {  
+         for (Int_t i = 1; i<=fQNBins; i++)
+           {
+             if ( fQSecondBackground->GetBinContent(i,j,k) < threshold) return;
+           }
+       }
+            
+   }
+  
+  
+  Info("SwapGeneratingHistograms","*******************************************");
+  Info("SwapGeneratingHistograms","*******************************************");
+  Info("SwapGeneratingHistograms","*******************************************");
+  Info("SwapGeneratingHistograms","****   SWAPPING HISTOGRAMS             ****");
+  Info("SwapGeneratingHistograms","*******************************************");
+  Info("SwapGeneratingHistograms","*******************************************");
+  Info("SwapGeneratingHistograms","*******************************************");
+  
+  
+  TH3D* h = fQSignal;
+  fQSignal = fQSecondSignal;
+  fQSecondSignal = h;
+  fQSecondSignal->Reset();
+  fQSecondSignal->SetDirectory(0x0);
+  
+  h = fQBackground;
+  fQBackground = fQSecondBackground;
+  fQSecondBackground = h;
+  fQSecondBackground->Reset();
+  fQSecondBackground->SetDirectory(0x0);
+  
+  fSwapped = kTRUE;
+  
+}
diff --git a/EVGEN/AliGenHBTosl.h b/EVGEN/AliGenHBTosl.h
new file mode 100644 (file)
index 0000000..c36ee94
--- /dev/null
@@ -0,0 +1,143 @@
+#ifndef ALIGENHBTOSL_H
+#define ALIGENHBTOSL_H
+//__________________________________________________________
+/////////////////////////////////////////////////////////////
+//                                                         //
+//  class AliGenHBTosl                                     //
+//                                                         //
+//  Genarator simulating particle correlations             //
+//                                                         //
+//  The main idea of the generator is to produce particles //
+//  according to some distribution of two particle         //
+//  property. In HBT they are qout,qsie and qlong.         //
+//  In order to be able to generate signal that produces   //
+//  given two particle correlation background must be      //
+//  known before in order to produce the shape of signal   //
+//  to randomize given distribution from.                  //
+//                                                         //
+//  The generator works as follows:                        //
+//  1. Coarse Background (fQCoarseBackground) is generated //
+//     ade  from the particles                             //
+//     given by the external generator (variable           //
+//     fGenerator) by the mixing technique.                //
+//  2. Coarse signal is prduced by multiplying Coarse      //
+//     background by a required function                   //
+//     See method FillCoarseSignal                         //
+//  3. Signal is randomized out of the coarse signal       //
+//     histogram (two particle property). First particle   //
+//     is taken from the external generator, and the       //
+//     second one is CALCULATED on the basis of the first  //
+//     one and the two particle property (qout,qside,qlong)//
+//     Background is made by the mixing out of the         //
+//     genereted signal events.                            //
+//     This step is cotinued up to the moment signal       //
+//     histogram has enough statistics (data member        //
+//     fMinFill)                                           //
+//     See method StartSignalPass1()                       //
+//  4. chi is calculated for each bin (chiarray variqable) // 
+//     (not the chi2 because sign is important)            //
+//     Two particle prioperty                              //
+//     (qout,qside,qlong) is chosen at the points that     //
+//     chi is the smallest. First particle is taken from   //
+//     the the external generator (fGenerator) and second's /
+//     momenta are caclulated out of their momenta and     //
+//     (qout,qside,qlong). Background is updated           //
+//     continuesely for all the events. This step is       //
+//     continued until stability conditions are fullfiled  //
+//     or maximum number of iteration is reached.          //
+//  5. The same as step 4 but events are stored.           //
+//                                                         //
+////////////////////////////////////////////////////////////
+
+#include "AliGenerator.h"
+
+class TH3D;
+class AliStack;
+class TParticle;
+class TVector;
+
+class AliGenHBTosl: public AliGenerator
+{
+ public:
+   AliGenHBTosl();
+   AliGenHBTosl(Int_t n,Int_t pid = 211);
+   virtual ~AliGenHBTosl();
+
+   void      Init();
+   void      Generate();
+   
+   void      SetGenerator(AliGenerator* gen){fGenerator = gen;}
+   void      SetDebug(Bool_t debug){fDebug = debug;}
+   Bool_t    GetDebug() const {return fDebug;}
+
+   void      Rotate(TVector3& relvector, TVector3& vector);
+   Double_t  Rotate(Double_t x,Double_t y,Double_t z);
+   void      SetSamplePhiRange(Float_t min,Float_t max){fSamplePhiMin = min; fSamplePhiMax = max;}
+
+   Int_t GetThreeD(TParticle* first,TParticle* second, Double_t qout, Double_t qside, Double_t qlong);
+      
+   
+ protected:
+
+   void GetOneD(TParticle* first, TParticle* second,Double_t qinv);
+   
+   void FillCoarse();
+   void FillCoarseSignal();
+   void StartSignal();
+   void StartSignalPass1();
+   void Mix(TList* eventbuffer,TH3D* denominator,TH3D* denominator2);
+   void Mix(AliStack* stack, TH3D* numerator, TH3D* numerator2);
+   Double_t GetQInv(TParticle* f, TParticle* s);
+   void     GetQOutQSideQLong(TParticle* f, TParticle* s,Double_t& out, Double_t& side, Double_t& lon);
+   Double_t GetQInvCorrTheorValue(Double_t qinv) const;
+   Double_t GetQOutQSideQLongCorrTheorValue(Double_t& out, Double_t& side, Double_t& lon) const;
+   
+   Double_t Scale(TH3D* num,TH3D* den);
+   void SetTrack(TParticle* p, Int_t& ntr);
+   void SetTrack(TParticle* p, Int_t& ntr, AliStack* stack);
+   
+   AliStack* RotateStack();
+   void SwapGeneratingHistograms();
+   void     TestCoarseSignal();
+   
+   Bool_t CheckParticle(TParticle* p, TParticle* aupair,AliStack* stack);
+ private:
+  TH3D*    fQCoarseBackground;//Initial Background
+  TH3D*    fQCoarseSignal;//
+  TH3D*    fQSignal;
+  TH3D*    fQBackground;
+
+  TH3D*    fQSecondSignal;
+  TH3D*    fQSecondBackground;
+  
+  Float_t  fQRange;
+  Int_t    fQNBins;
+  AliGenerator* fGenerator;
+  
+  TList*   fStackBuffer;
+  Int_t    fBufferSize;
+  Int_t    fNBinsToScale;
+  Bool_t   fDebug;
+  Bool_t   fSignalShapeCreated;
+  
+  Int_t    fMaxIterations;  //maximal nuber of iterations on startup
+  Float_t  fMaxChiSquereChange;//value of ChiSqr change in %, when sturtup process in stable
+  Float_t  fMaxChiSquerePerNDF;//value of the chi2 where generating histograms are considered as good
+  
+  
+  Double_t fQRadius;
+  
+  Int_t    fPID;
+  //we limit mixing to some finit phi range to make it faster
+  Float_t  fSamplePhiMin;//
+  Float_t  fSamplePhiMax;//
+  
+  Float_t  fSignalRegion;
+  
+  Int_t    fMinFill;
+  
+  Bool_t   fSwapped;
+  
+  ClassDef(AliGenHBTosl,1)
+};
+#endif
index ce645e7..ddab2da 100644 (file)
@@ -51,6 +51,7 @@
 #pragma link C++ class  AliGenSlowNucleons+;
 #pragma link C++ class  AliSlowNucleonModel+;
 #pragma link C++ class  AliSlowNucleonModelExp+;
+#pragma link C++ class  AliGenHBTosl+;
 #pragma link C++ class  AliGenMUONCocktail+;
 #endif
 
index ba9d705..fb82652 100644 (file)
@@ -1,12 +1,12 @@
 FSRCS= PartEloss/integralgauss.f  PartEloss/read.f PartEloss/readca3.f
 SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
-               AliGenCocktail.cxx AliGenCocktailAfterBurner.cxx AliGenFLUKAsource.cxx \
-               AliGenLib.cxx \
+                AliGenCocktail.cxx AliGenCocktailAfterBurner.cxx AliGenFLUKAsource.cxx \
+                AliGenLib.cxx \
                 AliGenMUONlib.cxx AliGenParam.cxx\
                 AliDimuCombinator.cxx AliGenHalo.cxx \
                 AliGenHaloProtvino.cxx \
                 AliGenExtFile.cxx AliGenScan.cxx AliGenPHOSlib.cxx \
-               AliGenDoubleScan.cxx AliGenCocktailEntry.cxx \
+                AliGenDoubleScan.cxx AliGenCocktailEntry.cxx \
                AliGenGSIlib.cxx AliGenPMDlib.cxx\
                AliGenMC.cxx\
                AliGenReader.cxx AliGenReaderCwn.cxx AliGenReaderTreeK.cxx \
@@ -18,7 +18,7 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                AliStructFuncType.cxx AliGenSlowNucleons.cxx \
                AliGenGeVSimEventHeader.cxx\
                AliSlowNucleonModel.cxx AliSlowNucleonModelExp.cxx \
-               AliGenMUONCocktail.cxx
+               AliGenMUONCocktail.cxx AliGenHBTosl.cxx 
 
 # Headerfiles for this particular package (Path respect to own directory)
 HDRS= $(SRCS:.cxx=.h)