00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "../EeTo4e/EeTo4e.h"
00013 #include "../EeTo4e/EeTo4eRandom.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 double p3[4][160001];
00041 } MOMSET_DEF;
00042
00043 #define MOMSET COMMON_BLOCK(MOMSET_DEF, momset)
00044 COMMON_BLOCK_DEF(MOMSET_DEF,MOMSET);
00045
00046
00047 typedef struct {
00048 double ecms;
00049 } BEAMENERGY_DEF;
00050 #define BEAMENERGY COMMON_BLOCK(BEAMENERGY_DEF, beamenergy)
00051 COMMON_BLOCK_DEF(BEAMENERGY_DEF,BEAMENERGY);
00052
00053
00054
00055
00056 typedef struct {
00057 double setcos;
00058 } COSEE_DEF;
00059 #define COSEE COMMON_BLOCK(COSEE_DEF, cosee)
00060 COMMON_BLOCK_DEF(COSEE_DEF,COSEE);
00061
00062
00063 typedef struct {
00064 int mccheck;
00065 } MCTRUTH_DEF;
00066 #define MCTRUTH COMMON_BLOCK(MCTRUTH_DEF, mctruth)
00067 COMMON_BLOCK_DEF(MCTRUTH_DEF,MCTRUTH);
00068
00069
00070 extern "C" {
00071 extern float flat_();
00072 }
00073
00074 float flat_(){
00075 float rr;
00076 double rd=EeTo4eRandom::random();
00077
00078 rr=(float)rd;
00079 return (float)EeTo4eRandom::random();
00080
00081 }
00082 extern "C" {
00083 extern void intxs_();
00084 }
00085
00086
00087 PROTOCCALLSFSUB1(GEVENT,gevent,INT)
00088 #define GEVENT(NEVENTS) CCALLSFSUB1(GEVENT,gevent,INT,NEVENTS)
00089
00090
00091 EeTo4e::EeTo4e(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
00092 {
00093 declareProperty("Ecms", m_Ecms = 3.65);
00094 declareProperty("cosee", m_cosee = 0.93);
00095 declareProperty("WriteMC", m_mctruth = 0);
00096 }
00097
00098 StatusCode EeTo4e::initialize(){
00099
00100 MsgStream log(messageService(), name());
00101
00102 log << MSG::INFO << "EeTo4e initialize" << endreq;
00103
00104
00105 static const bool CREATEIFNOTTHERE(true);
00106 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00107 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00108 {
00109 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00110 return RndmStatus;
00111 }
00112 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("EeTo4e");
00113 std::cout<<"==============================="<<engine<<endl;
00114 EeTo4eRandom::setRandomEngine(engine);
00115
00116 MCTRUTH.mccheck=m_mctruth;
00117 COSEE.setcos=m_cosee;
00118 BEAMENERGY.ecms=m_Ecms;
00119
00120
00121 getMaxEvent();
00122 std::cout<<"m_evtMax = "<<m_evtMax<<std::endl;
00123 intxs_();
00124 GEVENT(m_evtMax);
00125
00126 return StatusCode::SUCCESS;
00127 }
00128
00129 static int evtgen=1;
00130 StatusCode EeTo4e::execute()
00131 {
00132 MsgStream log(messageService(), name());
00133
00134
00135 int npart = 0;
00136 int pid1,pid2,pid3,pid4,pst1,pst2;
00137 pid1 = 11;
00138 pid2 =-11;
00139 pid3 = 11;
00140 pid4 =-11;
00141 pst1 = 1;
00142 pst2 = 1;
00143
00144
00145
00146
00147 GenEvent* evt = new GenEvent(1,1);
00148
00149 GenVertex* prod_vtx = new GenVertex();
00150
00151 evt->add_vertex( prod_vtx );
00152
00153
00154 GenParticle* p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p1[1][evtgen], MOMSET.p1[2][evtgen],
00155 MOMSET.p1[3][evtgen], MOMSET.p1[0][evtgen]),
00156 11, 3);
00157 p->suggest_barcode( ++npart );
00158 prod_vtx->add_particle_in(p);
00159
00160
00161
00162 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q1[1][evtgen], MOMSET.q1[2][evtgen],
00163 MOMSET.q1[3][evtgen], MOMSET.q1[4][evtgen]),
00164 -11, 3);
00165
00166 p->suggest_barcode( ++npart );
00167 prod_vtx->add_particle_in(p);
00168
00169
00170 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p2[1][evtgen], MOMSET.p2[2][evtgen],
00171 MOMSET.p2[3][evtgen], MOMSET.p2[0][evtgen]),
00172 pid1,pst1);
00173 p->suggest_barcode( ++npart );
00174 prod_vtx->add_particle_out(p);
00175
00176
00177 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q2[1][evtgen], MOMSET.q2[2][evtgen],
00178 MOMSET.q2[3][evtgen], MOMSET.q2[0][evtgen]),
00179 pid2, pst2);
00180 p->suggest_barcode( ++npart );
00181 prod_vtx->add_particle_out(p);
00182
00183
00184 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p3[1][evtgen], MOMSET.p3[2][evtgen],
00185 MOMSET.p3[3][evtgen], MOMSET.p3[0][evtgen]),
00186 pid3,pst1);
00187 p->suggest_barcode( ++npart );
00188 prod_vtx->add_particle_out(p);
00189
00190
00191 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q3[1][evtgen], MOMSET.q3[2][evtgen],
00192 MOMSET.q3[3][evtgen], MOMSET.q3[0][evtgen]),
00193 pid4, pst2);
00194 p->suggest_barcode( ++npart );
00195 prod_vtx->add_particle_out(p);
00196
00197
00198
00199 evtgen++;
00200
00201 if( log.level() < MSG::INFO )
00202 {
00203 evt->print();
00204 }
00205
00206
00207 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00208 if (anMcCol!=0)
00209 {
00210
00211 MsgStream log(messageService(), name());
00212 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00213 McGenEvent* mcEvent = new McGenEvent(evt);
00214 anMcCol->push_back(mcEvent);
00215 }
00216 else
00217 {
00218
00219 McGenEventCol *mcColl = new McGenEventCol;
00220 McGenEvent* mcEvent = new McGenEvent(evt);
00221 mcColl->push_back(mcEvent);
00222 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00223 if (sc != StatusCode::SUCCESS)
00224 {
00225 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00226 delete mcColl;
00227 delete evt;
00228 delete mcEvent;
00229 return StatusCode::FAILURE;
00230 }
00231 else
00232 {
00233
00234 }
00235 }
00236
00237 return StatusCode::SUCCESS;
00238 }
00239
00240 StatusCode EeTo4e::finalize()
00241 {
00242 MsgStream log(messageService(), name());
00243 char delcmd[300];
00244 strcpy(delcmd,"cat ");
00245 strcat(delcmd,"fresult.dat");
00246 system(delcmd);
00247
00248 std::cout<< "EeTo4e finalized" << endl;
00249 return StatusCode::SUCCESS;
00250 }
00251
00252 StatusCode EeTo4e::getMaxEvent() {
00253 IProperty* appPropMgr=0;
00254 StatusCode status =
00255 serviceLocator()->getService("ApplicationMgr", IProperty::interfaceID(),
00256 reinterpret_cast<IInterface*&>( appPropMgr ));
00257 if( status.isFailure() ) return status;
00258
00259 IntegerProperty evtMax("EvtMax",0);
00260 status = appPropMgr->getProperty( &evtMax );
00261 if (status.isFailure()) return status;
00262
00263 m_evtMax = evtMax.value();
00264 return status;
00265 }
00266