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