Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtD0gammaDalitz.cpp
1 //--------------------------------------------------------------------------
2 //
3 // Module: EvtD0gammaDalitz.cc
4 //
5 // Modification history:
6 //
7 //    JGT     February 13, 2012         Module created
8 //
9 //------------------------------------------------------------------------
10
11 #include <cstdlib>
12 #include <string>
13 #include <cmath>
14
15 #include "EvtGenBase/EvtPatches.hh"
16
17 #include "EvtGenBase/EvtParticle.hh"
18 #include "EvtGenBase/EvtGenKine.hh"
19 #include "EvtGenBase/EvtPDL.hh"
20 #include "EvtGenBase/EvtReport.hh"
21 #include "EvtGenBase/EvtResonance.hh"
22 #include "EvtGenBase/EvtResonance2.hh"
23 #include "EvtGenModels/EvtD0gammaDalitz.hh"
24 #include "EvtGenBase/EvtConst.hh"
25 #include "EvtGenBase/EvtFlatte.hh"
26 #include "EvtGenBase/EvtDecayTable.hh"
27
28
29 // Initialize the static variables.
30 const EvtSpinType::spintype& EvtD0gammaDalitz::_SCALAR = EvtSpinType::SCALAR;
31 const EvtSpinType::spintype& EvtD0gammaDalitz::_VECTOR = EvtSpinType::VECTOR;
32 const EvtSpinType::spintype& EvtD0gammaDalitz::_TENSOR = EvtSpinType::TENSOR;
33
34 const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_EtaPic   = EvtDalitzReso::EtaPic;
35 const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_PicPicKK = EvtDalitzReso::PicPicKK;
36
37 const EvtDalitzReso::NumType& EvtD0gammaDalitz::_RBW   = EvtDalitzReso::RBW_CLEO_ZEMACH;
38 const EvtDalitzReso::NumType& EvtD0gammaDalitz::_GS    = EvtDalitzReso::GS_CLEO_ZEMACH;
39 const EvtDalitzReso::NumType& EvtD0gammaDalitz::_KMAT  = EvtDalitzReso::K_MATRIX;
40
41 const EvtCyclic3::Pair& EvtD0gammaDalitz::_AB = EvtCyclic3::AB;
42 const EvtCyclic3::Pair& EvtD0gammaDalitz::_AC = EvtCyclic3::AC;
43 const EvtCyclic3::Pair& EvtD0gammaDalitz::_BC = EvtCyclic3::BC;
44
45
46 EvtD0gammaDalitz::EvtD0gammaDalitz()
47 {
48   /* Empty constructor. */
49 }
50
51
52 EvtD0gammaDalitz::~EvtD0gammaDalitz()
53 {
54   /* Empty destructor. */
55 }
56
57
58 std::string EvtD0gammaDalitz::getName()
59 {
60   return "D0GAMMADALITZ";
61 }
62
63
64 EvtDecayBase* EvtD0gammaDalitz::clone()
65 {
66   return new EvtD0gammaDalitz;
67 }
68
69
70 void EvtD0gammaDalitz::init()
71 {
72   // check that there are 0 arguments
73   checkNArg(0);
74
75   // Check that this model is valid for the specified decay.
76   checkNDaug( 3 );
77   checkSpinParent  (    _SCALAR );
78   checkSpinDaughter( 0, _SCALAR );
79   checkSpinDaughter( 1, _SCALAR );
80   checkSpinDaughter( 2, _SCALAR );
81
82   // Get the values of the EvtId objects from the data files.
83   readPDGValues();
84
85   // Get the EvtId of the D0 and its 3 daughters.
86   getParentId();
87
88   EvtId dau[ 3 ];
89   for ( int index = 0; index < 3; index++ )
90   {
91     dau[ index ] = getDaug( index );
92   }
93
94   // Look for K0bar h+ h-. The order will be K[0SL] h+ h-
95   for ( int index = 0; index < 3; index++ )
96   {
97     if      ( ( dau[ index ] == _K0B ) || ( dau[ index ] == _KS ) || ( dau[ index ] == _KL ) )
98     {
99       _d1 = index;
100     }
101     else if ( ( dau[ index ] == _PIP ) || ( dau[ index ] == _KP ) )
102     {
103       _d2 = index;
104     }
105     else if ( ( dau[ index ] == _PIM ) || ( dau[ index ] == _KM ) )
106     {
107       _d3 = index;
108     }
109     else
110     {
111       reportInvalidAndExit();
112     }
113   }
114
115   // Check if we're dealing with Ks pi pi or with Ks K K.
116   _isKsPiPi = false;
117   if ( dau[ _d2 ] == _PIP || dau[ _d2 ] == _PIM )
118   {
119     _isKsPiPi = true;
120   }
121
122 }
123
124
125 void EvtD0gammaDalitz::initProbMax()
126 {
127   setProbMax( 5200. );
128 }
129
130
131 void EvtD0gammaDalitz::decay( EvtParticle* part )
132 {
133   // Check if the D is from a B+- -> D0 K+- decay with the appropriate model.
134   EvtParticle* parent = part->getParent(); // If there are no mistakes, should be B+ or B-.
135   if (parent != 0 && EvtDecayTable::getInstance()->getDecayFunc( parent )->getName() == "BTODDALITZCPK" )
136   {
137     EvtId parId = parent->getId();
138     if ( ( parId == _BP ) || ( parId == _BM ) ||
139          ( parId == _B0 ) || ( parId == _B0B) )
140     {
141       _bFlavor = parId;
142     }
143     else
144     {
145       reportInvalidAndExit();
146     }
147   }
148   else
149   {
150     reportInvalidAndExit();
151   }
152
153   // Read the D decay parameters from the B decay model.
154   // Gamma angle in rad.
155   double gamma = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 0 );
156   // Strong phase in rad.
157   double delta = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 1 );
158   // Ratio between B->D0K and B->D0barK
159   double rB    = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 2 );
160
161   // Same structure for all of these decays.
162   part->initializePhaseSpace( getNDaug(), getDaugs() );
163   EvtVector4R pA = part->getDaug( _d1 )->getP4();
164   EvtVector4R pB = part->getDaug( _d2 )->getP4();
165   EvtVector4R pC = part->getDaug( _d3 )->getP4();
166
167   // Squared invariant masses.
168   double mSqAB = ( pA + pB ).mass2();
169   double mSqAC = ( pA + pC ).mass2();
170   double mSqBC = ( pB + pC ).mass2();
171
172   EvtComplex amp( 1.0, 0.0 );
173
174   // Direct and conjugated amplitudes.
175   EvtComplex ampDir;
176   EvtComplex ampCnj;
177
178   if ( _isKsPiPi )
179   {
180     // Direct and conjugated Dalitz points.
181     EvtDalitzPoint pointDir( _mKs, _mPi, _mPi, mSqAB, mSqBC, mSqAC );
182     EvtDalitzPoint pointCnj( _mKs, _mPi, _mPi, mSqAC, mSqBC, mSqAB );
183
184     // Direct and conjugated amplitudes.
185     ampDir = dalitzKsPiPi( pointDir );
186     ampCnj = dalitzKsPiPi( pointCnj );
187   }
188   else
189   {
190     // Direct and conjugated Dalitz points.
191     EvtDalitzPoint pointDir( _mKs, _mK, _mK, mSqAB, mSqBC, mSqAC );
192     EvtDalitzPoint pointCnj( _mKs, _mK, _mK, mSqAC, mSqBC, mSqAB );
193
194     // Direct and conjugated amplitudes.
195     ampDir = dalitzKsKK( pointDir );
196     ampCnj = dalitzKsKK( pointCnj );
197   }
198
199   if ( _bFlavor == _BP || _bFlavor == _B0 )
200   {
201     amp = ampCnj + rB * exp( EvtComplex( 0., delta + gamma ) ) * ampDir;
202   }
203   else
204   {
205     amp = ampDir + rB * exp( EvtComplex( 0., delta - gamma ) ) * ampCnj;
206   }
207
208   vertex( amp );
209
210   return;
211
212 }
213
214
215 EvtComplex EvtD0gammaDalitz::dalitzKsPiPi( const EvtDalitzPoint& point ) const
216 {
217   static const EvtDalitzPlot plot( _mKs, _mPi, _mPi, _mD0 );
218
219   EvtComplex amp = 0.;
220
221   // This corresponds to relativistic Breit-Wigner distributions. Not K-matrix.
222   // Defining resonances.
223   static EvtDalitzReso KStarm      ( plot, _BC, _AC, _VECTOR, 0.893606, 0.0463407, _RBW );
224   static EvtDalitzReso KStarp      ( plot, _BC, _AB, _VECTOR, 0.893606, 0.0463407, _RBW );
225   static EvtDalitzReso rho0        ( plot, _AC, _BC, _VECTOR, 0.7758  , 0.1464   , _GS  );
226   static EvtDalitzReso omega       ( plot, _AC, _BC, _VECTOR, 0.78259 , 0.00849  , _RBW );
227   static EvtDalitzReso f0_980      ( plot, _AC, _BC, _SCALAR, 0.975   , 0.044    , _RBW );
228   static EvtDalitzReso f0_1370     ( plot, _AC, _BC, _SCALAR, 1.434   , 0.173    , _RBW );
229   static EvtDalitzReso f2_1270     ( plot, _AC, _BC, _TENSOR, 1.2754  , 0.1851   , _RBW );
230   static EvtDalitzReso K0Starm_1430( plot, _BC, _AC, _SCALAR, 1.459   , 0.175    , _RBW );
231   static EvtDalitzReso K0Starp_1430( plot, _BC, _AB, _SCALAR, 1.459   , 0.175    , _RBW );
232   static EvtDalitzReso K2Starm_1430( plot, _BC, _AC, _TENSOR, 1.4256  , 0.0985   , _RBW );
233   static EvtDalitzReso K2Starp_1430( plot, _BC, _AB, _TENSOR, 1.4256  , 0.0985   , _RBW );
234   static EvtDalitzReso sigma       ( plot, _AC, _BC, _SCALAR, 0.527699, 0.511861 , _RBW );
235   static EvtDalitzReso sigma2      ( plot, _AC, _BC, _SCALAR, 1.03327 , 0.0987890, _RBW );
236   static EvtDalitzReso KStarm_1680 ( plot, _BC, _AC, _VECTOR, 1.677   , 0.205    , _RBW );
237
238   // Adding terms to the amplitude with their corresponding amplitude and phase terms.
239   amp += EvtComplex(   .848984 ,   .893618  );
240   amp += EvtComplex( -1.16356  ,  1.19933   ) * KStarm      .evaluate( point );
241   amp += EvtComplex(   .106051 , - .118513  ) * KStarp      .evaluate( point );
242   amp += EvtComplex(  1.0      ,  0.0       ) * rho0        .evaluate( point );
243   amp += EvtComplex( - .0249569,   .0388072 ) * omega       .evaluate( point );
244   amp += EvtComplex( - .423586 , - .236099  ) * f0_980      .evaluate( point );
245   amp += EvtComplex( -2.16486  ,  3.62385   ) * f0_1370     .evaluate( point );
246   amp += EvtComplex(   .217748 , - .133327  ) * f2_1270     .evaluate( point );
247   amp += EvtComplex(  1.62128  ,  1.06816   ) * K0Starm_1430.evaluate( point );
248   amp += EvtComplex(   .148802 ,   .0897144 ) * K0Starp_1430.evaluate( point );
249   amp += EvtComplex(  1.15489  , - .773363  ) * K2Starm_1430.evaluate( point );
250   amp += EvtComplex(   .140865 , - .165378  ) * K2Starp_1430.evaluate( point );
251   amp += EvtComplex( -1.55556  , - .931685  ) * sigma       .evaluate( point );
252   amp += EvtComplex( - .273791 , - .0535596 ) * sigma2      .evaluate( point );
253   amp += EvtComplex( -1.69720  ,   .128038  ) * KStarm_1680 .evaluate( point );
254
255   return amp;
256 }
257
258
259 EvtComplex EvtD0gammaDalitz::dalitzKsKK( const EvtDalitzPoint& point ) const
260 {
261   static const EvtDalitzPlot plot( _mKs, _mK, _mK, _mD0 );
262
263   // Defining resonances.
264   static EvtDalitzReso a00_980 ( plot, _AC, _BC, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
265   static EvtDalitzReso phi     ( plot, _AC, _BC, _VECTOR, 1.01943,       .00459319    , _RBW      );
266   static EvtDalitzReso a0p_980 ( plot, _AC, _AB, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
267   static EvtDalitzReso f0_1370 ( plot, _AC, _BC, _SCALAR, 1.350  ,       .265         , _RBW      );
268   static EvtDalitzReso a0m_980 ( plot, _AB, _AC, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
269   static EvtDalitzReso f0_980  ( plot, _AC, _BC, _SCALAR, 0.965  , _RBW, .695   , .165, _PicPicKK );
270   static EvtDalitzReso f2_1270 ( plot, _AC, _BC, _TENSOR, 1.2754 ,       .1851        , _RBW      );
271   static EvtDalitzReso a00_1450( plot, _AC, _BC, _SCALAR, 1.474  ,       .265         , _RBW      );
272   static EvtDalitzReso a0p_1450( plot, _AC, _AB, _SCALAR, 1.474  ,       .265         , _RBW      );
273   static EvtDalitzReso a0m_1450( plot, _AB, _AC, _SCALAR, 1.474  ,       .265         , _RBW      );
274
275   // Adding terms to the amplitude with their corresponding amplitude and phase terms.
276   EvtComplex amp( 0., 0. ); // Phase space amplitude.
277   amp += EvtComplex( 1.0          , 0.0        ) * a00_980 .evaluate( point );
278   amp += EvtComplex( -.126314     ,  .188701   ) * phi     .evaluate( point );
279   amp += EvtComplex( -.561428     ,  .0135338  ) * a0p_980 .evaluate( point );
280   amp += EvtComplex(  .035        , -.00110488 ) * f0_1370 .evaluate( point );
281   amp += EvtComplex( -.0872735    ,  .0791190  ) * a0m_980 .evaluate( point );
282   amp += EvtComplex( 0.           , 0.         ) * f0_980  .evaluate( point );
283   amp += EvtComplex(  .257341     , -.0408343  ) * f2_1270 .evaluate( point );
284   amp += EvtComplex( -.0614342    , -.649930   ) * a00_1450.evaluate( point );
285   amp += EvtComplex( -.104629     ,  .830120   ) * a0p_1450.evaluate( point );
286   amp += EvtComplex( 0.           , 0.         ) * a0m_1450.evaluate( point );
287
288   return 2.8 * amp; // Multiply by 2.8 in order to reuse the same probmax as Ks pi pi.
289 }
290
291
292 void EvtD0gammaDalitz::readPDGValues()
293 {
294   // Define the EvtIds.
295   _BP  = EvtPDL::getId( "B+"      );
296   _BM  = EvtPDL::getId( "B-"      );
297   _B0  = EvtPDL::getId( "B0"      );
298   _B0B = EvtPDL::getId( "anti-B0" );
299   _D0  = EvtPDL::getId( "D0"      );
300   _D0B = EvtPDL::getId( "anti-D0" );
301   _KM  = EvtPDL::getId( "K-"      );
302   _KP  = EvtPDL::getId( "K+"      );
303   _K0  = EvtPDL::getId( "K0"      );
304   _K0B = EvtPDL::getId( "anti-K0" );
305   _KL  = EvtPDL::getId( "K_L0"    );
306   _KS  = EvtPDL::getId( "K_S0"    );
307   _PIM = EvtPDL::getId( "pi-"     );
308   _PIP = EvtPDL::getId( "pi+"     );
309
310   // Read the relevant masses.
311   _mD0 = EvtPDL::getMass( _D0  );
312   _mKs = EvtPDL::getMass( _KS  );
313   _mPi = EvtPDL::getMass( _PIP );
314   _mK  = EvtPDL::getMass( _KP  );
315 }
316
317
318 void EvtD0gammaDalitz::reportInvalidAndExit() const
319 {
320   report( ERROR, "EvtD0gammaDalitz" ) << "EvtD0gammaDalitz: Invalid mode." << std::endl;
321   exit( 1 );
322 }