StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtThPair.cc
1 /***************************************************************************
2  *
3  *
4  *
5  * Author: Laurent Conin, Fabrice Retiere, Subatech, France
6  ***************************************************************************
7  *
8  * Description : implementation of StHbtThPair
9  *
10  ***************************************************************************
11  *
12  *
13  *
14  ***************************************************************************/
15 
16 #include <Stiostream.h>
17 #include <Stsstream.h>
18 #include <stdlib.h>
19 
20 #include "StHbtMaker/Base/StHbtThPair.hh"
21 #include "StHbtMaker/Base/StHbtFsiWeight.hh"
22 
23 #ifdef __ROOT__
24 ClassImp(StHbtThPair)
25 #endif
26 
27 StHbtThPair::StHbtThPair(){
28  mMomentum1=mMomentum2=mEmPoint1=mEmPoint2=0;
29  mPid1=mPid2=0;
30  mMeasPair=0;
31  mWeight=0;
32  mWeightNum=mWeightDen=1.;
33  mWeightOk=false;
34  mPairPurity=1;
35 };
36 
37 void StHbtThPair::setVariables(const StHbtPair*){
38  /* no-op */
39 }
40 
41 void StHbtThPair::setMomRes1(int aPid){
42  if(aPid==7 || aPid==8){// pion
43  mMomRes1=1;
44  mPLoss1[0]=-0.00275; // (PRec-Pmc)P = [0]+[1]*P^[2]
45  mPLoss1[1]=-0.000335;
46  mPLoss1[2]=-2.052;
47  mPtRes1[0]=0.0178; // Dpt/pt = [0] + [1]*Pt^[2]
48  mPtRes1[1]=7.95e-5;
49  mPtRes1[2]=-2.273;
50  mPhiRes1[0]=0.0642; // DPhi = [0] +[1]*P^[2]
51  mPhiRes1[1]=0.0322;
52  mPhiRes1[2]=-1.511;
53  mThetaRes1[0]=0.0779; // DTheta = [0] + P^[2]
54  mThetaRes1[1]=0.0445;
55  mThetaRes1[2]=-1.528;
56  }
57  if(aPid==11 || aPid==12){//kaon
58  mMomRes1=1;
59  }
60 }
61 void StHbtThPair::setMomRes2(int aPid){
62  if(aPid==7 || aPid==8){// pion
63  mMomRes2=1;
64  mPLoss2[0]=-0.00275; // (PRec-Pmc)P = [0]+[1]*P^[2]
65  mPLoss2[1]=-0.000335;
66  mPLoss2[2]=-2.052;
67  mPtRes2[0]=0.0178; // Dpt/pt = [0] + [1]*Pt^[2]
68  mPtRes2[1]=7.95e-5;
69  mPtRes2[2]=-2.273;
70  mPhiRes2[0]=0.0642; // DPhi = [0] +[1]*P^[2]
71  mPhiRes2[1]=0.0322;
72  mPhiRes2[2]=-1.511;
73  mThetaRes2[0]=0.0779; // DTheta = [0] + P^[2]
74  mThetaRes2[1]=0.0445;
75  mThetaRes2[2]=-1.528;
76  }
77  if(aPid==11 || aPid==12){//kaon
78  mMomRes2=1;
79  }
80 }
81 
82 
83 void StHbtThPair::UpdateWeight() {
84  if (mWeight) {
85  mWeightNum=mWeight->GetWeight(this);
86  mWeightNum=(mWeightNum-1.)*mPairPurity+1.;
87  mWeightDen=mWeight->GetWeightDen();
88  if(mWeightNum<=0. || mWeightNum>1000.){
89  std::ostringstream tCom;
90  tCom << "echo ---> " << mWeightNum << " "
91  << mEmPoint1 << " " << mEmPoint2 << " "
92  << mMomentum1 << " " << mMomentum2 << " >> Err.txt" << ends;
93  system(tCom.str().c_str());
94  mWeightNum=1;
95  }
96  } else {
97  cout << "StHbtThPair Error - No Weight Generator plugged - set to 1 " <<endl;
98  mWeightNum=mWeightDen=1.;
99  }
100  mWeightOk=true;
101 }
102 
103 StHbtString StHbtThPair::Report() {
104  std::ostringstream tStr;
105  tStr << "Default StHbtThPair Report" << endl;
106  if (mWeight) {
107  tStr << mWeight->Report() << endl;
108  } else {
109  tStr << "No Weight Generator plugged - Weight set to 1 " << endl;
110  }
111  StHbtString returnThis = tStr.str();
112  return returnThis;
113 }
114 
115 double StHbtThPair::RealqSideCMS() const {
116  double x1 = mMomentum1->x(); double y1 = mMomentum1->y();
117  double x2 = mMomentum2->x(); double y2 = mMomentum2->y();
118 
119  double xt = x1+x2; double yt = y1+y2;
120  double k1 = ::sqrt(xt*xt+yt*yt);
121 
122  double tmp = 2.0*(x1*y2-x2*y1)/k1;
123 
124  return (tmp);
125 }
126 double StHbtThPair::RealqOutCMS() const {
127  double dx = mMomentum1->x() - mMomentum2->x();
128  double xt = mMomentum1->x() + mMomentum2->x();
129 
130  double dy = mMomentum1->y() - mMomentum2->y();
131  double yt = mMomentum1->y() + mMomentum2->y();
132 
133  double k1 = (::sqrt(xt*xt+yt*yt));
134  double k2 = (dx*xt+dy*yt);
135  double tmp = k2/k1;
136  return (tmp);
137 }
138 double StHbtThPair::RealqLongCMS() const {
139  double dz = mMomentum1->z() - mMomentum2->z();
140  double zz = mMomentum1->z() + mMomentum2->z();
141 
142  double dt = mMomentum1->t() - mMomentum2->t();
143  double tt = mMomentum1->t() + mMomentum2->t();
144 
145  double beta = zz/tt;
146  double gamma = 1.0/::sqrt(1.0 - beta*beta);
147 
148  double temp = gamma*(dz - beta*dt);
149  return (temp);
150 }
151 
152 //________________________________
153 double StHbtThPair::RealqOutPf() const
154 {
155  double dt = mMomentum1->t() - mMomentum2->t();
156  double tt = mMomentum1->t() + mMomentum2->t();
157 
158  double xt = mMomentum1->x() + mMomentum2->x();
159  double yt = mMomentum1->y() + mMomentum2->y();
160 
161  double k1 = ::sqrt(xt*xt + yt*yt);
162  double bOut = k1/tt;
163  double gOut = 1.0/::sqrt(1.0 - bOut*bOut);
164 
165  double temp = gOut*(this->RealqOutCMS() - bOut*dt);
166  return (temp);
167 }
168 
169 //___________________________________
170 double StHbtThPair::RealqSidePf() const
171 {
172  return(this->RealqSideCMS());
173 }
174 
175 //___________________________________
176 
177 double StHbtThPair::RealqLongPf() const
178 {
179  return(this->RealqLongCMS());
180 }
181 
182 //___________________________________
183 
184 void StHbtThPair::calcMomParameters() const{ // fortran like function! faster?
185  mMomParCalculated=1;
186 
187  double tPx1(mMomentum1->px());
188  double tPy1(mMomentum1->py());
189  double tPz1(mMomentum1->pz());
190  double tE1(mMomentum1->e());
191  double tPx(tPx1+mMomentum2->px());
192  double tPy(tPy1+mMomentum2->py());
193  double tPz(tPz1+mMomentum2->pz());
194  double tE(tE1 +mMomentum2->e());
195  mPt = tPx*tPx + tPy*tPy;
196  double tP = tPz*tPz;
197  double tMt = tE*tE - tP;
198  tP += mPt;
199  tP = ::sqrt(tP);
200  double tM = ::sqrt(tMt - mPt);
201  tMt = ::sqrt(tMt);
202  mPt = ::sqrt(mPt);
203  mBetat = mPt/tMt;
204  mUt = mPt/tMt;
205 
206  // Boost to LCMS
207  double tBeta = tPz/tE;
208  double tGamma = tE/tMt;
209  mKStarLong = tGamma * (tPz1 - tBeta * tE1);
210  double tE1L = tGamma * (tE1 - tBeta * tPz1);
211 
212  // Rotate in transverse plane
213  mKStarOut = ( tPx1*tPx + tPy1*tPy)/mPt;
214  mKStarSide = (-tPx1*tPy + tPy1*tPx)/mPt;
215 
216  // Boost to pair cms
217  mKStarOut = tMt/tM * (mKStarOut - mPt/tMt * tE1L);
218 
219  mKStar = ::sqrt(mKStarOut*mKStarOut+mKStarSide*mKStarSide+
220  mKStarLong*mKStarLong);
221 
222  mCVK = (mKStarOut*mPt + mKStarLong*tPz)/mKStar/mCVK;
223  mKStarLong = (tPz>=0)* mKStarLong - (tPz<0)*mKStarLong;
224 
225 }
226 
227 
228 void StHbtThPair::calcPosParameters() const{ // fortran like function! faster?
229  mPosParCalculated=1;
230 
231  double tPx = mMomentum1->px()+mMomentum2->px();
232  double tPy = mMomentum1->py()+mMomentum2->py();
233  double tPz = mMomentum1->pz()+mMomentum2->pz();
234  double tE = mMomentum1->e()+mMomentum2->e();
235  double tPt = tPx*tPx + tPy*tPy;
236  double tMt = tE*tE - tPz*tPz;
237  double tM = ::sqrt(tMt - tPt);
238  tMt = ::sqrt(tMt);
239  tPt = ::sqrt(tPt);
240 
241  double tDX = mEmPoint1->x()-mEmPoint2->x();
242  double tDY = mEmPoint1->y()-mEmPoint2->y();
243  mRLong = mEmPoint1->z()-mEmPoint2->z();
244  mDTime = mEmPoint1->t()-mEmPoint2->t();
245 
246  mRTrans = tDX>0.? ::sqrt(tDX*tDX+tDY*tDY) : -1.*::sqrt(tDX*tDX+tDY*tDY);
247  mROut = (tDX*tPx + tDY*tPy)/tPt;
248  mRSide = (-tDX*tPy + tDY*tPx)/tPt;
249  mRSidePairCMS = mRSide;
250 
251  // Lab -> LCMS -> PRF method
252  double tBeta = tPz/tE;
253  double tGamma = tE/tMt;
254  mRLongPairCMS = tGamma*(mRLong - tBeta* mDTime);
255  mDTimePairLCMS = tGamma*(mDTime - tBeta* mRLong);
256  tBeta = tPt/tMt;
257  tGamma = tMt/tM;
258  mROutPairCMS = tGamma*(mROut - tBeta* mDTimePairLCMS);
259  mDTimePairCMS = tGamma*(mDTimePairLCMS - tBeta* mROut);
260  mRStar = ::sqrt(mROutPairCMS*mROutPairCMS + mRSidePairCMS*mRSidePairCMS +
261  mRLongPairCMS*mRLongPairCMS);
262 
263 
264  /*
265  // Lab -> PRF method
266  double tRS12 = DX*Px+DY*Py+mRLong*Pz;
267  double tH1 = (tRS12/(E+M) - mDTime) / M;
268  double DXPairCMS = DX + Px* tH1;
269  double DYPairCMS = DY + Py* tH1;
270  mRLongPairCMS = mRLong + Pz* tH1;
271  //mRLong = Pz*mRLong>0? mRLong : -mRLong;
272  //mRLongPairCMS = Pz*mRLongPairCMS>0? mRLongPairCMS : -mRLongPairCMS;
273  mDTimePairCMS = (E* mDTime - tRS12) / M;
274  mRStar = ::sqrt(DXPairCMS*DXPairCMS + DYPairCMS*DYPairCMS +
275  mRLongPairCMS * mRLongPairCMS);
276  mROutPairCMS = (DXPairCMS*Px + DYPairCMS*Py)/Ptrans;
277  mRSidePairCMS = (-DXPairCMS*Py + DYPairCMS*Px)/Ptrans;
278  */
279 }
Generated by   doxygen 1.8.5