StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtThPairDoubleGauss.cxx
1 /***************************************************************************
2  *
3  *
4  *
5  * Author: Laurent Conin, Fabrice Retiere, Subatech, France
6  * Adam Kisiel, Warsaw University of Technology
7  ***************************************************************************
8  *
9  * Description : implementation of StHbtThPAirGAuss
10  *
11  ***************************************************************************
12  *
13  *
14  *
15  ***************************************************************************/
16 
17 #include "StHbtMaker/ThCorrFctn/StHbtThPairDoubleGauss.h"
18 #include "StHbtMaker/Infrastructure/StHbtParticle.hh"
19 #include "StHbtMaker/Base/StHbtFsiWeight.hh"
20 
21 #ifdef __ROOT__
22 ClassImp(StHbtThPairDoubleGauss)
23 #endif
24 
25 StHbtThPairDoubleGauss::StHbtThPairDoubleGauss() : StHbtThPair(){
26  mEmPoint1=&mPos1;
27  mEmPoint2=&mPos2;
28  mUseHidMom=false;
29  mUseHidPid=false;
30  mMomentum1=&mMom1;
31  mMomentum2=&mMom2;
32  mRef=RCMSDG;
33  mSizeX1=mSizeY1=mSizeZ1=1.;
34  mTime1=0;
35  mSizeX2=mSizeY2=mSizeZ2=1.;
36  mTime2=0;
37  mProb1=0.5;
38  mBetaRCMS=0.;
39  mGammaRCMS=1.;
40  mCoreHalo=1;
41  mXShift = 0.0;
42  mYShift = 0.0;
43  mZShift = 0.0;
44  mTShift = 0.0;
45  mRand = *(new TRandom2(31811));
46 
47  Int_t nbins = 100;
48  Double_t posMin = -10.0;
49  Double_t posMax = 10.0;
50  Double_t ptMin = 0.0;
51  Double_t ptMax = 1.0;
52 
53  mPosDist1 = new StHbt3DHisto("PosDist1","PosDist1",nbins,posMin,posMax,nbins,posMin,posMax,nbins,posMin,posMax);
54  mPosDist2 = new StHbt3DHisto("PosDist2","PosDist2",nbins,posMin,posMax,nbins,posMin,posMax,nbins,posMin,posMax);
55  mTDist1 = new StHbt1DHisto("TDist1","PosDist1",nbins,posMin,posMax);
56  mTDist2 = new StHbt1DHisto("TDist2","PosDist2",nbins,posMin,posMax);
57  mPosPtDist1 = new StHbt2DHisto("PosPtDist1","PosPtDist1",nbins,posMin,posMax,nbins,ptMin,ptMax);
58  mPosPtDist2 = new StHbt2DHisto("PosPtDist2","PosPtDist2",nbins,posMin,posMax,nbins,ptMin,ptMax);
59 };
60 
61 StHbtThPairDoubleGauss::~StHbtThPairDoubleGauss()
62 {
63  delete mPosDist1;
64  delete mPosDist2;
65  delete mTDist1;
66  delete mTDist2;
67  delete mPosPtDist1;
68  delete mPosPtDist2;
69 };
70 
71 void StHbtThPairDoubleGauss::Set(const StHbtPair* aPair){
72  SetMomentum_PID (aPair);
73  SetPosition(aPair);
74  BoostPosition();
75 
76  mMomParCalculated=0;
77  mPosParCalculated=0;
78 
79  mMeasPair=aPair;
80  mWeightOk=false;
81 
82 }
83 
84 void StHbtThPairDoubleGauss::SetMomentum_PID( const StHbtPair* aPair ){
85 
86  StHbtSmearedHiddenInfo* tSmearedHidInf1=0;
87  StHbtSmearedHiddenInfo* tSmearedHidInf2=0;
88  const StHbtEvtGenHiddenInfo* tEvtGenHidInf1=0;
89  const StHbtEvtGenHiddenInfo* tEvtGenHidInf2=0;
90  StHbtShiftedHiddenInfo* tShiftedHidInf1=0;
91  StHbtShiftedHiddenInfo* tShiftedHidInf2=0;
92 
93  if (mUseHidMom||mUseHidPid) {
94  if (mHiddenInfoType == SMEAR) {
95  tSmearedHidInf1=dynamic_cast< StHbtSmearedHiddenInfo*>
96  (aPair->track1()->getHiddenInfo());
97  tSmearedHidInf2=dynamic_cast< StHbtSmearedHiddenInfo*>
98  (aPair->track2()->getHiddenInfo());
99  if ((tSmearedHidInf1==0)||(tSmearedHidInf2==0)) {
100  if (tSmearedHidInf1==0) {
101  StHbtMomRes *mMomRes1 = new StHbtMomRes(mPid1);
102  mMomRes1->setMult(mResMult);
103  StHbtLorentzVector *temp;
104  temp = GenerateFreezeOut(1);
105  aPair->track1()->SetHiddenInfo(new StHbtSmearedHiddenInfo(aPair->track1()->FourMomentum(), *temp, mPid1, &mRand, mMomRes1));
106  tSmearedHidInf1 = dynamic_cast< StHbtSmearedHiddenInfo*>(aPair->track1()->getHiddenInfo());
107  delete mMomRes1;
108  delete temp;
109  }
110  if (tSmearedHidInf2==0) {
111  StHbtMomRes *mMomRes2 = new StHbtMomRes(mPid2);
112  mMomRes2->setMult(mResMult);
113  StHbtLorentzVector *temp;
114  temp = GenerateFreezeOut(2);
115  aPair->track2()->SetHiddenInfo(new StHbtSmearedHiddenInfo(aPair->track2()->FourMomentum(), *temp, mPid2, &mRand, mMomRes2));
116  tSmearedHidInf2 = dynamic_cast< StHbtSmearedHiddenInfo*>(aPair->track2()->getHiddenInfo());
117  delete mMomRes2;
118  delete temp;
119  }
120  }
121  }
122  else if (mHiddenInfoType == SHIFT) {
123  tShiftedHidInf1=dynamic_cast< StHbtShiftedHiddenInfo*>
124  (aPair->track1()->getHiddenInfo());
125  tShiftedHidInf2=dynamic_cast< StHbtShiftedHiddenInfo*>
126  (aPair->track2()->getHiddenInfo());
127  if ((tShiftedHidInf1==0)||(tShiftedHidInf2==0)) {
128  if (tShiftedHidInf1==0) {
129  StHbtMomRes *mMomRes1 = new StHbtMomRes(mPid1);
130  mMomRes1->setMult(mResMult);
131  aPair->track1()->SetHiddenInfo(new StHbtShiftedHiddenInfo(aPair->track1()->FourMomentum(), mPid1, &mRand, mMomRes1, mShift, PHISHIFT));
132  tShiftedHidInf1 = dynamic_cast< StHbtShiftedHiddenInfo*>(aPair->track1()->getHiddenInfo());
133  delete mMomRes1;
134  }
135  if (tShiftedHidInf2==0) {
136  StHbtMomRes *mMomRes2 = new StHbtMomRes(mPid2);
137  mMomRes2->setMult(mResMult);
138  aPair->track2()->SetHiddenInfo(new StHbtShiftedHiddenInfo(aPair->track2()->FourMomentum(), mPid2, &mRand, mMomRes2, mShift, PHISHIFT));
139  tShiftedHidInf2 = dynamic_cast< StHbtShiftedHiddenInfo*>(aPair->track2()->getHiddenInfo());
140  delete mMomRes2;
141  }
142  }
143  }
144  else {
145  tEvtGenHidInf1=dynamic_cast<const StHbtEvtGenHiddenInfo*>
146  (aPair->track1()->HiddenInfo());
147  tEvtGenHidInf2=dynamic_cast<const StHbtEvtGenHiddenInfo*>
148  (aPair->track2()->HiddenInfo());
149  }
150  }
151 
152  if (mUseHidMom&&tSmearedHidInf2&&tSmearedHidInf1) {
153  // cout << "Setting Smeared HiddenMomentum " << endl;
154  mMomentum1=&(tSmearedHidInf1->mSmearedMom);
155  mMomentum2=&(tSmearedHidInf2->mSmearedMom);
156  }
157  else if (mUseHidMom&&tShiftedHidInf2&&tShiftedHidInf1) {
158  mMomentum1=&(tShiftedHidInf1->mShiftedMom);
159  mMomentum2=&(tShiftedHidInf2->mShiftedMom);
160  }
161  else if (mUseHidMom&&tEvtGenHidInf2&&tEvtGenHidInf1) {
162  // mMomentum1=&(tEvtGenHidInf1->getFreezeOutMomEn());
163  // mMomentum2=&(tEvtGenHidInf2->getFreezeOutMomEn());
164  mMomentum1= new StHbtLorentzVector(*tEvtGenHidInf1->getFreezeOutMomEn());
165  mMomentum2= new StHbtLorentzVector(*tEvtGenHidInf2->getFreezeOutMomEn());
166  }
167  else{
168  mMom1=aPair->track1()->FourMomentum();
169  mMom1.setE(::sqrt(mMassSq1+mMom1.vect().mag2()));
170  mMom2=aPair->track2()->FourMomentum();
171  mMom2.setE(::sqrt(mMassSq2+mMom2.vect().mag2()));
172  if ((!mUseHidMom)&&(mBetaRCMS>0)) {
173  double tE=mMom1.e();
174  mMom1.setE(mGammaRCMS*(tE-mBetaRCMS*mMom1.pz()));
175  mMom1.setPz(mGammaRCMS*(mMom1.pz()-mBetaRCMS*tE));
176  tE=mMom2.e();
177  mMom2.setE(mGammaRCMS*(tE-mBetaRCMS*mMom2.pz()));
178  mMom2.setPz(mGammaRCMS*(mMom2.pz()-mBetaRCMS*tE));
179  }
180  }
181 }
182 
183 void StHbtThPairDoubleGauss::SetPosition( const StHbtPair* aPair) {
184  /*
185  * Set the random freeze-out space coordinates
186  * according to either of two models:
187  * Core-Halo - the coordinates are randomly generated
188  * either with core (1st) or halo (2nd) parametrization
189  * TwoSources - The first particle coordinates are generated
190  * from the 1st set of parameters, coordinates for te second -
191  * from the second set.
192  */
193 
194  const StHbtSmearedHiddenInfo* tSmearedHidInf1=0;
195  const StHbtSmearedHiddenInfo* tSmearedHidInf2=0;
196 
197  if (mHiddenInfoType == SMEAR) {
198  tSmearedHidInf1 = dynamic_cast<const StHbtSmearedHiddenInfo*>(aPair->track1()->HiddenInfo());
199  tSmearedHidInf2 = dynamic_cast<const StHbtSmearedHiddenInfo*>(aPair->track2()->HiddenInfo());
200  }
201 
202  if ((mHiddenInfoType == SMEAR) && (tSmearedHidInf1 != 0) && (tSmearedHidInf2 != 0))
203  {
204  // cout << "Setting position for pair: " << this << endl;
205  mPos1 = tSmearedHidInf1->getFreezeOut();
206  mPos2 = tSmearedHidInf2->getFreezeOut();
207  // cout << "Set the position 1 to: X: " << mPos1.x() << " Y: " << mPos1.y() << " Z: " << mPos1.z() << " T: " << mPos1.t() << endl;
208  // cout << "Set the position 2 to: X: " << mPos2.x() << " Y: " << mPos2.y() << " Z: " << mPos2.z() << " T: " << mPos2.t() << endl;
209 
210  }
211  else
212  {
213  StHbtLorentzVector* temp;
214  temp = GenerateFreezeOut(1);
215  mPos1 = *temp;
216  delete temp;
217  temp = GenerateFreezeOut(2);
218  mPos2 = *temp;
219  delete temp;
220  // cout << mPos1.x() << " " << mPos1.y() << endl;
221  }
222 
223 };
224 
225 StHbtLorentzVector* StHbtThPairDoubleGauss::GenerateFreezeOut(int partno) {
226 
227  StHbtLorentzVector* mPos;
228  mPos = new StHbtLorentzVector();
229  float x,y,z,t;
230 
231  // cout << "Generating Freeze-out! " << partno << endl;
232  mRand.Rannor(x,y);
233  mRand.Rannor(z,t);
234  if (mCoreHalo == 1) {
235  if (mRand.Uniform() < mProb1) {
236  mPos->setX(x*mSizeX1);
237  mPos->setY(y*mSizeY1);
238  mPos->setZ(z*mSizeZ1);
239  mPos->setT(t*mTime1);
240  }
241  else {
242  mPos->setX(x*mSizeX2);
243  mPos->setY(y*mSizeY2);
244  mPos->setZ(z*mSizeZ2);
245  mPos->setT(t*mTime2);
246  }
247  }
248  else if (mCoreHalo ==2)
249  {
250  Double_t phi;
251  phi = mRand.Rndm();
252 
253  mPos->setX(x*mSizeX1*cos(phi*2*TMath::Pi()));
254  mPos->setY(x*mSizeX1*sin(phi*2*TMath::Pi()));
255  mPos->setZ(z*mSizeZ2);
256  mPos->setT(t*mTime2);
257  }
258  else if (mCoreHalo ==3)
259  {
260  if (partno == 1)
261  {
262  mPos->setX(0.0);
263  mPos->setY(0.0);
264  mPos->setZ(0.0);
265  mPos->setT(0.0);
266  }
267  else
268  {
269  double tROut = mRand.Gaus(mTime1,mSizeX1); // reuse of long
270  double tRSide = mRand.Gaus(0,mSizeX1);
271  double tRLong = mRand.Gaus(0,mSizeX1);
272  double tDTime = tROut;
273 
274  double tPx = mMomentum1->x()+mMomentum2->x();
275  double tPy = mMomentum1->y()+mMomentum2->y();
276  double tPz = mMomentum1->z()+mMomentum2->z();
277  double tE = mMomentum1->e()+mMomentum2->e();
278 
279  double tPt = tPx*tPx+tPy*tPy;
280  double tMt = tE*tE-tPz*tPz;
281  double tM = ::sqrt(tMt - tPt);
282  tPt = ::sqrt(tPt);
283  tMt = ::sqrt(tMt);
284 
285  tROut *= (tMt/tM); // Rout*gammaT
286  tDTime *= (tPt/tM); // Rout*betaT*gammaT
287  double ttDTime = tDTime;
288  tDTime += (tPz/tE*tRLong);
289  tDTime *= (tE/tMt);
290  tRLong += (tPz/tE*ttDTime);
291  tRLong *= (tE/tMt);
292 
293  tPx /= tPt;
294  tPy /= tPt;
295 
296  mPos->setX(tROut*tPx-tRSide*tPy);
297  mPos->setY(tROut*tPy+tRSide*tPx);
298  mPos->setZ(tRLong);
299  mPos->setT(tDTime);
300 
301  // mPos->setX(x*mSizeX1*cos(phi*2*TMath::Pi()));
302  // mPos->setY(x*mSizeX1*sin(phi*2*TMath::Pi()));
303  // mPos->setZ(z*mSizeZ2);
304  // mPos->setT(t*mTime2);
305  }
306  }
307  else
308  {
309  if (partno == 1)
310  {
311  mPos->setX(x*mSizeX1);
312  mPos->setY(y*mSizeY1);
313  mPos->setZ(z*mSizeZ1);
314  mPos->setT(t*mTime1);
315  // mPos->setT(mTime1);
316  }
317  else
318  {
319  mPos->setX(mXShift+x*mSizeX2);
320  mPos->setY(mYShift+y*mSizeY2);
321  mPos->setZ(mZShift+z*mSizeZ2);
322  mPos->setT(mTShift+t*mTime2);
323  // mPos->setT(mTime2);
324  }
325  }
326  if (partno == 1)
327  {
328  mTDist1->Fill(mPos->t());
329  mPosDist1->Fill(mPos->x(),mPos->y(),mPos->z());
330  // mPosPtDist1->Fill(hypot(mPos->x(),mPos->y()), hypot(partMom->x(),partMom->y()));
331  }
332  else
333  {
334  mTDist2->Fill(mPos->t());
335  mPosDist2->Fill(mPos->x(),mPos->y(),mPos->z());
336  // mPosPtDist2->Fill(hypot(mPos->x(),mPos->y()), hypot(partMom->y(),partMom->x()));
337  }
338  return mPos;
339 }
340 
341 void StHbtThPairDoubleGauss::BoostPosition(){
342  double tBeta;
343  double tGamma;
344  double tT;
345  switch (mRef) {
346  case RCMSDG: break;
347  case LCMSDG:
348  tBeta=(mMomentum1->pz()+mMomentum2->pz())/(mMomentum1->e()+mMomentum2->e());
349  tGamma=::sqrt(1/1-tBeta*tBeta);
350  tT=mPos1.t();
351  mPos1.setT(tGamma*(tT-tBeta*mPos1.z()));
352  mPos1.setZ(tGamma*(tBeta*mPos1.z()-tBeta*tT));
353  tT=mPos2.t();
354  mPos2.setT(tGamma*(tT-tBeta*mPos2.z()));
355  mPos2.setZ(tGamma*(mPos2.z()-tBeta*tT));
356  break;
357  case PRFDG:
358  StHbtLorentzVector tBoost=*mMomentum1+*mMomentum2;
359  mPos1=mPos1.boost(tBoost);
360  mPos2=mPos2.boost(tBoost);
361  break;
362  }
363 };
364 
365 inline void StHbtThPairDoubleGauss::SetResolutionMult(const double mult) {
366  mResMult = mult;
367 }
368 
369 inline void StHbtThPairDoubleGauss::SetMomentumShift(const double shift)
370 {
371  mShift = shift;
372 }
373 
374 inline void StHbtThPairDoubleGauss::UseSmearedHiddenInfo()
375 {
376  mHiddenInfoType = SMEAR;
377 }
378 
379 inline void StHbtThPairDoubleGauss::UseShiftedHiddenInfo()
380 {
381  mHiddenInfoType = SHIFT;
382 }
383 
384 inline void StHbtThPairDoubleGauss::UseEvtGenHiddenInfo()
385 {
386  mHiddenInfoType = EVTGEN;
387 }
388 
389 inline void StHbtThPairDoubleGauss::SetSizes(double aXYZ1, double aT1, double aXYZ2, double aT2)
390 {mSizeX1=mSizeY1=mSizeZ1=aXYZ1; mTime1=aT1; mSizeX2=mSizeY2=mSizeZ2=aXYZ2; mTime2=aT2; }
391 inline void StHbtThPairDoubleGauss::SetSizes(double aX1,double aY1,double aZ1, double aT1,double aX2,double aY2,double aZ2, double aT2)
392 {mSizeX1=aX1; mSizeY1=aY1; mSizeZ1=aZ1; mTime1=aT1; mSizeX2=aX2; mSizeY2=aY2; mSizeZ2=aZ2; mTime2=aT2; }
393 inline void StHbtThPairDoubleGauss::SetSize1(double aSize,double aTime) {mSizeX1=aSize;mSizeY1=aSize;mSizeZ1=aSize;mTime1=aTime;};
394 inline void StHbtThPairDoubleGauss::SetSize1(double aSizeX,double aSizeY, double aSizeZ,double aTime)
395 {mSizeX1=aSizeX;mSizeY1=aSizeY;mSizeZ1=aSizeZ;mTime1=aTime;};
396 inline void StHbtThPairDoubleGauss::SetSize2(double aSize,double aTime) {mSizeX2=aSize;mSizeY2=aSize;mSizeZ2=aSize;mTime2=aTime;};
397 inline void StHbtThPairDoubleGauss::SetSize2(double aSizeX,double aSizeY, double aSizeZ,double aTime)
398 {mSizeX2=aSizeX;mSizeY2=aSizeY;mSizeZ2=aSizeZ;mTime2=aTime;};
399 
400 inline void StHbtThPairDoubleGauss::UseHiddenMomentum(){mUseHidMom=1;};
401 inline void StHbtThPairDoubleGauss::UseParticleMomentum(){
402  mUseHidMom=0;
403  if (mUseHidPid){
404  mMomentum1=&mMom1;
405  mMomentum2=&mMom2;
406  }
407 };
408 
409 inline void StHbtThPairDoubleGauss::UseHiddenPid() {
410  mUseHidPid=true;
411  if (mUseHidMom){
412  mMomentum1=&mMom1;
413  mMomentum2=&mMom2;
414  }
415 };
416 
417 inline void StHbtThPairDoubleGauss::UseFixedPid( int const tPid1,double const tMass1,int const tPid2, double const tMass2) {
418  mUseHidPid=false;mPid1=tPid1;mPid2=tPid2;mMassSq1=tMass1*tMass1;mMassSq2=tMass2*tMass2; };
419 inline void StHbtThPairDoubleGauss::UseFixedPid( int const tPid1,double const tMass1) {
420  mUseHidPid=false;mPid1=tPid1;mPid2=tPid1;mMassSq1=tMass1*tMass1;mMassSq2=tMass1*tMass1; };
421 
422 inline void StHbtThPairDoubleGauss::SetRCMS() {mRef=RCMSDG;};
423 inline void StHbtThPairDoubleGauss::SetLCMS(){mRef=LCMSDG;};
424 inline void StHbtThPairDoubleGauss::SetPRF(){mRef=PRFDG;};
425 
426 inline void StHbtThPairDoubleGauss::SetBoostRCMS(double aPlab,double aMBeam, double aMTarget){
427  double tEBeamLab=::sqrt(aPlab*aPlab+aMBeam*aMBeam);
428  mGammaRCMS=(tEBeamLab+aMTarget)/::sqrt(aMBeam*aMBeam+aMTarget*aMTarget+2*tEBeamLab*aMTarget);
429  mBetaRCMS=::sqrt(1.-1/(mGammaRCMS*mGammaRCMS));
430 }
431 
432 inline void StHbtThPairDoubleGauss::SetFirstProb(double aProb1) {
433  if (aProb1<0.0) mProb1 = 0.0;
434  else if (aProb1>1.0) mProb1 = 1.0;
435  else mProb1 = aProb1;
436 }
437 
438 inline void StHbtThPairDoubleGauss::SetCoreHalo() { mCoreHalo=1; }
439 inline void StHbtThPairDoubleGauss::SetTwoSources() { mCoreHalo=0; }
440 inline void StHbtThPairDoubleGauss::SetRadialGaus() { mCoreHalo=2; }
441 inline void StHbtThPairDoubleGauss::SetPRFGaus() { mCoreHalo=3; }
442 inline void StHbtThPairDoubleGauss::SetPositionShift(double aX, double aY, double aZ, double aT)
443 {
444  mXShift = aX;
445  mYShift = aY;
446  mZShift = aZ;
447  mTShift = aT;
448 }
449 
450 void StHbtThPairDoubleGauss::Write()
451 {
452  mPosDist1->Write();
453  mPosDist2->Write();
454  mTDist1->Write();
455  mTDist2->Write();
456  mPosPtDist1->Write();
457  mPosPtDist2->Write();
458 
459 }
Generated by   doxygen 1.8.5