Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtAbsLineShape.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 //
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   _mass=mass;
47   _width=width;
48   _spin=sp;
49   _maxRange=maxRange;
50   double maxdelta=15.0*width;
51   //if ( width>0.001 ) {
52   //  if ( 5.0*width < 0.6 ) maxdelta = 0.6;
53   //}
54   if ( maxRange > 0.00001 ) {
55     _massMax=mass+maxdelta;
56     _massMin=mass-maxRange;
57   }
58   else{
59     _massMax=mass+maxdelta;
60     _massMin=mass-15.0*width;
61   }
62   if ( _massMin< 0. ) _massMin=0.;
63   _massMax=mass+maxdelta;
64 }
65
66 EvtAbsLineShape::EvtAbsLineShape(const EvtAbsLineShape& x){
67
68   _includeDecayFact = x._includeDecayFact;
69   _includeBirthFact = x._includeBirthFact;
70   _mass=x._mass;
71   _massMax=x._massMax;
72   _massMin=x._massMin;
73   _width=x._width;
74   _spin=x._spin;
75   _maxRange=x._maxRange;
76 }
77
78 EvtAbsLineShape& EvtAbsLineShape::operator=(const EvtAbsLineShape& x){
79
80   _includeDecayFact = x._includeDecayFact;
81   _includeBirthFact = x._includeBirthFact;
82   _mass=x._mass;
83   _massMax=x._massMax;
84   _massMin=x._massMin;
85   _width=x._width;
86   _spin=x._spin;
87   _maxRange=x._maxRange;
88   return *this; 
89 }
90
91 EvtAbsLineShape* EvtAbsLineShape::clone() {
92
93   return new EvtAbsLineShape(*this);
94 }
95
96
97 double EvtAbsLineShape::rollMass() {
98
99   double ymin, ymax;
100   double temp;
101
102   if ( _width < 0.0001 ) {
103     return _mass;
104   }
105   else{
106     ymin = atan( 2.0*(_massMin-_mass)/_width);
107     ymax = atan( 2.0*(_massMax-_mass)/_width);
108
109     temp= ( _mass + ((_width/2.0)*tan(EvtRandom::Flat(ymin,ymax))));
110
111     return temp;
112   }
113 }
114 double EvtAbsLineShape::getRandMass(EvtId *parId, int /* nDaug */, EvtId * /*dauId*/, EvtId */*othDaugId*/, double maxMass, double */*dauMasses*/) {
115
116   if ( _width< 0.0001) return _mass;
117   //its not flat - but generated according to a BW
118
119   if (maxMass>0&&maxMass<_massMin) {
120     report(DEBUG,"EvtGen") << "In EvtAbsLineShape::getRandMass:"<<endl;
121     report(DEBUG,"EvtGen") << "Cannot create a particle with a minimal mass of "
122                            << _massMin << " from a "<<EvtPDL::name(*parId)
123                            << " decay with available left-over mass-energy " << maxMass
124                            << ". Returning 0.0 mass. The rest of this decay chain will probably fail..." << endl;
125     return 0.0;
126   }
127
128   double mMin=_massMin;
129   double mMax=_massMax;
130   if ( maxMass>-0.5 && maxMass< mMax) mMax=maxMass;
131   double ymin = atan( 2.0*(mMin-_mass)/_width);
132   double ymax = atan( 2.0*(mMax-_mass)/_width);
133   
134   return ( _mass + ((_width/2.0)*tan(EvtRandom::Flat(ymin,ymax))));
135   //  return EvtRandom::Flat(_massMin,_massMax);
136 };
137
138 double EvtAbsLineShape::getMassProb(double mass, double massPar, int nDaug, double *massDau) {
139
140   double dTotMass=0.;
141   if ( nDaug>1) {
142     int i;
143     for (i=0; i<nDaug; i++) {
144       dTotMass+=massDau[i];
145     }
146     //report(INFO,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl;
147     //    if ( (mass-dTotMass)<0.0001 ) return 0.;
148     if ( (mass<dTotMass) ) return 0.;
149   }
150   if ( _width< 0.0001) return 1.;
151
152   // no parent - lets not panic
153   if ( massPar>0.0000000001 ) {
154     if ( mass > massPar) return 0.;
155   }
156   //Otherwise return the right value.
157   //Fortunately we have generated events according to a non-rel BW, so 
158   //just return..
159   //EvtPropBreitWigner bw(_mass,_width);
160   //EvtPropFactor<EvtTwoBodyVertex> f(bw);
161   //EvtComplex fm=f.eval(mass);
162   //EvtComplex fm0=f.eval(_mass);
163   //return (abs(fm)*abs(fm))/(abs(fm0)*abs(fm0));
164   return 1.0;
165 }
166
167