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