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