1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 #include "AliGenHBTosl.h"
21 //__________________________________________________________
22 /////////////////////////////////////////////////////////////
24 // class AliGenHBTosl //
26 // Genarator simulating particle correlations //
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. //
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 //
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. //
68 ////////////////////////////////////////////////////////////
76 #include <TParticle.h>
80 #include <TStopwatch.h>
83 #include "AliGenCocktailAfterBurner.h"
84 #include "AliGeVSimParticle.h"
85 #include "AliGenGeVSim.h"
86 #include "AliGenHIJINGpara.h"
89 /***********************************************************/
90 ClassImp(AliGenHBTosl)
92 AliGenHBTosl::AliGenHBTosl():
94 fQCoarseBackground(0x0),
99 fQSecondBackground(0x0),
104 fNBinsToScale(Int_t(fQNBins*0.1)),
106 fMaxIterations(1000),
107 fMaxChiSquereChange(0.01),
108 fMaxChiSquerePerNDF(1.5),
111 fSamplePhiMin(-0.01),
112 fSamplePhiMax(TMath::TwoPi()+0.01),
117 //default constructor
119 /***********************************************************/
121 AliGenHBTosl::AliGenHBTosl(Int_t n, Int_t pid):
123 fQCoarseBackground(0x0),
128 fQSecondBackground(0x0),
133 fNBinsToScale(Int_t(fQNBins*0.1)),
135 fSignalShapeCreated(kFALSE),
136 fMaxIterations(1000),
137 fMaxChiSquereChange(0.01),
138 fMaxChiSquerePerNDF(1.5),
141 fSamplePhiMin(-0.01),
142 fSamplePhiMax(TMath::TwoPi()+0.01),
147 //default constructor
149 /***********************************************************/
151 AliGenHBTosl::~AliGenHBTosl()
154 delete fQCoarseSignal;
155 delete fQCoarseBackground;
159 delete fQSecondSignal;
160 delete fQSecondBackground;
163 /***********************************************************/
165 void AliGenHBTosl::Init()
167 //Initializes generator
168 if (fGenerator == 0x0)
171 AliGenHIJINGpara* bkggen = new AliGenHIJINGpara(fNpart*4);
175 AliGenGeVSim * gevsim = new AliGenGeVSim(0.0);
176 AliGeVSimParticle* kplus = new AliGeVSimParticle(fPID,1,fNpart, 0.17, 0.9);
177 gevsim->AddParticleType(kplus);
184 AliMevSimConfig *c = new AliMevSimConfig(1);
185 c->SetRectPlane(1); // reaction plane control, model 4
188 AliGenMevSim *mevsim = new AliGenMevSim(c);
189 mevsim->SetPtRange(0.001, 3);
190 mevsim->SetMomentumRange(0.1, 3);
191 mevsim->SetTrackingFlag(0);
192 mevsim->SetOrigin(0.0, 0.0, 0.0);
193 mevsim->SetSigma(0.0, 0.0, 0.0);
194 AliMevSimParticle *kplus = new AliMevSimParticle(kKPlus, fNpart, 0, 0.25, 0.0, 2, 0.15, 0.0, 0.0 );
195 mevsim->AddParticleType(kplus);
199 fGenerator->SetOrigin(fOrigin[0],fOrigin[1],fOrigin[2]);
200 static const Double_t kDegToRadCF = 180./TMath::Pi();
201 fGenerator->SetMomentumRange(fPtMin,fPtMax);
202 fGenerator->SetPhiRange(kDegToRadCF*fPhiMin,kDegToRadCF*fPhiMax);
203 fGenerator->SetYRange(fYMin,fYMax);
204 fGenerator->SetThetaRange(kDegToRadCF*fThetaMin,kDegToRadCF*fThetaMax);
209 // fQCoarseBackground = new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
210 // fQCoarseSignal = new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
211 // fQCoarseBackground->Sumw2();
212 // fQCoarseSignal->Sumw2();
214 fQSignal = new TH3D("fQSignal1","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
215 fQBackground = new TH3D("fQBackground1","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
217 fQSecondSignal = new TH3D("fQSignal2","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
218 fQSecondBackground = new TH3D("fQBackground2","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
221 fQBackground->Sumw2();
222 fQSecondSignal->Sumw2();
223 fQSecondBackground->Sumw2();
225 fLogFile = new ofstream("BadEvent",ios::out);
228 /***********************************************************/
230 void AliGenHBTosl::Generate()
234 ofstream& logfile = *fLogFile;
235 logfile<<"Generate"<<"Attempts to generate "<<fNpart<<" particles.";
238 if (fStackBuffer == 0x0) fStackBuffer = new TList();
239 //Here is initialization level
240 if (fSignalShapeCreated == kFALSE)
242 TH3D *hs = 0x0, *hb = 0x0;
245 file = TFile::Open("QTSignal.root");
248 hs = (TH3D*)file->Get("fQSignal1");
249 if (hs) hs->SetDirectory(0x0);
253 file = TFile::Open("QTBackground.root");
256 hb = (TH3D*)file->Get("fQBackground1");
257 if (hb) hb->SetDirectory(0x0);
263 Info("Generate","**********************************");
264 Info("Generate","Found starting histograms in files");
265 Info("Generate","**********************************");
273 TH3D *cs = 0x0, *cb = 0x0;
274 file = TFile::Open("QTCoarseBackground.root");
277 cb = (TH3D*)file->Get("fQCoarseBackground");
278 if (cb) cb->SetDirectory(0x0);
282 file = TFile::Open("QTCoarseSignal.root");
285 cs = (TH3D*)file->Get("fQCoarseSignal");
286 if (cs) cs->SetDirectory(0x0);
293 Info("Generate","Got Coarse signal and bkg from files");
294 delete fQCoarseBackground;
295 delete fQCoarseSignal;
297 fQCoarseBackground = cb;
303 Info("Generate","Got Coarse bkg from file");
304 delete fQCoarseBackground;
305 fQCoarseBackground = cb;
309 fQCoarseBackground = new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
310 fQCoarseBackground->Sumw2();
312 FillCoarse(); //create coarse background - just to know the spectrum
315 fQCoarseSignal = new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
316 fQCoarseSignal->Sumw2();
317 FillCoarseSignal();//create first coarse signal by brutal multplication coarse background and required function shape
320 StartSignal(); //Initilizes the stack that is used for generation
322 fSignalShapeCreated = kTRUE;
325 AliStack* stack = RotateStack();
327 AliStack* genstack = fGenerator->GetStack();
330 genstack = new AliStack(fNpart);
331 fGenerator->SetStack(genstack);
338 fGenerator->Generate();
339 Int_t j = 0, ntr = 0;
340 if ( genstack->GetNtrack() < fNpart/2)
342 Warning("Generate","************************************************************");
343 Warning("Generate","Generator generated (%d) less particles then expected (%d).",
344 stack->GetNtrack(),fNpart/2);
345 Warning("Generate","************************************************************");
348 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
349 work->SetDirectory(0x0);
352 Double_t*** chiarray = new Double_t** [fQNBins+1];
353 Double_t*** sigarray = new Double_t** [fQNBins+1];
355 for (Int_t i = 1; i<=fQNBins; i++)
357 chiarray[i] = new Double_t* [fQNBins+1];
358 sigarray[i] = new Double_t* [fQNBins+1];
360 for (Int_t k = 1; k<=fQNBins; k++)
362 chiarray[i][k] = new Double_t [fQNBins+1];
363 sigarray[i][k] = new Double_t [fQNBins+1];
368 Double_t scale = Scale(fQSignal,fQBackground);
369 work->Divide(fQSignal,fQBackground,scale);
371 Double_t binwdh = work->GetBinWidth(1)/2.;
373 for (Int_t k = 1; k<=fQNBins; k++)
375 Double_t z = work->GetZaxis()->GetBinCenter(k);
376 for (Int_t j = 1; j<=fQNBins; j++)
378 Double_t y = work->GetYaxis()->GetBinCenter(j);
379 for (Int_t i = 1; i<=fQNBins; i++)
381 sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
382 Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qinv)
383 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
384 Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
385 chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
393 Int_t middlebin = fQNBins/2;
395 for (Int_t k = middlebin-5; k < middlebin+5; k++)
397 Double_t tx = work->GetXaxis()->GetBinCenter(30);
398 Double_t ty = work->GetYaxis()->GetBinCenter(30);
399 Double_t tz = work->GetZaxis()->GetBinCenter(k);
400 sprintf(msg,"% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
405 for (Int_t k = middlebin-5; k < middlebin+5; k++)
407 sprintf(msg,"% 6.5f ",work->GetBinContent(30,30,k));
412 for (Int_t k = middlebin-5; k < middlebin+5; k++)
414 sprintf(msg,"% 6.5f ",chiarray[30][30][k]);
419 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
420 TParticle* second = &particle;
422 Bool_t shortloop = kTRUE;
423 Int_t sc = 0;//security check against infinite loop
425 while ( (ntr+1) < fNpart)
447 loopmax = fQNBins/2+fQNBins/4;
448 loopmin = fQNBins/2-fQNBins/4;
452 for (Int_t k = loopmin; k <=loopmax; k++ )
454 qlong = work->GetZaxis()->GetBinCenter(k);
455 for (Int_t j = loopmin; j<=loopmax; j++)
457 qside = work->GetYaxis()->GetBinCenter(j);
458 for (Int_t i = loopmin; i<=loopmax; i++)
460 qout = work->GetXaxis()->GetBinCenter(i);
461 if (chiarray[xmax][ymax][zmax] < chiarray[i][j][k])
468 // Double_t qdist = TMath::Sqrt(qout*qout + qside*qside + qlong*qlong);
470 // Double_t fact = chiarray[i][j][k];//chiarray is chi2
471 // if (fact > work->GetBinError(i,j,k))//if differece between what we want and
472 // { //what we have is bigger than stat. error
473 // xmax = i; //we force to fill that bin
481 Double_t qlongc = work->GetZaxis()->GetBinCenter(zmax);
482 Double_t qsidec = work->GetYaxis()->GetBinCenter(ymax);
483 Double_t qoutc = work->GetXaxis()->GetBinCenter(xmax);
486 sprintf(msg,"Generate Fill bin chi2(%d,%d,%d)=%f",xmax,ymax,zmax,chiarray[xmax][ymax][zmax]);
489 qout = gRandom->Uniform(qoutc -binwdh, qoutc +binwdh);
490 qside = gRandom->Uniform(qsidec-binwdh, qsidec+binwdh);
491 qlong = gRandom->Uniform(qlongc-binwdh, qlongc+binwdh);
493 TParticle* first = 0;
494 while (j < genstack->GetNtrack())
496 TParticle* tmpp = genstack->Particle(j++);
497 if (tmpp->GetPdgCode() == fPID)
499 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
509 if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
513 //Here put the check if this particle do not fall into signal region with other partticle
515 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
518 //Info("StartSignal","Can not find momenta for this OSL and particle");
521 //in case this particle is falling into signal area with another
522 //particle we take a another pair
523 //it can intruduce artificial correlations
524 Bool_t checkresult = CheckParticle(second,first,stack);
525 if ( checkresult && (sc < 10) )
532 //Put on output stack
534 SetTrack(second,ntr);
536 //Put on internal stack
538 SetTrack(first,etmp,stack);
539 SetTrack(second,etmp,stack);
541 Double_t y = GetQOutQSideQLongCorrTheorValue(qoutc,qsidec,qlongc);
543 sigarray[xmax][ymax][zmax] ++;
544 chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
545 chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
549 Mix(fStackBuffer,fQBackground,fQSecondSignal); //upgrate background
550 Mix(stack,fQSignal,fQSecondBackground); //upgrate signal
554 for (Int_t i = 1; i<=fQNBins; i++)
556 for (Int_t k = 1; k<=fQNBins; k++)
558 delete [] chiarray[i][k];
559 delete [] sigarray[i][k];
561 delete [] chiarray[i];
562 delete [] sigarray[i];
567 /***********************************************************/
569 void AliGenHBTosl::GetOneD(TParticle* first, TParticle* second,Double_t qinv)
571 //deprecated method that caclulates momenta of the second particle
572 // out of qinv and the first particle
573 //first particle is rotated that only X is non-zero
576 Double_t m = first->GetMass();
577 Double_t msqrd = m*m;
578 Double_t fourmassSquered = 4.*msqrd;
580 //Condition that R must fullfill to be possible to have qinv less smaller then randomized
581 // Double_t rRange = qinv*TMath::Sqrt(qinv*qinv + fourmassSquered)/fourmassSquered;
582 // Double_t r = gRandom->Uniform(rRange);
584 Double_t r = gRandom->Uniform(qinv);
585 Double_t phi = gRandom->Uniform(TMath::TwoPi());
587 Double_t firstPx = first->P();//first particle is rotated that only X is non-zero thus P==Px
588 Double_t px = 2.*msqrd*firstPx + firstPx*qinv*qinv;
589 Double_t temp = qinv*qinv*qinv*qinv + fourmassSquered * (qinv*qinv - r*r );
592 Error("GetOneD","temp is less then 0: %f",temp);
595 temp = temp*(msqrd+firstPx*firstPx);
597 px = (px - TMath::Sqrt(temp))/(2.*msqrd);
599 Double_t py = r*TMath::Sin(phi);
600 Double_t pz = r*TMath::Cos(phi);
602 TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
603 TVector3 vector(px,py,pz);
604 Rotate(firstpvector,vector);
606 Double_t e = TMath::Sqrt(msqrd + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
607 second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
608 // TParticle* f = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, firstPx,0,0,e=TMath::Sqrt(msqrd+firstPx*firstPx),0.0,0.0,0.0,0.0);
609 // TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
610 // px, py, pz, e, vx, vy, vz, tof);
612 AliDebug(1,Form("Randomized qinv = %f, obtained = %f",qinv,GetQInv(first,second)));
615 /***********************************************************/
617 Int_t AliGenHBTosl::GetThreeD(TParticle* first,TParticle* second, Double_t qout, Double_t qside, Double_t qlong)
619 //deprecated method that caclulates momenta of the second particle
620 //out of qout qside and qlong and the first particle
621 Double_t m = first->GetMass();
624 Double_t px = first->P();//first particle is rotated that only X is non-zero thus P==Px
625 Double_t px2 = px*px;
628 Double_t qout2 = qout*qout;
629 Double_t qside2 = qside*qside;
630 Double_t qlong2 = qlong*qlong;
633 Double_t util1 = 4.*px2 - qside2;//4*P1x^2 - Y^2
636 Info("GetThreeD","4.*px2* - qside2 is negative px: %f, qside: %f",px,qside);
639 Double_t util2 = TMath::Sqrt(px2*qout2*util1);
642 Double_t p2x,p2y,p2z;
644 // if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))
648 Double_t tmp = px*(2.*px2 - qside2);
653 tmp = qout - TMath::Sqrt(util1);
654 p2y = - (tmp*qside)/(2.*px);
657 tmp = 4.*m2 + 2.*qout2+qlong2;
659 tmp -= 2.*util2;//!!!
663 Double_t m2px2 = m2+px2;
664 Double_t tmp2 = m2px2*tmp;
666 tmp = 4.*(m2px2+qout2) + qlong2;
669 tmp *= 4.*(m2px2) + qlong2;
670 tmp *= qlong2*qlong2;
675 Error("","Argument of sqrt is negative");
679 tmp2 += TMath::Sqrt(tmp);
681 tmp = 8.0*px*m2px2*m2px2;
682 p2z = -TMath::Sqrt(tmp2/tmp);
683 if (qlong < 0) p2z = -p2z;
688 Double_t tmp = px*(2.*px2 - qside2);
693 tmp = qout - TMath::Sqrt(util1);
694 p2y = - (tmp*qside)/(2.*px);
697 tmp = 4.*m2 + 2.*qout2+qlong2;
699 tmp += 2.*util2;//!!!
703 Double_t m2px2 = m2+px2;
704 Double_t tmp2 = m2px2*tmp;
706 tmp = 4.*(m2px2+qout2) + qlong2;
709 tmp *= 4.*(m2px2) + qlong2;
710 tmp *= qlong2*qlong2;
715 Error("","Argument of sqrt is negative");
719 tmp2 += TMath::Sqrt(tmp);
721 tmp = 8.0*px*m2px2*m2px2;
722 p2z = -TMath::Sqrt(tmp2/tmp);
723 if (qlong < 0) p2z = -p2z;
726 // if ( (qside >= 0) && (qout >= 0) && (qlong >= 0)) p2z = -p2z;
728 TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
729 TVector3 vector(p2x,p2y,p2z);
730 Rotate(firstpvector,vector);
732 Double_t e = TMath::Sqrt(m2 + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
733 second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
736 if ( AliDebugLevel() > 3 )
738 e=TMath::Sqrt(m2+px*px);
739 TParticle* f = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, px , 0.0, 0.0, e,0.0,0.0,0.0,0.0);
741 e = TMath::Sqrt(m2 + p2x*p2x + p2y*p2y + p2z*p2z);
742 TParticle* s = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, p2x, p2y, p2z, e, 0.0, 0.0, 0.0, 0.0);
745 GetQOutQSideQLong(f,s,qo, qs, ql);
747 Info("GetThreeD","TEST");
750 Info("GetThreeD","Required %f %f %f",qout,qside,qlong);
751 Info("GetThreeD","Got %f %f %f",qo,qs,ql);
755 Info("GetThreeD","!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
760 /***********************************************************/
762 void AliGenHBTosl::StartSignal()
764 //Starts the signal histograms
765 ofstream& logfile = *fLogFile;
767 logfile<<"************************************************"<<endl;
768 logfile<<"************************************************"<<endl;
769 logfile<<" StartSignal "<<endl;
770 logfile<<"************************************************"<<endl;
771 logfile<<"************************************************"<<endl;
777 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
778 TParticle* second = &particle;
780 TIter next(fStackBuffer);
781 while(( stack=(AliStack*)next() ))
786 AliStack* genstack = fGenerator->GetStack();
789 genstack = new AliStack(fNpart);
790 fGenerator->SetStack(genstack);
798 //We alread have detailed histograms and we do not need Coarse anymore
799 delete fQCoarseSignal;
800 delete fQCoarseBackground;
801 fQCoarseSignal = 0x0;
802 fQCoarseBackground = 0x0;
805 const Double_t kNDF = fQNBins*fQNBins*fQNBins;
807 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
809 work->SetDirectory(0x0);
812 Double_t binwdh = work->GetBinWidth(1)/2.;
814 Double_t*** chiarray = new Double_t** [fQNBins+1];
815 Double_t*** sigarray = new Double_t** [fQNBins+1];
817 for (Int_t i = 1; i<=fQNBins; i++)
819 chiarray[i] = new Double_t* [fQNBins+1];
820 sigarray[i] = new Double_t* [fQNBins+1];
822 for (Int_t k = 1; k<=fQNBins; k++)
824 chiarray[i][k] = new Double_t [fQNBins+1];
825 sigarray[i][k] = new Double_t [fQNBins+1];
830 Float_t chisqrchange = fMaxChiSquereChange + 1.;
831 Float_t chisqrPerDF = fMaxChiSquerePerNDF;
832 Float_t chisqrold = 0.0;
835 Int_t niterations = 1;
836 Int_t rotaxisorder = 1;//defines order of looping over 3D histo (X,Y,Z or Y,Z,X or Z,X,Y)
839 Bool_t shortloop = kTRUE;
840 TCanvas* c1 = new TCanvas();
844 Info("StartSignal","\n\n\n\nSecond Pass\n\n\n\n");
846 while ( ( (chisqrPerDF > fMaxChiSquereChange) || flag) && (niterations++ < fMaxIterations) )
849 logfile<<"StartSignal\n";
850 logfile<<" Row 1 Theory, 2 current value, 3 Chi2 \n";
852 Double_t chisqrnew = 0.0;
855 Double_t scale = Scale(fQSignal,fQBackground);
856 work->Divide(fQSignal,fQBackground,scale);
858 if ( (counter%100) == 0)
862 sprintf(buff,"QTWorkPass2.%3d.root",counter);
863 TFile* file = TFile::Open(buff,"update");
865 work->SetDirectory(0x0);
869 sprintf(buff,"QTBackgroundPass2.%3d.root",counter);
870 file = TFile::Open(buff,"update");
871 fQBackground->Write();
872 fQBackground->SetDirectory(0x0);
876 sprintf(buff,"QTSignalPass2.%3d.root",counter);
877 file = TFile::Open(buff,"update");
879 fQSignal->SetDirectory(0x0);
885 Int_t novertresh = 0;
886 for (Int_t k = 1; k<=fQNBins; k++)
888 Double_t z = work->GetZaxis()->GetBinCenter(k);
889 for (Int_t j = 1; j<=fQNBins; j++)
891 Double_t y = work->GetYaxis()->GetBinCenter(j);
892 for (Int_t i = 1; i<=fQNBins; i++)
894 Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qout)
895 sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
896 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
897 Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
898 chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
899 Double_t be = work->GetBinError(i,j,k);
900 chisqrnew += diff*diff/(be*be);//add up chisq
902 //even if algorithm is stable (chi sqr change less then threshold)
903 //and any bin differs more then 5% from expected value we continue
904 Double_t fact = diff;
905 if (TMath::Abs(fact) > 0.1)
919 for (Int_t k = 25; k < 36; k++)
921 Double_t tx = work->GetXaxis()->GetBinCenter(30);
922 Double_t ty = work->GetYaxis()->GetBinCenter(30);
923 Double_t tz = work->GetZaxis()->GetBinCenter(k);
924 sprintf(msg,"% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
929 for (Int_t k = 25; k < 36; k++)
931 sprintf(msg,"%6.5f ",work->GetBinContent(30,30,k));
936 for (Int_t k = 25; k < 36; k++)
938 sprintf(msg,"% 6.5f ",chiarray[30][30][k]);
943 chisqrchange = TMath::Abs(chisqrnew - chisqrold)/chisqrnew;
944 chisqrold = chisqrnew;
946 chisqrPerDF = chisqrnew/kNDF;
948 Info("StartSignal","Iteration %d Chi-sq change %f%%",niterations,chisqrchange*100.0);
949 Info("StartSignal","ChiSq = %f, NDF = %f, ChiSq/NDF = %f",chisqrnew, kNDF, chisqrPerDF );
950 Info("StartSignal","novertresh = %d",novertresh);
953 stack = RotateStack();
955 fGenerator->Generate();
956 Int_t ninputparticle = 0, ntr = 0;
957 if ( genstack->GetNtrack() < fNpart/2)
959 Warning("StartSignal","**********************************");
960 Warning("StartSignal","Generator generated (%d) less particles then expected (%d).",
961 genstack->GetNtrack(),fNpart/2);
962 Warning("StartSignal","**********************************");
965 Int_t sc = 0; //security check against infinite loop
966 while ( (ntr+1) < fNpart)//ntr is number of track generated up to now
981 switch (rotaxisorder)
1000 if (rotaxisorder > 3) rotaxisorder = 1;
1014 // Bool_t force = kFALSE;
1015 for ( k = 1; k <=nrange;k++ )
1017 for ( j = 1; j<=nrange; j++)
1019 for ( i = 1; i<=nrange; i++)
1021 if ( (chiarray[*cx][*cy][*cz]) > (chiarray[xmax][ymax][zmax]) )
1028 // Double_t fact = chiarray[*cx][*cy][*cz];//chiarray is chi2*qinv
1029 // if (fact > work->GetBinError(*cx,*cy,*cz))//if differece between what we want and
1030 // { //what we have is bigger than stat. error
1031 // //we force to fill that bin
1032 // qout = work->GetXaxis()->GetBinCenter(*cx);
1033 // qside = work->GetYaxis()->GetBinCenter(*cy);
1034 // qlong = work->GetZaxis()->GetBinCenter(*cz);
1036 // Info("StartSignal"," bin: (%d,%d,%d) loop status (%d,%d,%d) \nUsing Force: chiarray: %f \nq(o,s,l): (%f,%f,%f) signal: %d background: %d binerror: %f",
1037 // *cx,*cy,*cz,i,j,k,fact,qout,qside,qlong,
1038 // (Int_t)sigarray[*cx][*cy][*cz],(Int_t)fQBackground->GetBinContent(*cx,*cy,*cz),work->GetBinError(*cx,*cy,*cz));
1044 // if (force) break;
1046 // if (force) break;
1050 qout = work->GetXaxis()->GetBinCenter(xmax);
1051 qside = work->GetYaxis()->GetBinCenter(ymax);
1052 qlong = work->GetZaxis()->GetBinCenter(zmax);
1054 // Info("StartSignal"," bin: (%d,%d,%d) chiarray: %f \nq(o,s,l): (%f,%f,%f) signal: %d background: %d binerror: %f",
1055 // xmax,ymax,zmax,chiarray[xmax][ymax][zmax],qout,qside,qlong,
1056 // (Int_t)sigarray[xmax][ymax][zmax],
1057 // (Int_t)fQBackground->GetBinContent(xmax,ymax,zmax),
1058 // work->GetBinError(xmax,ymax,zmax));
1060 qout = gRandom->Uniform(qout-binwdh,qout+binwdh);
1061 qside = gRandom->Uniform(qside-binwdh,qside+binwdh);
1062 qlong = gRandom->Uniform(qlong-binwdh,qlong+binwdh);
1064 TParticle* first = 0;
1065 while (ninputparticle < genstack->GetNtrack())
1067 TParticle* tmpp = genstack->Particle(ninputparticle++);
1068 if (tmpp->GetPdgCode() == fPID)
1070 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1080 if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
1084 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1087 Info("StartSignal","Can not find momenta for this OSL and particle OSL = %f %f %f",qout,qside,qlong);
1093 //in case this particle is falling into signal area with another
1094 //particle we take a another pair
1095 //it can intruduce artificial correlations
1096 Bool_t checkresult = CheckParticle(second,first,stack);
1097 if ( checkresult && (sc < 10) )
1104 //Put on output stack
1105 SetTrack(first,ntr,stack);
1106 SetTrack(second,ntr,stack);
1108 Double_t y = GetQOutQSideQLongCorrTheorValue(qout,qside,qlong);
1110 sigarray[xmax][ymax][zmax] ++;
1111 chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
1112 chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
1115 Info("StartSignal","Mixing background...");
1116 Mix(fStackBuffer,fQBackground,fQSecondBackground); //upgrate background
1117 Info("StartSignal","Mixing signal...");
1118 Mix(stack,fQSignal,fQSecondSignal); //upgrate background
1121 // if ( (chisqrPerDF < 2.0) && (fSwapped == kFALSE) )
1123 // SwapGeneratingHistograms();
1127 TFile* filef = TFile::Open("QTBackground.root","recreate");
1128 fQBackground->Write();
1129 fQBackground->SetDirectory(0x0);
1133 filef = TFile::Open("QTSignal.root","recreate");
1135 fQSignal->SetDirectory(0x0);
1143 for (Int_t i = 1; i<=fQNBins; i++)
1145 for (Int_t k = 1; k<=fQNBins; k++)
1147 delete [] chiarray[i][k];
1148 delete [] sigarray[i][k];
1150 delete [] chiarray[i];
1151 delete [] sigarray[i];
1157 /***********************************************************/
1159 void AliGenHBTosl::StartSignalPass1()
1161 //This method makes first part of the initialization of working histograms
1162 //It randomizes qout, qside and qlong from the coarse signal histogram
1164 Bool_t flag = kTRUE;
1165 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
1166 TParticle* second = &particle;
1171 Info("StartSignalPass1","\n\nFirst Pass\n\n");
1175 Info("StartSignalPass1","NextEvent");
1176 AliStack* stack = RotateStack();
1177 AliStack* genstack = fGenerator->GetStack();
1179 fGenerator->Generate();
1180 Int_t j = 0, ntr = 0;
1181 if ( genstack->GetNtrack() < fNpart/2)
1183 Warning("StartSignalPass1","**********************************");
1184 Warning("StartSignalPass1","Generator generated (%d) less particles then expected (%d).",
1185 genstack->GetNtrack(),fNpart/2);
1186 Warning("StartSignalPass1","**********************************");
1189 Int_t sc = 0;//security check against infinite loop
1190 while ((ntr+1)<fNpart)
1193 // Info("StartSignal","Number of track on output stack: = %d", ntr);
1194 // Info("StartSignal","Number of track on input stack: = %d\n", j);
1196 TParticle* first = 0;
1197 while (j < genstack->GetNtrack())
1199 TParticle* tmpp = genstack->Particle(j++);
1200 if (tmpp->GetPdgCode() == fPID)
1202 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1209 Info("StartSignalPass1","Particle did not pass the safety check 1");
1217 if ( fDebug > 2 ) Info("StartSignalPass1","No more particles of that type");
1222 SetTrack(first,ntr,stack);
1224 fQCoarseSignal->GetRandom3(qout,qside,qlong);
1226 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1229 //Info("StartSignal","Can not find momenta for this OSL and particle");
1232 //in case this particle is falling into signal area with another
1233 //particle we take a another pair
1234 //it can intruduce artificial correlations
1235 Bool_t checkresult = CheckParticle(second,first,stack);
1236 if ( checkresult && (sc < 10) )
1239 Info("StartSignalPass1","Particle did not pass the safety check 2");
1246 SetTrack(second,ntr,stack);
1249 Mix(stack,fQSignal,fQSecondSignal);
1250 Mix(fStackBuffer,fQBackground,fQSecondBackground);
1254 for (Int_t k = 1; k<=fQNBins; k++)
1256 for (Int_t j = 1; j<=fQNBins; j++)
1258 for (Int_t i = 1; i<=fQNBins; i++)
1260 if ( (fQBackground->GetBinContent(i,j,k) < fMinFill) )
1262 //(fQSignal->GetBinContent(i,j,k) < fMinFill) ||
1263 Info("StartSignalPass1","bin (%d,%d,%d): signal=%f background=%f",i,j,k,
1264 fQSignal->GetBinContent(i,j,k),fQBackground->GetBinContent(i,j,k));
1265 flag = kTRUE;//continue while
1266 break;//breakes for not while
1269 if (flag == kTRUE) break;
1271 if (flag == kTRUE) break;
1278 /***********************************************************/
1280 void AliGenHBTosl::FillCoarseSignal()
1282 //Makes coarse signal by multiplying the coarse background and required function
1283 Info("FillCoarseSignal","START");
1284 for (Int_t k = 1; k <=fQNBins ;k++ )
1286 Double_t z = fQCoarseBackground->GetZaxis()->GetBinCenter(k);
1287 for (Int_t j = 1; j <=fQNBins; j++)
1289 Double_t y = fQCoarseBackground->GetYaxis()->GetBinCenter(j);
1290 for (Int_t i = 1; i <=fQNBins; i++)
1292 Double_t x = fQCoarseBackground->GetXaxis()->GetBinCenter(i);
1293 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);
1294 Info("FillCoarseSignal","Bin (%d,%d,%d): osl(%f,%f,%f)=%f",i,j,k,x,y,z,v);
1295 fQCoarseSignal->SetBinContent(i,j,k,v*fQCoarseBackground->GetBinContent(i,j,k));
1300 //if (AliDebugLevel())
1303 Info("FillCoarseSignal","DONE");
1305 /***********************************************************/
1307 void AliGenHBTosl::FillCoarse()
1309 //creates the statistical background histogram on the base of input from
1311 Info("FillCoarse","START");
1317 TH3D tmph("tmph","tmph",2,0,1,2,0,1,2,0,1);
1322 // if (niter > 20) break;
1324 cout<<niter++<<" bincont "<<fQCoarseBackground->GetBinContent(30,30,28)
1325 <<" "<<fQCoarseBackground->GetBinContent(30,30,29)
1326 <<" "<<fQCoarseBackground->GetBinContent(30,30,30)
1327 <<" "<<fQCoarseBackground->GetBinContent(30,30,31)
1328 <<" "<<fQCoarseBackground->GetBinContent(30,30,32)
1332 stack = RotateStack();
1333 fGenerator->SetStack(stack);
1335 fGenerator->Generate();
1337 Mix(fStackBuffer,fQCoarseBackground,&tmph);
1341 Info("FillCoarse","fMinFill = %d",fMinFill);
1342 for (Int_t k = 1; k<=fQNBins; k++)
1344 for (Int_t j = 1; j<=fQNBins; j++)
1346 for (Int_t i = 1; i<=fQNBins; i++)
1348 if ( fQCoarseBackground->GetBinContent(i,j,k) < fMinFill )
1351 Info("FillCoarse","bin (%d,%d,%d)=%f",i,j,k,fQCoarseBackground->GetBinContent(i,j,k));
1363 fGenerator->SetStack(0x0);
1364 Info("FillCoarse","DONE");
1367 /***********************************************************/
1369 void AliGenHBTosl::Mix(TList* eventbuffer,TH3D* denominator,TH3D* denominator2)
1371 //Fills denominators
1372 //Mixes events stored in the eventbuffer and fills the background histograms
1373 static TStopwatch stoper;
1375 if (eventbuffer == 0x0)
1377 Error("Mix","Buffer List is null.");
1381 if (denominator == 0x0)
1383 Error("Mix","Denominator histogram is null.");
1387 if (denominator2 == 0x0)
1389 Error("Mix","Denominator2 histogram is null.");
1393 Info("Mix","%s",denominator->GetName());
1396 TIter next(eventbuffer);
1397 AliStack* firstevent;
1398 AliStack* secondevent = 0x0;
1400 while(( firstevent=(AliStack*)next() ))
1402 if (secondevent == 0x0)
1404 secondevent = firstevent;
1407 // Info("Mix","Mixing %#x with %#x",firstevent,secondevent);
1408 for(Int_t j = 0; j < firstevent->GetNtrack(); j++ )
1410 TParticle* firstpart = firstevent->Particle(j);
1412 Float_t phi = firstpart->Phi();
1413 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1415 // Info("Mix","Mixing %d phi %f min %f max %f",j,phi,fSamplePhiMin,fSamplePhiMax);
1417 for(Int_t i = 0; i < secondevent->GetNtrack(); i++ )
1419 TParticle* secondpart = secondevent->Particle(i);
1420 phi = secondpart->Phi();
1421 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1426 GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1427 denominator->Fill(qout,qside,qlong);
1428 denominator2->Fill(qout,qside,qlong);
1432 secondevent = firstevent;
1438 /***********************************************************/
1440 void AliGenHBTosl::Mix(AliStack* stack, TH3D* numerator, TH3D* numerator2)
1442 //fils numerator with particles from stack
1443 static TStopwatch stoper;
1446 Error("Mix","Stack is null.");
1450 if ( (numerator == 0x0) || (numerator2 == 0x0) )
1452 Error("Mix","Numerator histogram is null.");
1456 Info("Mix","%s",numerator->GetName());
1459 for(Int_t j = 0; j < stack->GetNtrack(); j++ )
1461 TParticle* firstpart = stack->Particle(j);
1462 Float_t phi = firstpart->Phi();
1463 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1465 for(Int_t i = j+1; i < stack->GetNtrack(); i++ )
1467 TParticle* secondpart = stack->Particle(i);
1468 phi = secondpart->Phi();
1469 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1473 GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1474 numerator->Fill(qout,qside,qlong);
1475 numerator2->Fill(qout,qside,qlong);
1482 /***********************************************************/
1484 Double_t AliGenHBTosl::GetQInv(TParticle* f, TParticle* s)
1487 // cout<<f->Px()<<" "<<s->Px()<<endl;
1488 Double_t pxdiff = f->Px() - s->Px();
1489 Double_t pydiff = f->Py() - s->Py();
1490 Double_t pzdiff = f->Pz() - s->Pz();
1491 Double_t ediff = f->Energy() - s->Energy();
1493 Double_t qinvl = ediff*ediff - ( pxdiff*pxdiff + pydiff*pydiff + pzdiff*pzdiff );
1494 Double_t qinv = TMath::Sqrt(TMath::Abs(qinvl));
1497 /***********************************************************/
1499 void AliGenHBTosl::GetQOutQSideQLong(TParticle* f, TParticle* s,Double_t& out, Double_t& side, Double_t& lon)
1501 //returns qout,qside and qlong of the pair of particles
1502 out = side = lon = 10e5;
1504 Double_t pxsum = f->Px() + s->Px();
1505 Double_t pysum = f->Py() + s->Py();
1506 Double_t pzsum = f->Pz() + s->Pz();
1507 Double_t esum = f->Energy() + s->Energy();
1508 Double_t pxdiff = f->Px() - s->Px();
1509 Double_t pydiff = f->Py() - s->Py();
1510 Double_t pzdiff = f->Pz() - s->Pz();
1511 Double_t ediff = f->Energy() - s->Energy();
1512 Double_t kt = 0.5*TMath::Hypot(pxsum,pysum);
1514 Double_t k2 = pxsum*pxdiff+pysum*pydiff;
1525 side = (f->Px()*s->Py()-s->Px()*f->Py())/kt;
1528 Double_t beta = pzsum/esum;
1529 Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
1531 lon = gamma * ( pzdiff - beta*ediff );
1533 // out = TMath::Abs(out);
1534 // side = TMath::Abs(side);
1535 // lon = TMath::Abs(lon);
1538 /***********************************************************/
1540 Double_t AliGenHBTosl::Scale(TH3D* num, TH3D* den)
1542 //Calculates the factor that should be used to scale
1543 //quatience of num and den to 1 at tail
1545 AliDebug(1,"Entered");
1548 AliError("No numerator");
1553 AliError("No denominator");
1557 if(fNBinsToScale < 1)
1560 AliError("Number of bins for scaling is smaller thnan 1");
1562 Int_t fNBinsToScaleX = fNBinsToScale;
1563 Int_t fNBinsToScaleY = fNBinsToScale;
1564 Int_t fNBinsToScaleZ = fNBinsToScale;
1566 Int_t nbinsX = num->GetNbinsX();
1567 if (fNBinsToScaleX > nbinsX)
1569 AliError("Number of X bins for scaling is bigger thnan number of bins in histograms");
1573 Int_t nbinsY = num->GetNbinsX();
1574 if (fNBinsToScaleY > nbinsY)
1576 AliError("Number of Y bins for scaling is bigger thnan number of bins in histograms");
1580 Int_t nbinsZ = num->GetNbinsZ();
1581 if (fNBinsToScaleZ > nbinsZ)
1583 AliError("Number of Z bins for scaling is bigger thnan number of bins in histograms");
1587 AliDebug(1,"No errors detected");
1589 Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
1590 Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
1591 Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
1593 Double_t densum = 0.0;
1594 Double_t numsum = 0.0;
1596 for (Int_t k = offsetZ; k<nbinsZ; k++)
1597 for (Int_t j = offsetY; j<nbinsY; j++)
1598 for (Int_t i = offsetX; i<nbinsX; i++)
1600 if ( num->GetBinContent(i,j,k) > 0.0 )
1603 densum += den->GetBinContent(i,j,k);
1604 numsum += num->GetBinContent(i,j,k);
1608 AliDebug(1,Form("numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
1609 numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ));
1611 if (numsum == 0) return 0.0;
1612 Double_t ret = densum/numsum;
1614 AliDebug(1,Form("returning %f",ret));
1618 /***********************************************************/
1620 void AliGenHBTosl::TestCoarseSignal()
1622 //Tests how works filling from generated histogram shape
1623 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
1625 // for (Int_t i = 0; i < fQCoarseBackground->GetEntries() ;i++)
1628 // fQCoarseSignal->GetRandom3(x,y,z);
1629 // work->Fill(x,y,z);
1632 TCanvas* c1 = new TCanvas();
1635 c1->SaveAs("QTwork.root");
1636 TFile* file = TFile::Open("QTwork.root","update");
1638 work->SetDirectory(0x0);
1641 fQCoarseSignal->Draw();
1642 c1->SaveAs("QTCoarseSignal.root");
1643 file = TFile::Open("QTCoarseSignal.root","update");
1644 fQCoarseSignal->Write();
1645 fQCoarseSignal->SetDirectory(0x0);
1648 fQCoarseBackground->Draw();
1649 c1->SaveAs("QTCoarseBackground.root");
1650 file = TFile::Open("QTCoarseBackground.root","update");
1651 fQCoarseBackground->Write();
1652 fQCoarseBackground->SetDirectory(0x0);
1655 TH1 *result = (TH1*)fQCoarseBackground->Clone("ratio");
1656 result->SetTitle("ratio");
1657 Float_t normfactor = Scale(work,fQCoarseBackground);
1658 result->Divide(work,fQCoarseBackground,normfactor);//work
1663 c1->SaveAs("QTresult.root");
1664 file = TFile::Open("QTresult.root","update");
1666 result->SetDirectory(0x0);
1672 /***********************************************************/
1674 void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr)
1676 //Shortcut to PushTrack(bla,bla,bla,bla.............)
1679 Error("SetTrack(TParticle*,Int_t&)","Particle has zero momentum");
1684 Int_t pdg = p->GetPdgCode();
1685 Double_t px = p->Px();
1686 Double_t py = p->Py();
1687 Double_t pz = p->Pz();
1688 Double_t e = p->Energy();
1689 Double_t vx = p->Vx();
1690 Double_t vy = p->Vy();
1691 Double_t vz = p->Vz();
1692 Double_t tof = p->T();
1695 p->GetPolarisation(pol);
1697 Double_t polx = pol.X();
1698 Double_t poly = pol.Y();
1699 Double_t polz = pol.Z();
1700 TMCProcess mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
1701 Float_t weight = p->GetWeight();
1703 AliGenerator::PushTrack(fTrackIt, -1, pdg, px, py, pz, e, vx, vy, vz, tof,polx, poly, polz, mech, ntr, weight);
1705 /***********************************************************/
1707 void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr, AliStack* stack) const
1709 //Shortcut to SetTrack(bla,bla,bla,bla.............)
1712 Error("SetTrack(TParticle*,Int_t&,AliStack*)","Particle has zero momentum");
1716 Int_t pdg = p->GetPdgCode();
1717 Double_t px = p->Px();
1718 Double_t py = p->Py();
1719 Double_t pz = p->Pz();
1720 Double_t e = p->Energy();
1722 stack->PushTrack(fTrackIt, -1, pdg, px, py, pz, e, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, kPPrimary, ntr,1,0);
1724 /***********************************************************/
1726 void AliGenHBTosl::Rotate(TVector3& relvector, TVector3& vector)
1728 //This method rotates vector about the angeles that are needed to rotate
1729 //relvector from postion (firstPx,0,0) to its actual positon
1730 //In other words: To make equations easier
1732 static TVector3 first;
1733 if (AliDebugLevel()>=1)
1735 first.SetXYZ(relvector.x(),relvector.y(),relvector.z());
1738 Double_t firstPx = TMath::Sqrt( relvector.x()*relvector.x() +
1739 relvector.y()*relvector.y() +
1740 relvector.z()*relvector.z() );
1742 Double_t rotAngleZ = -TMath::ATan2(relvector.y(),relvector.x());//calculating rot angles
1743 relvector.RotateZ(rotAngleZ);
1744 rotAngleZ = -rotAngleZ;
1745 Double_t rotAngleY = -TMath::ATan2(relvector.z(),relvector.x());
1747 vector.RotateY(rotAngleY);
1748 vector.RotateZ(rotAngleZ);
1750 if (AliDebugLevel()>5)
1752 TVector3 test(firstPx,0.0,0.0);
1753 test.RotateY(rotAngleY);
1754 test.RotateZ(rotAngleZ);
1755 AliInfo(Form("Rotation test: px %f %f",first.x(),test.x()));
1756 AliInfo(Form("Rotation test: py %f %f",first.y(),test.y()));
1757 AliInfo(Form("Rotation test: pz %f %f",first.z(),test.z()));
1760 /***********************************************************/
1762 Double_t AliGenHBTosl::Rotate(Double_t x,Double_t y,Double_t z)
1764 //Rotates vector to base where only x - coordinate is no-zero, and returns that
1766 Double_t xylength = TMath::Hypot(x,y);
1767 Double_t sinphi = -y/xylength;
1768 Double_t cosphi = x/xylength;
1770 Double_t xprime = cosphi*x - sinphi*y;
1771 Double_t yprime = sinphi*x + cosphi*y;
1774 Double_t a1 = -TMath::ATan2(v.Y(),v.X());
1776 if (AliDebugLevel()>5)
1778 AliInfo(Form("Xpr = %f Ypr = %f",xprime,yprime));
1779 AliInfo(Form("Calc sin = %f, and %f",sinphi,TMath::Sin(a1)));
1780 AliInfo(Form("Calc cos = %f, and %f",cosphi,TMath::Cos(a1)));
1783 Double_t xprimezlength = TMath::Hypot(xprime,z);
1785 Double_t sintheta = z/xprimezlength;
1786 Double_t costheta = xprime/xprimezlength;
1789 Double_t xbis = sintheta*z + costheta*(cosphi*x - sinphi*y);
1791 AliInfo(Form("Calculated rot %f, modulus %f",xbis,TMath::Sqrt(x*x+y*y+z*z)));
1794 /***********************************************************/
1796 AliStack* AliGenHBTosl::RotateStack()
1798 //swaps to next stack last goes to first and is reseted
1801 if ( fStackBuffer->GetSize() >= fBufferSize )
1803 stack = (AliStack*)fStackBuffer->Remove(fStackBuffer->Last());
1807 stack = new AliStack(fNpart);
1810 fStackBuffer->AddFirst(stack);
1814 /***********************************************************/
1816 Double_t AliGenHBTosl::GetQInvCorrTheorValue(Double_t qinv) const
1818 //Function (deprecated)
1819 static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1821 return 1.0 + 0.5*TMath::Exp(-qinv*qinv*fQRadius*fQRadius/kFactorsqrd);
1823 /***********************************************************/
1825 Double_t AliGenHBTosl::GetQOutQSideQLongCorrTheorValue(Double_t& out, Double_t& side, Double_t& lon) const
1827 //Theoretical function. Wa want to get correlation of the shape of this function
1828 static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1829 return 1.0 + 0.7*TMath::Exp(-fQRadius*fQRadius*(out*out+side*side+lon*lon)/kFactorsqrd);
1831 /***********************************************************/
1833 Bool_t AliGenHBTosl::CheckParticle(TParticle* p, TParticle* aupair ,AliStack* stack)
1835 //Checks if a given particle is falling into signal region with any other particle
1836 //already existing on stack
1839 if (fSignalRegion <=0) return kFALSE;
1841 for (Int_t i = 0; i < stack->GetNtrack(); i++)
1843 TParticle* part = stack->Particle(i);
1844 if (part == aupair) continue;
1845 Double_t qout = 10e5;
1846 Double_t qside= 10e5;
1847 Double_t qlong= 10e5;
1848 GetQOutQSideQLong(p,part,qout,qside,qlong);
1850 if (TMath::Abs(qout) < fSignalRegion)
1851 if (TMath::Abs(qside) < fSignalRegion)
1852 if (TMath::Abs(qlong) < fSignalRegion)
1857 /***********************************************************/
1859 void AliGenHBTosl::SwapGeneratingHistograms()
1861 //Checks if it is time to swap signal and background histograms
1862 //if yes it swaps them
1863 Int_t threshold = fMinFill;
1864 for (Int_t k = 1; k<=fQNBins; k++)
1866 for (Int_t j = 1; j<=fQNBins; j++)
1868 for (Int_t i = 1; i<=fQNBins; i++)
1870 if ( fQSecondBackground->GetBinContent(i,j,k) < threshold) return;
1877 Info("SwapGeneratingHistograms","*******************************************");
1878 Info("SwapGeneratingHistograms","*******************************************");
1879 Info("SwapGeneratingHistograms","*******************************************");
1880 Info("SwapGeneratingHistograms","**** SWAPPING HISTOGRAMS ****");
1881 Info("SwapGeneratingHistograms","*******************************************");
1882 Info("SwapGeneratingHistograms","*******************************************");
1883 Info("SwapGeneratingHistograms","*******************************************");
1887 fQSignal = fQSecondSignal;
1889 fQSecondSignal->Reset();
1890 fQSecondSignal->SetDirectory(0x0);
1893 fQBackground = fQSecondBackground;
1894 fQSecondBackground = h;
1895 fQSecondBackground->Reset();
1896 fQSecondBackground->SetDirectory(0x0);