BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Generator/BesEvtGen/BesEvtGen-00-03-58/src/EvtGen/EvtGenModels/EvtVPHOtoVISRHi.cc Source File

00001 //--------------------------------------------------------------------------
00002 //
00003 // Environment:
00004 //      This software is part of the EvtGen package developed jointly
00005 //      for the BaBar and CLEO collaborations.  If you use all or part
00006 //      of it, please give an appropriate acknowledgement.
00007 //
00008 // Copyright Information: See EvtGen/COPYRIGHT
00009 //      Copyright (C) 2004      Cornell
00010 //
00011 // Module: EvtVPHOtoVISR.cc
00012 //
00013 // Description: Routine to decay vpho -> (DDx) + ISR photon from 3.9 to 4.3 GeV, using CLEO-c data (Brian Lang)
00014 //
00015 // Modification history:
00016 //
00017 //    Ryd       March 20, 2004       Module created
00018 //
00019 //------------------------------------------------------------------------
00020 // 
00021 #include <stdlib.h>
00022 #include "EvtGenBase/EvtParticle.hh"
00023 #include "EvtGenBase/EvtGenKine.hh"
00024 #include "EvtGenBase/EvtPDL.hh"
00025 #include "EvtGenBase/EvtVector4C.hh"
00026 #include "EvtGenBase/EvtVector4R.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenBase/EvtRandom.hh"
00029 #include "EvtGenModels/EvtVPHOtoVISRHi.hh"
00030 #include <string>
00031 
00032 using std::endl;
00033 
00034 EvtVPHOtoVISRHi::~EvtVPHOtoVISRHi() {}
00035 
00036 void EvtVPHOtoVISRHi::getName(std::string&  model_name){
00037 
00038    model_name="VPHOTOVISRHI"; 
00039     
00040 }
00041 
00042 
00043 EvtDecayBase* EvtVPHOtoVISRHi::clone(){
00044 
00045   return new EvtVPHOtoVISRHi;
00046 
00047 }
00048 
00049 void EvtVPHOtoVISRHi::init(){
00050 
00051   // check that there are 0 or 1 arguments
00052   checkNArg(0,1);
00053 
00054   // check that there are 2 daughters
00055   checkNDaug(2);
00056 
00057   // check the parent and daughter spins
00058   checkSpinParent(EvtSpinType::VECTOR);
00059   checkSpinDaughter(0,EvtSpinType::VECTOR);
00060   checkSpinDaughter(1,EvtSpinType::PHOTON);
00061 }
00062 
00063 void EvtVPHOtoVISRHi::initProbMax() {
00064 
00065    setProbMax(20.0);
00066 
00067 }      
00068 
00069 void EvtVPHOtoVISRHi::decay( EvtParticle *p){
00070   //take photon along z-axis, either forward or backward.
00071   //Implement this as generating the photon momentum along 
00072   //the z-axis uniformly 
00073    double power=1;
00074    if (getNArg()==1) power=getArg(0);
00075    // define particle names
00076   static EvtId D0=EvtPDL::getId("D0");
00077   static EvtId D0B=EvtPDL::getId("anti-D0");
00078   static EvtId DP=EvtPDL::getId("D+");
00079   static EvtId DM=EvtPDL::getId("D-");
00080   static EvtId DSM=EvtPDL::getId("D_s-");
00081   static EvtId DSP=EvtPDL::getId("D_s+");
00082   static EvtId DSMS=EvtPDL::getId("D_s*-");
00083   static EvtId DSPS=EvtPDL::getId("D_s*+");
00084   static EvtId D0S=EvtPDL::getId("D*0");
00085   static EvtId D0BS=EvtPDL::getId("anti-D*0");
00086   static EvtId DPS=EvtPDL::getId("D*+");
00087   static EvtId DMS=EvtPDL::getId("D*-");
00088   // setup some parameters
00089   double w=p->mass();
00090   double s=w*w;
00091   double L=2.0*log(w/0.000511);
00092   double alpha=1/137.0;
00093   double beta=(L-1)*2.0*alpha/EvtConst::pi;
00094   // make sure only 2 or 3 body are present
00095   assert (p->getDaug(0)->getNDaug() == 2 || p->getDaug(0)->getNDaug() == 3);
00096 
00097   // determine minimum rest mass of parent
00098   double md1 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(0)->getId());
00099   double md2 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(1)->getId());
00100   double minResMass = md1+md2;
00101   if (p->getDaug(0)->getNDaug() == 3) {
00102      double md3 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(2)->getId());
00103      minResMass = minResMass + md3;
00104   }
00105   
00106   // calculate the maximum energy of the ISR photon
00107   double pgmax=(s-minResMass*minResMass)/(2.0*w);
00108   double pgz=0.99*pgmax*exp(log(EvtRandom::Flat(1.0))/(beta*power));
00109   if (EvtRandom::Flat(1.0)<0.5) pgz=-pgz;
00110   
00111   double k=fabs(pgz);
00112   // print of ISR energy 
00113   // std::cout << "Energy ISR :"<< k <<std::endl;
00114   EvtVector4R p4g(k,0.0,0.0,pgz);
00115 
00116   EvtVector4R p4res=p->getP4Restframe()-p4g;
00117 
00118   double mres=p4res.mass();
00119 
00120   // set masses
00121   p->getDaug(0)->init(getDaug(0),p4res);
00122   p->getDaug(1)->init(getDaug(1),p4g);
00123 
00124   
00125   // determine XS - langbw
00126   // very crude way of determining XS just a simple straight line Approx.
00127   // this was determined by eye.
00128   // lots of cout statements to make plots to check that things are working as expected
00129   double sigma=9.0;
00130   if (mres<=3.9) sigma = 0.00001;
00131 
00132   bool sigmacomputed(false);
00133 
00134   // DETERMINE XS FOR D*D*
00135   if (p->getDaug(0)->getNDaug() == 2 
00136       &&((p->getDaug(0)->getDaug(0)->getId()==D0S 
00137           && p->getDaug(0)->getDaug(1)->getId()==D0BS)
00138          ||(p->getDaug(0)->getDaug(0)->getId()==DPS 
00139             && p->getDaug(0)->getDaug(1)->getId()==DMS))){
00140      if(mres>4.18) {
00141         sigma*=5./9.*(1.-1.*sqrt((4.18-mres)*(4.18-mres))/(4.3-4.18));
00142      }  
00143      else if(mres>4.07 && mres<=4.18) {
00144      sigma*=5./9.;
00145      }  
00146      else if (mres<=4.07&&mres>4.03)
00147      {
00148         sigma*=(5./9. - 1.5/9.*sqrt((4.07-mres)*(4.07-mres))/(4.07-4.03)); 
00149      }
00150      else if (mres<=4.03&& mres>=4.013)
00151      {
00152         sigma*=(3.5/9. - 3.5/9.*sqrt((4.03-mres)*(4.03-mres))/(4.03-4.013)); 
00153      }
00154      else{     
00155         sigma=0.00001; 
00156      }
00157      sigmacomputed = true;
00158 //     std::cout << "DSDSXS "<<sigma<< " " <<  mres<<std::endl;
00159   }
00160   
00161   // DETERMINE XS FOR D*D
00162   if(p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==D0S 
00163                                          && p->getDaug(0)->getDaug(1)->getId()==D0B)
00164                                         ||(p->getDaug(0)->getDaug(0)->getId()==DPS 
00165                                            && p->getDaug(0)->getDaug(1)->getId()==DM) 
00166                                         ||(p->getDaug(0)->getDaug(0)->getId()==D0BS
00167                                            && p->getDaug(0)->getDaug(1)->getId()==D0)
00168                                         ||(p->getDaug(0)->getDaug(0)->getId()==DMS 
00169                                            && p->getDaug(0)->getDaug(1)->getId()==DP)) )
00170   {
00171      if(mres>=4.2){
00172         sigma*=1.5/9.;
00173      }
00174      else if( mres>4.06 && mres<4.2){
00175         sigma*=((1.5/9.+2.5/9.*sqrt((4.2-mres)*(4.2-mres))/(4.2-4.06)));
00176      }  
00177      else if(mres>=4.015 && mres<4.06){
00178         sigma*=((4./9.+3./9.*sqrt((4.06-mres)*(4.06-mres))/(4.06-4.015)));
00179      }  
00180      else if (mres<4.015 && mres>=3.9){
00181         sigma*=((7./9.-7/9.*sqrt((4.015-mres)*(4.015-mres))/(4.015-3.9))); 
00182      } 
00183      else { 
00184         sigma = 0.00001; 
00185      }
00186      sigmacomputed = true;
00187 //     std::cout << "DSDXS "<<sigma<< " " <<  mres<<std::endl;
00188   }
00189      
00190   // DETERMINE XS FOR Ds*Ds*
00191   if (((p->getDaug(0)->getDaug(0)->getId()==DSPS && p->getDaug(0)->getDaug(1)->getId()==DSMS)))
00192   {
00193      if(mres>(2.112+2.112)){
00194       sigma=0.4; 
00195      }
00196      else  {
00197 //      sigma=0.4; 
00198 //      sigma = 0 surely below Ds*Ds* threshold? - ponyisi
00199         sigma=0.00001;
00200      }
00201      sigmacomputed = true;
00202 //     std::cout << "DsSDsSXS "<<sigma<< " " <<  mres<<std::endl;
00203   }
00204 
00205   // DETERMINE XS FOR Ds*Ds
00206   if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==DSPS 
00207                                           && p->getDaug(0)->getDaug(1)->getId()==DSM)
00208                                          || (p->getDaug(0)->getDaug(0)->getId()==DSMS
00209                                              && p->getDaug(0)->getDaug(1)->getId()==DSP)))
00210   {
00211      if(mres>4.26){
00212         sigma=0.05; 
00213      } 
00214      else if (mres>4.18 && mres<=4.26){
00215         sigma*=1./9.*(0.05+0.95*sqrt((4.26-mres)*(4.26-mres))/(4.26-4.18));
00216      } 
00217      else if (mres>4.16 && mres<=4.18){
00218         sigma*=1/9.; 
00219      } 
00220      else if (mres<=4.16 && mres>4.08){
00221         sigma*=1/9.*(1-sqrt((4.16-mres)*(4.16-mres))/(4.16-4.08)); 
00222      }
00223      else if (mres<=(4.08)){
00224         sigma=0.00001; 
00225      }
00226      sigmacomputed = true;
00227 //     std::cout << "DsSDsXS "<<sigma<< " " <<  mres<<std::endl;
00228   }
00229 
00230   // DETERMINE XS FOR DD
00231   if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==D0 
00232                                           && p->getDaug(0)->getDaug(1)->getId()==D0B)
00233                                          ||(p->getDaug(0)->getDaug(0)->getId()==DP 
00234                                             && p->getDaug(0)->getDaug(1)->getId()==DM))){ 
00235      sigma*=0.4/9.;  
00236      sigmacomputed = true;
00237 //     std::cout << "DDXS "<<sigma<< " " <<  mres<<std::endl;
00238   } 
00239   
00240   // DETERMINE XS FOR DsDs
00241   if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==DSP && p->getDaug(0)->getDaug(1)->getId()==DSM))){
00242      sigma*=0.2/9.;
00243      sigmacomputed = true;
00244 //     std::cout << "DsDsXS "<<sigma<< " " <<  mres<<std::endl;
00245   } 
00246 
00247   // DETERMINE XS FOR MULTIBODY
00248   if (p->getDaug(0)->getNDaug() == 3){
00249      if(mres>4.03){
00250         sigma*=0.5/9.;
00251      }
00252      else {
00253         sigma=0.00001; 
00254      }
00255      sigmacomputed = true;
00256 //     std::cout << "DSDpiXS "<<sigma<< " " <<  mres<<std::endl;
00257   }
00258 
00259   if (! sigmacomputed) {
00260     report(ERROR,"EvtGen") << "VPHOTOVISRHI: This model requires daughters to be listed in a particular order." << endl
00261                            << "The following are acceptable:" << endl
00262                            << "D0 anti-D0" << endl
00263                            << "D+ D-" << endl
00264                            << "D*0 anti-D0" << endl
00265                            << "anti-D*0 D0" << endl
00266                            << "D*+ D-" << endl
00267                            << "D*- D+" << endl
00268                            << "D*0 anti-D*0" << endl
00269                            << "D*+ D*-" << endl
00270                            << "D_s+ D_s-" << endl
00271                            << "D_s*+ D_s-" << endl
00272                            << "D_s*- D_s+" << endl
00273                            << "D_s*+ D_s*-" << endl
00274                            << "(D* D pi can be in any order)" << endl
00275                            << "Aborting..." << endl;
00276     assert(0);
00277   }
00278 
00279   if(sigma<0) sigma = 0.0;
00280 
00281 //   static double sigmax=sigma;
00282 //   if (sigma>sigmax){
00283 //      sigmax=sigma;
00284 //   }
00285   
00286   static int count=0;
00287   
00288   count++;
00289   
00290 //   if (count%10000==0){
00291 //      std::cout << "sigma :"<<sigma<<std::endl;
00292 //      std::cout << "sigmax:"<<sigmax<<std::endl;
00293 //   }
00294   
00295   double norm=sqrt(sigma);
00296   
00297 //  EvtParticle* d=p->getDaug(0);
00298   
00299   
00300   vertex(0,0,0,norm*p->eps(0)*p->epsParent(0).conj());
00301   vertex(1,0,0,norm*p->eps(1)*p->epsParent(0).conj());
00302   vertex(2,0,0,norm*p->eps(2)*p->epsParent(0).conj());
00303   
00304   vertex(0,1,0,norm*p->eps(0)*p->epsParent(1).conj());
00305   vertex(1,1,0,norm*p->eps(1)*p->epsParent(1).conj());
00306   vertex(2,1,0,norm*p->eps(2)*p->epsParent(1).conj());
00307   
00308   vertex(0,2,0,norm*p->eps(0)*p->epsParent(2).conj());
00309   vertex(1,2,0,norm*p->eps(1)*p->epsParent(2).conj());
00310   vertex(2,2,0,norm*p->eps(2)*p->epsParent(2).conj());
00311   
00312   vertex(0,0,1,norm*p->eps(0)*p->epsParent(0).conj());
00313   vertex(1,0,1,norm*p->eps(1)*p->epsParent(0).conj());
00314   vertex(2,0,1,norm*p->eps(2)*p->epsParent(0).conj());
00315   
00316   vertex(0,1,1,norm*p->eps(0)*p->epsParent(1).conj());
00317   vertex(1,1,1,norm*p->eps(1)*p->epsParent(1).conj());
00318   vertex(2,1,1,norm*p->eps(2)*p->epsParent(1).conj());
00319   
00320   vertex(0,2,1,norm*p->eps(0)*p->epsParent(2).conj());
00321   vertex(1,2,1,norm*p->eps(1)*p->epsParent(2).conj());
00322   vertex(2,2,1,norm*p->eps(2)*p->epsParent(2).conj());
00323   
00324   return;
00325 }