BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtVPHOtoVISRHi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 2004 Cornell
10//
11// Module: EvtVPHOtoVISR.cc
12//
13// Description: Routine to decay vpho -> (DDx) + ISR photon from 3.9 to 4.3 GeV, using CLEO-c
14// data (Brian Lang)
15//
16// Modification history:
17//
18// Ryd March 20, 2004 Module created
19//
20//------------------------------------------------------------------------
21//
22#include "EvtVPHOtoVISRHi.hh"
23#include "../EvtGenBase/EvtGenKine.hh"
24#include "../EvtGenBase/EvtPDL.hh"
25#include "../EvtGenBase/EvtParticle.hh"
26#include "../EvtGenBase/EvtRandom.hh"
27#include "../EvtGenBase/EvtReport.hh"
28#include "../EvtGenBase/EvtVector4C.hh"
29#include "../EvtGenBase/EvtVector4R.hh"
30#include <stdlib.h>
31#include <string>
32
33using std::endl;
34
35EvtVPHOtoVISRHi::~EvtVPHOtoVISRHi() {}
36
37void EvtVPHOtoVISRHi::getName( std::string& model_name ) { model_name = "VPHOTOVISRHI"; }
38
39EvtDecayBase* EvtVPHOtoVISRHi::clone() { return new EvtVPHOtoVISRHi; }
40
41void EvtVPHOtoVISRHi::init() {
42
43 // check that there are 0 or 1 arguments
44 checkNArg( 0, 1 );
45
46 // check that there are 2 daughters
47 checkNDaug( 2 );
48
49 // check the parent and daughter spins
50 checkSpinParent( EvtSpinType::VECTOR );
51 checkSpinDaughter( 0, EvtSpinType::VECTOR );
52 checkSpinDaughter( 1, EvtSpinType::PHOTON );
53}
54
55void EvtVPHOtoVISRHi::initProbMax() { setProbMax( 20.0 ); }
56
57void EvtVPHOtoVISRHi::decay( EvtParticle* p ) {
58 // take photon along z-axis, either forward or backward.
59 // Implement this as generating the photon momentum along
60 // the z-axis uniformly
61 double power = 1;
62 if ( getNArg() == 1 ) power = getArg( 0 );
63 // define particle names
64 static EvtId D0 = EvtPDL::getId( "D0" );
65 static EvtId D0B = EvtPDL::getId( "anti-D0" );
66 static EvtId DP = EvtPDL::getId( "D+" );
67 static EvtId DM = EvtPDL::getId( "D-" );
68 static EvtId DSM = EvtPDL::getId( "D_s-" );
69 static EvtId DSP = EvtPDL::getId( "D_s+" );
70 static EvtId DSMS = EvtPDL::getId( "D_s*-" );
71 static EvtId DSPS = EvtPDL::getId( "D_s*+" );
72 static EvtId D0S = EvtPDL::getId( "D*0" );
73 static EvtId D0BS = EvtPDL::getId( "anti-D*0" );
74 static EvtId DPS = EvtPDL::getId( "D*+" );
75 static EvtId DMS = EvtPDL::getId( "D*-" );
76 // setup some parameters
77 double w = p->mass();
78 double s = w * w;
79 double L = 2.0 * log( w / 0.000511 );
80 double alpha = 1 / 137.0;
81 double beta = ( L - 1 ) * 2.0 * alpha / EvtConst::pi;
82 // make sure only 2 or 3 body are present
83 assert( p->getDaug( 0 )->getNDaug() == 2 || p->getDaug( 0 )->getNDaug() == 3 );
84
85 // determine minimum rest mass of parent
86 double md1 = EvtPDL::getMeanMass( p->getDaug( 0 )->getDaug( 0 )->getId() );
87 double md2 = EvtPDL::getMeanMass( p->getDaug( 0 )->getDaug( 1 )->getId() );
88 double minResMass = md1 + md2;
89 if ( p->getDaug( 0 )->getNDaug() == 3 )
90 {
91 double md3 = EvtPDL::getMeanMass( p->getDaug( 0 )->getDaug( 2 )->getId() );
92 minResMass = minResMass + md3;
93 }
94
95 // calculate the maximum energy of the ISR photon
96 double pgmax = ( s - minResMass * minResMass ) / ( 2.0 * w );
97 double pgz = 0.99 * pgmax * exp( log( EvtRandom::Flat( 1.0 ) ) / ( beta * power ) );
98 if ( EvtRandom::Flat( 1.0 ) < 0.5 ) pgz = -pgz;
99
100 double k = fabs( pgz );
101 // print of ISR energy
102 // std::cout << "Energy ISR :"<< k <<std::endl;
103 EvtVector4R p4g( k, 0.0, 0.0, pgz );
104
105 EvtVector4R p4res = p->getP4Restframe() - p4g;
106
107 double mres = p4res.mass();
108
109 // set masses
110 p->getDaug( 0 )->init( getDaug( 0 ), p4res );
111 p->getDaug( 1 )->init( getDaug( 1 ), p4g );
112
113 // determine XS - langbw
114 // very crude way of determining XS just a simple straight line Approx.
115 // this was determined by eye.
116 // lots of cout statements to make plots to check that things are working as expected
117 double sigma = 9.0;
118 if ( mres <= 3.9 ) sigma = 0.00001;
119
120 bool sigmacomputed( false );
121
122 // DETERMINE XS FOR D*D*
123 if ( p->getDaug( 0 )->getNDaug() == 2 &&
124 ( ( p->getDaug( 0 )->getDaug( 0 )->getId() == D0S &&
125 p->getDaug( 0 )->getDaug( 1 )->getId() == D0BS ) ||
126 ( p->getDaug( 0 )->getDaug( 0 )->getId() == DPS &&
127 p->getDaug( 0 )->getDaug( 1 )->getId() == DMS ) ) )
128 {
129 if ( mres > 4.18 )
130 {
131 sigma *=
132 5. / 9. * ( 1. - 1. * sqrt( ( 4.18 - mres ) * ( 4.18 - mres ) ) / ( 4.3 - 4.18 ) );
133 }
134 else if ( mres > 4.07 && mres <= 4.18 ) { sigma *= 5. / 9.; }
135 else if ( mres <= 4.07 && mres > 4.03 )
136 {
137 sigma *=
138 ( 5. / 9. - 1.5 / 9. * sqrt( ( 4.07 - mres ) * ( 4.07 - mres ) ) / ( 4.07 - 4.03 ) );
139 }
140 else if ( mres <= 4.03 && mres >= 4.013 )
141 {
142 sigma *= ( 3.5 / 9. -
143 3.5 / 9. * sqrt( ( 4.03 - mres ) * ( 4.03 - mres ) ) / ( 4.03 - 4.013 ) );
144 }
145 else { sigma = 0.00001; }
146 sigmacomputed = true;
147 // std::cout << "DSDSXS "<<sigma<< " " << mres<<std::endl;
148 }
149
150 // DETERMINE XS FOR D*D
151 if ( p->getDaug( 0 )->getNDaug() == 2 &&
152 ( ( p->getDaug( 0 )->getDaug( 0 )->getId() == D0S &&
153 p->getDaug( 0 )->getDaug( 1 )->getId() == D0B ) ||
154 ( p->getDaug( 0 )->getDaug( 0 )->getId() == DPS &&
155 p->getDaug( 0 )->getDaug( 1 )->getId() == DM ) ||
156 ( p->getDaug( 0 )->getDaug( 0 )->getId() == D0BS &&
157 p->getDaug( 0 )->getDaug( 1 )->getId() == D0 ) ||
158 ( p->getDaug( 0 )->getDaug( 0 )->getId() == DMS &&
159 p->getDaug( 0 )->getDaug( 1 )->getId() == DP ) ) )
160 {
161 if ( mres >= 4.2 ) { sigma *= 1.5 / 9.; }
162 else if ( mres > 4.06 && mres < 4.2 )
163 {
164 sigma *= ( ( 1.5 / 9. +
165 2.5 / 9. * sqrt( ( 4.2 - mres ) * ( 4.2 - mres ) ) / ( 4.2 - 4.06 ) ) );
166 }
167 else if ( mres >= 4.015 && mres < 4.06 )
168 {
169 sigma *= ( ( 4. / 9. +
170 3. / 9. * sqrt( ( 4.06 - mres ) * ( 4.06 - mres ) ) / ( 4.06 - 4.015 ) ) );
171 }
172 else if ( mres < 4.015 && mres >= 3.9 )
173 {
174 sigma *= ( ( 7. / 9. -
175 7 / 9. * sqrt( ( 4.015 - mres ) * ( 4.015 - mres ) ) / ( 4.015 - 3.9 ) ) );
176 }
177 else { sigma = 0.00001; }
178 sigmacomputed = true;
179 // std::cout << "DSDXS "<<sigma<< " " << mres<<std::endl;
180 }
181
182 // DETERMINE XS FOR Ds*Ds*
183 if ( ( ( p->getDaug( 0 )->getDaug( 0 )->getId() == DSPS &&
184 p->getDaug( 0 )->getDaug( 1 )->getId() == DSMS ) ) )
185 {
186 if ( mres > ( 2.112 + 2.112 ) ) { sigma = 0.4; }
187 else
188 {
189 // sigma=0.4;
190 // sigma = 0 surely below Ds*Ds* threshold? - ponyisi
191 sigma = 0.00001;
192 }
193 sigmacomputed = true;
194 // std::cout << "DsSDsSXS "<<sigma<< " " << mres<<std::endl;
195 }
196
197 // DETERMINE XS FOR Ds*Ds
198 if ( p->getDaug( 0 )->getNDaug() == 2 &&
199 ( ( p->getDaug( 0 )->getDaug( 0 )->getId() == DSPS &&
200 p->getDaug( 0 )->getDaug( 1 )->getId() == DSM ) ||
201 ( p->getDaug( 0 )->getDaug( 0 )->getId() == DSMS &&
202 p->getDaug( 0 )->getDaug( 1 )->getId() == DSP ) ) )
203 {
204 if ( mres > 4.26 ) { sigma = 0.05; }
205 else if ( mres > 4.18 && mres <= 4.26 )
206 {
207 sigma *= 1. / 9. *
208 ( 0.05 + 0.95 * sqrt( ( 4.26 - mres ) * ( 4.26 - mres ) ) / ( 4.26 - 4.18 ) );
209 }
210 else if ( mres > 4.16 && mres <= 4.18 ) { sigma *= 1 / 9.; }
211 else if ( mres <= 4.16 && mres > 4.08 )
212 { sigma *= 1 / 9. * ( 1 - sqrt( ( 4.16 - mres ) * ( 4.16 - mres ) ) / ( 4.16 - 4.08 ) ); }
213 else if ( mres <= ( 4.08 ) ) { sigma = 0.00001; }
214 sigmacomputed = true;
215 // std::cout << "DsSDsXS "<<sigma<< " " << mres<<std::endl;
216 }
217
218 // DETERMINE XS FOR DD
219 if ( p->getDaug( 0 )->getNDaug() == 2 &&
220 ( ( p->getDaug( 0 )->getDaug( 0 )->getId() == D0 &&
221 p->getDaug( 0 )->getDaug( 1 )->getId() == D0B ) ||
222 ( p->getDaug( 0 )->getDaug( 0 )->getId() == DP &&
223 p->getDaug( 0 )->getDaug( 1 )->getId() == DM ) ) )
224 {
225 sigma *= 0.4 / 9.;
226 sigmacomputed = true;
227 // std::cout << "DDXS "<<sigma<< " " << mres<<std::endl;
228 }
229
230 // DETERMINE XS FOR DsDs
231 if ( p->getDaug( 0 )->getNDaug() == 2 &&
232 ( ( p->getDaug( 0 )->getDaug( 0 )->getId() == DSP &&
233 p->getDaug( 0 )->getDaug( 1 )->getId() == DSM ) ) )
234 {
235 sigma *= 0.2 / 9.;
236 sigmacomputed = true;
237 // std::cout << "DsDsXS "<<sigma<< " " << mres<<std::endl;
238 }
239
240 // DETERMINE XS FOR MULTIBODY
241 if ( p->getDaug( 0 )->getNDaug() == 3 )
242 {
243 if ( mres > 4.03 ) { sigma *= 0.5 / 9.; }
244 else { sigma = 0.00001; }
245 sigmacomputed = true;
246 // std::cout << "DSDpiXS "<<sigma<< " " << mres<<std::endl;
247 }
248
249 if ( !sigmacomputed )
250 {
251 report( ERROR, "EvtGen" )
252 << "VPHOTOVISRHI: This model requires daughters to be listed in a particular order."
253 << endl
254 << "The following are acceptable:" << endl
255 << "D0 anti-D0" << endl
256 << "D+ D-" << endl
257 << "D*0 anti-D0" << endl
258 << "anti-D*0 D0" << endl
259 << "D*+ D-" << endl
260 << "D*- D+" << endl
261 << "D*0 anti-D*0" << endl
262 << "D*+ D*-" << endl
263 << "D_s+ D_s-" << endl
264 << "D_s*+ D_s-" << endl
265 << "D_s*- D_s+" << endl
266 << "D_s*+ D_s*-" << endl
267 << "(D* D pi can be in any order)" << endl
268 << "Aborting..." << endl;
269 assert( 0 );
270 }
271
272 if ( sigma < 0 ) sigma = 0.0;
273
274 // static double sigmax=sigma;
275 // if (sigma>sigmax){
276 // sigmax=sigma;
277 // }
278
279 static int count = 0;
280
281 count++;
282
283 // if (count%10000==0){
284 // std::cout << "sigma :"<<sigma<<std::endl;
285 // std::cout << "sigmax:"<<sigmax<<std::endl;
286 // }
287
288 double norm = sqrt( sigma );
289
290 // EvtParticle* d=p->getDaug(0);
291
292 vertex( 0, 0, 0, norm * p->eps( 0 ) * p->epsParent( 0 ).conj() );
293 vertex( 1, 0, 0, norm * p->eps( 1 ) * p->epsParent( 0 ).conj() );
294 vertex( 2, 0, 0, norm * p->eps( 2 ) * p->epsParent( 0 ).conj() );
295
296 vertex( 0, 1, 0, norm * p->eps( 0 ) * p->epsParent( 1 ).conj() );
297 vertex( 1, 1, 0, norm * p->eps( 1 ) * p->epsParent( 1 ).conj() );
298 vertex( 2, 1, 0, norm * p->eps( 2 ) * p->epsParent( 1 ).conj() );
299
300 vertex( 0, 2, 0, norm * p->eps( 0 ) * p->epsParent( 2 ).conj() );
301 vertex( 1, 2, 0, norm * p->eps( 1 ) * p->epsParent( 2 ).conj() );
302 vertex( 2, 2, 0, norm * p->eps( 2 ) * p->epsParent( 2 ).conj() );
303
304 vertex( 0, 0, 1, norm * p->eps( 0 ) * p->epsParent( 0 ).conj() );
305 vertex( 1, 0, 1, norm * p->eps( 1 ) * p->epsParent( 0 ).conj() );
306 vertex( 2, 0, 1, norm * p->eps( 2 ) * p->epsParent( 0 ).conj() );
307
308 vertex( 0, 1, 1, norm * p->eps( 0 ) * p->epsParent( 1 ).conj() );
309 vertex( 1, 1, 1, norm * p->eps( 1 ) * p->epsParent( 1 ).conj() );
310 vertex( 2, 1, 1, norm * p->eps( 2 ) * p->epsParent( 1 ).conj() );
311
312 vertex( 0, 2, 1, norm * p->eps( 0 ) * p->epsParent( 2 ).conj() );
313 vertex( 1, 2, 1, norm * p->eps( 1 ) * p->epsParent( 2 ).conj() );
314 vertex( 2, 2, 1, norm * p->eps( 2 ) * p->epsParent( 2 ).conj() );
315
316 return;
317}
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
count
DOUBLE_PRECISION count[3]
Definition EvtPhokharaDef.hh:63
alpha
double alpha
Definition EvtPhokharaDef.hh:17
w
double w
Definition EvtPhokharaDef.hh:28
report
ostream & report(Severity severity, const char *facility)
Definition EvtReport.cc:34
ERROR
@ ERROR
Definition EvtReport.hh:49
s
XmlRpcServer s
Definition HelloServer.cpp:11
EvtConst::pi
static const double pi
Definition EvtConst.hh:27
EvtDecayAmp::vertex
void vertex(const EvtComplex &amp)
Definition EvtDecayAmp.hh:36
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:492
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:479
EvtDecayBase::getArg
double getArg(int j)
Definition EvtDecayBase.cc:520
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:261
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:465
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:449
EvtDecayBase::getDaug
EvtId getDaug(int i)
Definition EvtDecayBase.hh:69
EvtId
Definition EvtId.hh:27
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition EvtPDL.hh:43
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition EvtPDL.cc:272
EvtParticle::epsParent
virtual EvtVector4C epsParent(int i) const
Definition EvtParticle.cc:582
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:100
EvtParticle::getP4Restframe
EvtVector4R getP4Restframe()
Definition EvtParticle.cc:710
EvtParticle::getNDaug
int getNDaug() const
Definition EvtParticle.cc:112
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::mass
double mass() const
Definition EvtParticle.cc:114
EvtParticle::eps
virtual EvtVector4C eps(int i) const
Definition EvtParticle.cc:592
EvtRandom::Flat
static double Flat(double min, double max)
Definition EvtRandom.cc:55
EvtVPHOtoVISRHi::decay
void decay(EvtParticle *p)
Definition EvtVPHOtoVISRHi.cc:57
EvtVPHOtoVISRHi::~EvtVPHOtoVISRHi
virtual ~EvtVPHOtoVISRHi()
Definition EvtVPHOtoVISRHi.cc:35
EvtVPHOtoVISRHi::clone
EvtDecayBase * clone()
Definition EvtVPHOtoVISRHi.cc:39
EvtVPHOtoVISRHi::getName
void getName(std::string &name)
Definition EvtVPHOtoVISRHi.cc:37
EvtVector4C::conj
EvtVector4C conj() const
Definition EvtVector4C.hh:197
EvtVector4R::mass
double mass() const
Definition EvtVector4R.cc:40