00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "Bhlumi/Bhlumi.h"
00014 #include "Bhlumi/BhlumiRandom.h"
00015
00016
00017 #include "HepMC/GenEvent.h"
00018 #include "HepMC/GenVertex.h"
00019 #include "HepMC/GenParticle.h"
00020 #include "CLHEP/Vector/LorentzVector.h"
00021
00022 #include "GaudiKernel/MsgStream.h"
00023 #include "GaudiKernel/ISvcLocator.h"
00024 #include "GaudiKernel/AlgFactory.h"
00025 #include "GaudiKernel/DataSvc.h"
00026 #include "GaudiKernel/SmartDataPtr.h"
00027
00028 #include "GeneratorObject/McGenEvent.h"
00029 #include "BesKernel/IBesRndmGenSvc.h"
00030
00031 #include "cfortran/cfortran.h"
00032
00033 #include <stdlib.h>
00034
00035 using namespace std;
00036
00037
00038
00039
00040 typedef struct {
00041 double p1[4];
00042 double q1[4];
00043 double p2[4];
00044 double q2[4];
00045 double phot[4][100];
00046 int nphot; } MOMSET_DEF;
00047
00048 #define MOMSET COMMON_BLOCK(MOMSET_DEF, momset)
00049 COMMON_BLOCK_DEF(MOMSET_DEF,MOMSET);
00050
00051
00052
00053
00054 PROTOCCALLSFSUB1(GLIMIT,glimit,INT)
00055 #define GLIMIT(LENMX) CCALLSFSUB1(GLIMIT,glimit,INT,LENMX)
00056
00057 PROTOCCALLSFSUB1(DUMPS,dumps,INT)
00058 #define DUMPS(NOUT) CCALLSFSUB1(DUMPS,dumps,INT,NOUT)
00059
00060 PROTOCCALLSFSUB3(BHLUMI,bhlumi,INT,DOUBLEV,INTV)
00061 #define BHLUMI(MODE,XPAR,NPAR) CCALLSFSUB3(BHLUMI,bhlumi,INT,DOUBLEV,INTV,MODE,XPAR,NPAR)
00062
00063 PROTOCCALLSFSUB3(MARINI, marini, INT, INT, INT)
00064 #define MARINI(IJKLIN, NTOTIN, NTOT2N) CCALLSFSUB3(MARINI, marini, INT, INT, INT, IJKLIN, NTOTIN, NTOT2N)
00065
00066 extern "C" {
00067 extern double ranmarr_();
00068 extern void marran_(float* rvec, int* nma);
00069 extern void ecuran_(double* rvec, int* nma);
00070 extern void carran_(double* rvec, int* nma);
00071 }
00072
00073 double ranmarr_()
00074 {
00075 return BhlumiRandom::random();
00076 }
00077
00078 void marran_(float* rvec, int* nma)
00079 {
00080 int nmax = *nma;
00081 assert(nmax<100);
00082 double rvecd[100];
00083 BhlumiRandom::FlatArray(rvecd, nmax);
00084 for(int i=0; i<nmax; i++) rvec[i]=rvecd[i];
00085 }
00086
00087 void carran_(double* rvec, int* nma)
00088 {
00089 int nmax = *nma;
00090 assert(nmax<100);
00091 double rvecd[100];
00092 BhlumiRandom::FlatArray(rvecd, nmax);
00093 for(int i=0; i<nmax; i++) rvec[i]=rvecd[i];
00094 }
00095
00096 void ecuran_(double* rvec, int* nma)
00097 {
00098 int nmax = *nma;
00099 assert(nmax<100);
00100 double rvecd[100];
00101 BhlumiRandom::FlatArray(rvecd, nmax);
00102 for(int i=0; i<nmax; i++) rvec[i]=rvecd[i];
00103 }
00104
00105
00106 Bhlumi::Bhlumi(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
00107 {
00108 declareProperty("CMEnergy", m_cmEnergy = 3.097);
00109 declareProperty("AngleMode", m_angleMode = 0);
00110
00111 declareProperty("MinThetaAngle", m_minThetaAngle = 0.24);
00112 declareProperty("MaxThetaAngle", m_maxThetaAngle = 0.58);
00113 declareProperty("SoftPhotonCut", m_infraredCut = 1E-4);
00114 m_initSeed.clear();
00115
00116
00117 }
00118
00119 StatusCode Bhlumi::initialize(){
00120
00121 MsgStream log(messageService(), name());
00122
00123 log << MSG::INFO << "Bhlumi initialize" << endreq;
00124
00125 static const bool CREATEIFNOTTHERE(true);
00126 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00127 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00128 {
00129 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00130 return RndmStatus;
00131 }
00132 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("Bhlumi");
00133 engine->showStatus();
00134 BhlumiRandom::setRandomEngine(engine);
00135
00136 GLIMIT(50000);
00137
00147
00148
00149
00150
00151 int KeyGen = 3;
00152 int KeyRem = 1;
00153 KeyRem = 0;
00155 int KeyWgt = 2;
00156 KeyWgt = 0;
00157 int KeyRnd = 1;
00158 int KeyOpt =1000*KeyGen +100*KeyRem +10*KeyWgt +KeyRnd;
00161 int KeyPia = 0;
00162 int KeyMod = 2;
00163 KeyPia = 2;
00164 int KeyZet = 0;
00165 KeyZet = 1;
00166 int KeyRad =1000*KeyZet +10*KeyMod +KeyPia;
00168 npar[0]= KeyOpt;
00169 npar[1]= KeyRad;
00170 double CmsEne = m_cmEnergy;
00171 xpar[0]= CmsEne;
00172 double th1,th2,thmin,thmax;
00173 if(m_angleMode==0){
00174 th1 = m_minThetaAngle;
00175 th2 = m_maxThetaAngle;
00176 thmin = 0.7*th1;
00177 thmax = 2.0*th2;
00178 if(thmin<0.||thmax>3.1415926) {
00179 log << MSG::WARNING << "input angles exceed range (0~pi), so effect will be unprospectable" << endreq;
00180 return StatusCode::FAILURE;
00181 }
00182 }
00183 else if(m_angleMode==1){
00184 th1 = m_minThetaAngle*3.1415926/180.;
00185 th2 = m_maxThetaAngle*3.1415926/180.;
00186
00187 thmin = th1;
00188 thmax = th2;
00189 }
00190 else if(m_angleMode==2){
00191 th1 = acos(max(m_minThetaAngle,m_maxThetaAngle));
00192 th2 = acos(min(m_minThetaAngle,m_maxThetaAngle));
00193
00194 thmin = th1;
00195 thmax = th2;
00196 }
00197 else{
00198 log << MSG::FATAL << "unknown angle mode!" << endreq;
00199 return StatusCode::FAILURE;
00200 }
00201 if(thmin<0.||thmax>3.1415926) {
00202 log << MSG::FATAL << "input angles exceed range (0~pi), unprospectable" << endreq;
00203 return StatusCode::FAILURE;
00204 }
00205 else if(thmin>thmax) {
00206 log << MSG::FATAL << "thmin>thmax, unprospectable" << endreq;
00207 return StatusCode::FAILURE;
00208 }
00209 if(KeyWgt == 2) thmin=th1;
00210 xpar[1]= CmsEne*CmsEne*(1-cos(thmin))/2;
00211 xpar[2]= CmsEne*CmsEne*(1-cos(thmax))/2;
00212 xpar[3]= m_infraredCut;
00213
00214
00215
00216 BHLUMI(-1,xpar,npar);
00217
00218 return StatusCode::SUCCESS;
00219 }
00220
00221 StatusCode Bhlumi::execute()
00222 {
00223 MsgStream log(messageService(), name());
00224 log << MSG::INFO << "Bhlumi executing" << endreq;
00225
00226
00227 BHLUMI( 0,xpar,npar);
00228
00229 if( log.level() < MSG::INFO )
00230 {
00231 DUMPS(6);
00232
00233
00234 }
00235
00236 int npart = 0;
00237
00238
00239 GenEvent* evt = new GenEvent(1,1);
00240
00241 GenVertex* prod_vtx = new GenVertex();
00242
00243 evt->add_vertex( prod_vtx );
00244
00245
00246 GenParticle* p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p1[0], MOMSET.p1[1],
00247 MOMSET.p1[2], MOMSET.p1[3]),
00248 -11, 3);
00249 p->suggest_barcode( ++npart );
00250 prod_vtx->add_particle_in(p);
00251
00252
00253 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q1[0], MOMSET.q1[1], MOMSET.q1[2], MOMSET.q1[3]),
00254 11, 3);
00255 p->suggest_barcode( ++npart );
00256 prod_vtx->add_particle_in(p);
00257
00258
00259 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p2[0], MOMSET.p2[1],
00260 MOMSET.p2[2], MOMSET.p2[3]),
00261 -11, 1);
00262 p->suggest_barcode( ++npart );
00263 prod_vtx->add_particle_out(p);
00264
00265
00266 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q2[0], MOMSET.q2[1], MOMSET.q2[2], MOMSET.q2[3]),
00267 11, 1);
00268 p->suggest_barcode( ++npart );
00269 prod_vtx->add_particle_out(p);
00270
00271 int iphot=0;
00272 for (iphot=0; iphot<MOMSET.nphot; iphot++)
00273 {
00274
00275 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.phot[0][iphot], MOMSET.phot[1][iphot],
00276 MOMSET.phot[2][iphot], MOMSET.phot[3][iphot]),
00277 22, 1);
00278 p->suggest_barcode( ++npart );
00279 prod_vtx->add_particle_out(p);
00280 }
00281
00282 if( log.level() < MSG::INFO )
00283 {
00284 evt->print();
00285 }
00286
00287
00288 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00289 if (anMcCol!=0)
00290 {
00291
00292 MsgStream log(messageService(), name());
00293 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00294 McGenEvent* mcEvent = new McGenEvent(evt);
00295 anMcCol->push_back(mcEvent);
00296 }
00297 else
00298 {
00299
00300 McGenEventCol *mcColl = new McGenEventCol;
00301 McGenEvent* mcEvent = new McGenEvent(evt);
00302 mcColl->push_back(mcEvent);
00303 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00304 if (sc != StatusCode::SUCCESS)
00305 {
00306 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00307 delete mcColl;
00308 delete evt;
00309 delete mcEvent;
00310 return StatusCode::FAILURE;
00311 }
00312 else
00313 {
00314 log << MSG::INFO << "McGenEventCol created and " << npart <<" particles stored in McGenEvent" << endreq;
00315 }
00316 }
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336 return StatusCode::SUCCESS;
00337 }
00338
00339 StatusCode Bhlumi::finalize()
00340 {
00341 MsgStream log(messageService(), name());
00342
00343 BHLUMI( 2,xpar,npar);
00344
00345 log << MSG::INFO << "Bhlumi finalized" << endreq;
00346
00347 return StatusCode::SUCCESS;
00348 }