]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/EvtGenBase/EvtParticle.cpp
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates EvtGen Code
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtParticle.cpp
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"
0ca57c2f 48#include "EvtGenBase/EvtStatus.hh"
da0e9ce3 49
50using std::endl;
51
52
53
54EvtParticle::~EvtParticle() {
55 delete _decayProb;
56}
57
58EvtParticle::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
72void EvtParticle::setFirstOrNot() {
73 _first=0;
74}
75void EvtParticle::resetFirstOrNot() {
76 _first=1;
77}
78
79void EvtParticle::setChannel( int i ) {
80 _channel=i;
81}
82
83EvtParticle *EvtParticle::getDaug(int i) { return _daug[i]; }
84
85EvtParticle *EvtParticle::getParent() const { return _parent;}
86
87void EvtParticle::setLifetime(double tau){
88 _t=tau;
89}
90
91void EvtParticle::setLifetime(){
92 if (_genlifetime){
93 _t=-log(EvtRandom::Flat())*EvtPDL::getctau(getId());
94 }
95}
96
97double EvtParticle::getLifetime(){
98
99 return _t;
100}
101
102void EvtParticle::addDaug(EvtParticle *node) {
103 node->_daug[node->_ndaug++]=this;
104 _ndaug=0;
105 _parent=node;
106}
107
108
109int EvtParticle::firstornot() const { return _first;}
110
111EvtId EvtParticle::getId() const { return _id;}
112
0ca57c2f 113int EvtParticle::getPDGId() const {return EvtPDL::getStdHep(_id);}
114
da0e9ce3 115EvtSpinType::spintype EvtParticle::getSpinType() const
116 { return EvtPDL::getSpinType(_id);}
117
118int EvtParticle::getSpinStates() const
119 { return EvtSpinType::getSpinStates(EvtPDL::getSpinType(_id));}
120
121const EvtVector4R& EvtParticle::getP4() const { return _p;}
122
123int EvtParticle::getChannel() const { return _channel;}
124
125size_t EvtParticle::getNDaug() const { return _ndaug;}
126
127double EvtParticle::mass() const {
128
129 return _p.mass();
130}
131
132
133void EvtParticle::setDiagonalSpinDensity(){
134
135 _rhoForward.setDiag(getSpinStates());
136}
137
138void 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
158void 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
184void 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
213void 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
394void 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 464bool 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
508double 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
540void 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
552void EvtParticle::deleteTree(){
553
554 this->deleteDaughters();
555
556 delete this;
557
558}
559
560EvtVector4C 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
571EvtVector4C 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
582EvtVector4C 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
593EvtVector4C 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
604EvtDiracSpinor 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
615EvtDiracSpinor 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
626EvtDiracSpinor 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
637EvtDiracSpinor 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
648EvtTensor4C 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
659EvtTensor4C 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
671EvtRaritaSchwinger 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
682EvtRaritaSchwinger 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
695EvtVector4R 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
710EvtVector4R 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
727EvtVector4R EvtParticle::getP4Restframe() const {
728
729 return EvtVector4R(mass(),0.0,0.0,0.0);
730
731}
732
733EvtVector4R 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
756EvtParticle * 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
786void 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
820void 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
839void 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
875void 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
892void 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
919void 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
930std::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
953std::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
963void 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
998void init_vector( EvtParticle **part ){
999 *part = (EvtParticle *) new EvtVectorParticle;
1000}
1001
1002
1003void init_scalar( EvtParticle **part ){
1004 *part = (EvtParticle *) new EvtScalarParticle;
1005}
1006
1007void init_tensor( EvtParticle **part ){
1008 *part = (EvtParticle *) new EvtTensorParticle;
1009}
1010
1011void init_dirac( EvtParticle **part ){
1012 *part = (EvtParticle *) new EvtDiracParticle;
1013}
1014
1015void init_photon( EvtParticle **part ){
1016 *part = (EvtParticle *) new EvtPhotonParticle;
1017}
1018
1019void init_neutrino( EvtParticle **part ){
1020 *part = (EvtParticle *) new EvtNeutrinoParticle;
1021}
1022
1023void init_string( EvtParticle **part ){
1024 *part = (EvtParticle *) new EvtStringParticle;
1025}
1026
1027double 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 1133void 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 1156void 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
1189void EvtParticle::setDecayProb(double prob) {
1190
1191 if ( _decayProb == 0 ) _decayProb=new double;
1192 *_decayProb=prob;
1193}
0ca57c2f 1194
1195std::string EvtParticle::getName() {
1196
1197 std::string theName = _id.getName();
1198 return theName;
1199
1200}