00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "../EeToeeV/EeToeeV.h"
00013 #include "../EeToeeV/EeToeeVRandom.h"
00014 #include "BesKernel/IBesRndmGenSvc.h"
00015
00016
00017 #include "HepMC/GenEvent.h"
00018 #include "HepMC/GenVertex.h"
00019 #include "HepMC/GenParticle.h"
00020
00021 #include "GaudiKernel/MsgStream.h"
00022 #include "GaudiKernel/ISvcLocator.h"
00023 #include "GaudiKernel/AlgFactory.h"
00024 #include "GaudiKernel/DataSvc.h"
00025 #include "GaudiKernel/SmartDataPtr.h"
00026
00027 #include "GeneratorObject/McGenEvent.h"
00028 #include "CLHEP/Vector/LorentzVector.h"
00029 #include "cfortran/cfortran.h"
00030
00031 #include <stdlib.h>
00032 #include <time.h>
00033
00034 typedef struct {
00035 double q1[4][160001];
00036 double p1[4][160001];
00037 double q2[4][160001];
00038 double p2[4][160001];
00039 double q3[4][160001];
00040 } MOMSET_DEF;
00041
00042 #define MOMSET COMMON_BLOCK(MOMSET_DEF, momset)
00043 COMMON_BLOCK_DEF(MOMSET_DEF,MOMSET);
00044
00045
00046 typedef struct {
00047 double ecms;
00048 } BEAMENERGY_DEF;
00049 #define BEAMENERGY COMMON_BLOCK(BEAMENERGY_DEF, beamenergy)
00050 COMMON_BLOCK_DEF(BEAMENERGY_DEF,BEAMENERGY);
00051
00052
00053
00054 typedef struct {
00055 int vct;
00056 } VECTOR_DEF;
00057 #define VECTOR COMMON_BLOCK(VECTOR_DEF, vector)
00058 COMMON_BLOCK_DEF(VECTOR_DEF,VECTOR);
00059
00060
00061 typedef struct {
00062 int mccheck;
00063 } MCTRUTH_DEF;
00064 #define MCTRUTH COMMON_BLOCK(MCTRUTH_DEF, mctruth)
00065 COMMON_BLOCK_DEF(MCTRUTH_DEF,MCTRUTH);
00066
00067
00068
00069 extern "C" {
00070 extern float flat_();
00071 }
00072
00073 float flat_(){
00074 float rr;
00075 double rd=EeToeeVRandom::random();
00076
00077 rr=(float)rd;
00078 return (float)EeToeeVRandom::random();
00079
00080 }
00081 extern "C" {
00082 extern void intxs_();
00083 }
00084
00085
00086 PROTOCCALLSFSUB1(GEVENT,gevent,INT)
00087 #define GEVENT(NEVENTS) CCALLSFSUB1(GEVENT,gevent,INT,NEVENTS)
00088
00089
00090 EeToeeV::EeToeeV(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
00091 {
00092 declareProperty("Ecms", m_Ecms = 3.65);
00093 declareProperty("Vector", m_vect = "phi");
00094 declareProperty("WriteMC", m_mctruth = 0);
00095 }
00096
00097
00098
00099
00100 StatusCode EeToeeV::initialize(){
00101
00102 MsgStream log(messageService(), name());
00103
00104 log << MSG::INFO << "EeToeeV initialize" << endreq;
00105
00106
00107 static const bool CREATEIFNOTTHERE(true);
00108 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00109 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00110 {
00111 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00112 return RndmStatus;
00113 }
00114 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("EeToeeV");
00115 std::cout<<"==============================="<<engine<<endl;
00116 EeToeeVRandom::setRandomEngine(engine);
00117
00118 MCTRUTH.mccheck=m_mctruth;
00119 BEAMENERGY.ecms=m_Ecms;
00120 double mpar=10.;
00121 if(m_vect=="omega"){
00122 VECTOR.vct=1;
00123 mpar=0.782;
00124 }else if(m_vect=="phi"){
00125 VECTOR.vct=2;
00126 mpar=1.019;
00127 }else if(m_vect=="J/psi"){
00128 VECTOR.vct=3;
00129 mpar=3.097;
00130 }else if(m_vect=="psi(2S)"){
00131 VECTOR.vct=4;
00132 mpar=3.686;
00133 }else if(m_vect=="psi(3770)"){
00134 VECTOR.vct=5;
00135 mpar=3.773;
00136 }else if(m_vect=="psi(4040)"){
00137 VECTOR.vct=6;
00138 mpar=4.039;
00139 }else if(m_vect=="psi(4160)"){
00140 VECTOR.vct=7;
00141 mpar=4.153;
00142 }else if(m_vect=="psi(4415)"){
00143 VECTOR.vct=8;
00144 mpar=4.421;
00145 }else{
00146 std::cout<<"EeToeeV::initialize() Bad vector "<<std::endl; abort();
00147 }
00148 if(m_Ecms<mpar){std::cout<<"EeToeeV:initialize: the Ecms less than the vector mass"<<std::endl;abort();}
00149
00150
00151 getMaxEvent();
00152 std::cout<<"m_evtMax = "<<m_evtMax<<std::endl;
00153 intxs_();
00154 GEVENT(m_evtMax);
00155
00156 return StatusCode::SUCCESS;
00157 }
00158
00159 static int evtgen=1;
00160 StatusCode EeToeeV::execute()
00161 {
00162 MsgStream log(messageService(), name());
00163
00164
00165 int npart = 0;
00166 int pid1,pid2,pid3,pid4,pst1,pst2;
00167 pid1 = 11;
00168 pid2 =-11;
00169 pst1 = 1;
00170 pst2 = 1;
00171
00172 if(m_vect=="omega"){
00173 pid3=223;
00174 }else if(m_vect=="phi"){
00175 pid3=333;
00176 }else if(m_vect=="J/psi"){
00177 pid3=443;
00178 }else if(m_vect=="psi(2S)"){
00179 pid3=100443;
00180 }else if(m_vect=="psi(3770)"){
00181 pid3=30443;
00182 }else if(m_vect=="psi(4040)"){
00183 pid3=9000443;
00184 }else if(m_vect=="psi(4160)"){
00185 pid3=9010443;
00186 }else if(m_vect=="psi(4415)"){
00187 pid3=9020443;
00188 }
00189
00190
00191 GenEvent* evt = new GenEvent(1,1);
00192
00193 GenVertex* prod_vtx = new GenVertex();
00194
00195 evt->add_vertex( prod_vtx );
00196
00197
00198 GenParticle* p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p1[1][evtgen], MOMSET.p1[2][evtgen],
00199 MOMSET.p1[3][evtgen], MOMSET.p1[0][evtgen]),
00200 11, 3);
00201 p->suggest_barcode( ++npart );
00202 prod_vtx->add_particle_in(p);
00203
00204
00205
00206 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q1[1][evtgen], MOMSET.q1[2][evtgen],
00207 MOMSET.q1[3][evtgen], MOMSET.q1[4][evtgen]),
00208 -11, 3);
00209
00210 p->suggest_barcode( ++npart );
00211 prod_vtx->add_particle_in(p);
00212
00213
00214 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p2[1][evtgen], MOMSET.p2[2][evtgen],
00215 MOMSET.p2[3][evtgen], MOMSET.p2[0][evtgen]),
00216 pid1,pst1);
00217 p->suggest_barcode( ++npart );
00218 prod_vtx->add_particle_out(p);
00219
00220
00221 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q2[1][evtgen], MOMSET.q2[2][evtgen],
00222 MOMSET.q2[3][evtgen], MOMSET.q2[0][evtgen]),
00223 pid2, pst1);
00224 p->suggest_barcode( ++npart );
00225 prod_vtx->add_particle_out(p);
00226
00227
00228 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q3[1][evtgen], MOMSET.q3[2][evtgen],
00229 MOMSET.q3[3][evtgen], MOMSET.q3[0][evtgen]),
00230 pid3,pst2);
00231 p->suggest_barcode( ++npart );
00232 prod_vtx->add_particle_out(p);
00233
00234 evtgen++;
00235
00236 if( log.level() < MSG::INFO )
00237 {
00238 evt->print();
00239 }
00240
00241
00242 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00243 if (anMcCol!=0)
00244 {
00245
00246 MsgStream log(messageService(), name());
00247 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00248 McGenEvent* mcEvent = new McGenEvent(evt);
00249 anMcCol->push_back(mcEvent);
00250 }
00251 else
00252 {
00253
00254
00255 McGenEventCol *mcColl = new McGenEventCol;
00256 McGenEvent* mcEvent = new McGenEvent(evt);
00257 mcColl->push_back(mcEvent);
00258 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00259 if (sc != StatusCode::SUCCESS)
00260 {
00261 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00262 delete mcColl;
00263 delete evt;
00264 delete mcEvent;
00265 return StatusCode::FAILURE;
00266 }
00267 else
00268 {
00269
00270 }
00271 }
00272
00273 return StatusCode::SUCCESS;
00274 }
00275
00276 StatusCode EeToeeV::finalize()
00277 {
00278 MsgStream log(messageService(), name());
00279 char delcmd[300];
00280 strcpy(delcmd,"cat ");
00281 strcat(delcmd,"fresult.dat");
00282 system(delcmd);
00283
00284 std::cout<< "EeToeeV finalized" << endl;
00285 return StatusCode::SUCCESS;
00286 }
00287
00288 StatusCode EeToeeV::getMaxEvent() {
00289 IProperty* appPropMgr=0;
00290 StatusCode status =
00291 serviceLocator()->getService("ApplicationMgr", IProperty::interfaceID(),
00292 reinterpret_cast<IInterface*&>( appPropMgr ));
00293 if( status.isFailure() ) return status;
00294
00295 IntegerProperty evtMax("EvtMax",0);
00296 status = appPropMgr->getProperty( &evtMax );
00297 if (status.isFailure()) return status;
00298
00299 m_evtMax = evtMax.value();
00300 return status;
00301 }
00302