]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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) 1998 Caltech, UCSB | |
10 | // | |
11 | // Module: EvtParticle.cc | |
12 | // | |
13 | // Description: Class to describe all particles | |
14 | // | |
15 | // Modification history: | |
16 | // | |
17 | // DJL/RYD September 25, 1996 Module created | |
18 | // | |
19 | //------------------------------------------------------------------------ | |
20 | // | |
21 | #include "EvtGenBase/EvtPatches.hh" | |
22 | #include <iostream> | |
23 | #include <math.h> | |
24 | #include <stdio.h> | |
25 | #include <stdlib.h> | |
26 | #include <sys/stat.h> | |
27 | #include "EvtGenBase/EvtParticle.hh" | |
28 | #include "EvtGenBase/EvtId.hh" | |
29 | #include "EvtGenBase/EvtRandom.hh" | |
30 | #include "EvtGenBase/EvtRadCorr.hh" | |
31 | #include "EvtGenBase/EvtPDL.hh" | |
32 | #include "EvtGenBase/EvtDecayTable.hh" | |
33 | #include "EvtGenBase/EvtDiracParticle.hh" | |
34 | #include "EvtGenBase/EvtScalarParticle.hh" | |
35 | #include "EvtGenBase/EvtVectorParticle.hh" | |
36 | #include "EvtGenBase/EvtTensorParticle.hh" | |
37 | #include "EvtGenBase/EvtPhotonParticle.hh" | |
38 | #include "EvtGenBase/EvtNeutrinoParticle.hh" | |
39 | #include "EvtGenBase/EvtRaritaSchwingerParticle.hh" | |
40 | #include "EvtGenBase/EvtStringParticle.hh" | |
41 | #include "EvtGenBase/EvtStdHep.hh" | |
42 | #include "EvtGenBase/EvtSecondary.hh" | |
43 | #include "EvtGenBase/EvtReport.hh" | |
44 | #include "EvtGenBase/EvtGenKine.hh" | |
45 | #include "EvtGenBase/EvtCPUtil.hh" | |
46 | #include "EvtGenBase/EvtParticleFactory.hh" | |
47 | #include "EvtGenBase/EvtIdSet.hh" | |
0ca57c2f | 48 | #include "EvtGenBase/EvtStatus.hh" |
da0e9ce3 | 49 | |
50 | using std::endl; | |
51 | ||
52 | ||
53 | ||
54 | EvtParticle::~EvtParticle() { | |
55 | delete _decayProb; | |
56 | } | |
57 | ||
58 | EvtParticle::EvtParticle() { | |
59 | _ndaug=0; | |
60 | _parent=0; | |
61 | _channel=-10; | |
62 | _t=0.0; | |
63 | _genlifetime=1; | |
64 | _first=1; | |
65 | _isInit=false; | |
66 | _validP4=false; | |
67 | _isDecayed=false; | |
68 | _decayProb=0; | |
69 | // _mix=false; | |
70 | } | |
71 | ||
72 | void EvtParticle::setFirstOrNot() { | |
73 | _first=0; | |
74 | } | |
75 | void EvtParticle::resetFirstOrNot() { | |
76 | _first=1; | |
77 | } | |
78 | ||
79 | void EvtParticle::setChannel( int i ) { | |
80 | _channel=i; | |
81 | } | |
82 | ||
83 | EvtParticle *EvtParticle::getDaug(int i) { return _daug[i]; } | |
84 | ||
85 | EvtParticle *EvtParticle::getParent() const { return _parent;} | |
86 | ||
87 | void EvtParticle::setLifetime(double tau){ | |
88 | _t=tau; | |
89 | } | |
90 | ||
91 | void EvtParticle::setLifetime(){ | |
92 | if (_genlifetime){ | |
93 | _t=-log(EvtRandom::Flat())*EvtPDL::getctau(getId()); | |
94 | } | |
95 | } | |
96 | ||
97 | double EvtParticle::getLifetime(){ | |
98 | ||
99 | return _t; | |
100 | } | |
101 | ||
102 | void EvtParticle::addDaug(EvtParticle *node) { | |
103 | node->_daug[node->_ndaug++]=this; | |
104 | _ndaug=0; | |
105 | _parent=node; | |
106 | } | |
107 | ||
108 | ||
109 | int EvtParticle::firstornot() const { return _first;} | |
110 | ||
111 | EvtId EvtParticle::getId() const { return _id;} | |
112 | ||
0ca57c2f | 113 | int EvtParticle::getPDGId() const {return EvtPDL::getStdHep(_id);} |
114 | ||
da0e9ce3 | 115 | EvtSpinType::spintype EvtParticle::getSpinType() const |
116 | { return EvtPDL::getSpinType(_id);} | |
117 | ||
118 | int EvtParticle::getSpinStates() const | |
119 | { return EvtSpinType::getSpinStates(EvtPDL::getSpinType(_id));} | |
120 | ||
121 | const EvtVector4R& EvtParticle::getP4() const { return _p;} | |
122 | ||
123 | int EvtParticle::getChannel() const { return _channel;} | |
124 | ||
125 | size_t EvtParticle::getNDaug() const { return _ndaug;} | |
126 | ||
127 | double EvtParticle::mass() const { | |
128 | ||
129 | return _p.mass(); | |
130 | } | |
131 | ||
132 | ||
133 | void EvtParticle::setDiagonalSpinDensity(){ | |
134 | ||
135 | _rhoForward.setDiag(getSpinStates()); | |
136 | } | |
137 | ||
138 | void EvtParticle::setVectorSpinDensity(){ | |
139 | ||
140 | if (getSpinStates()!=3) { | |
141 | report(ERROR,"EvtGen")<<"Error in EvtParticle::setVectorSpinDensity"<<endl; | |
142 | report(ERROR,"EvtGen")<<"spin_states:"<<getSpinStates()<<endl; | |
143 | report(ERROR,"EvtGen")<<"particle:"<<EvtPDL::name(_id).c_str()<<endl; | |
144 | ::abort(); | |
145 | } | |
146 | ||
147 | EvtSpinDensity rho; | |
148 | ||
149 | //Set helicity +1 and -1 to 1. | |
150 | rho.setDiag(getSpinStates()); | |
151 | rho.set(1,1,EvtComplex(0.0,0.0)); | |
152 | ||
153 | setSpinDensityForwardHelicityBasis(rho); | |
154 | ||
155 | } | |
156 | ||
157 | ||
158 | void EvtParticle::setSpinDensityForwardHelicityBasis(const EvtSpinDensity& rho){ | |
159 | ||
160 | EvtSpinDensity R=rotateToHelicityBasis(); | |
161 | ||
162 | assert(R.getDim()==rho.getDim()); | |
163 | ||
164 | int n=rho.getDim(); | |
165 | ||
166 | _rhoForward.setDim(n); | |
167 | ||
168 | int i,j,k,l; | |
169 | ||
170 | for(i=0;i<n;i++){ | |
171 | for(j=0;j<n;j++){ | |
172 | EvtComplex tmp=0.0; | |
173 | for(k=0;k<n;k++){ | |
174 | for(l=0;l<n;l++){ | |
175 | tmp+=R.get(l,i)*rho.get(l,k)*conj(R.get(k,j)); | |
176 | } | |
177 | } | |
178 | _rhoForward.set(i,j,tmp); | |
179 | } | |
180 | } | |
181 | ||
182 | } | |
183 | ||
184 | void EvtParticle::setSpinDensityForwardHelicityBasis(const EvtSpinDensity& rho, | |
185 | double alpha, | |
186 | double beta, | |
187 | double gamma){ | |
188 | ||
189 | EvtSpinDensity R=rotateToHelicityBasis(alpha,beta,gamma); | |
190 | ||
191 | assert(R.getDim()==rho.getDim()); | |
192 | ||
193 | int n=rho.getDim(); | |
194 | ||
195 | _rhoForward.setDim(n); | |
196 | ||
197 | int i,j,k,l; | |
198 | ||
199 | for(i=0;i<n;i++){ | |
200 | for(j=0;j<n;j++){ | |
201 | EvtComplex tmp=0.0; | |
202 | for(k=0;k<n;k++){ | |
203 | for(l=0;l<n;l++){ | |
204 | tmp+=R.get(l,i)*rho.get(l,k)*conj(R.get(k,j)); | |
205 | } | |
206 | } | |
207 | _rhoForward.set(i,j,tmp); | |
208 | } | |
209 | } | |
210 | ||
211 | } | |
212 | ||
213 | void EvtParticle::initDecay(bool useMinMass) { | |
214 | ||
215 | EvtParticle* p=this; | |
216 | // carefull - the parent mass might be fixed in stone.. | |
217 | EvtParticle* par=p->getParent(); | |
218 | double parMass=-1.; | |
219 | if ( par != 0 ) { | |
220 | if ( par->hasValidP4() ) parMass=par->mass(); | |
221 | for (size_t i=0;i<par->getNDaug();i++) { | |
222 | EvtParticle *tDaug=par->getDaug(i); | |
223 | if ( p != tDaug ) | |
0ca57c2f | 224 | parMass-=EvtPDL::getMinMass(tDaug->getId()); |
da0e9ce3 | 225 | } |
226 | } | |
0ca57c2f | 227 | |
da0e9ce3 | 228 | if ( _isInit ) { |
229 | //we have already been here - just reroll the masses! | |
230 | if ( _ndaug>0) { | |
231 | for(size_t ii=0;ii<_ndaug;ii++){ | |
0ca57c2f | 232 | if ( _ndaug==1 || EvtPDL::getWidth(p->getDaug(ii)->getId()) > 0.0000001) |
233 | p->getDaug(ii)->initDecay(useMinMass); | |
234 | else p->getDaug(ii)->setMass(EvtPDL::getMeanMass(p->getDaug(ii)->getId())); | |
da0e9ce3 | 235 | } |
236 | } | |
0ca57c2f | 237 | |
da0e9ce3 | 238 | EvtId *dauId=0; |
239 | double *dauMasses=0; | |
240 | if ( _ndaug > 0) { | |
241 | dauId=new EvtId[_ndaug]; | |
242 | dauMasses=new double[_ndaug]; | |
243 | for (size_t j=0;j<_ndaug;j++) { | |
0ca57c2f | 244 | dauId[j]=p->getDaug(j)->getId(); |
245 | dauMasses[j]=p->getDaug(j)->mass(); | |
da0e9ce3 | 246 | } |
247 | } | |
248 | EvtId *parId=0; | |
249 | EvtId *othDauId=0; | |
250 | EvtParticle *tempPar=p->getParent(); | |
251 | if (tempPar) { | |
252 | parId=new EvtId(tempPar->getId()); | |
253 | if ( tempPar->getNDaug()==2 ) { | |
0ca57c2f | 254 | if ( tempPar->getDaug(0) == this ) othDauId=new EvtId(tempPar->getDaug(1)->getId()); |
255 | else othDauId=new EvtId(tempPar->getDaug(0)->getId()); | |
da0e9ce3 | 256 | } |
257 | } | |
258 | if ( p->getParent() && _validP4==false ) { | |
259 | if ( !useMinMass ) { | |
0ca57c2f | 260 | p->setMass(EvtPDL::getRandMass(p->getId(),parId,_ndaug,dauId,othDauId,parMass,dauMasses)); |
da0e9ce3 | 261 | } |
262 | else p->setMass(EvtPDL::getMinMass(p->getId())); | |
263 | } | |
264 | if ( parId) delete parId; | |
265 | if ( othDauId) delete othDauId; | |
266 | if ( dauId) delete [] dauId; | |
267 | if ( dauMasses) delete [] dauMasses; | |
268 | return; | |
269 | } | |
0ca57c2f | 270 | |
da0e9ce3 | 271 | |
272 | //Will include effects of mixing here | |
273 | //added by Lange Jan4,2000 | |
274 | static EvtId BS0=EvtPDL::getId("B_s0"); | |
275 | static EvtId BSB=EvtPDL::getId("anti-B_s0"); | |
276 | static EvtId BD0=EvtPDL::getId("B0"); | |
277 | static EvtId BDB=EvtPDL::getId("anti-B0"); | |
278 | static EvtId D0=EvtPDL::getId("D0"); | |
279 | static EvtId D0B=EvtPDL::getId("anti-D0"); | |
280 | static EvtId U4S=EvtPDL::getId("Upsilon(4S)"); | |
281 | static EvtIdSet borUps(BS0,BSB,BD0,BDB,U4S); | |
0ca57c2f | 282 | |
da0e9ce3 | 283 | //only makes sense if there is no parent particle which is a B or an Upsilon |
284 | bool hasBorUps=false; | |
285 | if ( getParent() && borUps.contains(getParent()->getId()) ) hasBorUps=true; | |
286 | // if ( (getNDaug()==0)&&(getParent()==0) && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB)){ | |
287 | EvtId thisId=getId(); | |
288 | // remove D0 mixing for now. | |
289 | // if ( (getNDaug()==0 && !hasBorUps) && (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB||thisId==D0||thisId==D0B)){ | |
290 | if ( (getNDaug()==0 && !hasBorUps) && (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB)){ | |
291 | double t; | |
292 | int mix; | |
0ca57c2f | 293 | EvtCPUtil::getInstance()->incoherentMix(getId(), t, mix); |
da0e9ce3 | 294 | setLifetime(t); |
295 | ||
296 | if (mix) { | |
297 | ||
298 | EvtScalarParticle* scalar_part; | |
299 | ||
300 | scalar_part=new EvtScalarParticle; | |
301 | if (getId()==BS0) { | |
0ca57c2f | 302 | EvtVector4R p_init(EvtPDL::getMass(BSB),0.0,0.0,0.0); |
303 | scalar_part->init(EvtPDL::chargeConj(getId()),p_init); | |
da0e9ce3 | 304 | } |
305 | else if (getId()==BSB) { | |
0ca57c2f | 306 | EvtVector4R p_init(EvtPDL::getMass(BS0),0.0,0.0,0.0); |
307 | scalar_part->init(EvtPDL::chargeConj(getId()),p_init); | |
da0e9ce3 | 308 | } |
309 | else if (getId()==BD0) { | |
0ca57c2f | 310 | EvtVector4R p_init(EvtPDL::getMass(BDB),0.0,0.0,0.0); |
311 | scalar_part->init(EvtPDL::chargeConj(getId()),p_init); | |
da0e9ce3 | 312 | } |
313 | else if (getId()==BDB) { | |
0ca57c2f | 314 | EvtVector4R p_init(EvtPDL::getMass(BD0),0.0,0.0,0.0); |
315 | scalar_part->init(EvtPDL::chargeConj(getId()),p_init); | |
da0e9ce3 | 316 | } |
317 | else if (getId()==D0) { | |
0ca57c2f | 318 | EvtVector4R p_init(EvtPDL::getMass(D0B),0.0,0.0,0.0); |
319 | scalar_part->init(EvtPDL::chargeConj(getId()),p_init); | |
da0e9ce3 | 320 | } |
321 | else if (getId()==D0B) { | |
0ca57c2f | 322 | EvtVector4R p_init(EvtPDL::getMass(D0),0.0,0.0,0.0); |
323 | scalar_part->init(EvtPDL::chargeConj(getId()),p_init); | |
da0e9ce3 | 324 | } |
0ca57c2f | 325 | |
da0e9ce3 | 326 | scalar_part->setLifetime(0); |
327 | scalar_part->setDiagonalSpinDensity(); | |
0ca57c2f | 328 | |
da0e9ce3 | 329 | insertDaugPtr(0,scalar_part); |
330 | ||
331 | _ndaug=1; | |
332 | _isInit=true; | |
333 | p=scalar_part; | |
334 | p->initDecay(useMinMass); | |
335 | return; | |
336 | ||
337 | ||
338 | } | |
339 | } | |
da0e9ce3 | 340 | |
341 | EvtDecayBase *decayer; | |
0ca57c2f | 342 | decayer = EvtDecayTable::getInstance()->getDecayFunc(p); |
da0e9ce3 | 343 | |
344 | if ( decayer ) { | |
345 | p->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs()); | |
346 | //then loop over the daughters and init their decay | |
347 | for(size_t i=0;i<p->getNDaug();i++){ | |
348 | // std::cout << EvtPDL::name(p->getDaug(i)->getId()) << " " << i << " " << p->getDaug(i)->getSpinType() << " " << EvtPDL::name(p->getId()) << std::endl; | |
349 | if ( EvtPDL::getWidth(p->getDaug(i)->getId()) > 0.0000001) | |
0ca57c2f | 350 | p->getDaug(i)->initDecay(useMinMass); |
da0e9ce3 | 351 | else p->getDaug(i)->setMass(EvtPDL::getMeanMass(p->getDaug(i)->getId())); |
352 | } | |
353 | } | |
354 | ||
355 | int j; | |
356 | EvtId *dauId=0; | |
357 | double *dauMasses=0; | |
358 | int nDaugT=p->getNDaug(); | |
359 | if ( nDaugT > 0) { | |
360 | dauId=new EvtId[nDaugT]; | |
361 | dauMasses=new double[nDaugT]; | |
362 | for (j=0;j<nDaugT;j++) { | |
363 | dauId[j]=p->getDaug(j)->getId(); | |
364 | dauMasses[j]=p->getDaug(j)->mass(); | |
365 | } | |
366 | } | |
367 | ||
368 | EvtId *parId=0; | |
369 | EvtId *othDauId=0; | |
370 | EvtParticle *tempPar=p->getParent(); | |
371 | if (tempPar) { | |
372 | parId=new EvtId(tempPar->getId()); | |
373 | if ( tempPar->getNDaug()==2 ) { | |
374 | if ( tempPar->getDaug(0) == this ) othDauId=new EvtId(tempPar->getDaug(1)->getId()); | |
375 | else othDauId=new EvtId(tempPar->getDaug(0)->getId()); | |
376 | } | |
377 | } | |
378 | if ( p->getParent() && p->hasValidP4()==false ) { | |
379 | if ( !useMinMass ) { | |
380 | p->setMass(EvtPDL::getRandMass(p->getId(),parId,p->getNDaug(),dauId,othDauId,parMass,dauMasses)); | |
381 | } | |
382 | else { | |
383 | p->setMass(EvtPDL::getMinMass(p->getId())); | |
384 | } | |
385 | } | |
386 | if ( parId) delete parId; | |
387 | if ( othDauId) delete othDauId; | |
388 | if ( dauId) delete [] dauId; | |
389 | if ( dauMasses) delete [] dauMasses; | |
390 | _isInit=true; | |
391 | } | |
392 | ||
393 | ||
394 | void EvtParticle::decay(){ | |
395 | //P is particle to decay, typically 'this' but sometime | |
396 | //modified by mixing | |
397 | EvtParticle* p=this; | |
398 | //Did it mix? | |
399 | //if ( p->getMixed() ) { | |
400 | //should take C(p) - this should only | |
401 | //happen the first time we call decay for this | |
402 | //particle | |
403 | //p->takeCConj(); | |
404 | // p->setUnMixed(); | |
405 | //} | |
406 | ||
407 | EvtDecayBase *decayer; | |
0ca57c2f | 408 | decayer = EvtDecayTable::getInstance()->getDecayFunc(p); |
da0e9ce3 | 409 | // if ( decayer ) { |
410 | // report(INFO,"EvtGen") << "calling decay for " << EvtPDL::name(p->getId()) << " " << p->mass() << " " << p->getP4() << " " << p->getNDaug() << " " << p << endl; | |
411 | // report(INFO,"EvtGen") << "NDaug= " << decayer->getNDaug() << endl; | |
412 | // int ti; | |
413 | // for ( ti=0; ti<decayer->getNDaug(); ti++) | |
414 | // report(INFO,"EvtGen") << "Daug " << ti << " " << EvtPDL::name(decayer->getDaug(ti)) << endl; | |
415 | // } | |
416 | //if (p->_ndaug>0) { | |
417 | // report(INFO,"EvtGen") <<"Is decaying particle with daughters!!!!!"<<endl; | |
418 | // ::abort(); | |
419 | //return; | |
420 | //call initdecay first - April 29,2002 - Lange | |
421 | //} | |
422 | ||
423 | //if there are already daughters, then this step is already done! | |
424 | // figure out the masses | |
0ca57c2f | 425 | bool massTreeOK(true); |
426 | if ( _ndaug == 0 ) { | |
427 | massTreeOK = generateMassTree(); | |
428 | } | |
429 | ||
430 | if (massTreeOK == false) { | |
431 | report(INFO,"EvtGen")<<"Could not decay "<<EvtPDL::name(p->getId()) | |
432 | <<" with mass "<<p->mass() | |
433 | <<" to decay channel number "<<_channel<<endl; | |
434 | _isDecayed = false; | |
435 | return; | |
436 | } | |
da0e9ce3 | 437 | |
438 | static EvtId BS0=EvtPDL::getId("B_s0"); | |
439 | static EvtId BSB=EvtPDL::getId("anti-B_s0"); | |
440 | static EvtId BD0=EvtPDL::getId("B0"); | |
441 | static EvtId BDB=EvtPDL::getId("anti-B0"); | |
0ca57c2f | 442 | // static EvtId D0=EvtPDL::getId("D0"); |
443 | // static EvtId D0B=EvtPDL::getId("anti-D0"); | |
da0e9ce3 | 444 | |
445 | EvtId thisId=getId(); | |
446 | // remove D0 mixing for now.. | |
447 | // if ( _ndaug==1 && (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB||thisId==D0||thisId==D0B) ) { | |
448 | if ( _ndaug==1 && (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB) ) { | |
449 | p=p->getDaug(0); | |
0ca57c2f | 450 | decayer = EvtDecayTable::getInstance()->getDecayFunc(p); |
da0e9ce3 | 451 | } |
452 | //now we have accepted a set of masses - time | |
453 | if ( decayer != 0) { | |
454 | decayer->makeDecay(p); | |
455 | } | |
456 | else{ | |
457 | p->_rhoBackward.setDiag(p->getSpinStates()); | |
458 | } | |
459 | ||
460 | _isDecayed=true; | |
461 | return; | |
462 | } | |
463 | ||
0ca57c2f | 464 | bool EvtParticle::generateMassTree() { |
465 | ||
466 | bool isOK(true); | |
467 | ||
da0e9ce3 | 468 | double massProb=1.; |
469 | double ranNum=2.; | |
470 | int counter=0; | |
471 | EvtParticle *p=this; | |
472 | while (massProb<ranNum) { | |
473 | //check it out the first time. | |
474 | p->initDecay(); | |
475 | massProb=p->compMassProb(); | |
476 | ranNum=EvtRandom::Flat(); | |
477 | counter++; | |
478 | ||
479 | if ( counter > 10000 ) { | |
480 | if ( counter == 10001 ) { | |
481 | report(INFO,"EvtGen") << "Too many iterations to determine the mass tree. Parent mass= "<< p->mass() << " " << massProb <<endl; | |
482 | p->printTree(); | |
483 | report(INFO,"EvtGen") << "will take next combo with non-zero likelihood\n"; | |
484 | } | |
485 | if ( massProb>0. ) massProb=2.0; | |
486 | if ( counter > 20000 ) { | |
487 | // one last try - take the minimum masses | |
488 | p->initDecay(true); | |
489 | p->printTree(); | |
490 | massProb=p->compMassProb(); | |
491 | if ( massProb>0. ) { | |
492 | massProb=2.0; | |
493 | report(INFO,"EvtGen") << "Taking the minimum mass of all particles in the chain\n"; | |
494 | } | |
495 | else { | |
496 | report(INFO,"EvtGen") << "Sorry, no luck finding a valid set of masses. This may be a pathological combo\n"; | |
0ca57c2f | 497 | isOK = false; |
498 | break; | |
da0e9ce3 | 499 | } |
500 | } | |
501 | } | |
502 | } | |
0ca57c2f | 503 | |
504 | return isOK; | |
505 | ||
da0e9ce3 | 506 | } |
507 | ||
508 | double EvtParticle::compMassProb() { | |
509 | ||
510 | EvtParticle *p=this; | |
511 | double mass=p->mass(); | |
512 | double parMass=0.; | |
513 | if ( p->getParent()) { | |
514 | parMass=p->getParent()->mass(); | |
515 | } | |
516 | ||
517 | int nDaug=p->getNDaug(); | |
518 | double *dMasses=0; | |
519 | ||
520 | int i; | |
521 | if ( nDaug>0 ) { | |
522 | dMasses=new double[nDaug]; | |
523 | for (i=0; i<nDaug; i++) dMasses[i]=p->getDaug(i)->mass(); | |
524 | } | |
525 | ||
526 | double temp=1.0; | |
527 | temp=EvtPDL::getMassProb(p->getId(), mass, parMass, nDaug, dMasses); | |
528 | ||
529 | //If the particle already has a mass, we dont need to include | |
530 | //it in the probability calculation | |
531 | if ( (!p->getParent() || _validP4 ) && temp>0.0 ) temp=1.; | |
532 | ||
533 | delete [] dMasses; | |
534 | for (i=0; i<nDaug; i++) { | |
535 | temp*=p->getDaug(i)->compMassProb(); | |
536 | } | |
537 | return temp; | |
538 | } | |
539 | ||
540 | void EvtParticle::deleteDaughters(bool keepChannel){ | |
541 | ||
542 | for(size_t i=0;i<_ndaug;i++){ | |
543 | _daug[i]->deleteTree(); | |
544 | } | |
545 | ||
546 | _ndaug=0; | |
547 | if ( !keepChannel) _channel=-10; | |
548 | _first=1; | |
549 | _isInit=false; | |
550 | } | |
551 | ||
552 | void EvtParticle::deleteTree(){ | |
553 | ||
554 | this->deleteDaughters(); | |
555 | ||
556 | delete this; | |
557 | ||
558 | } | |
559 | ||
560 | EvtVector4C EvtParticle::epsParent(int i) const { | |
561 | EvtVector4C temp; | |
562 | printParticle(); | |
563 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
564 | <<"th polarization vector." | |
565 | <<" I.e. you thought it was a" | |
566 | <<" vector particle!" << endl; | |
567 | ::abort(); | |
568 | return temp; | |
569 | } | |
570 | ||
571 | EvtVector4C EvtParticle::eps(int i) const { | |
572 | EvtVector4C temp; | |
573 | printParticle(); | |
574 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
575 | <<"th polarization vector." | |
576 | <<" I.e. you thought it was a" | |
577 | <<" vector particle!" << endl; | |
578 | ::abort(); | |
579 | return temp; | |
580 | } | |
581 | ||
582 | EvtVector4C EvtParticle::epsParentPhoton(int i){ | |
583 | EvtVector4C temp; | |
584 | printParticle(); | |
585 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
586 | <<"th polarization vector of photon." | |
587 | <<" I.e. you thought it was a" | |
588 | <<" photon particle!" << endl; | |
589 | ::abort(); | |
590 | return temp; | |
591 | } | |
592 | ||
593 | EvtVector4C EvtParticle::epsPhoton(int i){ | |
594 | EvtVector4C temp; | |
595 | printParticle(); | |
596 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
597 | <<"th polarization vector of a photon." | |
598 | <<" I.e. you thought it was a" | |
599 | <<" photon particle!" << endl; | |
600 | ::abort(); | |
601 | return temp; | |
602 | } | |
603 | ||
604 | EvtDiracSpinor EvtParticle::spParent(int i) const { | |
605 | EvtDiracSpinor tempD; | |
da0e9ce3 | 606 | printParticle(); |
607 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
608 | <<"th dirac spinor." | |
609 | <<" I.e. you thought it was a" | |
610 | <<" Dirac particle!" << endl; | |
611 | ::abort(); | |
612 | return tempD; | |
613 | } | |
614 | ||
615 | EvtDiracSpinor EvtParticle::sp(int i) const { | |
616 | EvtDiracSpinor tempD; | |
da0e9ce3 | 617 | printParticle(); |
618 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
619 | <<"th dirac spinor." | |
620 | <<" I.e. you thought it was a" | |
621 | <<" Dirac particle!" << endl; | |
622 | ::abort(); | |
623 | return tempD; | |
624 | } | |
625 | ||
626 | EvtDiracSpinor EvtParticle::spParentNeutrino() const { | |
627 | EvtDiracSpinor tempD; | |
628 | printParticle(); | |
629 | report(ERROR,"EvtGen") << "and you have asked for the " | |
630 | <<"dirac spinor." | |
631 | <<" I.e. you thought it was a" | |
632 | <<" neutrino particle!" << endl; | |
633 | ::abort(); | |
634 | return tempD; | |
635 | } | |
636 | ||
637 | EvtDiracSpinor EvtParticle::spNeutrino() const { | |
638 | EvtDiracSpinor tempD; | |
639 | printParticle(); | |
640 | report(ERROR,"EvtGen") << "and you have asked for the " | |
641 | <<"dirac spinor." | |
642 | <<" I.e. you thought it was a" | |
643 | <<" neutrino particle!" << endl; | |
644 | ::abort(); | |
645 | return tempD; | |
646 | } | |
647 | ||
648 | EvtTensor4C EvtParticle::epsTensorParent(int i) const { | |
da0e9ce3 | 649 | EvtTensor4C tempC; |
650 | printParticle(); | |
651 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
652 | <<"th tensor." | |
653 | <<" I.e. you thought it was a" | |
654 | <<" Tensor particle!" << endl; | |
655 | ::abort(); | |
656 | return tempC; | |
657 | } | |
658 | ||
659 | EvtTensor4C EvtParticle::epsTensor(int i) const { | |
da0e9ce3 | 660 | EvtTensor4C tempC; |
661 | printParticle(); | |
662 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
663 | <<"th tensor." | |
664 | <<" I.e. you thought it was a" | |
665 | <<" Tensor particle!" << endl; | |
666 | ::abort(); | |
667 | return tempC; | |
668 | } | |
669 | ||
670 | ||
671 | EvtRaritaSchwinger EvtParticle::spRSParent(int i) const { | |
672 | EvtRaritaSchwinger tempD; | |
da0e9ce3 | 673 | printParticle(); |
674 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
675 | <<"th Rarita-Schwinger spinor." | |
676 | <<" I.e. you thought it was a" | |
677 | <<" RaritaSchwinger particle!" << std::endl; | |
678 | ::abort(); | |
679 | return tempD; | |
680 | } | |
681 | ||
682 | EvtRaritaSchwinger EvtParticle::spRS(int i) const { | |
683 | EvtRaritaSchwinger tempD; | |
da0e9ce3 | 684 | printParticle(); |
685 | report(ERROR,"EvtGen") << "and you have asked for the:"<<i | |
686 | <<"th Rarita-Schwinger spinor." | |
687 | <<" I.e. you thought it was a" | |
688 | <<" RaritaSchwinger particle!" << std::endl; | |
689 | ::abort(); | |
690 | return tempD; | |
691 | } | |
692 | ||
693 | ||
694 | ||
695 | EvtVector4R EvtParticle::getP4Lab() const { | |
696 | EvtVector4R temp,mom; | |
697 | const EvtParticle *ptemp; | |
698 | ||
699 | temp=this->getP4(); | |
700 | ptemp=this; | |
701 | ||
702 | while (ptemp->getParent()!=0) { | |
703 | ptemp=ptemp->getParent(); | |
704 | mom=ptemp->getP4(); | |
705 | temp=boostTo(temp,mom); | |
706 | } | |
707 | return temp; | |
708 | } | |
709 | ||
710 | EvtVector4R EvtParticle::getP4LabBeforeFSR() { | |
711 | EvtVector4R temp,mom; | |
712 | EvtParticle *ptemp; | |
713 | ||
714 | temp=this->_pBeforeFSR; | |
715 | ptemp=this; | |
716 | ||
717 | while (ptemp->getParent()!=0) { | |
718 | ptemp=ptemp->getParent(); | |
719 | mom=ptemp->getP4(); | |
720 | temp=boostTo(temp,mom); | |
721 | } | |
722 | return temp; | |
723 | } | |
724 | ||
725 | ||
726 | ||
727 | EvtVector4R EvtParticle::getP4Restframe() const { | |
728 | ||
729 | return EvtVector4R(mass(),0.0,0.0,0.0); | |
730 | ||
731 | } | |
732 | ||
733 | EvtVector4R EvtParticle::get4Pos() const { | |
734 | ||
735 | EvtVector4R temp,mom; | |
736 | EvtParticle *ptemp; | |
737 | ||
738 | temp.set(0.0,0.0,0.0,0.0); | |
739 | ptemp=getParent(); | |
740 | ||
741 | if (ptemp==0) return temp; | |
742 | ||
743 | temp=(ptemp->_t/ptemp->mass())*(ptemp->getP4()); | |
744 | ||
745 | while (ptemp->getParent()!=0) { | |
746 | ptemp=ptemp->getParent(); | |
747 | mom=ptemp->getP4(); | |
748 | temp=boostTo(temp,mom); | |
749 | temp=temp+(ptemp->_t/ptemp->mass())*(ptemp->getP4()); | |
750 | } | |
751 | ||
752 | return temp; | |
753 | } | |
754 | ||
755 | ||
756 | EvtParticle * EvtParticle::nextIter(EvtParticle *rootOfTree) { | |
757 | ||
758 | EvtParticle *bpart; | |
759 | EvtParticle *current; | |
760 | ||
761 | current=this; | |
762 | size_t i; | |
763 | ||
764 | if (_ndaug!=0) return _daug[0]; | |
765 | ||
766 | do{ | |
767 | bpart=current->_parent; | |
768 | if (bpart==0) return 0; | |
769 | i=0; | |
770 | while (bpart->_daug[i]!=current) {i++;} | |
771 | ||
772 | if ( bpart==rootOfTree ) { | |
773 | if ( i+1 == bpart->_ndaug ) return 0; | |
774 | } | |
775 | ||
776 | i++; | |
777 | current=bpart; | |
778 | ||
779 | }while(i>=bpart->_ndaug); | |
780 | ||
781 | return bpart->_daug[i]; | |
782 | ||
783 | } | |
784 | ||
785 | ||
786 | void EvtParticle::makeStdHep(EvtStdHep& stdhep,EvtSecondary& secondary, | |
787 | EvtId *list_of_stable){ | |
788 | ||
789 | //first add particle to the stdhep list; | |
790 | stdhep.createParticle(getP4Lab(),get4Pos(),-1,-1, | |
791 | EvtPDL::getStdHep(getId())); | |
792 | ||
793 | int ii=0; | |
794 | ||
795 | //lets see if this is a longlived particle and terminate the | |
796 | //list building! | |
797 | ||
798 | while (list_of_stable[ii]!=EvtId(-1,-1)) { | |
799 | if (getId()==list_of_stable[ii]){ | |
800 | secondary.createSecondary(0,this); | |
801 | return; | |
802 | } | |
803 | ii++; | |
804 | } | |
805 | ||
806 | ||
807 | ||
808 | for(size_t i=0;i<_ndaug;i++){ | |
809 | stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),0,0, | |
810 | EvtPDL::getStdHep(_daug[i]->getId())); | |
811 | } | |
812 | ||
813 | for(size_t i=0;i<_ndaug;i++){ | |
814 | _daug[i]->makeStdHepRec(1+i,1+i,stdhep,secondary,list_of_stable); | |
815 | } | |
816 | return; | |
817 | ||
818 | } | |
819 | ||
820 | void EvtParticle::makeStdHep(EvtStdHep& stdhep){ | |
821 | ||
822 | //first add particle to the stdhep list; | |
823 | stdhep.createParticle(getP4Lab(),get4Pos(),-1,-1, | |
824 | EvtPDL::getStdHep(getId())); | |
825 | ||
826 | for(size_t i=0;i<_ndaug;i++){ | |
827 | stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(),0,0, | |
828 | EvtPDL::getStdHep(_daug[i]->getId())); | |
829 | } | |
830 | ||
831 | for(size_t i=0;i<_ndaug;i++){ | |
832 | _daug[i]->makeStdHepRec(1+i,1+i,stdhep); | |
833 | } | |
834 | return; | |
835 | ||
836 | } | |
837 | ||
838 | ||
839 | void EvtParticle::makeStdHepRec(int firstparent,int lastparent, | |
840 | EvtStdHep& stdhep, | |
841 | EvtSecondary& secondary, | |
842 | EvtId *list_of_stable){ | |
843 | ||
844 | ||
845 | int ii=0; | |
846 | ||
847 | //lets see if this is a longlived particle and terminate the | |
848 | //list building! | |
849 | ||
850 | while (list_of_stable[ii]!=EvtId(-1,-1)) { | |
851 | if (getId()==list_of_stable[ii]){ | |
852 | secondary.createSecondary(firstparent,this); | |
853 | return; | |
854 | } | |
855 | ii++; | |
856 | } | |
857 | ||
858 | ||
859 | ||
860 | int parent_num=stdhep.getNPart(); | |
861 | for(size_t i=0;i<_ndaug;i++){ | |
862 | stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(), | |
863 | firstparent,lastparent, | |
864 | EvtPDL::getStdHep(_daug[i]->getId())); | |
865 | } | |
866 | ||
867 | for(size_t i=0;i<_ndaug;i++){ | |
868 | _daug[i]->makeStdHepRec(parent_num+i,parent_num+i,stdhep, | |
869 | secondary,list_of_stable); | |
870 | } | |
871 | return; | |
872 | ||
873 | } | |
874 | ||
875 | void EvtParticle::makeStdHepRec(int firstparent,int lastparent, | |
876 | EvtStdHep& stdhep){ | |
877 | ||
878 | int parent_num=stdhep.getNPart(); | |
879 | for(size_t i=0;i<_ndaug;i++){ | |
880 | stdhep.createParticle(_daug[i]->getP4Lab(),_daug[i]->get4Pos(), | |
881 | firstparent,lastparent, | |
882 | EvtPDL::getStdHep(_daug[i]->getId())); | |
883 | } | |
884 | ||
885 | for(size_t i=0;i<_ndaug;i++){ | |
886 | _daug[i]->makeStdHepRec(parent_num+i,parent_num+i,stdhep); | |
887 | } | |
888 | return; | |
889 | ||
890 | } | |
891 | ||
892 | void EvtParticle::printTreeRec(unsigned int level) const { | |
893 | ||
894 | size_t newlevel,i; | |
895 | newlevel = level +1; | |
896 | ||
897 | ||
898 | if (_ndaug!=0) { | |
899 | if ( level > 0 ) { | |
900 | for (i=0;i<(5*level);i++) { | |
901 | report(INFO,"") <<" "; | |
902 | } | |
903 | } | |
904 | report(INFO,"") << EvtPDL::name(_id).c_str(); | |
905 | report(INFO,"") << " -> "; | |
906 | for(i=0;i<_ndaug;i++){ | |
907 | report(INFO,"") << EvtPDL::name(_daug[i]->getId()).c_str()<<" "; | |
908 | } | |
909 | for(i=0;i<_ndaug;i++){ | |
0ca57c2f | 910 | report(INFO,"") << _daug[i]->mass()<< " " << _daug[i]->getP4() << " " <<_daug[i]->getSpinStates() << "; "; |
da0e9ce3 | 911 | } |
912 | report(INFO,"")<<endl; | |
913 | for(i=0;i<_ndaug;i++){ | |
914 | _daug[i]->printTreeRec(newlevel); | |
915 | } | |
916 | } | |
917 | } | |
918 | ||
919 | void EvtParticle::printTree() const { | |
920 | ||
921 | report(INFO,"EvtGen") << "This is the current decay chain"<<endl; | |
922 | report(INFO,"") << "This top particle is "<< | |
923 | EvtPDL::name(_id).c_str()<<" " << this->mass() << " " << this->getP4() << endl; | |
924 | ||
925 | this->printTreeRec(0); | |
926 | report(INFO,"EvtGen") << "End of decay chain."<<endl; | |
927 | ||
928 | } | |
929 | ||
930 | std::string EvtParticle::treeStrRec(unsigned int level) const { | |
931 | ||
932 | size_t newlevel,i; | |
933 | newlevel = level +1; | |
934 | ||
935 | std::string retval=""; | |
936 | ||
937 | for(i=0;i<_ndaug;i++){ | |
938 | retval+=EvtPDL::name(_daug[i]->getId()); | |
939 | if ( _daug[i]->getNDaug() > 0 ) { | |
940 | retval+= " ("; | |
941 | retval+= _daug[i]->treeStrRec(newlevel); | |
942 | retval+= ") "; | |
943 | } | |
944 | else{ | |
945 | if ( i+1 !=_ndaug) retval+=" "; | |
946 | } | |
947 | } | |
948 | ||
949 | return retval; | |
950 | } | |
951 | ||
952 | ||
953 | std::string EvtParticle::treeStr() const { | |
954 | ||
955 | std::string retval=EvtPDL::name(_id); | |
956 | retval+=" -> "; | |
957 | ||
958 | retval+=treeStrRec(0); | |
959 | ||
960 | return retval; | |
961 | } | |
962 | ||
963 | void EvtParticle::printParticle() const { | |
964 | ||
965 | switch (EvtPDL::getSpinType(_id)){ | |
966 | case EvtSpinType::SCALAR: | |
967 | report(INFO,"EvtGen") << "This is a scalar particle:"<<EvtPDL::name(_id).c_str()<<"\n"; | |
968 | break; | |
969 | case EvtSpinType::VECTOR: | |
970 | report(INFO,"EvtGen") << "This is a vector particle:"<<EvtPDL::name(_id).c_str()<<"\n"; | |
971 | break; | |
972 | case EvtSpinType::TENSOR: | |
973 | report(INFO,"EvtGen") << "This is a tensor particle:"<<EvtPDL::name(_id).c_str()<<"\n"; | |
974 | break; | |
975 | case EvtSpinType::DIRAC: | |
976 | report(INFO,"EvtGen") << "This is a dirac particle:"<<EvtPDL::name(_id).c_str()<<"\n"; | |
977 | break; | |
978 | case EvtSpinType::PHOTON: | |
979 | report(INFO,"EvtGen") << "This is a photon:"<<EvtPDL::name(_id).c_str()<<"\n"; | |
980 | break; | |
981 | case EvtSpinType::NEUTRINO: | |
982 | report(INFO,"EvtGen") << "This is a neutrino:"<<EvtPDL::name(_id).c_str()<<"\n"; | |
983 | break; | |
984 | case EvtSpinType::STRING: | |
985 | report(INFO,"EvtGen") << "This is a string:"<<EvtPDL::name(_id).c_str()<<"\n"; | |
986 | break; | |
987 | default: | |
988 | report(INFO,"EvtGen") <<"Unknown particle type in EvtParticle::printParticle()"<<endl; | |
989 | break; | |
990 | } | |
991 | report(INFO,"EvtGen") << "Number of daughters:"<<_ndaug<<"\n"; | |
992 | ||
993 | ||
994 | } | |
995 | ||
996 | ||
997 | ||
998 | void init_vector( EvtParticle **part ){ | |
999 | *part = (EvtParticle *) new EvtVectorParticle; | |
1000 | } | |
1001 | ||
1002 | ||
1003 | void init_scalar( EvtParticle **part ){ | |
1004 | *part = (EvtParticle *) new EvtScalarParticle; | |
1005 | } | |
1006 | ||
1007 | void init_tensor( EvtParticle **part ){ | |
1008 | *part = (EvtParticle *) new EvtTensorParticle; | |
1009 | } | |
1010 | ||
1011 | void init_dirac( EvtParticle **part ){ | |
1012 | *part = (EvtParticle *) new EvtDiracParticle; | |
1013 | } | |
1014 | ||
1015 | void init_photon( EvtParticle **part ){ | |
1016 | *part = (EvtParticle *) new EvtPhotonParticle; | |
1017 | } | |
1018 | ||
1019 | void init_neutrino( EvtParticle **part ){ | |
1020 | *part = (EvtParticle *) new EvtNeutrinoParticle; | |
1021 | } | |
1022 | ||
1023 | void init_string( EvtParticle **part ){ | |
1024 | *part = (EvtParticle *) new EvtStringParticle; | |
1025 | } | |
1026 | ||
1027 | double EvtParticle::initializePhaseSpace( | |
1028 | unsigned int numdaughter,EvtId *daughters, | |
1029 | bool forceDaugMassReset, double poleSize, | |
1030 | int whichTwo1, int whichTwo2) { | |
1031 | ||
1032 | double m_b; | |
1033 | unsigned int i; | |
1034 | //lange | |
1035 | // this->makeDaughters(numdaughter,daughters); | |
1036 | ||
1037 | static EvtVector4R p4[100]; | |
1038 | static double mass[100]; | |
1039 | ||
1040 | m_b = this->mass(); | |
1041 | ||
1042 | //lange - Jan2,2002 - Need to check to see if the daughters of the parent | |
1043 | // have changed. If so, delete them and start over. | |
1044 | //report(INFO,"EvtGen") << "the parent is\n"; | |
1045 | //if ( this->getParent() ) { | |
1046 | // if ( this->getParent()->getParent() ) this->getParent()->getParent()->printTree(); | |
1047 | // this->getParent()->printTree(); | |
1048 | //} | |
1049 | //report(INFO,"EvtGen") << "and this is\n"; | |
1050 | //if ( this) this->printTree(); | |
1051 | bool resetDaughters=false; | |
1052 | ||
1053 | if ( numdaughter != this->getNDaug() && this->getNDaug() > 0 ) resetDaughters=true; | |
1054 | if ( numdaughter == this->getNDaug() ) | |
1055 | for (i=0; i<numdaughter;i++) { | |
1056 | if ( this->getDaug(i)->getId() != daughters[i] ) resetDaughters=true; | |
1057 | //report(INFO,"EvtGen") << EvtPDL::name(this->getDaug(i)->getId()) | |
1058 | // << " " << EvtPDL::name(daughters[i]) << endl; | |
1059 | } | |
1060 | ||
1061 | if ( resetDaughters || forceDaugMassReset) { | |
1062 | bool t1=true; | |
1063 | //but keep the decay channel of the parent. | |
1064 | this->deleteDaughters(t1); | |
1065 | this->makeDaughters(numdaughter,daughters); | |
0ca57c2f | 1066 | bool massTreeOK = this->generateMassTree(); |
1067 | if (massTreeOK == false) {return 0.0;} | |
da0e9ce3 | 1068 | } |
1069 | ||
1070 | double weight=0.; | |
1071 | for (i=0; i<numdaughter;i++) { | |
1072 | mass[i]=this->getDaug(i)->mass(); | |
1073 | } | |
1074 | ||
1075 | if ( poleSize<-0.1) { | |
1076 | //special case to enforce 4-momentum conservation in 1->1 decays | |
1077 | if (numdaughter==1) { | |
1078 | this->getDaug(0)->init(daughters[0],EvtVector4R(m_b,0.0,0.0,0.0)); | |
1079 | } | |
1080 | else{ | |
1081 | EvtGenKine::PhaseSpace( numdaughter, mass, p4, m_b ); | |
1082 | for(i=0;i<numdaughter;i++){ | |
1083 | this->getDaug(i)->init(daughters[i],p4[i]); | |
1084 | } | |
1085 | } | |
1086 | } | |
1087 | else { | |
1088 | if ( numdaughter != 3 ) { | |
1089 | report(ERROR,"EvtGen") << "Only can generate pole phase space " | |
1090 | << "distributions for 3 body final states" | |
1091 | << endl<<"Will terminate."<<endl; | |
1092 | ::abort(); | |
1093 | } | |
1094 | bool ok=false; | |
1095 | if ( (whichTwo1 == 1 && whichTwo2 == 0 ) || | |
1096 | (whichTwo1 == 0 && whichTwo2 == 1 ) ) { | |
1097 | weight=EvtGenKine::PhaseSpacePole( m_b, mass[0], mass[1], mass[2], | |
1098 | poleSize, p4); | |
1099 | this->getDaug(0)->init(daughters[0],p4[0]); | |
1100 | this->getDaug(1)->init(daughters[1],p4[1]); | |
1101 | this->getDaug(2)->init(daughters[2],p4[2]); | |
1102 | ok=true; | |
1103 | } | |
1104 | if ( (whichTwo1 == 1 && whichTwo2 == 2 ) || | |
1105 | (whichTwo1 == 2 && whichTwo2 == 1 ) ) { | |
1106 | weight=EvtGenKine::PhaseSpacePole( m_b, mass[2], mass[1], mass[0], | |
1107 | poleSize, p4); | |
1108 | this->getDaug(0)->init(daughters[0],p4[2]); | |
1109 | this->getDaug(1)->init(daughters[1],p4[1]); | |
1110 | this->getDaug(2)->init(daughters[2],p4[0]); | |
1111 | ok=true; | |
1112 | } | |
1113 | if ( (whichTwo1 == 0 && whichTwo2 == 2 ) || | |
1114 | (whichTwo1 == 2 && whichTwo2 == 0 ) ) { | |
1115 | weight=EvtGenKine::PhaseSpacePole( m_b, mass[1], mass[0], mass[2], | |
1116 | poleSize, p4); | |
1117 | this->getDaug(0)->init(daughters[0],p4[1]); | |
1118 | this->getDaug(1)->init(daughters[1],p4[0]); | |
1119 | this->getDaug(2)->init(daughters[2],p4[2]); | |
1120 | ok=true; | |
1121 | } | |
1122 | if ( !ok) { | |
1123 | report(ERROR,"EvtGen") << "Invalid pair of particle to generate a pole dist " | |
1124 | << whichTwo1 << " " << whichTwo2 | |
1125 | << endl<<"Will terminate."<<endl; | |
1126 | ::abort(); | |
1127 | } | |
1128 | } | |
1129 | ||
1130 | return weight; | |
1131 | } | |
1132 | ||
0ca57c2f | 1133 | void EvtParticle::makeDaughters(unsigned int ndaugstore, std::vector<EvtId> idVector) { |
1134 | ||
1135 | // Convert the STL vector method to use the array method for now, since the | |
1136 | // array method pervades most of the EvtGen code... | |
1137 | ||
1138 | unsigned int nVector = idVector.size(); | |
1139 | if (nVector < ndaugstore) { | |
1140 | report(ERROR,"EvtGen") << "Asking to make "<<ndaugstore<<" daughters when there " | |
1141 | << "are only "<<nVector<<" EvtId values available"<<endl; | |
1142 | return; | |
1143 | } | |
1144 | ||
1145 | EvtId *idArray=new EvtId[ndaugstore]; | |
1146 | unsigned int i; | |
1147 | for (i = 0; i < ndaugstore; i++) { | |
1148 | idArray[i] = idVector[i]; | |
1149 | } | |
1150 | ||
1151 | this->makeDaughters(ndaugstore, idArray); | |
1152 | ||
1153 | delete[] idArray; | |
1154 | } | |
1155 | ||
da0e9ce3 | 1156 | void EvtParticle::makeDaughters( unsigned int ndaugstore, EvtId *id){ |
1157 | ||
1158 | unsigned int i; | |
1159 | if ( _channel < 0 ) { | |
1160 | setChannel(0); | |
1161 | } | |
1162 | EvtParticle* pdaug; | |
1163 | if (_ndaug!=0 ){ | |
1164 | if (_ndaug!=ndaugstore){ | |
1165 | report(ERROR,"EvtGen") << "Asking to make a different number of " | |
1166 | << "daughters than what was previously created."<<endl; | |
1167 | report(ERROR,"EvtGen") << "Original parent:"<<EvtPDL::name(_id)<<endl; | |
1168 | for (size_t i=0;i<_ndaug;i++){ | |
1169 | report(ERROR,"EvtGen") << "Original daugther:"<<EvtPDL::name(getDaug(i)->getId())<<endl; | |
1170 | } | |
1171 | for (size_t i=0;i<ndaugstore;i++){ | |
1172 | report(ERROR,"EvtGen") << "New Daug:"<<EvtPDL::name(id[i])<<endl; | |
1173 | } | |
1174 | report(ERROR,"EvtGen") << "Will terminate."<<endl; | |
1175 | ::abort(); | |
1176 | } | |
1177 | } | |
1178 | else{ | |
1179 | for(i=0;i<ndaugstore;i++){ | |
1180 | pdaug=EvtParticleFactory::particleFactory(EvtPDL::getSpinType(id[i])); | |
1181 | pdaug->setId(id[i]); | |
1182 | pdaug->addDaug(this); | |
1183 | } | |
1184 | ||
1185 | } //else | |
1186 | } //makeDaughters | |
1187 | ||
1188 | ||
1189 | void EvtParticle::setDecayProb(double prob) { | |
1190 | ||
1191 | if ( _decayProb == 0 ) _decayProb=new double; | |
1192 | *_decayProb=prob; | |
1193 | } | |
0ca57c2f | 1194 | |
1195 | std::string EvtParticle::getName() { | |
1196 | ||
1197 | std::string theName = _id.getName(); | |
1198 | return theName; | |
1199 | ||
1200 | } |