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