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