Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtDMix.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: EvtDMix.cc
12 //
13 // Description: Routine to decay a particle according th phase space
14 //
15 // Modification history:
16 //
17 //    RYD       January 8, 1997       Module created
18 //
19 //------------------------------------------------------------------------
20 //
21 #include "EvtGenBase/EvtPatches.hh"
22 #include <stdlib.h>
23 #include "EvtGenBase/EvtParticle.hh"
24 #include "EvtGenBase/EvtGenKine.hh"
25 #include "EvtGenBase/EvtPDL.hh"
26 #include "EvtGenBase/EvtReport.hh"
27 #include "EvtGenModels/EvtDMix.hh"
28 #include "EvtGenBase/EvtRandom.hh"
29 #include <string>
30
31 EvtDMix::~EvtDMix() {}
32
33 std::string EvtDMix::getName(){
34
35   return "DMIX";     
36
37 }
38
39 EvtDecayBase* EvtDMix::clone(){
40
41   return new EvtDMix;
42
43 }
44
45
46 void EvtDMix::init(){
47
48   // check that there are 0 arguments
49   checkNArg(3);
50   _rd=getArg(0);
51   _xpr=getArg(1);
52   _ypr=getArg(2);
53 }
54
55 void EvtDMix::initProbMax(){
56
57   noProbMax();
58
59 }
60
61 void EvtDMix::decay( EvtParticle *p ){
62
63   //unneeded - lange - may13-02
64   //if ( p->getNDaug() != 0 ) {
65     //Will end up here because maxrate multiplies by 1.2
66   //  report(DEBUG,"EvtGen") << "In EvtDMix: has "
67   //                       <<" daugthers should not be here!"<<endl;
68   //  return;
69   //}
70
71   p->initializePhaseSpace(getNDaug(),getDaugs());
72
73   double ctau=EvtPDL::getctau(p->getId());
74   if ( ctau==0. ) return;
75
76
77   double pdf, random, gt, weight;
78
79   double maxPdf=_rd + sqrt(_rd)*_ypr*50. + 2500.0*(_xpr*_xpr+_ypr*_ypr)/4.0;
80   bool keepGoing=true;
81   while ( keepGoing ) {
82     random=EvtRandom::Flat();
83     gt=-log(random);
84     weight=random;
85     pdf=_rd + sqrt(_rd)*_ypr*gt + gt*gt*(_xpr*_xpr+_ypr*_ypr)/4.0;
86     pdf*=exp(-1.0*gt);
87     pdf/=weight;
88     if ( pdf > maxPdf )
89       std::cout << pdf << " " << weight << " " << maxPdf << " " << gt << std::endl;
90     if ( pdf > maxPdf*EvtRandom::Flat()) keepGoing=false;
91   }
92
93
94   p->setLifetime(gt*ctau);
95
96   return ;
97 }
98
99