TEvtGen/EvtGenModels/EvtPHOTOS.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) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtPHOTOS.cc
  12 //
  13 // Description: This routine takes the particle *p and applies
  14 //              the PHOTOS package to generate final state radiation
  15 //              on the produced mesons.
  16 //
  17 // Modification history:
  18 //
  19 //    RYD     October 1, 1997        Module created
  20 //
  21 //------------------------------------------------------------------------
  22 //
  23 #include "EvtGenBase/EvtPatches.hh"
  24 #include "EvtGenBase/EvtPatches.hh"
  25 #include "EvtGenBase/EvtIdSet.hh"
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtPhotonParticle.hh"
  28 #include "EvtGenBase/EvtPDL.hh"
  29 #include "EvtGenModels/EvtPHOTOS.hh"
  30 #include "EvtGenBase/EvtReport.hh"
  31 #include <stdlib.h>
  32
  33 extern "C" void begevtgenstorex_(int *,int *,int *,int *,
  34                                 int *,int *,int *,int *,
  35                                 double *,double *,double *,
  36                                 double *,double *,double *,
  37                                 double *,double *,double *);
  38
  39 extern "C" void begevtgengetx_(int *,int *,int *,int *,
  40                               int *,int *,int *,int *,
  41                               double *,double *,double *,
  42                               double *,double *,double *,
  43                               double *,double *,double *);
  44
  45 extern "C" void heplst_(int *);
  46
  47 extern "C" void photos_(int *);
  48
  49 extern "C" void phoini_();
  50
  51
  52 EvtPHOTOS::EvtPHOTOS(std::string photontype){
  53
  54     _photontype=photontype;
  55
  56 }
  57
  58 void EvtPHOTOS::doRadCorr( EvtParticle *p){
  59
  60   static int first=1;
  61
  62   //added by Lange Jan4,2000
  63   //allow to set photon tupe
  64   static EvtId GAMM=EvtPDL::getId(_photontype);
  65
  66   if (GAMM==EvtId(-1,-1)) {
  67     report(ERROR,"EvtGen") << "In EvtPHOTOS::doRadCorr():Particle:"<<
  68       _photontype<<" is not in EvtPDL"<<std::endl;
  69      ::abort();
  70   }
  71
  72   if (first) {
  73     first=0;
  74     phoini_();
  75   }
  76
  77   double mpho=EvtPDL::getMeanMass(GAMM);
  78
  79   int entry,eventnum,numparticle,istat,partnum,mother;
  80   int daugfirst,dauglast;
  81
  82   int numparticlephotos;
  83
  84   double px,py,pz,e,m,x,y,z,t;
  85
  86   static EvtId dq=EvtPDL::getId("d");
  87   static EvtId adq=EvtPDL::getId("anti-d");
  88   static EvtId uq=EvtPDL::getId("u");
  89   static EvtId auq=EvtPDL::getId("anti-u");
  90   static EvtId sq=EvtPDL::getId("s");
  91   static EvtId asq=EvtPDL::getId("anti-s");
  92   static EvtId cq=EvtPDL::getId("c");
  93   static EvtId acq=EvtPDL::getId("anti-c");
  94   static EvtId bq=EvtPDL::getId("b");
  95   static EvtId abq=EvtPDL::getId("anti-b");
  96   static EvtId tq=EvtPDL::getId("t");
  97   static EvtId atq=EvtPDL::getId("anti-t");
  98   static EvtId vpho=EvtPDL::getId("vpho");
  99
 100   static EvtIdSet quarks(dq,adq,uq,auq,sq,asq,cq,acq,bq,abq,tq,atq);
 101
 102   if ( p->getId() == vpho ) return;
 103   if ( p->getNDaug() > 10 ) return;
 104
 105   px=0.0;
 106   py=0.0;
 107   pz=0.0;
 108   e=p->mass();
 109   m=p->mass();
 110   x=0.0;
 111   y=0.0;
 112   z=0.0;
 113   t=0.0;
 114
 115   entry=1;
 116   eventnum=1;
 117   numparticle=1;
 118   istat=2;
 119   partnum=EvtPDL::getStdHep(p->getId());
 120   mother=0;
 121   daugfirst=2;
 122   dauglast=1+p->getNDaug();
 123
 124   begevtgenstorex_(&entry,&eventnum,&numparticle,&istat,&partnum,
 125                   &mother,&daugfirst,&dauglast,
 126                   &px,&py,&pz,&e,&m,&x,&y,&z,&t);
 127
 128   //  std::cout << EvtPDL::name(p->getId()) <<" " ;
 129   for(size_t i=0;i<p->getNDaug();i++){
 130
 131     //No quarks to photos
 132     if (quarks.contains(p->getDaug(i)->getId())==1) continue;
 133
 134     px=p->getDaug(i)->getP4().get(1);
 135     py=p->getDaug(i)->getP4().get(2);
 136     pz=p->getDaug(i)->getP4().get(3);
 137     e=p->getDaug(i)->getP4().get(0);
 138     m=p->getDaug(i)->mass();
 139     x=0.0;
 140     y=0.0;
 141     z=0.0;
 142     t=0.0;
 143
 144     //    std::cout << EvtPDL::name(p->getDaug(i)->getId()) << " " ;
 145     entry+=1;
 146     eventnum=1;
 147     numparticle+=1;
 148     istat=1;
 149     partnum=EvtPDL::getStdHep(p->getDaug(i)->getId());
 150     mother=1;
 151     daugfirst=0;
 152     dauglast=0;
 153
 154     begevtgenstorex_(&entry,&eventnum,&numparticle,&istat,&partnum,
 155                     &mother,&daugfirst,&dauglast,
 156                     &px,&py,&pz,&e,&m,&x,&y,&z,&t);
 157
 158
 159   }
 160   //  std::cout << std::endl;
 161   //can't use heplst since the common block used by the BaBar
 162   //implementation of PHOTOS  is renamed due to real*4 vs real*8
 163   //problems.
 164
 165   //int mlst=1;
 166
 167   //heplst_(&mlst);
 168
 169   entry=1;
 170
 171   //  report(INFO,"EvtGen") << "Doing photos " << EvtPDL::name(p->getId()) << endl;
 172   photos_(&entry);
 173   //  report(INFO,"EvtGen") << "done\n";
 174   begevtgengetx_(&entry,&eventnum,&numparticlephotos,&istat,&partnum,
 175                     &mother,&daugfirst,&dauglast,
 176                     &px,&py,&pz,&e,&m,&x,&y,&z,&t);
 177
 178
 179   //report(INFO,"EvtGen") << "numparticlephotos:"<<numparticlephotos<<endl;
 180
 181   if (numparticle==numparticlephotos) return;
 182
 183   EvtVector4R new4mom;
 184
 185   int np;
 186
 187   for(size_t i=0;i<p->getNDaug();i++){
 188
 189     entry=i+2;
 190
 191     begevtgengetx_(&entry,&eventnum,&np,&istat,&partnum,
 192                     &mother,&daugfirst,&dauglast,
 193                     &px,&py,&pz,&e,&m,&x,&y,&z,&t);
 194
 195     //this is needed to ensure that photos does not
 196     //change the masses. But it will violate energy conservation!
 197     double mp=p->getDaug(i)->mass();
 198     e=sqrt(mp*mp+px*px+py*py+pz*pz);
 199
 200     new4mom.set(e,px,py,pz);
 201
 202     p->getDaug(i)->setP4WithFSR(new4mom);
 203
 204
 205
 206   }
 207
 208   for(entry=numparticle+1;entry<=numparticlephotos;entry++){
 209
 210     begevtgengetx_(&entry,&eventnum,&np,&istat,&partnum,
 211                     &mother,&daugfirst,&dauglast,
 212                     &px,&py,&pz,&e,&m,&x,&y,&z,&t);
 213
 214     //Hack here to give the photon the mass of the generated particle
 215     e=sqrt(mpho*mpho+px*px+py*py+pz*pz);
 216     new4mom.set(e,px,py,pz);
 217
 218     //new4mom.dump();
 219
 220     EvtPhotonParticle* gamma;
 221     gamma=new EvtPhotonParticle;
 222     gamma->init(GAMM,new4mom);
 223     gamma->setFSRP4toZero();
 224     //    report(INFO,"EvtGen") << gamma << " " << p << " "<< px << " " << py << " " << pz << " " << p->getNDaug() << " " << EvtPDL::name(p->getId())<<" " << entry << " " <<numparticlephotos<<endl;
 225     gamma->addDaug(p);
 226
 227 //    p->getDaug(i)->set_type(EvtSpinType::PHOTON);
 228
 229   }
 230   return ;
 231 }
 232