TEvtGen/EvtGenBase/EvtRelBreitWignerBarrierFact.cxx

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