00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <iostream>
00023 #include <fstream>
00024 #include <stdlib.h>
00025 #include <ctype.h>
00026 #include <math.h>
00027 #include <algorithm>
00028 #include "EvtGenBase/EvtAbsLineShape.hh"
00029 #include "EvtGenBase/EvtRandom.hh"
00030 #include "EvtGenBase/EvtComplex.hh"
00031 #include "EvtGenBase/EvtTwoBodyVertex.hh"
00032 #include "EvtGenBase/EvtPropBreitWigner.hh"
00033 #include "EvtGenBase/EvtPDL.hh"
00034 #include "EvtGenBase/EvtReport.hh"
00035 using std::fstream;
00036 using namespace std;
00037
00038 EvtAbsLineShape::EvtAbsLineShape() {
00039 }
00040
00041 EvtAbsLineShape::~EvtAbsLineShape() {
00042 }
00043
00044 EvtAbsLineShape::EvtAbsLineShape(double mass, double width, double maxRange, EvtSpinType::spintype sp) {
00045
00046 _includeDecayFact = false;
00047 _includeBirthFact = false;
00048 _applyFixForSP8 = false;
00049 _mass=mass;
00050 _width=width;
00051 _spin=sp;
00052 _maxRange=maxRange;
00053 double maxdelta=15.0*width;
00054
00055
00056
00057 if ( maxRange > 0.00001 ) {
00058 _massMax=mass+maxdelta;
00059 _massMin=mass-maxRange;
00060 }
00061 else{
00062 _massMax=mass+maxdelta;
00063 _massMin=mass-15.0*width;
00064 }
00065 if ( _massMin< 0. ) _massMin=0.;
00066 _massMax=mass+maxdelta;
00067 }
00068
00069 EvtAbsLineShape::EvtAbsLineShape(const EvtAbsLineShape& x){
00070
00071 _includeDecayFact = x._includeDecayFact;
00072 _includeBirthFact = x._includeBirthFact;
00073 _mass=x._mass;
00074 _massMax=x._massMax;
00075 _massMin=x._massMin;
00076 _width=x._width;
00077 _spin=x._spin;
00078 _maxRange=x._maxRange;
00079 _applyFixForSP8 = x._applyFixForSP8;
00080 }
00081
00082 EvtAbsLineShape& EvtAbsLineShape::operator=(const EvtAbsLineShape& x){
00083
00084 _includeDecayFact = x._includeDecayFact;
00085 _includeBirthFact = x._includeBirthFact;
00086 _mass=x._mass;
00087 _massMax=x._massMax;
00088 _massMin=x._massMin;
00089 _width=x._width;
00090 _spin=x._spin;
00091 _maxRange=x._maxRange;
00092 _applyFixForSP8 = x._applyFixForSP8;
00093 return *this;
00094 }
00095
00096 EvtAbsLineShape* EvtAbsLineShape::clone() {
00097
00098 return new EvtAbsLineShape(*this);
00099 }
00100
00101
00102 double EvtAbsLineShape::rollMass() {
00103
00104 double ymin, ymax;
00105 double temp;
00106
00107 if ( _width < 0.0001 ) {
00108 return _mass;
00109 }
00110 else{
00111 ymin = atan( 2.0*(_massMin-_mass)/_width);
00112 ymax = atan( 2.0*(_massMax-_mass)/_width);
00113
00114 temp= ( _mass + ((_width/2.0)*tan(EvtRandom::Flat(ymin,ymax))));
00115
00116 return temp;
00117 }
00118 }
00119 double EvtAbsLineShape::getRandMass(EvtId *parId, int nDaug, EvtId *dauId, EvtId *othDaugId, double maxMass, double *dauMasses) {
00120
00121 if ( _width< 0.0001) return _mass;
00122
00123
00124 if (maxMass>0&&maxMass<_massMin) {
00125 report(ERROR,"EvtGen") << "In EvtAbsLineShape::getRandMass"<<endl;
00126 report(ERROR,"EvtGen") << "Decaying particle "<<EvtPDL::name(*parId)
00127 << " with mass "<<maxMass<<endl;
00128 report(ERROR,"EvtGen") << " to particle"
00129 << " with a minimal mass of "<< _massMin<<endl;
00130 }
00131
00132 double mMin=_massMin;
00133 double mMax=_massMax;
00134 if ( maxMass>-0.5 && maxMass< mMax) mMax=maxMass;
00135 double ymin = atan( 2.0*(mMin-_mass)/_width);
00136 double ymax = atan( 2.0*(mMax-_mass)/_width);
00137
00138 loop:
00139 double themass = ( _mass + ((_width/2.0)*tan(EvtRandom::Flat(ymin,ymax))));
00140
00141 if(fabs(_addFactorPn)>0.00000001){
00142 double massOthD=EvtPDL::getMeanMass(*othDaugId);
00143 double massParent=EvtPDL::getMeanMass(*parId);
00144 double phsp,maxp,maxp1,maxp2;
00145 if(themass+massOthD <massParent ){
00146 EvtTwoBodyVertex vb(themass,massOthD,massParent,1.0);
00147 phsp = vb.pD();
00148 } else {return themass;}
00149
00150 if( (massOthD + mMax)< massParent){
00151 EvtTwoBodyVertex vb1(massOthD,mMax,massParent,1);
00152 EvtTwoBodyVertex vb2(massOthD,mMin,massParent,1);
00153 maxp = vb1.pD();
00154 maxp1 = pow(maxp,_addFactorPn*2.0);
00155 maxp = vb2.pD();
00156 maxp2= pow(maxp,_addFactorPn*2.0);
00157 maxp = max(maxp1,maxp2);
00158 }else {return themass;}
00159
00160 double wt = pow(phsp,_addFactorPn*2.0)/maxp;
00161 double rdm = EvtRandom::Flat(0.0,1.0);
00162
00163 if(rdm> wt ) goto loop;
00164
00165 }
00166
00167 return themass;
00168
00169 }
00170
00171 double EvtAbsLineShape::getMassProb(double mass, double massPar, int nDaug, double *massDau) {
00172
00173 double dTotMass=0.;
00174 if ( nDaug>1) {
00175 int i;
00176 for (i=0; i<nDaug; i++) {
00177 dTotMass+=massDau[i];
00178 }
00179
00180
00181 if ( (mass<dTotMass) ) return 0.;
00182 }
00183 if ( _width< 0.0001) return 1.;
00184
00185
00186 if ( massPar>0.0000000001 ) {
00187 if ( mass > massPar) return 0.;
00188 }
00189
00190
00191
00192
00193
00194
00195
00196
00197 return 1.0;
00198 }
00199
00200