1 //--------------------------------------------------------------------------
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.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtLineShape.cc
13 // Description: Store particle properties for one particle.
15 // Modification history:
17 // Lange March 10, 2001 Module created
18 // Dvoretskii June 03, 2002 Reimplemented rollMass()
20 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenBase/EvtPatches.hh"
25 #include "EvtGenBase/EvtPredGen.hh"
27 #include "EvtGenBase/EvtRelBreitWignerBarrierFact.hh"
28 #include "EvtGenBase/EvtTwoBodyVertex.hh"
29 #include "EvtGenBase/EvtPropBreitWignerRel.hh"
30 #include "EvtGenBase/EvtPDL.hh"
31 #include "EvtGenBase/EvtAmpPdf.hh"
32 #include "EvtGenBase/EvtMassAmp.hh"
33 #include "EvtGenBase/EvtSpinType.hh"
34 #include "EvtGenBase/EvtIntervalFlatPdf.hh"
36 EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact() {
40 EvtRelBreitWignerBarrierFact::~EvtRelBreitWignerBarrierFact() {
43 EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(double mass, double width, double maxRange, EvtSpinType::spintype sp) :
44 EvtAbsLineShape(mass,width,maxRange,sp)
45 { // double mDaug1, double mDaug2, int l) {
47 _includeDecayFact=true;
48 _includeBirthFact=true;
57 double maxdelta = 15.0*width;
59 if ( maxRange > 0.00001 ) {
60 _massMax=mass+maxdelta;
61 _massMin=mass-maxRange;
64 _massMax=mass+maxdelta;
65 _massMin=mass-15.0*width;
68 _massMax=mass+maxdelta;
69 if ( _massMin< 0. ) _massMin=0.;
72 EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(const EvtRelBreitWignerBarrierFact& x) :
77 _blattDecay=x._blattDecay;
78 _blattBirth=x._blattBirth;
79 _maxRange=x._maxRange;
80 _includeDecayFact=x._includeDecayFact;
81 _includeBirthFact=x._includeBirthFact;
82 _errorCond=x._errorCond;
86 EvtRelBreitWignerBarrierFact& EvtRelBreitWignerBarrierFact::operator=(const EvtRelBreitWignerBarrierFact& x) {
92 _blattDecay=x._blattDecay;
93 _blattBirth=x._blattBirth;
94 _maxRange=x._maxRange;
95 _includeDecayFact=x._includeDecayFact;
96 _includeBirthFact=x._includeBirthFact;
97 _errorCond=x._errorCond;
102 EvtAbsLineShape* EvtRelBreitWignerBarrierFact::clone() {
104 return new EvtRelBreitWignerBarrierFact(*this);
108 double EvtRelBreitWignerBarrierFact::getMassProb(double mass, double massPar,int nDaug, double *massDau) {
111 //return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
112 if (nDaug!=2) return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
117 for (i=0; i<nDaug; i++) {
118 dTotMass+=massDau[i];
120 //report(INFO,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl;
121 // if ( (mass-dTotMass)<0.0001 ) return 0.;
122 //report(INFO,"EvtGen") << mass << " " << dTotMass << endl;
123 if ( (mass<dTotMass) ) return 0.;
125 if ( _width< 0.0001) return 1.;
127 if ( massPar>0.0000000001 ) {
128 if ( mass > massPar) return 0.;
131 if ( _errorCond ) return 0.;
133 // we did all the work in getRandMass
137 double EvtRelBreitWignerBarrierFact::getRandMass(EvtId *parId,int nDaug, EvtId *dauId, EvtId *othDaugId, double maxMass, double *dauMasses) {
138 if ( nDaug!=2) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
140 if ( _width< 0.00001) return _mass;
142 //first figure out L - take the lowest allowed.
144 EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
145 EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);
147 int t1=EvtSpinType::getSpin2(spinD1);
148 int t2=EvtSpinType::getSpin2(spinD2);
149 int t3=EvtSpinType::getSpin2(_spin);
154 // the user has overridden the partial wave to use.
155 for (unsigned int vC=0; vC<_userSetPW.size(); vC++) {
156 if ( dauId[0]==_userSetPWD1[vC] && dauId[1]==_userSetPWD2[vC] ) {
157 Lmin=2*_userSetPW[vC];
159 if ( dauId[0]==_userSetPWD2[vC] && dauId[1]==_userSetPWD1[vC] ) {
160 Lmin=2*_userSetPW[vC];
164 // allow for special cases.
167 //There are some things I don't know how to deal with
168 if ( t3>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
169 if ( t1>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
170 if ( t2>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
172 //figure the min and max allowwed "spins" for the daughters state
173 Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2));
175 assert(Lmin==0||Lmin==2||Lmin==4);
178 //double massD1=EvtPDL::getMeanMass(dauId[0]);
179 //double massD2=EvtPDL::getMeanMass(dauId[1]);
180 double massD1=dauMasses[0];
181 double massD2=dauMasses[1];
183 // I'm not sure how to define the vertex factor here - so retreat to nonRel code.
184 if ( (massD1+massD2)> _mass ) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
186 //parent vertex factor not yet implemented
187 double massOthD=-10.;
188 double massParent=-10.;
191 EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
192 EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
194 int tt1=EvtSpinType::getSpin2(spinOth);
195 int tt2=EvtSpinType::getSpin2(spinPar);
196 int tt3=EvtSpinType::getSpin2(_spin);
199 //figure the min and max allowwed "spins" for the daughters state
200 if ( (tt1<=4) && ( tt2<=4) ) {
201 birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2));
202 if (birthl<0) birthl=0;
204 massOthD=EvtPDL::getMeanMass(*othDaugId);
205 massParent=EvtPDL::getMeanMass(*parId);
209 // allow user to override
210 for (size_t vC=0; vC<_userSetBirthPW.size(); vC++) {
211 if ( *othDaugId==_userSetBirthOthD[vC] && *parId==_userSetBirthPar[vC]){
212 birthl=2*_userSetBirthPW[vC];
217 double massM=_massMax;
218 if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass;
220 //special case... if the parent mass is _fixed_ we can do a little better
221 //and only for a two body decay as that seems to be where we have problems
223 // Define relativistic propagator amplitude
225 EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
226 vd.set_f(_blattDecay);
227 EvtPropBreitWignerRel bw(_mass,_width);
228 EvtMassAmp amp(bw,vd);
231 if ( _includeDecayFact) {
233 amp.addDeathFactFF();
235 if ( massParent>-1.) {
236 if ( _includeBirthFact ) {
238 EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
239 vb.set_f(_blattBirth);
242 amp.addBirthFactFF();
247 EvtAmpPdf<EvtPoint1D> pdf(amp);
249 // Estimate maximum and create predicate for accept reject
252 double tempMaxLoc=_mass;
253 if ( maxMass>-0.5 && maxMass<_mass) tempMaxLoc=maxMass;
254 double tempMax=_massMax;
255 if ( maxMass>-0.5 && maxMass<_massMax) tempMax=maxMass;
256 double tempMinMass=_massMin;
257 if ( massD1+massD2 > _massMin) tempMinMass=massD1+massD2;
259 //redo sanity check - is there a solution to our problem.
260 //if not return an error condition that is caught by the
261 //mass prob calculation above.
262 if ( tempMinMass > tempMax ) {
267 if ( tempMaxLoc < tempMinMass) tempMaxLoc=tempMinMass;
269 double safetyFactor=1.2;
271 EvtPdfMax<EvtPoint1D> max(safetyFactor*pdf.evaluate(EvtPoint1D(tempMinMass,tempMax,tempMaxLoc)));
273 EvtPdfPred<EvtPoint1D> pred(pdf);
276 EvtIntervalFlatPdf flat(tempMinMass,tempMax);
277 EvtPdfGen<EvtPoint1D> gen(flat);
278 EvtPredGen<EvtPdfGen<EvtPoint1D>,EvtPdfPred<EvtPoint1D> > predgen(gen,pred);
280 EvtPoint1D point = predgen();
281 return point.value();