TEvtGen/EvtGen/EvtGen.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: EvtGen.cc
  12 //
  13 // Description: Main class to provide user interface to EvtGen.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     March 24, 1998        Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdio.h>
  23 #include <fstream>
  24 #include <math.h>
  25 #include "EvtGenBase/EvtComplex.hh"
  26 #include <stdlib.h>
  27 #include "EvtGen/EvtGen.hh"
  28 #include "EvtGenBase/EvtVector4R.hh"
  29 #include "EvtGenBase/EvtVectorParticle.hh"
  30 #include "EvtGenBase/EvtParticle.hh"
  31 #include "EvtGenBase/EvtScalarParticle.hh"
  32 #include "EvtGenBase/EvtDecayTable.hh"
  33 #include "EvtGenBase/EvtPDL.hh"
  34 #include "EvtGenBase/EvtStdHep.hh"
  35 #include "EvtGenBase/EvtSecondary.hh"
  36 #include "EvtGenBase/EvtReport.hh"
  37 #include "EvtGenBase/EvtId.hh"
  38 #include "EvtGenBase/EvtRandom.hh"
  39 #include "EvtGenBase/EvtRandomEngine.hh"
  40 #include "EvtGenBase/EvtSimpleRandomEngine.hh"
  41 #include "EvtGenBase/EvtParticleFactory.hh"
  42 //#include "CLHEP/Vector/LorentzVector.h"
  43 #include "TLorentzVector.h"
  44 #include "EvtGenModels/EvtModelReg.hh"
  45 #include "EvtGenBase/EvtStatus.hh"
  46 #include "EvtGenBase/EvtAbsRadCorr.hh"
  47 #include "EvtGenBase/EvtRadCorr.hh"
  48 #include "EvtGenModels/EvtPHOTOS.hh"
  49 using std::endl;
  50 using std::fstream;
  51 using std::ifstream;
  52
  53
  54 extern "C" void begevtgenstore_(int *,int *,int *,int *,
  55                                 int *,int *,int *,int *,int *,
  56                                 double *,double *,double *,
  57                                 double *,double *,double *,
  58                                 double *,double *,double *);
  59
  60 extern "C" {
  61 extern void evtgen_(float svertex[3],float *e_cms,float *beta_zs,
  62                     int *mode);
  63 }
  64
  65 EvtGen::~EvtGen(){
  66
  67   //This is a bit uggly, should not do anything
  68   //in a destructor. This will fail if EvtGen is made a static
  69   //because then this destructor might be called _after_
  70   //the destructoin of objects that it depends on, e.g., EvtPDL.
  71
  72   if (getenv("EVTINFO")){
  73     EvtDecayTable::printSummary();
  74   }
  75
  76 }
  77
  78 EvtGen::EvtGen(const char* const decayName,
  79                const char* const pdtTableName,
  80                EvtRandomEngine* randomEngine,
  81                EvtAbsRadCorr* isrEngine,
  82                const std::list<EvtDecayBase*>* extraModels){
  83
  84
  85   report(INFO,"EvtGen") << "Initializing EvtGen"<<endl;
  86
  87   report(INFO,"EvtGen") << "Storing known decay models"<<endl;
  88   EvtModelReg dummy(extraModels);
  89
  90   report(INFO,"EvtGen") << "Main decay file name  :"<<decayName<<endl;
  91   report(INFO,"EvtGen") << "PDT table file name   :"<<pdtTableName<<endl;
  92
  93   report(INFO,"EvtGen") << "Initializing RadCorr=PHOTOS"<<endl;
  94   if (isrEngine==0){
  95     static EvtPHOTOS defaultRadCorrEngine;
  96     EvtRadCorr::setRadCorrEngine(&defaultRadCorrEngine);
  97   }
  98   else{
  99     EvtRadCorr::setRadCorrEngine(isrEngine);
 100   }
 101
 102   _pdl.readPDT(pdtTableName);
 103
 104   EvtDecayTable::readDecayFile(decayName,false);
 105
 106   if (randomEngine==0){
 107     static EvtSimpleRandomEngine defaultRandomEngine;
 108     EvtRandom::setRandomEngine((EvtRandomEngine*)&defaultRandomEngine);
 109     report(INFO,"EvtGen") <<"No random engine given in "
 110                           <<"EvtGen::EvtGen constructor, "
 111                           <<"will use default EvtSimpleRandomEngine."<<endl;
 112   }
 113   else{
 114     EvtRandom::setRandomEngine(randomEngine);
 115   }
 116
 117   report(INFO,"EvtGen") << "Done initializing EvtGen"<<endl;
 118
 119 }
 120
 121
 122 void EvtGen::readUDecay(const char* const uDecayName){
 123
 124   ifstream indec;
 125
 126   if ( uDecayName[0] == 0) {
 127     report(INFO,"EvtGen") << "Is not reading a user decay file!"<<endl;
 128   }
 129   else{
 130     indec.open(uDecayName);
 131     if (indec) {
 132       EvtDecayTable::readDecayFile(uDecayName,true);
 133     }
 134     else{
 135
 136       report(INFO,"EvtGen") << "Can not find UDECAY file '"
 137                             <<uDecayName<<"'.  Stopping"<<endl;
 138       ::abort();
 139     }
 140   }
 141
 142 }
 143
 144
 145 void EvtGen::generateDecay(int stdhepid,
 146                            EvtVector4R P,
 147                            EvtVector4R D,
 148                            EvtStdHep *evtStdHep,
 149                            EvtSpinDensity *spinDensity ){
 150
 151
 152   EvtParticle *p;
 153
 154   if(spinDensity==0){
 155     p=EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(stdhepid),P);
 156   }
 157   else{
 158     p=EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(stdhepid),
 159                                           P,*spinDensity);
 160   }
 161
 162   generateDecay(p);
 163   //  p->Decay();
 164
 165   evtStdHep->init();
 166
 167   p->makeStdHep(*evtStdHep);
 168
 169   evtStdHep->translate(D);
 170
 171   p->deleteTree();
 172
 173
 174 }
 175
 176 void EvtGen::generateDecay(EvtParticle *p){
 177
 178   int times=0;
 179   do{
 180     times+=1;
 181     EvtStatus::initRejectFlag();
 182
 183     p->decay();
 184     //ok then finish.
 185     if ( EvtStatus::getRejectFlag()==0 ) {
 186       times=0;
 187     }
 188     else{
 189       for (size_t ii=0;ii<p->getNDaug();ii++){
 190         EvtParticle *temp=p->getDaug(ii);
 191         temp->deleteTree();
 192       }
 193       p->resetFirstOrNot();
 194       p->resetNDaug();
 195
 196     }
 197
 198     if ( times==10000) {
 199       report(ERROR,"EvtGen") << "Your event has been rejected 10000 times!"<<endl;
 200       report(ERROR,"EvtGen") << "Will now abort."<<endl;
 201       ::abort();
 202       times=0;
 203     }
 204   } while (times);
 205
 206 }
 207
 208
 209
 210 //void EvtGen::generateEvent(int stdhepid,CLHEP::HepLorentzVector P,CLHEP::HepLorentzVector D){
 211 void EvtGen::generateEvent(int stdhepid,TLorentzVector P,TLorentzVector D){
 212   EvtParticle *root_part;
 213   EvtVectorParticle *vector_part;
 214
 215   vector_part=new EvtVectorParticle;
 216   EvtVector4R p_init;
 217
 218   //p_init.set(P.t(),P.x(),P.y(),P.z());
 219   p_init.set(P.E(),P.Px(),P.Py(),P.Pz());
 220   vector_part->init(EvtPDL::evtIdFromStdHep(stdhepid),p_init);
 221
 222   root_part=(EvtParticle *)vector_part;
 223
 224   root_part->setVectorSpinDensity();
 225
 226   generateEvent(root_part,D);
 227
 228   root_part->deleteTree();
 229
 230 }
 231
 232 //void EvtGen::generateEvent(EvtParticle *root_part,CLHEP::HepLorentzVector D){
 233 void EvtGen::generateEvent(EvtParticle *root_part,TLorentzVector D){
 234   int i;
 235
 236   static int nevent=0;
 237
 238   nevent++;
 239
 240   static EvtStdHep evtstdhep;
 241   //  static EvtSecondary evtsecondary;
 242
 243   int j;
 244   int istat;
 245   int partnum;
 246   double px,py,pz,e,m;
 247   double x,y,z,t;
 248
 249   EvtVector4R p4,x4;
 250
 251   generateDecay(root_part);
 252   //  root_part->Decay();
 253
 254   int          npart=0;
 255
 256   EvtId        list_of_stable[10];
 257   EvtParticle* stable_parent[10];
 258
 259
 260   list_of_stable[0]=EvtId(-1,-1);
 261   stable_parent[0]=0;
 262
 263   evtstdhep.init();
 264   //  evtsecondary.init();
 265   // root_part->makeStdHep(evtstdhep,evtsecondary,list_of_stable);
 266   root_part->makeStdHep(evtstdhep);
 267
 268   //report(INFO,"EvtGen") << evtstdhep;
 269   //report(INFO,"EvtGen") << evtsecondary;
 270
 271   npart=evtstdhep.getNPart();
 272
 273   for(i=0;i<evtstdhep.getNPart();i++){
 274
 275     j=i+1;
 276
 277     int jmotherfirst=evtstdhep.getFirstMother(i)+1;
 278     int jmotherlast=evtstdhep.getLastMother(i)+1;
 279     int jdaugfirst=evtstdhep.getFirstDaughter(i)+1;
 280     int jdauglast=evtstdhep.getLastDaughter(i)+1;
 281
 282     partnum=evtstdhep.getStdHepID(i);
 283
 284     istat=evtstdhep.getIStat(i);
 285
 286     p4=evtstdhep.getP4(i);
 287     x4=evtstdhep.getX4(i);
 288
 289     px=p4.get(1);
 290     py=p4.get(2);
 291     pz=p4.get(3);
 292     e=p4.get(0);
 293
 294     x=x4.get(1)+D.X();
 295     y=x4.get(2)+D.Y();
 296     z=x4.get(3)+D.Z();
 297     t=x4.get(0)+D.T();
 298
 299     m=p4.mass();
 300
 301     begevtgenstore_(&j,&nevent,&npart,&istat,
 302                     &partnum,&jmotherfirst,&jmotherlast,
 303                     &jdaugfirst,&jdauglast,
 304                     &px,&py,&pz,&e,
 305                     &m,&x,&y,&z,&t);
 306   }
 307
 308 }
 309
 310
 311
 312
 313
 314
 315
 316
 317