TEvtGen/EvtGen/EvtGenBase/EvtRelBreitWignerBarrierFact.cpp

   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   _blattDecay=3.0;
  53   _blattBirth=1.0;
  54   _maxRange=maxRange;
  55   _errorCond=false;
  56
  57   double maxdelta = 15.0*width;
  58
  59   if ( maxRange > 0.00001 ) {
  60     _massMax=mass+maxdelta;
  61     _massMin=mass-maxRange;
  62   }
  63   else{
  64     _massMax=mass+maxdelta;
  65     _massMin=mass-15.0*width;
  66   }
  67
  68   _massMax=mass+maxdelta;
  69   if ( _massMin< 0. ) _massMin=0.;
  70 }
  71
  72 EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(const EvtRelBreitWignerBarrierFact& x) :
  73   EvtAbsLineShape(x)
  74 {
  75   _massMax=x._massMax;
  76   _massMin=x._massMin;
  77   _blattDecay=x._blattDecay;
  78   _blattBirth=x._blattBirth;
  79   _maxRange=x._maxRange;
  80   _includeDecayFact=x._includeDecayFact;
  81   _includeBirthFact=x._includeBirthFact;
  82   _errorCond=x._errorCond;
  83
  84 }
  85
  86 EvtRelBreitWignerBarrierFact& EvtRelBreitWignerBarrierFact::operator=(const EvtRelBreitWignerBarrierFact& x) {
  87   _mass=x._mass;
  88   _width=x._width;
  89   _spin=x._spin;
  90   _massMax=x._massMax;
  91   _massMin=x._massMin;
  92   _blattDecay=x._blattDecay;
  93   _blattBirth=x._blattBirth;
  94   _maxRange=x._maxRange;
  95   _includeDecayFact=x._includeDecayFact;
  96   _includeBirthFact=x._includeBirthFact;
  97   _errorCond=x._errorCond;
  98
  99   return *this;
 100 }
 101
 102 EvtAbsLineShape* EvtRelBreitWignerBarrierFact::clone() {
 103
 104   return new EvtRelBreitWignerBarrierFact(*this);
 105 }
 106
 107
 108 double EvtRelBreitWignerBarrierFact::getMassProb(double mass, double massPar,int nDaug, double *massDau) {
 109
 110   _errorCond=false;
 111   //return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
 112   if (nDaug!=2) return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
 113
 114   double dTotMass=0.;
 115
 116   int i;
 117   for (i=0; i<nDaug; i++) {
 118     dTotMass+=massDau[i];
 119   }
 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.;
 124
 125   if ( _width< 0.0001) return 1.;
 126
 127   if ( massPar>0.0000000001 ) {
 128     if ( mass > massPar) return 0.;
 129   }
 130
 131   if ( _errorCond ) return 0.;
 132
 133   // we did all the work in getRandMass
 134   return 1.;
 135 }
 136
 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);
 139
 140   if ( _width< 0.00001) return _mass;
 141
 142   //first figure out L - take the lowest allowed.
 143
 144   EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
 145   EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);
 146
 147   int t1=EvtSpinType::getSpin2(spinD1);
 148   int t2=EvtSpinType::getSpin2(spinD2);
 149   int t3=EvtSpinType::getSpin2(_spin);
 150
 151   int Lmin=-10;
 152
 153
 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];
 158     }
 159     if ( dauId[0]==_userSetPWD2[vC] &&  dauId[1]==_userSetPWD1[vC] ) {
 160       Lmin=2*_userSetPW[vC];
 161     }
 162   }
 163
 164   // allow for special cases.
 165   if (Lmin<-1 ) {
 166
 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);
 171
 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));
 174     if (Lmin<0) Lmin=0;
 175     assert(Lmin==0||Lmin==2||Lmin==4);
 176   }
 177
 178   //double massD1=EvtPDL::getMeanMass(dauId[0]);
 179   //double massD2=EvtPDL::getMeanMass(dauId[1]);
 180   double massD1=dauMasses[0];
 181   double massD2=dauMasses[1];
 182
 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);
 185
 186   //parent vertex factor not yet implemented
 187   double massOthD=-10.;
 188   double massParent=-10.;
 189   int birthl=-10;
 190   if ( othDaugId) {
 191     EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
 192     EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
 193
 194     int tt1=EvtSpinType::getSpin2(spinOth);
 195     int tt2=EvtSpinType::getSpin2(spinPar);
 196     int tt3=EvtSpinType::getSpin2(_spin);
 197
 198
 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;
 203
 204       massOthD=EvtPDL::getMeanMass(*othDaugId);
 205       massParent=EvtPDL::getMeanMass(*parId);
 206
 207     }
 208
 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];
 213       }
 214     }
 215
 216   }
 217   double massM=_massMax;
 218   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   EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
 226   vd.set_f(_blattDecay);
 227   EvtPropBreitWignerRel bw(_mass,_width);
 228   EvtMassAmp amp(bw,vd);
 229
 230
 231   if ( _includeDecayFact) {
 232     amp.addDeathFact();
 233     amp.addDeathFactFF();
 234   }
 235   if ( massParent>-1.) {
 236     if ( _includeBirthFact ) {
 237
 238       EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
 239       vb.set_f(_blattBirth);
 240       amp.setBirthVtx(vb);
 241       amp.addBirthFact();
 242       amp.addBirthFactFF();
 243     }
 244   }
 245
 246
 247   EvtAmpPdf<EvtPoint1D> pdf(amp);
 248
 249   // Estimate maximum and create predicate for accept reject
 250
 251
 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;
 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   if ( tempMinMass > tempMax ) {
 263     _errorCond=true;
 264     return tempMinMass;
 265   }
 266
 267   if ( tempMaxLoc < tempMinMass) tempMaxLoc=tempMinMass;
 268
 269   double safetyFactor=1.2;
 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