TEvtGen/EvtGenBase/EvtAbsLineShape.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 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <iostream>
  23 #include <stdlib.h>
  24 #include <ctype.h>
  25 #include <math.h>
  26 #include "EvtGenBase/EvtAbsLineShape.hh"
  27 #include "EvtGenBase/EvtRandom.hh"
  28 #include "EvtGenBase/EvtComplex.hh"
  29 #include "EvtGenBase/EvtTwoBodyVertex.hh"
  30 #include "EvtGenBase/EvtPropBreitWigner.hh"
  31 #include "EvtGenBase/EvtPDL.hh"
  32 #include "EvtGenBase/EvtReport.hh"
  33
  34 using namespace std;
  35
  36 EvtAbsLineShape::EvtAbsLineShape() {
  37 }
  38
  39 EvtAbsLineShape::~EvtAbsLineShape() {
  40 }
  41
  42 EvtAbsLineShape::EvtAbsLineShape(double mass, double width, double maxRange, EvtSpinType::spintype sp) {
  43
  44   _includeDecayFact = false;
  45   _includeBirthFact = false;
  46   _applyFixForSP8 = false;
  47   _mass=mass;
  48   _width=width;
  49   _spin=sp;
  50   _maxRange=maxRange;
  51   double maxdelta=15.0*width;
  52   //if ( width>0.001 ) {
  53   //  if ( 5.0*width < 0.6 ) maxdelta = 0.6;
  54   //}
  55   if ( maxRange > 0.00001 ) {
  56     _massMax=mass+maxdelta;
  57     _massMin=mass-maxRange;
  58   }
  59   else{
  60     _massMax=mass+maxdelta;
  61     _massMin=mass-15.0*width;
  62   }
  63   if ( _massMin< 0. ) _massMin=0.;
  64   _massMax=mass+maxdelta;
  65 }
  66
  67 EvtAbsLineShape::EvtAbsLineShape(const EvtAbsLineShape& x){
  68
  69   _includeDecayFact = x._includeDecayFact;
  70   _includeBirthFact = x._includeBirthFact;
  71   _mass=x._mass;
  72   _massMax=x._massMax;
  73   _massMin=x._massMin;
  74   _width=x._width;
  75   _spin=x._spin;
  76   _maxRange=x._maxRange;
  77   _applyFixForSP8 = x._applyFixForSP8;
  78 }
  79
  80 EvtAbsLineShape& EvtAbsLineShape::operator=(const EvtAbsLineShape& x){
  81
  82   _includeDecayFact = x._includeDecayFact;
  83   _includeBirthFact = x._includeBirthFact;
  84   _mass=x._mass;
  85   _massMax=x._massMax;
  86   _massMin=x._massMin;
  87   _width=x._width;
  88   _spin=x._spin;
  89   _maxRange=x._maxRange;
  90   _applyFixForSP8 = x._applyFixForSP8;
  91   return *this;
  92 }
  93
  94 EvtAbsLineShape* EvtAbsLineShape::clone() {
  95
  96   return new EvtAbsLineShape(*this);
  97 }
  98
  99
 100 double EvtAbsLineShape::rollMass() {
 101
 102   double ymin, ymax;
 103   double temp;
 104
 105   if ( _width < 0.0001 ) {
 106     return _mass;
 107   }
 108   else{
 109     ymin = atan( 2.0*(_massMin-_mass)/_width);
 110     ymax = atan( 2.0*(_massMax-_mass)/_width);
 111
 112     temp= ( _mass + ((_width/2.0)*tan(EvtRandom::Flat(ymin,ymax))));
 113
 114     return temp;
 115   }
 116 }
 117 double EvtAbsLineShape::getRandMass(EvtId *parId, int nDaug, EvtId *dauId, EvtId *othDaugId, double maxMass, double *dauMasses) {
 118
 119   if ( _width< 0.0001) return _mass;
 120   //its not flat - but generated according to a BW
 121
 122   if (maxMass>0&&maxMass<_massMin) {
 123     report(ERROR,"EvtGen") << "In EvtAbsLineShape::getRandMass"<<endl;
 124     report(ERROR,"EvtGen") << "Decaying particle "<<EvtPDL::name(*parId)
 125                            << " with mass "<<maxMass<<endl;
 126     report(ERROR,"EvtGen") << " to particle"
 127                            << " with a minimal mass of "<< _massMin<<endl;
 128   }
 129
 130   double mMin=_massMin;
 131   double mMax=_massMax;
 132   if ( maxMass>-0.5 && maxMass< mMax) mMax=maxMass;
 133   double ymin = atan( 2.0*(mMin-_mass)/_width);
 134   double ymax = atan( 2.0*(mMax-_mass)/_width);
 135
 136   return ( _mass + ((_width/2.0)*tan(EvtRandom::Flat(ymin,ymax))));
 137   //  return EvtRandom::Flat(_massMin,_massMax);
 138 }
 139
 140 double EvtAbsLineShape::getMassProb(double mass, double massPar, int nDaug, double *massDau) {
 141
 142   double dTotMass=0.;
 143   if ( nDaug>1) {
 144     int i;
 145     for (i=0; i<nDaug; i++) {
 146       dTotMass+=massDau[i];
 147     }
 148     //report(INFO,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl;
 149     //    if ( (mass-dTotMass)<0.0001 ) return 0.;
 150     if ( (mass<dTotMass) ) return 0.;
 151   }
 152   if ( _width< 0.0001) return 1.;
 153
 154   // no parent - lets not panic
 155   if ( massPar>0.0000000001 ) {
 156     if ( mass > massPar) return 0.;
 157   }
 158   //Otherwise return the right value.
 159   //Fortunately we have generated events according to a non-rel BW, so
 160   //just return..
 161   //EvtPropBreitWigner bw(_mass,_width);
 162   //EvtPropFactor<EvtTwoBodyVertex> f(bw);
 163   //EvtComplex fm=f.eval(mass);
 164   //EvtComplex fm0=f.eval(_mass);
 165   //return (abs(fm)*abs(fm))/(abs(fm0)*abs(fm0));
 166   return 1.0;
 167 }
 168
 169