TEvtGen/EvtGenModels/EvtVPHOtoVISRHi.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) 2004      Cornell
  10 //
  11 // Module: EvtVPHOtoVISR.cc
  12 //
  13 // Description: Routine to decay vpho -> (DDx) + ISR photon from 3.9 to 4.3 GeV, using CLEO-c data (Brian Lang)
  14 //
  15 // Modification history:
  16 //
  17 //    Ryd       March 20, 2004       Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include <stdlib.h>
  22 #include "EvtGenBase/EvtParticle.hh"
  23 #include "EvtGenBase/EvtGenKine.hh"
  24 #include "EvtGenBase/EvtPDL.hh"
  25 #include "EvtGenBase/EvtVector4C.hh"
  26 #include "EvtGenBase/EvtVector4R.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenBase/EvtRandom.hh"
  29 #include "EvtGenModels/EvtVPHOtoVISRHi.hh"
  30 #include <string>
  31
  32 using std::endl;
  33
  34 EvtVPHOtoVISRHi::~EvtVPHOtoVISRHi() {}
  35
  36 std::string EvtVPHOtoVISRHi::getName(){
  37
  38   return "VPHOTOVISRHI";
  39
  40 }
  41
  42
  43 EvtDecayBase* EvtVPHOtoVISRHi::clone(){
  44
  45   return new EvtVPHOtoVISRHi;
  46
  47 }
  48
  49 void EvtVPHOtoVISRHi::init(){
  50
  51   // check that there are 0 or 1 arguments
  52   checkNArg(0,1);
  53
  54   // check that there are 2 daughters
  55   checkNDaug(2);
  56
  57   // check the parent and daughter spins
  58   checkSpinParent(EvtSpinType::VECTOR);
  59   checkSpinDaughter(0,EvtSpinType::VECTOR);
  60   checkSpinDaughter(1,EvtSpinType::PHOTON);
  61 }
  62
  63 void EvtVPHOtoVISRHi::initProbMax() {
  64
  65    setProbMax(20.0);
  66
  67 }
  68
  69 void EvtVPHOtoVISRHi::decay( EvtParticle *p){
  70   //take photon along z-axis, either forward or backward.
  71   //Implement this as generating the photon momentum along
  72   //the z-axis uniformly
  73    double power=1;
  74    if (getNArg()==1) power=getArg(0);
  75    // define particle names
  76   static EvtId D0=EvtPDL::getId("D0");
  77   static EvtId D0B=EvtPDL::getId("anti-D0");
  78   static EvtId DP=EvtPDL::getId("D+");
  79   static EvtId DM=EvtPDL::getId("D-");
  80   static EvtId DSM=EvtPDL::getId("D_s-");
  81   static EvtId DSP=EvtPDL::getId("D_s+");
  82   static EvtId DSMS=EvtPDL::getId("D_s*-");
  83   static EvtId DSPS=EvtPDL::getId("D_s*+");
  84   static EvtId D0S=EvtPDL::getId("D*0");
  85   static EvtId D0BS=EvtPDL::getId("anti-D*0");
  86   static EvtId DPS=EvtPDL::getId("D*+");
  87   static EvtId DMS=EvtPDL::getId("D*-");
  88   // setup some parameters
  89   double w=p->mass();
  90   double s=w*w;
  91   double L=2.0*log(w/0.000511);
  92   double alpha=1/137.0;
  93   double beta=(L-1)*2.0*alpha/EvtConst::pi;
  94   // make sure only 2 or 3 body are present
  95   assert (p->getDaug(0)->getNDaug() == 2 || p->getDaug(0)->getNDaug() == 3);
  96
  97   // determine minimum rest mass of parent
  98   double md1 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(0)->getId());
  99   double md2 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(1)->getId());
 100   double minResMass = md1+md2;
 101   if (p->getDaug(0)->getNDaug() == 3) {
 102      double md3 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(2)->getId());
 103      minResMass = minResMass + md3;
 104   }
 105
 106   // calculate the maximum energy of the ISR photon
 107   double pgmax=(s-minResMass*minResMass)/(2.0*w);
 108   double pgz=0.99*pgmax*exp(log(EvtRandom::Flat(1.0))/(beta*power));
 109   if (EvtRandom::Flat(1.0)<0.5) pgz=-pgz;
 110
 111   double k=fabs(pgz);
 112   // print of ISR energy
 113   // std::cout << "Energy ISR :"<< k <<std::endl;
 114   EvtVector4R p4g(k,0.0,0.0,pgz);
 115
 116   EvtVector4R p4res=p->getP4Restframe()-p4g;
 117
 118   double mres=p4res.mass();
 119
 120   // set masses
 121   p->getDaug(0)->init(getDaug(0),p4res);
 122   p->getDaug(1)->init(getDaug(1),p4g);
 123
 124
 125   // determine XS - langbw
 126   // very crude way of determining XS just a simple straight line Approx.
 127   // this was determined by eye.
 128   // lots of cout statements to make plots to check that things are working as expected
 129   double sigma=9.0;
 130   if (mres<=3.9) sigma = 0.00001;
 131
 132   bool sigmacomputed(false);
 133
 134   // DETERMINE XS FOR D*D*
 135   if (p->getDaug(0)->getNDaug() == 2
 136       &&((p->getDaug(0)->getDaug(0)->getId()==D0S
 137           && p->getDaug(0)->getDaug(1)->getId()==D0BS)
 138          ||(p->getDaug(0)->getDaug(0)->getId()==DPS
 139             && p->getDaug(0)->getDaug(1)->getId()==DMS))){
 140      if(mres>4.18) {
 141         sigma*=5./9.*(1.-1.*sqrt((4.18-mres)*(4.18-mres))/(4.3-4.18));
 142      }
 143      else if(mres>4.07 && mres<=4.18) {
 144      sigma*=5./9.;
 145      }
 146      else if (mres<=4.07&&mres>4.03)
 147      {
 148         sigma*=(5./9. - 1.5/9.*sqrt((4.07-mres)*(4.07-mres))/(4.07-4.03));
 149      }
 150      else if (mres<=4.03&& mres>=4.013)
 151      {
 152         sigma*=(3.5/9. - 3.5/9.*sqrt((4.03-mres)*(4.03-mres))/(4.03-4.013));
 153      }
 154      else{
 155         sigma=0.00001;
 156      }
 157      sigmacomputed = true;
 158 //     std::cout << "DSDSXS "<<sigma<< " " <<  mres<<std::endl;
 159   }
 160
 161   // DETERMINE XS FOR D*D
 162   if(p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==D0S
 163                                          && p->getDaug(0)->getDaug(1)->getId()==D0B)
 164                                         ||(p->getDaug(0)->getDaug(0)->getId()==DPS
 165                                            && p->getDaug(0)->getDaug(1)->getId()==DM)
 166                                         ||(p->getDaug(0)->getDaug(0)->getId()==D0BS
 167                                            && p->getDaug(0)->getDaug(1)->getId()==D0)
 168                                         ||(p->getDaug(0)->getDaug(0)->getId()==DMS
 169                                            && p->getDaug(0)->getDaug(1)->getId()==DP)) )
 170   {
 171      if(mres>=4.2){
 172         sigma*=1.5/9.;
 173      }
 174      else if( mres>4.06 && mres<4.2){
 175         sigma*=((1.5/9.+2.5/9.*sqrt((4.2-mres)*(4.2-mres))/(4.2-4.06)));
 176      }
 177      else if(mres>=4.015 && mres<4.06){
 178         sigma*=((4./9.+3./9.*sqrt((4.06-mres)*(4.06-mres))/(4.06-4.015)));
 179      }
 180      else if (mres<4.015 && mres>=3.9){
 181         sigma*=((7./9.-7/9.*sqrt((4.015-mres)*(4.015-mres))/(4.015-3.9)));
 182      }
 183      else {
 184         sigma = 0.00001;
 185      }
 186      sigmacomputed = true;
 187 //     std::cout << "DSDXS "<<sigma<< " " <<  mres<<std::endl;
 188   }
 189
 190   // DETERMINE XS FOR Ds*Ds*
 191   if (((p->getDaug(0)->getDaug(0)->getId()==DSPS && p->getDaug(0)->getDaug(1)->getId()==DSMS)))
 192   {
 193      if(mres>(2.112+2.112)){
 194       sigma=0.4;
 195      }
 196      else  {
 197 //      sigma=0.4;
 198 //      sigma = 0 surely below Ds*Ds* threshold? - ponyisi
 199         sigma=0.00001;
 200      }
 201      sigmacomputed = true;
 202 //     std::cout << "DsSDsSXS "<<sigma<< " " <<  mres<<std::endl;
 203   }
 204
 205   // DETERMINE XS FOR Ds*Ds
 206   if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==DSPS
 207                                           && p->getDaug(0)->getDaug(1)->getId()==DSM)
 208                                          || (p->getDaug(0)->getDaug(0)->getId()==DSMS
 209                                              && p->getDaug(0)->getDaug(1)->getId()==DSP)))
 210   {
 211      if(mres>4.26){
 212         sigma=0.05;
 213      }
 214      else if (mres>4.18 && mres<=4.26){
 215         sigma*=1./9.*(0.05+0.95*sqrt((4.26-mres)*(4.26-mres))/(4.26-4.18));
 216      }
 217      else if (mres>4.16 && mres<=4.18){
 218         sigma*=1/9.;
 219      }
 220      else if (mres<=4.16 && mres>4.08){
 221         sigma*=1/9.*(1-sqrt((4.16-mres)*(4.16-mres))/(4.16-4.08));
 222      }
 223      else if (mres<=(4.08)){
 224         sigma=0.00001;
 225      }
 226      sigmacomputed = true;
 227 //     std::cout << "DsSDsXS "<<sigma<< " " <<  mres<<std::endl;
 228   }
 229
 230   // DETERMINE XS FOR DD
 231   if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==D0
 232                                           && p->getDaug(0)->getDaug(1)->getId()==D0B)
 233                                          ||(p->getDaug(0)->getDaug(0)->getId()==DP
 234                                             && p->getDaug(0)->getDaug(1)->getId()==DM))){
 235      sigma*=0.4/9.;
 236      sigmacomputed = true;
 237 //     std::cout << "DDXS "<<sigma<< " " <<  mres<<std::endl;
 238   }
 239
 240   // DETERMINE XS FOR DsDs
 241   if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==DSP && p->getDaug(0)->getDaug(1)->getId()==DSM))){
 242      sigma*=0.2/9.;
 243      sigmacomputed = true;
 244 //     std::cout << "DsDsXS "<<sigma<< " " <<  mres<<std::endl;
 245   }
 246
 247   // DETERMINE XS FOR MULTIBODY
 248   if (p->getDaug(0)->getNDaug() == 3){
 249      if(mres>4.03){
 250         sigma*=0.5/9.;
 251      }
 252      else {
 253         sigma=0.00001;
 254      }
 255      sigmacomputed = true;
 256 //     std::cout << "DSDpiXS "<<sigma<< " " <<  mres<<std::endl;
 257   }
 258
 259   if (! sigmacomputed) {
 260     report(ERROR,"EvtGen") << "VPHOTOVISRHI: This model requires daughters to be listed in a particular order." << endl
 261                            << "The following are acceptable:" << endl
 262                            << "D0 anti-D0" << endl
 263                            << "D+ D-" << endl
 264                            << "D*0 anti-D0" << endl
 265                            << "anti-D*0 D0" << endl
 266                            << "D*+ D-" << endl
 267                            << "D*- D+" << endl
 268                            << "D*0 anti-D*0" << endl
 269                            << "D*+ D*-" << endl
 270                            << "D_s+ D_s-" << endl
 271                            << "D_s*+ D_s-" << endl
 272                            << "D_s*- D_s+" << endl
 273                            << "D_s*+ D_s*-" << endl
 274                            << "(D* D pi can be in any order)" << endl
 275                            << "Aborting..." << endl;
 276     assert(0);
 277   }
 278
 279   if(sigma<0) sigma = 0.0;
 280
 281 //   static double sigmax=sigma;
 282 //   if (sigma>sigmax){
 283 //      sigmax=sigma;
 284 //   }
 285
 286   static int count=0;
 287
 288   count++;
 289
 290 //   if (count%10000==0){
 291 //      std::cout << "sigma :"<<sigma<<std::endl;
 292 //      std::cout << "sigmax:"<<sigmax<<std::endl;
 293 //   }
 294
 295   double norm=sqrt(sigma);
 296
 297 //  EvtParticle* d=p->getDaug(0);
 298
 299
 300   vertex(0,0,0,norm*p->eps(0)*p->epsParent(0).conj());
 301   vertex(1,0,0,norm*p->eps(1)*p->epsParent(0).conj());
 302   vertex(2,0,0,norm*p->eps(2)*p->epsParent(0).conj());
 303
 304   vertex(0,1,0,norm*p->eps(0)*p->epsParent(1).conj());
 305   vertex(1,1,0,norm*p->eps(1)*p->epsParent(1).conj());
 306   vertex(2,1,0,norm*p->eps(2)*p->epsParent(1).conj());
 307
 308   vertex(0,2,0,norm*p->eps(0)*p->epsParent(2).conj());
 309   vertex(1,2,0,norm*p->eps(1)*p->epsParent(2).conj());
 310   vertex(2,2,0,norm*p->eps(2)*p->epsParent(2).conj());
 311
 312   vertex(0,0,1,norm*p->eps(0)*p->epsParent(0).conj());
 313   vertex(1,0,1,norm*p->eps(1)*p->epsParent(0).conj());
 314   vertex(2,0,1,norm*p->eps(2)*p->epsParent(0).conj());
 315
 316   vertex(0,1,1,norm*p->eps(0)*p->epsParent(1).conj());
 317   vertex(1,1,1,norm*p->eps(1)*p->epsParent(1).conj());
 318   vertex(2,1,1,norm*p->eps(2)*p->epsParent(1).conj());
 319
 320   vertex(0,2,1,norm*p->eps(0)*p->epsParent(2).conj());
 321   vertex(1,2,1,norm*p->eps(1)*p->epsParent(2).conj());
 322   vertex(2,2,1,norm*p->eps(2)*p->epsParent(2).conj());
 323
 324   return;
 325 }