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