00001
00002
00003
00004
00005
00006
00007
00008 #include "BesBdkRc/BesBdkRc.h"
00009 #include "BesBdkRc/BesBdkRcRandom.h"
00010
00011
00012 #include "HepMC/GenEvent.h"
00013
00014 #include "HepMC/HEPEVT_Wrapper.h"
00015 #include "HepMC/IO_HEPEVT.h"
00016 #include "HepMC/GenEvent.h"
00017
00018 #include "GaudiKernel/MsgStream.h"
00019 #include "GaudiKernel/ISvcLocator.h"
00020 #include "GaudiKernel/AlgFactory.h"
00021 #include "GaudiKernel/DataSvc.h"
00022 #include "GaudiKernel/SmartDataPtr.h"
00023 #include "BesKernel/IBesRndmGenSvc.h"
00024
00025 #include "GeneratorObject/McGenEvent.h"
00026 #include "cfortran/cfortran.h"
00027 #include <cstdlib>
00028 #include <cassert>
00029
00030 using namespace std;
00031
00032 typedef struct{
00033 double EB,W2MIN,K0,KS,EWE;
00034 }INPUT_DEF;
00035 #define INPUT COMMON_BLOCK(INPUT_DEF,input)
00036 COMMON_BLOCK_DEF(INPUT_DEF,INPUT);
00037
00038 typedef struct{
00039 double ME,MU,PI,ALFA,BARN;
00040 } PCONST_DEF;
00041 #define PCONST COMMON_BLOCK(PCONST_DEF, pconst)
00042 COMMON_BLOCK_DEF(PCONST_DEF,PCONST);
00043
00044 typedef struct{
00045 float TCHMIN, PCHMIN;
00046 } FINCUT_DEF;
00047 #define FINCUT COMMON_BLOCK(FINCUT_DEF, fincut)
00048 COMMON_BLOCK_DEF(FINCUT_DEF,FINCUT);
00049
00050 typedef struct{
00051 float W2MINR, EBEAM, MAXESW, KZERO;
00052 int ISEED;
00053 float QMASS[6];
00054 }LOCALC_DEF;
00055 #define LOCALC COMMON_BLOCK(LOCALC_DEF, localc)
00056 COMMON_BLOCK_DEF(LOCALC_DEF,LOCALC);
00057
00058 typedef struct{
00059 double COLCHA;
00060 int IFINAL;
00061 }COLCHC_DEF;
00062 #define COLCHC COMMON_BLOCK(COLCHC_DEF, colchc)
00063 COMMON_BLOCK_DEF(COLCHC_DEF,COLCHC);
00064
00065 typedef struct{
00066 int NEVCUT, NEVTOT, MAXNTRY;
00067 }DELPHC_DEF;
00068 #define DELPHC COMMON_BLOCK(DELPHC_DEF, delphc)
00069 COMMON_BLOCK_DEF(DELPHC_DEF,DELPHC);
00070
00071 typedef struct{
00072 int ITOT, NEVUNW;
00073 double CSTOT,CSERR;
00074 float CSCUT, CSCUTERR;
00075 }XSECTC_DEF;
00076 #define XSECTC COMMON_BLOCK(XSECTC_DEF, xsectc)
00077 COMMON_BLOCK_DEF(XSECTC_DEF,XSECTC);
00078
00079 typedef struct{
00080 int GOODEVT;
00081 }FLAGS_DEF;
00082 #define FLAGS COMMON_BLOCK(FLAGS_DEF, flags)
00083 COMMON_BLOCK_DEF(FLAGS_DEF,FLAGS);
00084
00085 typedef struct{
00086 int N;
00087 int K[5][4000];
00088 float P[5][4000],V[5][4000];
00089 }LUJETS_DEF;
00090 #define LUJETS COMMON_BLOCK(LUJETS_DEF,lujets)
00091 COMMON_BLOCK_DEF(LUJETS_DEF,LUJETS);
00092
00093 PROTOCCALLSFSUB1(LUHEPC, luhepc, INT)
00094 #define LUHEPC(ICONV) CCALLSFSUB1(LUHEPC, luhepc, INT, ICONV)
00095
00096 PROTOCCALLSFSUB1(LULIST, lulist, INT)
00097 #define LULIST(ICONV) CCALLSFSUB1(LULIST, lulist, INT, ICONV)
00098
00099 PROTOCCALLSFSUB0(HEPEVT_CLEAN, hepevt_clean)
00100 #define HEPEVT_CLEAN() CCALLSFSUB0(HEPEVT_CLEAN, hepevt_clean)
00101
00102 PROTOCCALLSFSUB0(HEPEVT_PRINT, hepevt_print)
00103 #define HEPEVT_PRINT() CCALLSFSUB0(HEPEVT_PRINT, hepevt_print)
00104
00105 PROTOCCALLSFSUB1(BDKRC,bdkrc,INT)
00106 #define BDKRC(IDUMP) CCALLSFSUB1(BDKRC, bdkrc, INT, IDUMP)
00107
00108 PROTOCCALLSFSUB3(FINISH,finish,PINT,PDOUBLE,PDOUBLE)
00109 #define FINISH(IN,SIGT,ER) CCALLSFSUB3(FINISH,finish,PINT,PDOUBLE,PDOUBLE,IN,SIGT,ER)
00110
00111 PROTOCCALLSFSUB4(RLUXGO,rluxgo,INT,INT,INT,INT)
00112 #define RLUXGO(lux,ISEED,K1,K2) CCALLSFSUB4(RLUXGO,rluxgo,INT,INT,INT,INT,lux,ISEED,K1,K2)
00113
00114 PROTOCCALLSFSUB0(GEN_INIT, gen_init)
00115 #define GEN_INIT() CCALLSFSUB0(GEN_INIT, gen_init)
00116
00117 PROTOCCALLSFSUB0(GEN1EVT, gen1evt)
00118 #define GEN1EVT() CCALLSFSUB0(GEN1EVT, gen1evt)
00119
00120 extern "C" {
00121 extern float flat_();
00122 }
00123
00124 float flat_()
00125 {
00126 return (float)BesBdkRcRandom::random();
00127 }
00128
00129 BesBdkRc::BesBdkRc(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
00130 {
00131
00132 declareProperty("CMEnergy", m_CMEnergy = 3.097);
00133 declareProperty("W2min",m_w2min=0.02);
00134 declareProperty("EstimatedMaxWeight", m_ewe=2.5);
00135 declareProperty("SoftPhotonMaxEnergy", m_kzero=0.002);
00136 qmass[0] = 0.2;
00137 qmass[1] = 0.2;
00138 qmass[2] = 0.5;
00139 qmass[3] = 1.5;
00140 qmass[4] = 4.5;
00141 qmass[5] = 180.;
00142 declareProperty("MaxNTry", m_maxNTry);
00143 declareProperty("FinalState", m_ifinal);
00144 declareProperty("MinTheta",m_tcmin);
00145 declareProperty("MinMomentum",m_pcmin);
00146
00147 m_numberEvent=0;
00148 toRad=M_PI/180.0;
00149 toDeg=180.0/M_PI;
00150 }
00151
00152 StatusCode BesBdkRc::initialize(){
00153 MsgStream log(messageService(), name());
00154 log << MSG::WARNING << "BesBdkRc initialize" << endreq;
00155
00156 static const bool CREATEIFNOTTHERE(true);
00157 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00158 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00159 {
00160 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00161 return RndmStatus;
00162 }
00163 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("BesBdkRc");
00164 engine->showStatus();
00165 BesBdkRcRandom::setRandomEngine(engine);
00166
00167
00168
00169
00170 FINCUT.TCHMIN=m_tcmin*toRad;
00171 FINCUT.PCHMIN=m_pcmin;
00172
00173
00174 DELPHC.MAXNTRY=m_maxNTry;
00175
00176 COLCHC.IFINAL=m_ifinal;
00177
00178 INPUT.EB=0.5*m_CMEnergy;
00179
00180 LOCALC.W2MINR=m_w2min;
00181 LOCALC.EBEAM=0.5*m_CMEnergy;
00182 LOCALC.MAXESW=m_ewe;
00183 LOCALC.KZERO=m_kzero;
00184
00185 for(int i=0; i<6;i++)LOCALC.QMASS[i]=qmass[i];
00186 GEN_INIT();
00187 XSECTC.NEVUNW=0;
00188
00189 return StatusCode::SUCCESS;
00190 }
00191
00192 StatusCode BesBdkRc::execute()
00193 {
00194 MsgStream log(messageService(), name());
00195 log << MSG::INFO << "BesBdkRc executing" << endreq;
00196 HepMC::HEPEVT_Wrapper::set_max_number_entries(2000);
00197 HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
00198 HepMC::IO_HEPEVT HepEvtIO;
00199
00200 HEPEVT_CLEAN();
00201 GEN1EVT();
00202
00203 if(FLAGS.GOODEVT!=1){
00204 log << MSG::ERROR<<" BesBdkRc: fail to generate good event"<<endl;
00205 return StatusCode::FAILURE;
00206 }
00207
00208 m_numberEvent++;
00209 if( log.level() < MSG::INFO )LULIST(1);
00210
00211 LUHEPC(1);
00212 HepMC::GenEvent* evt = HepEvtIO.read_next_event();
00213 evt->set_event_number(m_numberEvent);
00214 evt->set_signal_process_id(1);
00215
00216
00217 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00218 if (anMcCol!=0) {
00219
00220 MsgStream log(messageService(), name());
00221 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00222 McGenEvent* mcEvent = new McGenEvent(evt);
00223 anMcCol->push_back(mcEvent);
00224 } else {
00225
00226 McGenEventCol *mcColl = new McGenEventCol;
00227 McGenEvent* mcEvent = new McGenEvent(evt);
00228 mcColl->push_back(mcEvent);
00229 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00230 if (sc != StatusCode::SUCCESS) {
00231 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00232 delete mcColl;
00233 delete evt;
00234 delete mcEvent;
00235 return StatusCode::FAILURE;
00236 }
00237 }
00238
00239
00240
00241 return StatusCode::SUCCESS;
00242 }
00243
00244 StatusCode BesBdkRc::finalize()
00245 {
00246 MsgStream log(messageService(), name());
00247 log << MSG::INFO << "BesBdkRc finalized" << endreq;
00248 int itot;
00249 itot=XSECTC.ITOT;
00250 double cstot;
00251 double cserr;
00252 FINISH(itot,cstot,cserr);
00253 float effcut=0;
00254 float cscut=0;
00255 float efferr=0;
00256 float cscuterr=0;
00257 if(XSECTC.NEVUNW){
00258 effcut=float(DELPHC.NEVCUT)/float(XSECTC.NEVUNW);
00259 cscut=effcut*cstot;
00260 efferr=sqrt(effcut*(1.0-effcut)/float(XSECTC.NEVUNW));
00261 cscuterr = sqrt(cstot*efferr*cstot*efferr + effcut*effcut*cserr*cserr);
00262 }
00263 printf("BDKRC SUMMARY: Cross section after user cuts= %G +- %G nb\n",cscut,cscuterr);
00264 printf(" Cut acceptance = %G +- %G \n",effcut,efferr);
00265 return StatusCode::SUCCESS;
00266 }