00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include <stdio.h>
00024 #include <fstream>
00025 #include <math.h>
00026 #include "EvtGenBase/EvtComplex.hh"
00027 #include <stdlib.h>
00028 #include "EvtGen.hh"
00029 #include "EvtGenBase/EvtVector4R.hh"
00030 #include "EvtGenBase/EvtVectorParticle.hh"
00031 #include "EvtGenBase/EvtParticle.hh"
00032 #include "EvtGenBase/EvtScalarParticle.hh"
00033 #include "EvtGenBase/EvtDecayTable.hh"
00034 #include "EvtGenBase/EvtPDL.hh"
00035 #include "EvtGenBase/EvtStdHep.hh"
00036 #include "EvtGenBase/EvtSecondary.hh"
00037 #include "EvtGenBase/EvtReport.hh"
00038 #include "EvtGenBase/EvtId.hh"
00039 #include "EvtGenBase/EvtRandom.hh"
00040 #include "EvtGenBase/EvtRandomEngine.hh"
00041 #include "EvtGenBase/EvtParticleFactory.hh"
00042 #include "CLHEP/Vector/LorentzVector.h"
00043 #include "EvtGenModels/EvtModelReg.hh"
00044 #include "EvtGenBase/EvtStatus.hh"
00045 #include "EvtGenBase/EvtAbsRadCorr.hh"
00046 #include "EvtGenBase/EvtRadCorr.hh"
00047 #include "EvtGenModels/EvtPHOTOS.hh"
00048
00049 using std::endl;
00050 using std::fstream;
00051 using std::ifstream;
00052
00053 extern "C" void begevtgenstore_(int *,int *,int *,int *,
00054 int *,int *,int *,int *,int *,
00055 double *,double *,double *,
00056 double *,double *,double *,
00057 double *,double *,double *);
00058
00059 extern "C" {
00060 extern void evtgen_(float svertex[3],float *e_cms,float *beta_zs,
00061 int *mode);
00062 }
00063
00064 EvtGen::~EvtGen(){
00065
00066
00067
00068
00069
00070
00071 if (getenv("EVTINFO")){
00072 EvtDecayTable::printSummary();
00073 }
00074
00075 }
00076
00077 EvtGen::EvtGen(const char* const decayName,
00078 const char* const pdtTableName,
00079 EvtRandomEngine* randomEngine,
00080 EvtAbsRadCorr* isrEngine){
00081
00082
00083 report(INFO,"EvtGen") << "Initializing EvtGen"<<endl;
00084
00085 report(INFO,"EvtGen") << "Storing known decay models"<<endl;
00086 EvtModelReg dummy;
00087
00088 report(INFO,"EvtGen") << "Main decay file name :"<<decayName<<endl;
00089 report(INFO,"EvtGen") << "PDT table file name :"<<pdtTableName<<endl;
00090
00091 report(INFO,"EvtGen") << "Initializing RadCorr=PHOTOS"<<endl;
00092 if (isrEngine==0){
00093 static EvtPHOTOS defaultRadCorrEngine;
00094 EvtRadCorr::setRadCorrEngine(&defaultRadCorrEngine);
00095 report(INFO,"EvtGen") <<"No RadCorr engine given in "
00096 <<"EvtGen::EvtGen constructor, "
00097 <<"will use default EvtPHOTOS."<<endl;
00098 }
00099 else{
00100 EvtRadCorr::setRadCorrEngine(isrEngine);
00101 }
00102
00103 _pdl.readPDT(pdtTableName);
00104
00105
00106 ifstream indec;
00107
00108 EvtDecayTable::readDecayFile(decayName);
00109
00110 if (randomEngine==0){
00111 static EvtRandomEngine defaultRandomEngine;
00112 EvtRandom::setRandomEngine(&defaultRandomEngine);
00113 report(INFO,"EvtGen") <<"No random engine given in "
00114 <<"EvtGen::EvtGen constructor, "
00115 <<"will use default EvtRandomEngine."<<endl;
00116 }
00117 else{
00118 EvtRandom::setRandomEngine(randomEngine);
00119 }
00120
00121 report(INFO,"EvtGen") << "Done initializing EvtGen"<<endl;
00122
00123 }
00124
00125
00126 void EvtGen::readUDecay(const char* const uDecayName){
00127
00128 ifstream indec;
00129
00130 if ( uDecayName[0] == 0) {
00131 report(INFO,"EvtGen") << "Is not reading a user decay file!"<<endl;
00132 }
00133 else{
00134 indec.open(uDecayName);
00135 if (indec) {
00136 EvtDecayTable::readDecayFile(uDecayName);
00137
00138 report(INFO,"EvtGen") << "Reading "<<uDecayName
00139 <<" to override decay table."<<endl;
00140 }
00141 else{
00142
00143 report(INFO,"EvtGen") << "Can not find UDECAY file '"
00144 <<uDecayName<<"'. Stopping"<<endl;
00145 ::abort();
00146 }
00147 }
00148
00149 }
00150
00151
00152 void EvtGen::generateDecay(int stdhepid,
00153 EvtVector4R P,
00154 EvtVector4R D,
00155 EvtStdHep *evtStdHep,
00156 EvtSpinDensity *spinDensity ){
00157
00158
00159 EvtParticle *p;
00160
00161 if(spinDensity==0){
00162 p=EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(stdhepid),P);
00163 }
00164 else{
00165 p=EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(stdhepid),
00166 P,*spinDensity);
00167 }
00168
00169 generateDecay(p);
00170
00171
00172 evtStdHep->init();
00173
00174 p->makeStdHep(*evtStdHep);
00175
00176 evtStdHep->translate(D);
00177
00178 p->deleteTree();
00179
00180
00181 }
00182
00183 void EvtGen::generateDecay(EvtParticle *p){
00184
00185 int times=0;
00186 do{
00187 times+=1;
00188 EvtStatus::initRejectFlag();
00189 p->decay();
00190
00191 if ( EvtStatus::getRejectFlag()==0 ) {
00192 times=0;
00193 }
00194 else{
00195
00196 int ii;
00197 for (ii=0;ii<p->getNDaug();ii++){
00198 EvtParticle *temp=p->getDaug(ii);
00199 temp->deleteTree();
00200 }
00201 p->resetFirstOrNot();
00202 p->resetNDaug();
00203
00204 }
00205
00206 if ( times==10000) {
00207 report(ERROR,"EvtGen") << "Your event has been rejected 10000 times!"<<endl;
00208 report(ERROR,"EvtGen") << "Will now abort."<<endl;
00209 ::abort();
00210 times=0;
00211 }
00212 } while (times);
00213
00214 }
00215
00216
00217
00218 void EvtGen::generateEvent(int stdhepid,HepLorentzVector P,HepLorentzVector D){
00219
00220 EvtParticle *root_part;
00221 EvtVectorParticle *vector_part;
00222
00223 vector_part=new EvtVectorParticle;
00224 EvtVector4R p_init;
00225
00226 p_init.set(P.t(),P.x(),P.y(),P.z());
00227
00228 vector_part->init(EvtPDL::evtIdFromStdHep(stdhepid),p_init);
00229
00230 root_part=(EvtParticle *)vector_part;
00231
00232 root_part->setVectorSpinDensity();
00233
00234 generateEvent(root_part,D);
00235
00236 root_part->deleteTree();
00237
00238 }
00239
00240 void EvtGen::generateEvent(EvtParticle *root_part,HepLorentzVector D){
00241 int i;
00242 static int nevent=0;
00243 nevent++;
00244
00245 static EvtStdHep evtstdhep;
00246
00247
00248 int j;
00249 int istat;
00250 int partnum;
00251 double px,py,pz,e,m;
00252 double x,y,z,t;
00253
00254 EvtVector4R p4,x4;
00255 generateDecay(root_part);
00256
00257 int npart=0;
00258 EvtId list_of_stable[10];
00259 EvtParticle* stable_parent[10];
00260
00261
00262 list_of_stable[0]=EvtId(-1,-1);
00263 stable_parent[0]=0;
00264
00265 evtstdhep.init();
00266
00267
00268 root_part->makeStdHep(evtstdhep);
00269
00270
00271
00272
00273 npart=evtstdhep.getNPart();
00274 for(i=0;i<evtstdhep.getNPart();i++){
00275
00276 j=i+1;
00277
00278 int jmotherfirst=evtstdhep.getFirstMother(i)+1;
00279 int jmotherlast=evtstdhep.getLastMother(i)+1;
00280 int jdaugfirst=evtstdhep.getFirstDaughter(i)+1;
00281 int jdauglast=evtstdhep.getLastDaughter(i)+1;
00282
00283 partnum=evtstdhep.getStdHepID(i);
00284
00285 istat=evtstdhep.getIStat(i);
00286
00287 p4=evtstdhep.getP4(i);
00288 x4=evtstdhep.getX4(i);
00289 px=p4.get(1);
00290 py=p4.get(2);
00291 pz=p4.get(3);
00292 e=p4.get(0);
00293
00294 x=x4.get(1)+D.x();
00295 y=x4.get(2)+D.y();
00296 z=x4.get(3)+D.z();
00297 t=x4.get(0)+D.t();
00298
00299 m=p4.mass();
00300
00301 begevtgenstore_(&j,&nevent,&npart,&istat,
00302 &partnum,&jmotherfirst,&jmotherlast,
00303 &jdaugfirst,&jdauglast,
00304 &px,&py,&pz,&e,
00305 &m,&x,&y,&z,&t);
00306 }
00307
00308 }
00309
00310
00311
00312
00313
00314
00315
00316
00317