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