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
00023 #include "EvtGenBase/EvtPatches.hh"
00024
00025 #include "EvtGenBase/EvtPredGen.hh"
00026
00027 #include "EvtGenBase/EvtRelBreitWignerBarrierFact.hh"
00028 #include "EvtGenBase/EvtTwoBodyVertex.hh"
00029 #include "EvtGenBase/EvtPropBreitWignerRel.hh"
00030 #include "EvtGenBase/EvtPDL.hh"
00031 #include "EvtGenBase/EvtAmpPdf.hh"
00032 #include "EvtGenBase/EvtMassAmp.hh"
00033 #include "EvtGenBase/EvtSpinType.hh"
00034 #include "EvtGenBase/EvtIntervalFlatPdf.hh"
00035 #include <algorithm>
00036 EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact() {
00037
00038 }
00039
00040 EvtRelBreitWignerBarrierFact::~EvtRelBreitWignerBarrierFact() {
00041 }
00042
00043 EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(double mass, double width, double maxRange, EvtSpinType::spintype sp) :
00044 EvtAbsLineShape(mass,width,maxRange,sp)
00045 {
00046
00047 _includeDecayFact=true;
00048 _includeBirthFact=true;
00049 _mass=mass;
00050 _width=width;
00051 _spin=sp;
00052 _blatt=3.0;
00053 _maxRange=maxRange;
00054 _errorCond=false;
00055
00056 double maxdelta = 15.0*width;
00057
00058 if ( maxRange > 0.00001 ) {
00059 _massMax=mass+maxdelta;
00060 _massMin=mass-maxRange;
00061 }
00062 else{
00063 _massMax=mass+maxdelta;
00064 _massMin=mass-15.0*width;
00065 }
00066
00067 _massMax=mass+maxdelta;
00068 if ( _massMin< 0. ) _massMin=0.;
00069 }
00070
00071 EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(const EvtRelBreitWignerBarrierFact& x) :
00072 EvtAbsLineShape(x)
00073 {
00074 _massMax=x._massMax;
00075 _massMin=x._massMin;
00076 _blatt=x._blatt;
00077 _maxRange=x._maxRange;
00078 _includeDecayFact=x._includeDecayFact;
00079 _includeBirthFact=x._includeBirthFact;
00080 _errorCond=x._errorCond;
00081
00082 }
00083
00084 EvtRelBreitWignerBarrierFact& EvtRelBreitWignerBarrierFact::operator=(const EvtRelBreitWignerBarrierFact& x) {
00085 _mass=x._mass;
00086 _width=x._width;
00087 _spin=x._spin;
00088 _massMax=x._massMax;
00089 _massMin=x._massMin;
00090 _blatt=x._blatt;
00091 _maxRange=x._maxRange;
00092 _includeDecayFact=x._includeDecayFact;
00093 _includeBirthFact=x._includeBirthFact;
00094 _errorCond=x._errorCond;
00095
00096 return *this;
00097 }
00098
00099 EvtAbsLineShape* EvtRelBreitWignerBarrierFact::clone() {
00100
00101 return new EvtRelBreitWignerBarrierFact(*this);
00102 }
00103
00104
00105 double EvtRelBreitWignerBarrierFact::getMassProb(double mass, double massPar,int nDaug, double *massDau) {
00106
00107 _errorCond=false;
00108
00109 if (nDaug!=2) return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
00110
00111 double dTotMass=0.;
00112
00113 int i;
00114 for (i=0; i<nDaug; i++) {
00115 dTotMass+=massDau[i];
00116 }
00117
00118
00119
00120 if ( (mass<dTotMass) ) return 0.;
00121
00122 if ( _width< 0.0001) return 1.;
00123
00124 if ( massPar>0.0000000001 ) {
00125 if ( mass > massPar) return 0.;
00126 }
00127
00128 if ( _errorCond ) return 0.;
00129
00130
00131 return 1.;
00132 }
00133
00134 double EvtRelBreitWignerBarrierFact::getRandMass(EvtId *parId,int nDaug, EvtId *dauId, EvtId *othDaugId, double maxMass, double *dauMasses) {
00135 if ( nDaug!=2) {
00136 if(fabs(_addFactorPn) >0.00000001) EvtAbsLineShape::addFactorPn(_addFactorPn);
00137 double themass = EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
00138 return themass;
00139 }
00140
00141 if ( _width< 0.00001) return _mass;
00142
00143
00144
00145 EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
00146 EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);
00147
00148 int t1=EvtSpinType::getSpin2(spinD1);
00149 int t2=EvtSpinType::getSpin2(spinD2);
00150 int t3=EvtSpinType::getSpin2(_spin);
00151
00152 int Lmin=-10;
00153
00154
00155
00156 for ( int vC=0; vC<_userSetPW.size(); vC++) {
00157 if ( dauId[0]==_userSetPWD1[vC] && dauId[1]==_userSetPWD2[vC] ) Lmin=2*_userSetPW[vC];
00158 if ( dauId[0]==_userSetPWD2[vC] && dauId[1]==_userSetPWD1[vC] ) Lmin=2*_userSetPW[vC];
00159 }
00160
00161
00162 if (Lmin<-1 ) {
00163
00164
00165 if ( t3>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
00166 if ( t1>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
00167 if ( t2>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
00168
00169
00170 Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2));
00171 if (Lmin<0) Lmin=0;
00172 assert(Lmin==0||Lmin==2||Lmin==4);
00173 }
00174
00175
00176
00177 double massD1=dauMasses[0];
00178 double massD2=dauMasses[1];
00179
00180
00181 if ( (massD1+massD2)> _mass ) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
00182
00183
00184 double massOthD=-10.;
00185 double massParent=-10.;
00186 int birthl=-10;
00187 if ( othDaugId) {
00188 EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
00189 EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
00190
00191 int tt1=EvtSpinType::getSpin2(spinOth);
00192 int tt2=EvtSpinType::getSpin2(spinPar);
00193 int tt3=EvtSpinType::getSpin2(_spin);
00194
00195
00196
00197 if ( (tt1<=4) && ( tt2<=4) ) {
00198 birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2));
00199 if (birthl<0) birthl=0;
00200
00201 massOthD=EvtPDL::getMeanMass(*othDaugId);
00202 massParent=EvtPDL::getMeanMass(*parId);
00203
00204 }
00205
00206
00207
00208 for (size_t vC=0; vC<_userSetBirthPW.size(); vC++) {
00209 if ( *othDaugId==_userSetBirthOthD[vC] && *parId==_userSetBirthPar[vC]){
00210 birthl=2*_userSetBirthPW[vC];
00211 }
00212 }
00213
00214 }
00215 double massM=_massMax;
00216 if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass;
00217
00218
00219
00220
00221
00222
00223 EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
00224 vd.set_f(_blatt);
00225 EvtPropBreitWignerRel bw(_mass,_width);
00226 EvtMassAmp amp(bw,vd);
00227 if ( _includeDecayFact) {
00228 amp.addDeathFact();
00229 amp.addDeathFactFF();
00230 }
00231
00232 if(fabs(_addFactorPn) >0.00000001){
00233
00234 amp.addFactorPn( _addFactorPn);
00235 }
00236 if ( massParent>-1.) {
00237 if ( _includeBirthFact ) {
00238
00239 EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
00240 if ( _applyFixForSP8 ) vb.set_f(_blatt);
00241 amp.setBirthVtx(vb);
00242 amp.addBirthFact();
00243 amp.addBirthFactFF();
00244 }
00245 }
00246
00247
00248
00249 EvtAmpPdf<EvtPoint1D> pdf(amp);
00250
00251
00252
00253
00254 double tempMaxLoc=_mass;
00255 if ( maxMass>-0.5 && maxMass<_mass) tempMaxLoc=maxMass;
00256 double tempMax=_massMax;
00257 if ( maxMass>-0.5 && maxMass<_massMax) tempMax=maxMass;
00258 double tempMinMass=_massMin;
00259 if ( massD1+massD2 > _massMin) tempMinMass=massD1+massD2;
00260
00261
00262
00263
00264 if ( tempMinMass > tempMax ) {
00265 _errorCond=true;
00266 return tempMinMass;
00267 }
00268
00269 if ( tempMaxLoc < tempMinMass) tempMaxLoc=tempMinMass;
00270
00271 double safetyFactor=1.2;
00272 if ( _applyFixForSP8 ) safetyFactor=1.4;
00273
00274 EvtPdfMax<EvtPoint1D> max(safetyFactor*pdf.evaluate(EvtPoint1D(tempMinMass,tempMax,tempMaxLoc)));
00275
00276 EvtPdfPred<EvtPoint1D> pred(pdf);
00277 pred.setMax(max);
00278
00279 EvtIntervalFlatPdf flat(tempMinMass,tempMax);
00280 EvtPdfGen<EvtPoint1D> gen(flat);
00281 EvtPredGen<EvtPdfGen<EvtPoint1D>,EvtPdfPred<EvtPoint1D> > predgen(gen,pred);
00282
00283 EvtPoint1D point = predgen();
00284 return point.value();
00285
00286 }
00287
00288
00289
00290
00291
00292
00293
00294
00295