1 #include "AliGenHBTosl.h"
2 //__________________________________________________________
3 /////////////////////////////////////////////////////////////
5 // class AliGenHBTosl //
7 // Genarator simulating particle correlations //
9 // The main idea of the generator is to produce particles //
10 // according to some distribution of two particle //
11 // property. In HBT they are qout,qsie and qlong. //
12 // In order to be able to generate signal that produces //
13 // given two particle correlation background must be //
14 // known before in order to produce the shape of signal //
15 // to randomize given distribution from. //
17 // The generator works as follows: //
18 // 1. Coarse Background (fQCoarseBackground) is generated //
19 // ade from the particles //
20 // given by the external generator (variable //
21 // fGenerator) by the mixing technique. //
22 // 2. Coarse signal is prduced by multiplying Coarse //
23 // background by a required function //
24 // See method FillCoarseSignal //
25 // 3. Signal is randomized out of the coarse signal //
26 // histogram (two particle property). First particle //
27 // is taken from the external generator, and the //
28 // second one is CALCULATED on the basis of the first //
29 // one and the two particle property (qout,qside,qlong)//
30 // Background is made by the mixing out of the //
31 // genereted signal events. //
32 // This step is cotinued up to the moment signal //
33 // histogram has enough statistics (data member //
35 // See method StartSignalPass1() //
36 // 4. chi is calculated for each bin (chiarray variqable) //
37 // (not the chi2 because sign is important) //
38 // Two particle prioperty //
39 // (qout,qside,qlong) is chosen at the points that //
40 // chi is the smallest. First particle is taken from //
41 // the the external generator (fGenerator) and second's /
42 // momenta are caclulated out of their momenta and //
43 // (qout,qside,qlong). Background is updated //
44 // continuesely for all the events. This step is //
45 // continued until stability conditions are fullfiled //
46 // or maximum number of iteration is reached. //
47 // 5. The same as step 4 but events are stored. //
49 ////////////////////////////////////////////////////////////
57 #include <TParticle.h>
61 #include <TStopwatch.h>
64 #include "AliGenCocktailAfterBurner.h"
65 #include "AliGeVSimParticle.h"
66 #include "AliGenGeVSim.h"
67 #include "AliGenHIJINGpara.h"
70 /***********************************************************/
71 ClassImp(AliGenHBTosl)
73 AliGenHBTosl::AliGenHBTosl():
75 fQCoarseBackground(0x0),
80 fQSecondBackground(0x0),
85 fNBinsToScale(Int_t(fQNBins*0.1)),
88 fMaxChiSquereChange(0.01),
89 fMaxChiSquerePerNDF(1.5),
93 fSamplePhiMax(TMath::TwoPi()+0.01),
100 /***********************************************************/
102 AliGenHBTosl::AliGenHBTosl(Int_t n, Int_t pid):
104 fQCoarseBackground(0x0),
109 fQSecondBackground(0x0),
114 fNBinsToScale(Int_t(fQNBins*0.1)),
116 fSignalShapeCreated(kFALSE),
117 fMaxIterations(1000),
118 fMaxChiSquereChange(0.01),
119 fMaxChiSquerePerNDF(1.5),
122 fSamplePhiMin(-0.01),
123 fSamplePhiMax(TMath::TwoPi()+0.01),
128 //default constructor
130 /***********************************************************/
132 AliGenHBTosl::~AliGenHBTosl()
135 delete fQCoarseSignal;
136 delete fQCoarseBackground;
140 delete fQSecondSignal;
141 delete fQSecondBackground;
144 /***********************************************************/
146 void AliGenHBTosl::Init()
148 //Initializes generator
149 if (fGenerator == 0x0)
152 AliGenHIJINGpara* bkggen = new AliGenHIJINGpara(fNpart*4);
156 AliGenGeVSim * gevsim = new AliGenGeVSim(0.0);
157 AliGeVSimParticle* kplus = new AliGeVSimParticle(fPID,1,fNpart, 0.17, 0.9);
158 gevsim->AddParticleType(kplus);
165 AliMevSimConfig *c = new AliMevSimConfig(1);
166 c->SetRectPlane(1); // reaction plane control, model 4
169 AliGenMevSim *mevsim = new AliGenMevSim(c);
170 mevsim->SetPtRange(0.001, 3);
171 mevsim->SetMomentumRange(0.1, 3);
172 mevsim->SetTrackingFlag(0);
173 mevsim->SetOrigin(0.0, 0.0, 0.0);
174 mevsim->SetSigma(0.0, 0.0, 0.0);
175 AliMevSimParticle *kplus = new AliMevSimParticle(kKPlus, fNpart, 0, 0.25, 0.0, 2, 0.15, 0.0, 0.0 );
176 mevsim->AddParticleType(kplus);
180 fGenerator->SetOrigin(fOrigin[0],fOrigin[1],fOrigin[2]);
181 static const Double_t kDegToRadCF = 180./TMath::Pi();
182 fGenerator->SetMomentumRange(fPtMin,fPtMax);
183 fGenerator->SetPhiRange(kDegToRadCF*fPhiMin,kDegToRadCF*fPhiMax);
184 fGenerator->SetYRange(fYMin,fYMax);
185 fGenerator->SetThetaRange(kDegToRadCF*fThetaMin,kDegToRadCF*fThetaMax);
190 // fQCoarseBackground = new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
191 // fQCoarseSignal = new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
192 // fQCoarseBackground->Sumw2();
193 // fQCoarseSignal->Sumw2();
195 fQSignal = new TH3D("fQSignal1","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
196 fQBackground = new TH3D("fQBackground1","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
198 fQSecondSignal = new TH3D("fQSignal2","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
199 fQSecondBackground = new TH3D("fQBackground2","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
202 fQBackground->Sumw2();
203 fQSecondSignal->Sumw2();
204 fQSecondBackground->Sumw2();
206 fLogFile = new ofstream("BadEvent",ios::out);
209 /***********************************************************/
211 void AliGenHBTosl::Generate()
215 ofstream& logfile = *fLogFile;
216 logfile<<"Generate"<<"Attempts to generate "<<fNpart<<" particles.";
219 if (fStackBuffer == 0x0) fStackBuffer = new TList();
220 //Here is initialization level
221 if (fSignalShapeCreated == kFALSE)
223 TH3D *hs = 0x0, *hb = 0x0;
226 file = TFile::Open("QTSignal.root");
229 hs = (TH3D*)file->Get("fQSignal1");
230 if (hs) hs->SetDirectory(0x0);
234 file = TFile::Open("QTBackground.root");
237 hb = (TH3D*)file->Get("fQBackground1");
238 if (hb) hb->SetDirectory(0x0);
244 Info("Generate","**********************************");
245 Info("Generate","Found starting histograms in files");
246 Info("Generate","**********************************");
254 TH3D *cs = 0x0, *cb = 0x0;
255 file = TFile::Open("QTCoarseBackground.root");
258 cb = (TH3D*)file->Get("fQCoarseBackground");
259 if (cb) cb->SetDirectory(0x0);
263 file = TFile::Open("QTCoarseSignal.root");
266 cs = (TH3D*)file->Get("fQCoarseSignal");
267 if (cs) cs->SetDirectory(0x0);
274 Info("Generate","Got Coarse signal and bkg from files");
275 delete fQCoarseBackground;
276 delete fQCoarseSignal;
278 fQCoarseBackground = cb;
284 Info("Generate","Got Coarse bkg from file");
285 delete fQCoarseBackground;
286 fQCoarseBackground = cb;
290 fQCoarseBackground = new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
291 fQCoarseBackground->Sumw2();
293 FillCoarse(); //create coarse background - just to know the spectrum
296 fQCoarseSignal = new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
297 fQCoarseSignal->Sumw2();
298 FillCoarseSignal();//create first coarse signal by brutal multplication coarse background and required function shape
301 StartSignal(); //Initilizes the stack that is used for generation
303 fSignalShapeCreated = kTRUE;
306 AliStack* stack = RotateStack();
308 AliStack* genstack = fGenerator->GetStack();
311 genstack = new AliStack(fNpart);
312 fGenerator->SetStack(genstack);
319 fGenerator->Generate();
320 Int_t j = 0, ntr = 0;
321 if ( genstack->GetNtrack() < fNpart/2)
323 Warning("Generate","************************************************************");
324 Warning("Generate","Generator generated (%d) less particles then expected (%d).",
325 stack->GetNtrack(),fNpart/2);
326 Warning("Generate","************************************************************");
329 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
330 work->SetDirectory(0x0);
333 Double_t*** chiarray = new Double_t** [fQNBins+1];
334 Double_t*** sigarray = new Double_t** [fQNBins+1];
336 for (Int_t i = 1; i<=fQNBins; i++)
338 chiarray[i] = new Double_t* [fQNBins+1];
339 sigarray[i] = new Double_t* [fQNBins+1];
341 for (Int_t k = 1; k<=fQNBins; k++)
343 chiarray[i][k] = new Double_t [fQNBins+1];
344 sigarray[i][k] = new Double_t [fQNBins+1];
349 Double_t scale = Scale(fQSignal,fQBackground);
350 work->Divide(fQSignal,fQBackground,scale);
352 Double_t binwdh = work->GetBinWidth(1)/2.;
354 for (Int_t k = 1; k<=fQNBins; k++)
356 Double_t z = work->GetZaxis()->GetBinCenter(k);
357 for (Int_t j = 1; j<=fQNBins; j++)
359 Double_t y = work->GetYaxis()->GetBinCenter(j);
360 for (Int_t i = 1; i<=fQNBins; i++)
362 sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
363 Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qinv)
364 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
365 Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
366 chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
374 Int_t middlebin = fQNBins/2;
376 for (Int_t k = middlebin-5; k < middlebin+5; k++)
378 Double_t tx = work->GetXaxis()->GetBinCenter(30);
379 Double_t ty = work->GetYaxis()->GetBinCenter(30);
380 Double_t tz = work->GetZaxis()->GetBinCenter(k);
381 sprintf(msg,"% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
386 for (Int_t k = middlebin-5; k < middlebin+5; k++)
388 sprintf(msg,"% 6.5f ",work->GetBinContent(30,30,k));
393 for (Int_t k = middlebin-5; k < middlebin+5; k++)
395 sprintf(msg,"% 6.5f ",chiarray[30][30][k]);
400 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
401 TParticle* second = &particle;
403 Bool_t shortloop = kTRUE;
404 Int_t sc = 0;//security check against infinite loop
406 while ( (ntr+1) < fNpart)
428 loopmax = fQNBins/2+fQNBins/4;
429 loopmin = fQNBins/2-fQNBins/4;
433 for (Int_t k = loopmin; k <=loopmax; k++ )
435 qlong = work->GetZaxis()->GetBinCenter(k);
436 for (Int_t j = loopmin; j<=loopmax; j++)
438 qside = work->GetYaxis()->GetBinCenter(j);
439 for (Int_t i = loopmin; i<=loopmax; i++)
441 qout = work->GetXaxis()->GetBinCenter(i);
442 if (chiarray[xmax][ymax][zmax] < chiarray[i][j][k])
449 // Double_t qdist = TMath::Sqrt(qout*qout + qside*qside + qlong*qlong);
451 // Double_t fact = chiarray[i][j][k];//chiarray is chi2
452 // if (fact > work->GetBinError(i,j,k))//if differece between what we want and
453 // { //what we have is bigger than stat. error
454 // xmax = i; //we force to fill that bin
462 Double_t qlongc = work->GetZaxis()->GetBinCenter(zmax);
463 Double_t qsidec = work->GetYaxis()->GetBinCenter(ymax);
464 Double_t qoutc = work->GetXaxis()->GetBinCenter(xmax);
467 sprintf(msg,"Generate Fill bin chi2(%d,%d,%d)=%f",xmax,ymax,zmax,chiarray[xmax][ymax][zmax]);
470 qout = gRandom->Uniform(qoutc -binwdh, qoutc +binwdh);
471 qside = gRandom->Uniform(qsidec-binwdh, qsidec+binwdh);
472 qlong = gRandom->Uniform(qlongc-binwdh, qlongc+binwdh);
474 TParticle* first = 0;
475 while (j < genstack->GetNtrack())
477 TParticle* tmpp = genstack->Particle(j++);
478 if (tmpp->GetPdgCode() == fPID)
480 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
490 if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
494 //Here put the check if this particle do not fall into signal region with other partticle
496 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
499 //Info("StartSignal","Can not find momenta for this OSL and particle");
502 //in case this particle is falling into signal area with another
503 //particle we take a another pair
504 //it can intruduce artificial correlations
505 Bool_t checkresult = CheckParticle(second,first,stack);
506 if ( checkresult && (sc < 10) )
513 //Put on output stack
515 SetTrack(second,ntr);
517 //Put on internal stack
519 SetTrack(first,etmp,stack);
520 SetTrack(second,etmp,stack);
522 Double_t y = GetQOutQSideQLongCorrTheorValue(qoutc,qsidec,qlongc);
524 sigarray[xmax][ymax][zmax] ++;
525 chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
526 chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
530 Mix(fStackBuffer,fQBackground,fQSecondSignal); //upgrate background
531 Mix(stack,fQSignal,fQSecondBackground); //upgrate signal
535 for (Int_t i = 1; i<=fQNBins; i++)
537 for (Int_t k = 1; k<=fQNBins; k++)
539 delete [] chiarray[i][k];
540 delete [] sigarray[i][k];
542 delete [] chiarray[i];
543 delete [] sigarray[i];
548 /***********************************************************/
550 void AliGenHBTosl::GetOneD(TParticle* first, TParticle* second,Double_t qinv)
552 //deprecated method that caclulates momenta of the second particle
553 // out of qinv and the first particle
554 //first particle is rotated that only X is non-zero
557 Double_t m = first->GetMass();
558 Double_t msqrd = m*m;
559 Double_t fourmassSquered = 4.*msqrd;
561 //Condition that R must fullfill to be possible to have qinv less smaller then randomized
562 // Double_t rRange = qinv*TMath::Sqrt(qinv*qinv + fourmassSquered)/fourmassSquered;
563 // Double_t r = gRandom->Uniform(rRange);
565 Double_t r = gRandom->Uniform(qinv);
566 Double_t phi = gRandom->Uniform(TMath::TwoPi());
568 Double_t firstPx = first->P();//first particle is rotated that only X is non-zero thus P==Px
569 Double_t px = 2.*msqrd*firstPx + firstPx*qinv*qinv;
570 Double_t temp = qinv*qinv*qinv*qinv + fourmassSquered * (qinv*qinv - r*r );
573 Error("GetOneD","temp is less then 0: %f",temp);
576 temp = temp*(msqrd+firstPx*firstPx);
578 px = (px - TMath::Sqrt(temp))/(2.*msqrd);
580 Double_t py = r*TMath::Sin(phi);
581 Double_t pz = r*TMath::Cos(phi);
583 TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
584 TVector3 vector(px,py,pz);
585 Rotate(firstpvector,vector);
587 Double_t e = TMath::Sqrt(msqrd + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
588 second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
589 // 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);
590 // TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
591 // px, py, pz, e, vx, vy, vz, tof);
593 if (GetDebug()) Info("GetOneD","Randomized qinv = %f, obtained = %f",qinv,GetQInv(first,second));
596 /***********************************************************/
598 Int_t AliGenHBTosl::GetThreeD(TParticle* first,TParticle* second, Double_t qout, Double_t qside, Double_t qlong)
600 //deprecated method that caclulates momenta of the second particle
601 //out of qout qside and qlong and the first particle
602 Double_t m = first->GetMass();
605 Double_t px = first->P();//first particle is rotated that only X is non-zero thus P==Px
606 Double_t px2 = px*px;
609 Double_t qout2 = qout*qout;
610 Double_t qside2 = qside*qside;
611 Double_t qlong2 = qlong*qlong;
614 Double_t util1 = 4.*px2 - qside2;//4*P1x^2 - Y^2
617 Info("GetThreeD","4.*px2* - qside2 is negative px: %f, qside: %f",px,qside);
620 Double_t util2 = TMath::Sqrt(px2*qout2*util1);
623 Double_t p2x,p2y,p2z;
625 // if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))
629 Double_t tmp = px*(2.*px2 - qside2);
634 tmp = qout - TMath::Sqrt(util1);
635 p2y = - (tmp*qside)/(2.*px);
638 tmp = 4.*m2 + 2.*qout2+qlong2;
640 tmp -= 2.*util2;//!!!
644 Double_t m2px2 = m2+px2;
645 Double_t tmp2 = m2px2*tmp;
647 tmp = 4.*(m2px2+qout2) + qlong2;
650 tmp *= 4.*(m2px2) + qlong2;
651 tmp *= qlong2*qlong2;
656 Error("","Argument of sqrt is negative");
660 tmp2 += TMath::Sqrt(tmp);
662 tmp = 8.0*px*m2px2*m2px2;
663 p2z = -TMath::Sqrt(tmp2/tmp);
664 if (qlong < 0) p2z = -p2z;
669 Double_t tmp = px*(2.*px2 - qside2);
674 tmp = qout - TMath::Sqrt(util1);
675 p2y = - (tmp*qside)/(2.*px);
678 tmp = 4.*m2 + 2.*qout2+qlong2;
680 tmp += 2.*util2;//!!!
684 Double_t m2px2 = m2+px2;
685 Double_t tmp2 = m2px2*tmp;
687 tmp = 4.*(m2px2+qout2) + qlong2;
690 tmp *= 4.*(m2px2) + qlong2;
691 tmp *= qlong2*qlong2;
696 Error("","Argument of sqrt is negative");
700 tmp2 += TMath::Sqrt(tmp);
702 tmp = 8.0*px*m2px2*m2px2;
703 p2z = -TMath::Sqrt(tmp2/tmp);
704 if (qlong < 0) p2z = -p2z;
707 // if ( (qside >= 0) && (qout >= 0) && (qlong >= 0)) p2z = -p2z;
709 TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
710 TVector3 vector(p2x,p2y,p2z);
711 Rotate(firstpvector,vector);
713 Double_t e = TMath::Sqrt(m2 + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
714 second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
717 if ( GetDebug() > 3 )
719 e=TMath::Sqrt(m2+px*px);
720 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);
722 e = TMath::Sqrt(m2 + p2x*p2x + p2y*p2y + p2z*p2z);
723 TParticle* s = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, p2x, p2y, p2z, e, 0.0, 0.0, 0.0, 0.0);
726 GetQOutQSideQLong(f,s,qo, qs, ql);
728 Info("GetThreeD","TEST");
731 Info("GetThreeD","Required %f %f %f",qout,qside,qlong);
732 Info("GetThreeD","Got %f %f %f",qo,qs,ql);
736 Info("GetThreeD","!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
741 /***********************************************************/
743 void AliGenHBTosl::StartSignal()
745 //Starts the signal histograms
746 ofstream& logfile = *fLogFile;
748 logfile<<"************************************************"<<endl;
749 logfile<<"************************************************"<<endl;
750 logfile<<" StartSignal "<<endl;
751 logfile<<"************************************************"<<endl;
752 logfile<<"************************************************"<<endl;
758 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
759 TParticle* second = &particle;
761 TIter next(fStackBuffer);
762 while(( stack=(AliStack*)next() ))
767 AliStack* genstack = fGenerator->GetStack();
770 genstack = new AliStack(fNpart);
771 fGenerator->SetStack(genstack);
779 //We alread have detailed histograms and we do not need Coarse anymore
780 delete fQCoarseSignal;
781 delete fQCoarseBackground;
782 fQCoarseSignal = 0x0;
783 fQCoarseBackground = 0x0;
786 const Double_t kNDF = fQNBins*fQNBins*fQNBins;
788 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
790 work->SetDirectory(0x0);
793 Double_t binwdh = work->GetBinWidth(1)/2.;
795 Double_t*** chiarray = new Double_t** [fQNBins+1];
796 Double_t*** sigarray = new Double_t** [fQNBins+1];
798 for (Int_t i = 1; i<=fQNBins; i++)
800 chiarray[i] = new Double_t* [fQNBins+1];
801 sigarray[i] = new Double_t* [fQNBins+1];
803 for (Int_t k = 1; k<=fQNBins; k++)
805 chiarray[i][k] = new Double_t [fQNBins+1];
806 sigarray[i][k] = new Double_t [fQNBins+1];
811 Float_t chisqrchange = fMaxChiSquereChange + 1.;
812 Float_t chisqrPerDF = fMaxChiSquerePerNDF;
813 Float_t chisqrold = 0.0;
816 Int_t niterations = 1;
817 Int_t rotaxisorder = 1;//defines order of looping over 3D histo (X,Y,Z or Y,Z,X or Z,X,Y)
820 Bool_t shortloop = kTRUE;
821 TCanvas* c1 = new TCanvas();
825 Info("StartSignal","\n\n\n\nSecond Pass\n\n\n\n");
827 while ( ( (chisqrPerDF > fMaxChiSquereChange) || flag) && (niterations++ < fMaxIterations) )
830 logfile<<"StartSignal\n";
831 logfile<<" Row 1 Theory, 2 current value, 3 Chi2 \n";
833 Double_t chisqrnew = 0.0;
836 Double_t scale = Scale(fQSignal,fQBackground);
837 work->Divide(fQSignal,fQBackground,scale);
839 if ( (counter%100) == 0)
843 sprintf(buff,"QTWorkPass2.%3d.root",counter);
844 TFile* file = TFile::Open(buff,"update");
846 work->SetDirectory(0x0);
850 sprintf(buff,"QTBackgroundPass2.%3d.root",counter);
851 file = TFile::Open(buff,"update");
852 fQBackground->Write();
853 fQBackground->SetDirectory(0x0);
857 sprintf(buff,"QTSignalPass2.%3d.root",counter);
858 file = TFile::Open(buff,"update");
860 fQSignal->SetDirectory(0x0);
866 Int_t novertresh = 0;
867 for (Int_t k = 1; k<=fQNBins; k++)
869 Double_t z = work->GetZaxis()->GetBinCenter(k);
870 for (Int_t j = 1; j<=fQNBins; j++)
872 Double_t y = work->GetYaxis()->GetBinCenter(j);
873 for (Int_t i = 1; i<=fQNBins; i++)
875 Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qout)
876 sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
877 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
878 Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
879 chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
880 Double_t be = work->GetBinError(i,j,k);
881 chisqrnew += diff*diff/(be*be);//add up chisq
883 //even if algorithm is stable (chi sqr change less then threshold)
884 //and any bin differs more then 5% from expected value we continue
885 Double_t fact = diff;
886 if (TMath::Abs(fact) > 0.1)
900 for (Int_t k = 25; k < 36; k++)
902 Double_t tx = work->GetXaxis()->GetBinCenter(30);
903 Double_t ty = work->GetYaxis()->GetBinCenter(30);
904 Double_t tz = work->GetZaxis()->GetBinCenter(k);
905 sprintf(msg,"% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
910 for (Int_t k = 25; k < 36; k++)
912 sprintf(msg,"%6.5f ",work->GetBinContent(30,30,k));
917 for (Int_t k = 25; k < 36; k++)
919 sprintf(msg,"% 6.5f ",chiarray[30][30][k]);
924 chisqrchange = TMath::Abs(chisqrnew - chisqrold)/chisqrnew;
925 chisqrold = chisqrnew;
927 chisqrPerDF = chisqrnew/kNDF;
929 Info("StartSignal","Iteration %d Chi-sq change %f%%",niterations,chisqrchange*100.0);
930 Info("StartSignal","ChiSq = %f, NDF = %f, ChiSq/NDF = %f",chisqrnew, kNDF, chisqrPerDF );
931 Info("StartSignal","novertresh = %d",novertresh);
934 stack = RotateStack();
936 fGenerator->Generate();
937 Int_t ninputparticle = 0, ntr = 0;
938 if ( genstack->GetNtrack() < fNpart/2)
940 Warning("StartSignal","**********************************");
941 Warning("StartSignal","Generator generated (%d) less particles then expected (%d).",
942 genstack->GetNtrack(),fNpart/2);
943 Warning("StartSignal","**********************************");
946 Int_t sc = 0; //security check against infinite loop
947 while ( (ntr+1) < fNpart)//ntr is number of track generated up to now
962 switch (rotaxisorder)
981 if (rotaxisorder > 3) rotaxisorder = 1;
995 // Bool_t force = kFALSE;
996 for ( k = 1; k <=nrange;k++ )
998 for ( j = 1; j<=nrange; j++)
1000 for ( i = 1; i<=nrange; i++)
1002 if ( (chiarray[*cx][*cy][*cz]) > (chiarray[xmax][ymax][zmax]) )
1009 // Double_t fact = chiarray[*cx][*cy][*cz];//chiarray is chi2*qinv
1010 // if (fact > work->GetBinError(*cx,*cy,*cz))//if differece between what we want and
1011 // { //what we have is bigger than stat. error
1012 // //we force to fill that bin
1013 // qout = work->GetXaxis()->GetBinCenter(*cx);
1014 // qside = work->GetYaxis()->GetBinCenter(*cy);
1015 // qlong = work->GetZaxis()->GetBinCenter(*cz);
1017 // 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",
1018 // *cx,*cy,*cz,i,j,k,fact,qout,qside,qlong,
1019 // (Int_t)sigarray[*cx][*cy][*cz],(Int_t)fQBackground->GetBinContent(*cx,*cy,*cz),work->GetBinError(*cx,*cy,*cz));
1025 // if (force) break;
1027 // if (force) break;
1031 qout = work->GetXaxis()->GetBinCenter(xmax);
1032 qside = work->GetYaxis()->GetBinCenter(ymax);
1033 qlong = work->GetZaxis()->GetBinCenter(zmax);
1035 // Info("StartSignal"," bin: (%d,%d,%d) chiarray: %f \nq(o,s,l): (%f,%f,%f) signal: %d background: %d binerror: %f",
1036 // xmax,ymax,zmax,chiarray[xmax][ymax][zmax],qout,qside,qlong,
1037 // (Int_t)sigarray[xmax][ymax][zmax],
1038 // (Int_t)fQBackground->GetBinContent(xmax,ymax,zmax),
1039 // work->GetBinError(xmax,ymax,zmax));
1041 qout = gRandom->Uniform(qout-binwdh,qout+binwdh);
1042 qside = gRandom->Uniform(qside-binwdh,qside+binwdh);
1043 qlong = gRandom->Uniform(qlong-binwdh,qlong+binwdh);
1045 TParticle* first = 0;
1046 while (ninputparticle < genstack->GetNtrack())
1048 TParticle* tmpp = genstack->Particle(ninputparticle++);
1049 if (tmpp->GetPdgCode() == fPID)
1051 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1061 if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
1065 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1068 Info("StartSignal","Can not find momenta for this OSL and particle OSL = %f %f %f",qout,qside,qlong);
1074 //in case this particle is falling into signal area with another
1075 //particle we take a another pair
1076 //it can intruduce artificial correlations
1077 Bool_t checkresult = CheckParticle(second,first,stack);
1078 if ( checkresult && (sc < 10) )
1085 //Put on output stack
1086 SetTrack(first,ntr,stack);
1087 SetTrack(second,ntr,stack);
1089 Double_t y = GetQOutQSideQLongCorrTheorValue(qout,qside,qlong);
1091 sigarray[xmax][ymax][zmax] ++;
1092 chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
1093 chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
1096 Info("StartSignal","Mixing background...");
1097 Mix(fStackBuffer,fQBackground,fQSecondBackground); //upgrate background
1098 Info("StartSignal","Mixing signal...");
1099 Mix(stack,fQSignal,fQSecondSignal); //upgrate background
1102 // if ( (chisqrPerDF < 2.0) && (fSwapped == kFALSE) )
1104 // SwapGeneratingHistograms();
1108 TFile* filef = TFile::Open("QTBackground.root","recreate");
1109 fQBackground->Write();
1110 fQBackground->SetDirectory(0x0);
1114 filef = TFile::Open("QTSignal.root","recreate");
1116 fQSignal->SetDirectory(0x0);
1124 for (Int_t i = 1; i<=fQNBins; i++)
1126 for (Int_t k = 1; k<=fQNBins; k++)
1128 delete [] chiarray[i][k];
1129 delete [] sigarray[i][k];
1131 delete [] chiarray[i];
1132 delete [] sigarray[i];
1138 /***********************************************************/
1140 void AliGenHBTosl::StartSignalPass1()
1142 //This method makes first part of the initialization of working histograms
1143 //It randomizes qout, qside and qlong from the coarse signal histogram
1145 Bool_t flag = kTRUE;
1146 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
1147 TParticle* second = &particle;
1152 Info("StartSignalPass1","\n\nFirst Pass\n\n");
1156 Info("StartSignalPass1","NextEvent");
1157 AliStack* stack = RotateStack();
1158 AliStack* genstack = fGenerator->GetStack();
1160 fGenerator->Generate();
1161 Int_t j = 0, ntr = 0;
1162 if ( genstack->GetNtrack() < fNpart/2)
1164 Warning("StartSignalPass1","**********************************");
1165 Warning("StartSignalPass1","Generator generated (%d) less particles then expected (%d).",
1166 genstack->GetNtrack(),fNpart/2);
1167 Warning("StartSignalPass1","**********************************");
1170 Int_t sc = 0;//security check against infinite loop
1171 while ((ntr+1)<fNpart)
1174 // Info("StartSignal","Number of track on output stack: = %d", ntr);
1175 // Info("StartSignal","Number of track on input stack: = %d\n", j);
1177 TParticle* first = 0;
1178 while (j < genstack->GetNtrack())
1180 TParticle* tmpp = genstack->Particle(j++);
1181 if (tmpp->GetPdgCode() == fPID)
1183 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1190 Info("StartSignalPass1","Particle did not pass the safety check 1");
1198 if ( fDebug > 2 ) Info("StartSignalPass1","No more particles of that type");
1203 SetTrack(first,ntr,stack);
1205 fQCoarseSignal->GetRandom3(qout,qside,qlong);
1207 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1210 //Info("StartSignal","Can not find momenta for this OSL and particle");
1213 //in case this particle is falling into signal area with another
1214 //particle we take a another pair
1215 //it can intruduce artificial correlations
1216 Bool_t checkresult = CheckParticle(second,first,stack);
1217 if ( checkresult && (sc < 10) )
1220 Info("StartSignalPass1","Particle did not pass the safety check 2");
1227 SetTrack(second,ntr,stack);
1230 Mix(stack,fQSignal,fQSecondSignal);
1231 Mix(fStackBuffer,fQBackground,fQSecondBackground);
1235 for (Int_t k = 1; k<=fQNBins; k++)
1237 for (Int_t j = 1; j<=fQNBins; j++)
1239 for (Int_t i = 1; i<=fQNBins; i++)
1241 if ( (fQBackground->GetBinContent(i,j,k) < fMinFill) )
1243 //(fQSignal->GetBinContent(i,j,k) < fMinFill) ||
1244 Info("StartSignalPass1","bin (%d,%d,%d): signal=%f background=%f",i,j,k,
1245 fQSignal->GetBinContent(i,j,k),fQBackground->GetBinContent(i,j,k));
1246 flag = kTRUE;//continue while
1247 break;//breakes for not while
1250 if (flag == kTRUE) break;
1252 if (flag == kTRUE) break;
1259 /***********************************************************/
1261 void AliGenHBTosl::FillCoarseSignal()
1263 //Makes coarse signal by multiplying the coarse background and required function
1264 Info("FillCoarseSignal","START");
1265 for (Int_t k = 1; k <=fQNBins ;k++ )
1267 Double_t z = fQCoarseBackground->GetZaxis()->GetBinCenter(k);
1268 for (Int_t j = 1; j <=fQNBins; j++)
1270 Double_t y = fQCoarseBackground->GetYaxis()->GetBinCenter(j);
1271 for (Int_t i = 1; i <=fQNBins; i++)
1273 Double_t x = fQCoarseBackground->GetXaxis()->GetBinCenter(i);
1274 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);
1275 Info("FillCoarseSignal","Bin (%d,%d,%d): osl(%f,%f,%f)=%f",i,j,k,x,y,z,v);
1276 fQCoarseSignal->SetBinContent(i,j,k,v*fQCoarseBackground->GetBinContent(i,j,k));
1284 Info("FillCoarseSignal","DONE");
1286 /***********************************************************/
1288 void AliGenHBTosl::FillCoarse()
1290 //creates the statistical background histogram on the base of input from
1292 Info("FillCoarse","START");
1298 TH3D tmph("tmph","tmph",2,0,1,2,0,1,2,0,1);
1303 // if (niter > 20) break;
1305 cout<<niter++<<" bincont "<<fQCoarseBackground->GetBinContent(30,30,28)
1306 <<" "<<fQCoarseBackground->GetBinContent(30,30,29)
1307 <<" "<<fQCoarseBackground->GetBinContent(30,30,30)
1308 <<" "<<fQCoarseBackground->GetBinContent(30,30,31)
1309 <<" "<<fQCoarseBackground->GetBinContent(30,30,32)
1313 stack = RotateStack();
1314 fGenerator->SetStack(stack);
1316 fGenerator->Generate();
1318 Mix(fStackBuffer,fQCoarseBackground,&tmph);
1322 Info("FillCoarse","fMinFill = %d",fMinFill);
1323 for (Int_t k = 1; k<=fQNBins; k++)
1325 for (Int_t j = 1; j<=fQNBins; j++)
1327 for (Int_t i = 1; i<=fQNBins; i++)
1329 if ( fQCoarseBackground->GetBinContent(i,j,k) < fMinFill )
1332 Info("FillCoarse","bin (%d,%d,%d)=%f",i,j,k,fQCoarseBackground->GetBinContent(i,j,k));
1344 fGenerator->SetStack(0x0);
1345 Info("FillCoarse","DONE");
1348 /***********************************************************/
1350 void AliGenHBTosl::Mix(TList* eventbuffer,TH3D* denominator,TH3D* denominator2)
1352 //Fills denominators
1353 //Mixes events stored in the eventbuffer and fills the background histograms
1354 static TStopwatch stoper;
1356 if (eventbuffer == 0x0)
1358 Error("Mix","Buffer List is null.");
1362 if (denominator == 0x0)
1364 Error("Mix","Denominator histogram is null.");
1368 if (denominator2 == 0x0)
1370 Error("Mix","Denominator2 histogram is null.");
1374 Info("Mix","%s",denominator->GetName());
1377 TIter next(eventbuffer);
1378 AliStack* firstevent;
1379 AliStack* secondevent = 0x0;
1381 while(( firstevent=(AliStack*)next() ))
1383 if (secondevent == 0x0)
1385 secondevent = firstevent;
1388 // Info("Mix","Mixing %#x with %#x",firstevent,secondevent);
1389 for(Int_t j = 0; j < firstevent->GetNtrack(); j++ )
1391 TParticle* firstpart = firstevent->Particle(j);
1393 Float_t phi = firstpart->Phi();
1394 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1396 // Info("Mix","Mixing %d phi %f min %f max %f",j,phi,fSamplePhiMin,fSamplePhiMax);
1398 for(Int_t i = 0; i < secondevent->GetNtrack(); i++ )
1400 TParticle* secondpart = secondevent->Particle(i);
1401 phi = secondpart->Phi();
1402 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1407 GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1408 denominator->Fill(qout,qside,qlong);
1409 denominator2->Fill(qout,qside,qlong);
1413 secondevent = firstevent;
1419 /***********************************************************/
1421 void AliGenHBTosl::Mix(AliStack* stack, TH3D* numerator, TH3D* numerator2)
1423 //fils numerator with particles from stack
1424 static TStopwatch stoper;
1427 Error("Mix","Stack is null.");
1431 if ( (numerator == 0x0) || (numerator2 == 0x0) )
1433 Error("Mix","Numerator histogram is null.");
1437 Info("Mix","%s",numerator->GetName());
1440 for(Int_t j = 0; j < stack->GetNtrack(); j++ )
1442 TParticle* firstpart = stack->Particle(j);
1443 Float_t phi = firstpart->Phi();
1444 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1446 for(Int_t i = j+1; i < stack->GetNtrack(); i++ )
1448 TParticle* secondpart = stack->Particle(i);
1449 phi = secondpart->Phi();
1450 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1454 GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1455 numerator->Fill(qout,qside,qlong);
1456 numerator2->Fill(qout,qside,qlong);
1463 /***********************************************************/
1465 Double_t AliGenHBTosl::GetQInv(TParticle* f, TParticle* s)
1468 // cout<<f->Px()<<" "<<s->Px()<<endl;
1469 Double_t pxdiff = f->Px() - s->Px();
1470 Double_t pydiff = f->Py() - s->Py();
1471 Double_t pzdiff = f->Pz() - s->Pz();
1472 Double_t ediff = f->Energy() - s->Energy();
1474 Double_t qinvl = ediff*ediff - ( pxdiff*pxdiff + pydiff*pydiff + pzdiff*pzdiff );
1475 Double_t qinv = TMath::Sqrt(TMath::Abs(qinvl));
1478 /***********************************************************/
1480 void AliGenHBTosl::GetQOutQSideQLong(TParticle* f, TParticle* s,Double_t& out, Double_t& side, Double_t& lon)
1482 //returns qout,qside and qlong of the pair of particles
1483 out = side = lon = 10e5;
1485 Double_t pxsum = f->Px() + s->Px();
1486 Double_t pysum = f->Py() + s->Py();
1487 Double_t pzsum = f->Pz() + s->Pz();
1488 Double_t esum = f->Energy() + s->Energy();
1489 Double_t pxdiff = f->Px() - s->Px();
1490 Double_t pydiff = f->Py() - s->Py();
1491 Double_t pzdiff = f->Pz() - s->Pz();
1492 Double_t ediff = f->Energy() - s->Energy();
1493 Double_t kt = 0.5*TMath::Hypot(pxsum,pysum);
1495 Double_t k2 = pxsum*pxdiff+pysum*pydiff;
1506 side = (f->Px()*s->Py()-s->Px()*f->Py())/kt;
1509 Double_t beta = pzsum/esum;
1510 Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
1512 lon = gamma * ( pzdiff - beta*ediff );
1514 // out = TMath::Abs(out);
1515 // side = TMath::Abs(side);
1516 // lon = TMath::Abs(lon);
1519 /***********************************************************/
1521 Double_t AliGenHBTosl::Scale(TH3D* num, TH3D* den)
1523 //Calculates the factor that should be used to scale
1524 //quatience of num and den to 1 at tail
1526 if (GetDebug()) Info("Scale","Entered Scale()");
1529 Error("Scale","No numerator");
1534 Error("Scale","No denominator");
1538 if(fNBinsToScale < 1)
1541 Error("Scale","Number of bins for scaling is smaller thnan 1");
1543 Int_t fNBinsToScaleX = fNBinsToScale;
1544 Int_t fNBinsToScaleY = fNBinsToScale;
1545 Int_t fNBinsToScaleZ = fNBinsToScale;
1547 Int_t nbinsX = num->GetNbinsX();
1548 if (fNBinsToScaleX > nbinsX)
1550 Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
1554 Int_t nbinsY = num->GetNbinsX();
1555 if (fNBinsToScaleY > nbinsY)
1557 Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
1561 Int_t nbinsZ = num->GetNbinsZ();
1562 if (fNBinsToScaleZ > nbinsZ)
1564 Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
1568 if (GetDebug()>0) Info("Scale","No errors detected");
1570 Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
1571 Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
1572 Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
1574 Double_t densum = 0.0;
1575 Double_t numsum = 0.0;
1577 for (Int_t k = offsetZ; k<nbinsZ; k++)
1578 for (Int_t j = offsetY; j<nbinsY; j++)
1579 for (Int_t i = offsetX; i<nbinsX; i++)
1581 if ( num->GetBinContent(i,j,k) > 0.0 )
1584 densum += den->GetBinContent(i,j,k);
1585 numsum += num->GetBinContent(i,j,k);
1590 Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
1591 numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ);
1593 if (numsum == 0) return 0.0;
1594 Double_t ret = densum/numsum;
1596 if(GetDebug() > 0) Info("Scale","returning %f",ret);
1600 /***********************************************************/
1602 void AliGenHBTosl::TestCoarseSignal()
1604 //Tests how works filling from generated histogram shape
1605 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
1607 // for (Int_t i = 0; i < fQCoarseBackground->GetEntries() ;i++)
1610 // fQCoarseSignal->GetRandom3(x,y,z);
1611 // work->Fill(x,y,z);
1614 TCanvas* c1 = new TCanvas();
1617 c1->SaveAs("QTwork.root");
1618 TFile* file = TFile::Open("QTwork.root","update");
1620 work->SetDirectory(0x0);
1623 fQCoarseSignal->Draw();
1624 c1->SaveAs("QTCoarseSignal.root");
1625 file = TFile::Open("QTCoarseSignal.root","update");
1626 fQCoarseSignal->Write();
1627 fQCoarseSignal->SetDirectory(0x0);
1630 fQCoarseBackground->Draw();
1631 c1->SaveAs("QTCoarseBackground.root");
1632 file = TFile::Open("QTCoarseBackground.root","update");
1633 fQCoarseBackground->Write();
1634 fQCoarseBackground->SetDirectory(0x0);
1637 TH1 *result = (TH1*)fQCoarseBackground->Clone("ratio");
1638 result->SetTitle("ratio");
1639 Float_t normfactor = Scale(work,fQCoarseBackground);
1640 result->Divide(work,fQCoarseBackground,normfactor);//work
1645 c1->SaveAs("QTresult.root");
1646 file = TFile::Open("QTresult.root","update");
1648 result->SetDirectory(0x0);
1654 /***********************************************************/
1656 void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr)
1658 //Shortcut to PushTrack(bla,bla,bla,bla.............)
1661 Error("SetTrack(TParticle*,Int_t&)","Particle has zero momentum");
1666 Int_t pdg = p->GetPdgCode();
1667 Double_t px = p->Px();
1668 Double_t py = p->Py();
1669 Double_t pz = p->Pz();
1670 Double_t e = p->Energy();
1671 Double_t vx = p->Vx();
1672 Double_t vy = p->Vy();
1673 Double_t vz = p->Vz();
1674 Double_t tof = p->T();
1677 p->GetPolarisation(pol);
1679 Double_t polx = pol.X();
1680 Double_t poly = pol.Y();
1681 Double_t polz = pol.Z();
1682 TMCProcess mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
1683 Float_t weight = p->GetWeight();
1685 AliGenerator::PushTrack(fTrackIt, -1, pdg, px, py, pz, e, vx, vy, vz, tof,polx, poly, polz, mech, ntr, weight);
1687 /***********************************************************/
1689 void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr, AliStack* stack) const
1691 //Shortcut to SetTrack(bla,bla,bla,bla.............)
1694 Error("SetTrack(TParticle*,Int_t&,AliStack*)","Particle has zero momentum");
1698 Int_t pdg = p->GetPdgCode();
1699 Double_t px = p->Px();
1700 Double_t py = p->Py();
1701 Double_t pz = p->Pz();
1702 Double_t e = p->Energy();
1704 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);
1706 /***********************************************************/
1708 void AliGenHBTosl::Rotate(TVector3& relvector, TVector3& vector)
1710 //This method rotates vector about the angeles that are needed to rotate
1711 //relvector from postion (firstPx,0,0) to its actual positon
1712 //In other words: To make equations easier
1714 static TVector3 first;
1717 first.SetXYZ(relvector.x(),relvector.y(),relvector.z());
1720 Double_t firstPx = TMath::Sqrt( relvector.x()*relvector.x() +
1721 relvector.y()*relvector.y() +
1722 relvector.z()*relvector.z() );
1724 Double_t rotAngleZ = -TMath::ATan2(relvector.y(),relvector.x());//calculating rot angles
1725 relvector.RotateZ(rotAngleZ);
1726 rotAngleZ = -rotAngleZ;
1727 Double_t rotAngleY = -TMath::ATan2(relvector.z(),relvector.x());
1729 vector.RotateY(rotAngleY);
1730 vector.RotateZ(rotAngleZ);
1734 TVector3 test(firstPx,0.0,0.0);
1735 test.RotateY(rotAngleY);
1736 test.RotateZ(rotAngleZ);
1737 ::Info("Rotate","Rotation test: px %f %f",first.x(),test.x());
1738 ::Info("Rotate","Rotation test: py %f %f",first.y(),test.y());
1739 ::Info("Rotate","Rotation test: pz %f %f",first.z(),test.z());
1742 /***********************************************************/
1744 Double_t AliGenHBTosl::Rotate(Double_t x,Double_t y,Double_t z)
1746 //Rotates vector to base where only x - coordinate is no-zero, and returns that
1748 Double_t xylength = TMath::Hypot(x,y);
1749 Double_t sinphi = -y/xylength;
1750 Double_t cosphi = x/xylength;
1752 Double_t xprime = cosphi*x - sinphi*y;
1753 Double_t yprime = sinphi*x + cosphi*y;
1756 Double_t a1 = -TMath::ATan2(v.Y(),v.X());
1760 ::Info("Rotate","Xpr = %f Ypr = %f",xprime,yprime);
1761 ::Info("Rotate","Calc sin = %f, and %f",sinphi,TMath::Sin(a1));
1762 ::Info("Rotate","Calc cos = %f, and %f",cosphi,TMath::Cos(a1));
1765 Double_t xprimezlength = TMath::Hypot(xprime,z);
1767 Double_t sintheta = z/xprimezlength;
1768 Double_t costheta = xprime/xprimezlength;
1771 Double_t xbis = sintheta*z + costheta*(cosphi*x - sinphi*y);
1773 ::Info("Rotate","Calculated rot %f, modulus %f",xbis,TMath::Sqrt(x*x+y*y+z*z));
1776 /***********************************************************/
1778 AliStack* AliGenHBTosl::RotateStack()
1780 //swaps to next stack last goes to first and is reseted
1783 if ( fStackBuffer->GetSize() >= fBufferSize )
1785 stack = (AliStack*)fStackBuffer->Remove(fStackBuffer->Last());
1789 stack = new AliStack(fNpart);
1792 fStackBuffer->AddFirst(stack);
1796 /***********************************************************/
1798 Double_t AliGenHBTosl::GetQInvCorrTheorValue(Double_t qinv) const
1800 //Function (deprecated)
1801 static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1803 return 1.0 + 0.5*TMath::Exp(-qinv*qinv*fQRadius*fQRadius/kFactorsqrd);
1805 /***********************************************************/
1807 Double_t AliGenHBTosl::GetQOutQSideQLongCorrTheorValue(Double_t& out, Double_t& side, Double_t& lon) const
1809 //Theoretical function. Wa want to get correlation of the shape of this function
1810 static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1811 return 1.0 + 0.7*TMath::Exp(-fQRadius*fQRadius*(out*out+side*side+lon*lon)/kFactorsqrd);
1813 /***********************************************************/
1815 Bool_t AliGenHBTosl::CheckParticle(TParticle* p, TParticle* aupair ,AliStack* stack)
1817 //Checks if a given particle is falling into signal region with any other particle
1818 //already existing on stack
1821 if (fSignalRegion <=0) return kFALSE;
1823 for (Int_t i = 0; i < stack->GetNtrack(); i++)
1825 TParticle* part = stack->Particle(i);
1826 if (part == aupair) continue;
1827 Double_t qout = 10e5;
1828 Double_t qside= 10e5;
1829 Double_t qlong= 10e5;
1830 GetQOutQSideQLong(p,part,qout,qside,qlong);
1832 if (TMath::Abs(qout) < fSignalRegion)
1833 if (TMath::Abs(qside) < fSignalRegion)
1834 if (TMath::Abs(qlong) < fSignalRegion)
1839 /***********************************************************/
1841 void AliGenHBTosl::SwapGeneratingHistograms()
1843 //Checks if it is time to swap signal and background histograms
1844 //if yes it swaps them
1845 Int_t threshold = fMinFill;
1846 for (Int_t k = 1; k<=fQNBins; k++)
1848 for (Int_t j = 1; j<=fQNBins; j++)
1850 for (Int_t i = 1; i<=fQNBins; i++)
1852 if ( fQSecondBackground->GetBinContent(i,j,k) < threshold) return;
1859 Info("SwapGeneratingHistograms","*******************************************");
1860 Info("SwapGeneratingHistograms","*******************************************");
1861 Info("SwapGeneratingHistograms","*******************************************");
1862 Info("SwapGeneratingHistograms","**** SWAPPING HISTOGRAMS ****");
1863 Info("SwapGeneratingHistograms","*******************************************");
1864 Info("SwapGeneratingHistograms","*******************************************");
1865 Info("SwapGeneratingHistograms","*******************************************");
1869 fQSignal = fQSecondSignal;
1871 fQSecondSignal->Reset();
1872 fQSecondSignal->SetDirectory(0x0);
1875 fQBackground = fQSecondBackground;
1876 fQSecondBackground = h;
1877 fQSecondBackground->Reset();
1878 fQSecondBackground->SetDirectory(0x0);