Line data Source code
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: EvtLineShape.cc
12 : //
13 : // Description: Store particle properties for one particle.
14 : //
15 : // Modification history:
16 : //
17 : // Lange March 10, 2001 Module created
18 : // Dvoretskii June 03, 2002 Reimplemented rollMass()
19 : //
20 : //------------------------------------------------------------------------
21 : #include "EvtGenBase/EvtPatches.hh"
22 :
23 : #include "EvtGenBase/EvtPatches.hh"
24 :
25 : #include "EvtGenBase/EvtPredGen.hh"
26 :
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"
35 : #include <algorithm>
36 0 : EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact() {
37 :
38 0 : }
39 :
40 0 : EvtRelBreitWignerBarrierFact::~EvtRelBreitWignerBarrierFact() {
41 0 : }
42 :
43 : EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(double mass, double width, double maxRange, EvtSpinType::spintype sp) :
44 0 : EvtAbsLineShape(mass,width,maxRange,sp)
45 0 : { // double mDaug1, double mDaug2, int l) {
46 :
47 0 : _includeDecayFact=true;
48 0 : _includeBirthFact=true;
49 0 : _mass=mass;
50 0 : _width=width;
51 0 : _spin=sp;
52 0 : _blattDecay=3.0;
53 0 : _blattBirth=1.0;
54 0 : _maxRange=maxRange;
55 0 : _errorCond=false;
56 :
57 0 : double maxdelta = 15.0*width;
58 :
59 0 : if ( maxRange > 0.00001 ) {
60 0 : _massMax=mass+maxdelta;
61 0 : _massMin=mass-maxRange;
62 0 : }
63 : else{
64 : _massMax=mass+maxdelta;
65 0 : _massMin=mass-15.0*width;
66 : }
67 :
68 0 : _massMax=mass+maxdelta;
69 0 : if ( _massMin< 0. ) _massMin=0.;
70 0 : }
71 :
72 : EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(const EvtRelBreitWignerBarrierFact& x) :
73 0 : EvtAbsLineShape(x)
74 0 : {
75 0 : _massMax=x._massMax;
76 0 : _massMin=x._massMin;
77 0 : _blattDecay=x._blattDecay;
78 0 : _blattBirth=x._blattBirth;
79 0 : _maxRange=x._maxRange;
80 0 : _includeDecayFact=x._includeDecayFact;
81 0 : _includeBirthFact=x._includeBirthFact;
82 0 : _errorCond=x._errorCond;
83 :
84 0 : }
85 :
86 : EvtRelBreitWignerBarrierFact& EvtRelBreitWignerBarrierFact::operator=(const EvtRelBreitWignerBarrierFact& x) {
87 0 : _mass=x._mass;
88 0 : _width=x._width;
89 0 : _spin=x._spin;
90 0 : _massMax=x._massMax;
91 0 : _massMin=x._massMin;
92 0 : _blattDecay=x._blattDecay;
93 0 : _blattBirth=x._blattBirth;
94 0 : _maxRange=x._maxRange;
95 0 : _includeDecayFact=x._includeDecayFact;
96 0 : _includeBirthFact=x._includeBirthFact;
97 0 : _errorCond=x._errorCond;
98 :
99 0 : return *this;
100 : }
101 :
102 : EvtAbsLineShape* EvtRelBreitWignerBarrierFact::clone() {
103 :
104 0 : return new EvtRelBreitWignerBarrierFact(*this);
105 0 : }
106 :
107 :
108 : double EvtRelBreitWignerBarrierFact::getMassProb(double mass, double massPar,int nDaug, double *massDau) {
109 :
110 0 : _errorCond=false;
111 : //return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
112 0 : if (nDaug!=2) return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
113 :
114 : double dTotMass=0.;
115 :
116 : int i;
117 0 : for (i=0; i<nDaug; i++) {
118 0 : dTotMass+=massDau[i];
119 : }
120 : //report(Severity::Info,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl;
121 : // if ( (mass-dTotMass)<0.0001 ) return 0.;
122 : //report(Severity::Info,"EvtGen") << mass << " " << dTotMass << endl;
123 0 : if ( (mass<dTotMass) ) return 0.;
124 :
125 0 : if ( _width< 0.0001) return 1.;
126 :
127 0 : if ( massPar>0.0000000001 ) {
128 0 : if ( mass > massPar) return 0.;
129 : }
130 :
131 0 : if ( _errorCond ) return 0.;
132 :
133 : // we did all the work in getRandMass
134 0 : return 1.;
135 0 : }
136 :
137 : double EvtRelBreitWignerBarrierFact::getRandMass(EvtId *parId,int nDaug, EvtId *dauId, EvtId *othDaugId, double maxMass, double *dauMasses) {
138 0 : if ( nDaug!=2) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
139 :
140 0 : if ( _width< 0.00001) return _mass;
141 :
142 : //first figure out L - take the lowest allowed.
143 :
144 0 : EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
145 0 : EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);
146 :
147 0 : int t1=EvtSpinType::getSpin2(spinD1);
148 0 : int t2=EvtSpinType::getSpin2(spinD2);
149 0 : int t3=EvtSpinType::getSpin2(_spin);
150 :
151 : int Lmin=-10;
152 :
153 :
154 : // the user has overridden the partial wave to use.
155 0 : for (unsigned int vC=0; vC<_userSetPW.size(); vC++) {
156 0 : if ( dauId[0]==_userSetPWD1[vC] && dauId[1]==_userSetPWD2[vC] ) {
157 0 : Lmin=2*_userSetPW[vC];
158 0 : }
159 0 : if ( dauId[0]==_userSetPWD2[vC] && dauId[1]==_userSetPWD1[vC] ) {
160 0 : Lmin=2*_userSetPW[vC];
161 0 : }
162 : }
163 :
164 : // allow for special cases.
165 0 : if (Lmin<-1 ) {
166 :
167 : //There are some things I don't know how to deal with
168 0 : if ( t3>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
169 0 : if ( t1>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
170 0 : if ( t2>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
171 :
172 : //figure the min and max allowwed "spins" for the daughters state
173 0 : Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2));
174 0 : if (Lmin<0) Lmin=0;
175 0 : assert(Lmin==0||Lmin==2||Lmin==4);
176 : }
177 :
178 : //double massD1=EvtPDL::getMeanMass(dauId[0]);
179 : //double massD2=EvtPDL::getMeanMass(dauId[1]);
180 0 : double massD1=dauMasses[0];
181 0 : double massD2=dauMasses[1];
182 :
183 : // I'm not sure how to define the vertex factor here - so retreat to nonRel code.
184 0 : if ( (massD1+massD2)> _mass ) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
185 :
186 : //parent vertex factor not yet implemented
187 : double massOthD=-10.;
188 : double massParent=-10.;
189 : int birthl=-10;
190 0 : if ( othDaugId) {
191 0 : EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
192 0 : EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
193 :
194 0 : int tt1=EvtSpinType::getSpin2(spinOth);
195 0 : int tt2=EvtSpinType::getSpin2(spinPar);
196 0 : int tt3=EvtSpinType::getSpin2(_spin);
197 :
198 :
199 : //figure the min and max allowwed "spins" for the daughters state
200 0 : if ( (tt1<=4) && ( tt2<=4) ) {
201 0 : birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2));
202 0 : if (birthl<0) birthl=0;
203 :
204 0 : massOthD=EvtPDL::getMeanMass(*othDaugId);
205 0 : massParent=EvtPDL::getMeanMass(*parId);
206 :
207 0 : }
208 :
209 : // allow user to override
210 0 : for (size_t vC=0; vC<_userSetBirthPW.size(); vC++) {
211 0 : if ( *othDaugId==_userSetBirthOthD[vC] && *parId==_userSetBirthPar[vC]){
212 0 : birthl=2*_userSetBirthPW[vC];
213 0 : }
214 : }
215 :
216 0 : }
217 0 : double massM=_massMax;
218 0 : if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass;
219 :
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
222 :
223 : // Define relativistic propagator amplitude
224 :
225 0 : EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
226 0 : vd.set_f(_blattDecay);
227 0 : EvtPropBreitWignerRel bw(_mass,_width);
228 0 : EvtMassAmp amp(bw,vd);
229 :
230 :
231 0 : if ( _includeDecayFact) {
232 0 : amp.addDeathFact();
233 0 : amp.addDeathFactFF();
234 0 : }
235 0 : if ( massParent>-1.) {
236 0 : if ( _includeBirthFact ) {
237 :
238 0 : EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
239 0 : vb.set_f(_blattBirth);
240 0 : amp.setBirthVtx(vb);
241 0 : amp.addBirthFact();
242 0 : amp.addBirthFactFF();
243 0 : }
244 : }
245 :
246 :
247 0 : EvtAmpPdf<EvtPoint1D> pdf(amp);
248 :
249 : // Estimate maximum and create predicate for accept reject
250 :
251 :
252 0 : double tempMaxLoc=_mass;
253 0 : if ( maxMass>-0.5 && maxMass<_mass) tempMaxLoc=maxMass;
254 0 : double tempMax=_massMax;
255 0 : if ( maxMass>-0.5 && maxMass<_massMax) tempMax=maxMass;
256 0 : double tempMinMass=_massMin;
257 0 : if ( massD1+massD2 > _massMin) tempMinMass=massD1+massD2;
258 :
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 0 : if ( tempMinMass > tempMax ) {
263 0 : _errorCond=true;
264 0 : return tempMinMass;
265 : }
266 :
267 0 : if ( tempMaxLoc < tempMinMass) tempMaxLoc=tempMinMass;
268 :
269 : double safetyFactor=1.2;
270 :
271 0 : EvtPdfMax<EvtPoint1D> max(safetyFactor*pdf.evaluate(EvtPoint1D(tempMinMass,tempMax,tempMaxLoc)));
272 :
273 0 : EvtPdfPred<EvtPoint1D> pred(pdf);
274 0 : pred.setMax(max);
275 :
276 0 : EvtIntervalFlatPdf flat(tempMinMass,tempMax);
277 0 : EvtPdfGen<EvtPoint1D> gen(flat);
278 0 : EvtPredGen<EvtPdfGen<EvtPoint1D>,EvtPdfPred<EvtPoint1D> > predgen(gen,pred);
279 :
280 0 : EvtPoint1D point = predgen();
281 0 : return point.value();
282 :
283 0 : };
284 :
285 :
286 :
287 :
288 :
289 :
290 :
291 :
292 :
|
