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: EvtPropSLPole.cc
12 : //
13 : // Description: Routine to implement semileptonic decays according
14 : // to light cone sum rules
15 : //
16 : // Modification history:
17 : //
18 : // DJL April 23, 1998 Module created
19 : //
20 : //------------------------------------------------------------------------
21 : //
22 : #include "EvtGenBase/EvtPatches.hh"
23 : #include <stdlib.h>
24 : #include "EvtGenBase/EvtParticle.hh"
25 : #include "EvtGenBase/EvtGenKine.hh"
26 : #include "EvtGenBase/EvtPDL.hh"
27 : #include "EvtGenBase/EvtReport.hh"
28 : #include "EvtGenModels/EvtPropSLPole.hh"
29 : #include "EvtGenModels/EvtSLPoleFF.hh"
30 : #include "EvtGenBase/EvtSemiLeptonicScalarAmp.hh"
31 : #include "EvtGenBase/EvtSemiLeptonicVectorAmp.hh"
32 : #include "EvtGenBase/EvtSemiLeptonicTensorAmp.hh"
33 : #include "EvtGenBase/EvtIntervalFlatPdf.hh"
34 : #include "EvtGenBase/EvtScalarParticle.hh"
35 : #include "EvtGenBase/EvtVectorParticle.hh"
36 : #include "EvtGenBase/EvtTensorParticle.hh"
37 : #include "EvtGenBase/EvtTwoBodyVertex.hh"
38 : #include "EvtGenBase/EvtPropBreitWignerRel.hh"
39 : #include "EvtGenBase/EvtPDL.hh"
40 : #include "EvtGenBase/EvtAmpPdf.hh"
41 : #include "EvtGenBase/EvtMassAmp.hh"
42 : #include "EvtGenBase/EvtSpinType.hh"
43 : #include "EvtGenBase/EvtDecayTable.hh"
44 : #include <string>
45 :
46 0 : EvtPropSLPole::~EvtPropSLPole() {}
47 :
48 : std::string EvtPropSLPole::getName(){
49 :
50 0 : return "PROPSLPOLE";
51 :
52 : }
53 :
54 :
55 : EvtDecayBase* EvtPropSLPole::clone(){
56 :
57 0 : return new EvtPropSLPole;
58 :
59 0 : }
60 :
61 : void EvtPropSLPole::decay( EvtParticle *p ){
62 :
63 0 : if(! _isProbMaxSet){
64 :
65 0 : EvtId parnum,mesnum,lnum,nunum;
66 :
67 0 : parnum = getParentId();
68 0 : mesnum = getDaug(0);
69 0 : lnum = getDaug(1);
70 0 : nunum = getDaug(2);
71 :
72 0 : double mymaxprob = calcMaxProb(parnum,mesnum,
73 0 : lnum,nunum,SLPoleffmodel);
74 :
75 0 : setProbMax(mymaxprob);
76 :
77 0 : _isProbMaxSet = true;
78 :
79 0 : }
80 :
81 0 : double minKstMass = EvtPDL::getMinMass(p->getDaug(0)->getId());
82 0 : double maxKstMass = EvtPDL::getMaxMass(p->getDaug(0)->getId());
83 :
84 0 : EvtIntervalFlatPdf flat(minKstMass, maxKstMass);
85 0 : EvtPdfGen<EvtPoint1D> gen(flat);
86 0 : EvtPoint1D point = gen();
87 :
88 0 : double massKst = point.value();
89 :
90 0 : p->getDaug(0)->setMass(massKst);
91 0 : p->initializePhaseSpace(getNDaug(),getDaugs());
92 :
93 : // EvtVector4R p4meson = p->getDaug(0)->getP4();
94 :
95 0 : calcamp->CalcAmp(p,_amp2,SLPoleffmodel);
96 :
97 0 : EvtParticle *mesonPart = p->getDaug(0);
98 :
99 0 : double meson_BWAmp = calBreitWigner(mesonPart, point);
100 :
101 0 : int list[2];
102 0 : list[0]=0; list[1]=0;
103 0 : _amp2.vertex(0,0,_amp2.getAmp(list)*meson_BWAmp);
104 0 : list[0]=0; list[1]=1;
105 0 : _amp2.vertex(0,1,_amp2.getAmp(list)*meson_BWAmp);
106 :
107 0 : list[0]=1; list[1]=0;
108 0 : _amp2.vertex(1,0,_amp2.getAmp(list)*meson_BWAmp);
109 0 : list[0]=1; list[1]=1;
110 0 : _amp2.vertex(1,1,_amp2.getAmp(list)*meson_BWAmp);
111 :
112 0 : list[0]=2; list[1]=0;
113 0 : _amp2.vertex(2,0,_amp2.getAmp(list)*meson_BWAmp);
114 0 : list[0]=2; list[1]=1;
115 0 : _amp2.vertex(2,1,_amp2.getAmp(list)*meson_BWAmp);
116 :
117 :
118 : return;
119 :
120 0 : }
121 :
122 : void EvtPropSLPole::initProbMax(){
123 :
124 0 : _isProbMaxSet = false;
125 :
126 0 : return;
127 :
128 : }
129 :
130 :
131 : void EvtPropSLPole::init(){
132 :
133 0 : checkNDaug(3);
134 :
135 : //We expect the parent to be a scalar
136 : //and the daughters to be X lepton neutrino
137 :
138 0 : checkSpinParent(EvtSpinType::SCALAR);
139 0 : checkSpinDaughter(1,EvtSpinType::DIRAC);
140 0 : checkSpinDaughter(2,EvtSpinType::NEUTRINO);
141 :
142 0 : EvtSpinType::spintype mesontype=EvtPDL::getSpinType(getDaug(0));
143 :
144 0 : SLPoleffmodel = new EvtSLPoleFF(getNArg(),getArgs());
145 :
146 0 : if ( mesontype==EvtSpinType::SCALAR ) {
147 0 : calcamp = new EvtSemiLeptonicScalarAmp;
148 0 : }
149 0 : if ( mesontype==EvtSpinType::VECTOR ) {
150 0 : calcamp = new EvtSemiLeptonicVectorAmp;
151 0 : }
152 0 : if ( mesontype==EvtSpinType::TENSOR ) {
153 0 : calcamp = new EvtSemiLeptonicTensorAmp;
154 0 : }
155 :
156 0 : }
157 :
158 :
159 : double EvtPropSLPole::calBreitWignerBasic(double maxMass){
160 :
161 0 : if ( _width< 0.0001) return 1.0;
162 : //its not flat - but generated according to a BW
163 :
164 0 : double mMin=_massMin;
165 0 : double mMax=_massMax;
166 0 : if ( maxMass>-0.5 && maxMass< mMax) mMax=maxMass;
167 :
168 0 : double massGood = EvtRandom::Flat(mMin, mMax);
169 :
170 0 : double ampVal = sqrt(1.0/(pow(massGood-_mass, 2.0) + pow(_width, 2.0)/4.0));
171 :
172 : return ampVal;
173 :
174 0 : }
175 :
176 :
177 : double EvtPropSLPole::calBreitWigner(EvtParticle *pmeson, EvtPoint1D point){
178 :
179 0 : EvtId mesnum = pmeson->getId();
180 0 : double _mass = EvtPDL::getMeanMass(mesnum);
181 0 : double _width = EvtPDL::getWidth(mesnum);
182 0 : double _maxRange = EvtPDL::getMaxRange(mesnum);
183 0 : EvtSpinType::spintype mesontype=EvtPDL::getSpinType(mesnum);
184 0 : _includeDecayFact=true;
185 0 : _includeBirthFact=true;
186 0 : _spin = mesontype;
187 0 : _blatt = 3.0;
188 :
189 0 : double maxdelta = 15.0*_width;
190 :
191 0 : if ( _maxRange > 0.00001 ) {
192 0 : _massMax=_mass+maxdelta;
193 0 : _massMin=_mass-_maxRange;
194 0 : }
195 : else{
196 : _massMax=_mass+maxdelta;
197 0 : _massMin=_mass-15.0*_width;
198 : }
199 :
200 0 : _massMax=_mass+maxdelta;
201 0 : if ( _massMin< 0. ) _massMin=0.;
202 :
203 :
204 0 : EvtParticle* par=pmeson->getParent();
205 : double maxMass=-1.;
206 0 : if ( par != 0 ) {
207 0 : if ( par->hasValidP4() ) maxMass=par->mass();
208 0 : for ( size_t i=0;i<par->getNDaug();i++) {
209 0 : EvtParticle *tDaug=par->getDaug(i);
210 0 : if ( pmeson != tDaug )
211 0 : maxMass-=EvtPDL::getMinMass(tDaug->getId());
212 : }
213 0 : }
214 :
215 : EvtId *dauId=0;
216 : double *dauMasses=0;
217 0 : size_t nDaug = pmeson->getNDaug();
218 0 : if ( nDaug > 0) {
219 0 : dauId=new EvtId[nDaug];
220 0 : dauMasses=new double[nDaug];
221 0 : for (size_t j=0;j<nDaug;j++) {
222 0 : dauId[j]=pmeson->getDaug(j)->getId();
223 0 : dauMasses[j]=pmeson->getDaug(j)->mass();
224 : }
225 0 : }
226 : EvtId *parId=0;
227 : EvtId *othDaugId=0;
228 0 : EvtParticle *tempPar=pmeson->getParent();
229 0 : if (tempPar) {
230 0 : parId=new EvtId(tempPar->getId());
231 0 : if ( tempPar->getNDaug()==2 ) {
232 0 : if ( tempPar->getDaug(0) == pmeson ) othDaugId=new EvtId(tempPar->getDaug(1)->getId());
233 0 : else othDaugId=new EvtId(tempPar->getDaug(0)->getId());
234 : }
235 : }
236 :
237 0 : if ( nDaug!=2) return calBreitWignerBasic(maxMass);
238 :
239 0 : if ( _width< 0.00001) return 1.0;
240 :
241 : //first figure out L - take the lowest allowed.
242 :
243 0 : EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
244 0 : EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);
245 :
246 0 : int t1=EvtSpinType::getSpin2(spinD1);
247 0 : int t2=EvtSpinType::getSpin2(spinD2);
248 0 : int t3=EvtSpinType::getSpin2(_spin);
249 :
250 : int Lmin=-10;
251 :
252 : // allow for special cases.
253 0 : if (Lmin<-1 ) {
254 :
255 : //There are some things I don't know how to deal with
256 0 : if ( t3>4) return calBreitWignerBasic(maxMass);
257 0 : if ( t1>4) return calBreitWignerBasic(maxMass);
258 0 : if ( t2>4) return calBreitWignerBasic(maxMass);
259 :
260 : //figure the min and max allowwed "spins" for the daughters state
261 0 : Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2));
262 0 : if (Lmin<0) Lmin=0;
263 0 : assert(Lmin==0||Lmin==2||Lmin==4);
264 : }
265 :
266 : //double massD1=EvtPDL::getMeanMass(dauId[0]);
267 : //double massD2=EvtPDL::getMeanMass(dauId[1]);
268 0 : double massD1=dauMasses[0];
269 0 : double massD2=dauMasses[1];
270 :
271 : // I'm not sure how to define the vertex factor here - so retreat to nonRel code.
272 0 : if ( (massD1+massD2)> _mass ) return calBreitWignerBasic(maxMass);
273 :
274 : //parent vertex factor not yet implemented
275 : double massOthD=-10.;
276 : double massParent=-10.;
277 : int birthl=-10;
278 0 : if ( othDaugId) {
279 0 : EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
280 0 : EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
281 :
282 0 : int tt1=EvtSpinType::getSpin2(spinOth);
283 0 : int tt2=EvtSpinType::getSpin2(spinPar);
284 0 : int tt3=EvtSpinType::getSpin2(_spin);
285 :
286 : //figure the min and max allowwed "spins" for the daughters state
287 0 : if ( (tt1<=4) && ( tt2<=4) ) {
288 0 : birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2));
289 0 : if (birthl<0) birthl=0;
290 :
291 0 : massOthD=EvtPDL::getMeanMass(*othDaugId);
292 0 : massParent=EvtPDL::getMeanMass(*parId);
293 :
294 0 : }
295 :
296 0 : }
297 0 : double massM=_massMax;
298 0 : if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass;
299 :
300 : //special case... if the parent mass is _fixed_ we can do a little better
301 : //and only for a two body decay as that seems to be where we have problems
302 :
303 : // Define relativistic propagator amplitude
304 :
305 0 : EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
306 0 : vd.set_f(_blatt);
307 0 : EvtPropBreitWignerRel bw(_mass,_width);
308 0 : EvtMassAmp amp(bw,vd);
309 : // if ( _fixMassForMax) amp.fixUpMassForMax();
310 : // else std::cout << "problem problem\n";
311 0 : if ( _includeDecayFact) {
312 0 : amp.addDeathFact();
313 0 : amp.addDeathFactFF();
314 0 : }
315 0 : if ( massParent>-1.) {
316 0 : if ( _includeBirthFact ) {
317 :
318 0 : EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
319 0 : amp.setBirthVtx(vb);
320 0 : amp.addBirthFact();
321 0 : amp.addBirthFactFF();
322 0 : }
323 : }
324 :
325 0 : EvtAmpPdf<EvtPoint1D> pdf(amp);
326 :
327 0 : double ampVal = sqrt(pdf.evaluate(point));
328 :
329 0 : if ( parId) delete parId;
330 0 : if ( othDaugId) delete othDaugId;
331 0 : if ( dauId) delete [] dauId;
332 0 : if ( dauMasses) delete [] dauMasses;
333 :
334 : return ampVal;
335 :
336 0 : }
337 :
338 :
339 : double EvtPropSLPole::calcMaxProb( EvtId parent, EvtId meson,
340 : EvtId lepton, EvtId nudaug,
341 : EvtSemiLeptonicFF *FormFactors ) {
342 :
343 : //This routine takes the arguements parent, meson, and lepton
344 : //number, and a form factor model, and returns a maximum
345 : //probability for this semileptonic form factor model. A
346 : //brute force method is used. The 2D cos theta lepton and
347 : //q2 phase space is probed.
348 :
349 : //Start by declaring a particle at rest.
350 :
351 : //It only makes sense to have a scalar parent. For now.
352 : //This should be generalized later.
353 :
354 : EvtScalarParticle *scalar_part;
355 : EvtParticle *root_part;
356 :
357 0 : scalar_part=new EvtScalarParticle;
358 :
359 : //cludge to avoid generating random numbers!
360 0 : scalar_part->noLifeTime();
361 :
362 0 : EvtVector4R p_init;
363 :
364 0 : p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
365 0 : scalar_part->init(parent,p_init);
366 : root_part=(EvtParticle *)scalar_part;
367 : // root_part->set_type(EvtSpinType::SCALAR);
368 0 : root_part->setDiagonalSpinDensity();
369 :
370 : EvtParticle *daughter, *lep, *trino;
371 :
372 0 : EvtAmp amp;
373 :
374 0 : EvtId listdaug[3];
375 0 : listdaug[0] = meson;
376 0 : listdaug[1] = lepton;
377 0 : listdaug[2] = nudaug;
378 :
379 0 : amp.init(parent,3,listdaug);
380 :
381 0 : root_part->makeDaughters(3,listdaug);
382 0 : daughter=root_part->getDaug(0);
383 0 : lep=root_part->getDaug(1);
384 0 : trino=root_part->getDaug(2);
385 :
386 : EvtDecayBase *decayer;
387 0 : decayer = EvtDecayTable::getInstance()->getDecayFunc(daughter);
388 0 : if ( decayer ) {
389 0 : daughter->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs());
390 0 : for(int ii=0; ii<decayer->nRealDaughters(); ii++){
391 0 : daughter->getDaug(ii)->setMass(EvtPDL::getMeanMass(daughter->getDaug(ii)->getId()));
392 : }
393 0 : }
394 :
395 : //cludge to avoid generating random numbers!
396 0 : daughter->noLifeTime();
397 0 : lep->noLifeTime();
398 0 : trino->noLifeTime();
399 :
400 : //Initial particle is unpolarized, well it is a scalar so it is
401 : //trivial
402 0 : EvtSpinDensity rho;
403 0 : rho.setDiag(root_part->getSpinStates());
404 :
405 : double mass[3];
406 :
407 0 : double m = root_part->mass();
408 :
409 0 : EvtVector4R p4meson, p4lepton, p4nu, p4w;
410 : double q2max;
411 :
412 : double q2, elepton, plepton;
413 : int i,j;
414 : double erho,prho,costl;
415 :
416 : double maxfoundprob = 0.0;
417 : double prob = -10.0;
418 : int massiter;
419 :
420 0 : for (massiter=0;massiter<3;massiter++){
421 :
422 0 : mass[0] = EvtPDL::getMeanMass(meson);
423 0 : mass[1] = EvtPDL::getMeanMass(lepton);
424 0 : mass[2] = EvtPDL::getMeanMass(nudaug);
425 0 : if ( massiter==1 ) {
426 0 : mass[0] = EvtPDL::getMinMass(meson);
427 0 : }
428 0 : if ( massiter==2 ) {
429 0 : mass[0] = EvtPDL::getMaxMass(meson);
430 0 : if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
431 : }
432 :
433 0 : q2max = (m-mass[0])*(m-mass[0]);
434 :
435 : //loop over q2
436 :
437 0 : for (i=0;i<25;i++) {
438 0 : q2 = ((i+0.5)*q2max)/25.0;
439 :
440 0 : erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
441 :
442 0 : prho = sqrt(erho*erho-mass[0]*mass[0]);
443 :
444 0 : p4meson.set(erho,0.0,0.0,-1.0*prho);
445 0 : p4w.set(m-erho,0.0,0.0,prho);
446 :
447 : //This is in the W rest frame
448 0 : elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
449 0 : plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
450 :
451 0 : double probctl[3];
452 :
453 0 : for (j=0;j<3;j++) {
454 :
455 0 : costl = 0.99*(j - 1.0);
456 :
457 : //These are in the W rest frame. Need to boost out into
458 : //the B frame.
459 0 : p4lepton.set(elepton,0.0,
460 0 : plepton*sqrt(1.0-costl*costl),plepton*costl);
461 0 : p4nu.set(plepton,0.0,
462 0 : -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
463 :
464 0 : EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
465 0 : p4lepton=boostTo(p4lepton,boost);
466 0 : p4nu=boostTo(p4nu,boost);
467 :
468 : //Now initialize the daughters...
469 :
470 0 : daughter->init(meson,p4meson);
471 0 : lep->init(lepton,p4lepton);
472 0 : trino->init(nudaug,p4nu);
473 :
474 0 : calcamp->CalcAmp(root_part,amp,FormFactors);
475 :
476 0 : EvtPoint1D *point = new EvtPoint1D(mass[0]);
477 :
478 0 : double meson_BWAmp = calBreitWigner(daughter, *point);
479 :
480 0 : int list[2];
481 0 : list[0]=0; list[1]=0;
482 0 : amp.vertex(0,0,amp.getAmp(list)*meson_BWAmp);
483 0 : list[0]=0; list[1]=1;
484 0 : amp.vertex(0,1,amp.getAmp(list)*meson_BWAmp);
485 :
486 0 : list[0]=1; list[1]=0;
487 0 : amp.vertex(1,0,amp.getAmp(list)*meson_BWAmp);
488 0 : list[0]=1; list[1]=1;
489 0 : amp.vertex(1,1,amp.getAmp(list)*meson_BWAmp);
490 :
491 0 : list[0]=2; list[1]=0;
492 0 : amp.vertex(2,0,amp.getAmp(list)*meson_BWAmp);
493 0 : list[0]=2; list[1]=1;
494 0 : amp.vertex(2,1,amp.getAmp(list)*meson_BWAmp);
495 :
496 : //Now find the probability at this q2 and cos theta lepton point
497 : //and compare to maxfoundprob.
498 :
499 : //Do a little magic to get the probability!!
500 0 : prob = rho.normalizedProb(amp.getSpinDensity());
501 :
502 0 : probctl[j]=prob;
503 0 : }
504 :
505 : //probclt contains prob at ctl=-1,0,1.
506 : //prob=a+b*ctl+c*ctl^2
507 :
508 0 : double a=probctl[1];
509 0 : double b=0.5*(probctl[2]-probctl[0]);
510 0 : double c=0.5*(probctl[2]+probctl[0])-probctl[1];
511 :
512 : prob=probctl[0];
513 0 : if (probctl[1]>prob) prob=probctl[1];
514 0 : if (probctl[2]>prob) prob=probctl[2];
515 :
516 0 : if (fabs(c)>1e-20){
517 0 : double ctlx=-0.5*b/c;
518 0 : if (fabs(ctlx)<1.0){
519 0 : double probtmp=a+b*ctlx+c*ctlx*ctlx;
520 0 : if (probtmp>prob) prob=probtmp;
521 0 : }
522 :
523 0 : }
524 :
525 : //report(Severity::Debug,"EvtGen") << "prob,probctl:"<<prob<<" "
526 : // << probctl[0]<<" "
527 : // << probctl[1]<<" "
528 : // << probctl[2]<<endl;
529 :
530 0 : if ( prob > maxfoundprob ) {
531 : maxfoundprob = prob;
532 0 : }
533 :
534 0 : }
535 0 : if ( EvtPDL::getWidth(meson) <= 0.0 ) {
536 : //if the particle is narrow dont bother with changing the mass.
537 : massiter = 4;
538 0 : }
539 :
540 : }
541 0 : root_part->deleteTree();
542 :
543 0 : maxfoundprob *=1.1;
544 : return maxfoundprob;
545 :
546 0 : }
547 :
|
