00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "../Babayaga/Babayaga.h"
00013 #include "../Babayaga/BabayagaRandom.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
00033
00034
00035
00036
00037 typedef struct {
00038 double q1[4][160001];
00039 double p1[4][160001];
00040 double q2[4][160001];
00041 double p2[4][160001];
00042 double phot[4][10][160001];
00043 } MOMSET_DEF;
00044
00045 #define MOMSET COMMON_BLOCK(MOMSET_DEF, momset)
00046 COMMON_BLOCK_DEF(MOMSET_DEF,MOMSET);
00047 typedef struct {
00048 double CutNgam,CutEgam;
00049 } RADBB_DEF;
00050
00051 #define RADBB COMMON_BLOCK(RADBB_DEF, radbb)
00052 COMMON_BLOCK_DEF(RADBB_DEF,RADBB);
00053
00054
00055 typedef struct{
00056 int ncqph[160000];
00057 }ISRPHOTONS_DEF;
00058 #define ISRPHOTONS COMMON_BLOCK(ISRPHOTONS_DEF, isrphotons)
00059 COMMON_BLOCK_DEF(ISRPHOTONS_DEF,ISRPHOTONS);
00060
00061
00062 typedef struct {
00063 double ebeam;
00064 } BEAMENERGY_DEF;
00065 #define BEAMENERGY COMMON_BLOCK(BEAMENERGY_DEF, beamenergy)
00066 COMMON_BLOCK_DEF(BEAMENERGY_DEF,BEAMENERGY);
00067
00068
00069
00070 typedef struct {
00071 double thmin,thmax,emin,zmax,egmin,thgmin,thgmax;
00072 } EXPCUTS_DEF;
00073 #define EXPCUTS COMMON_BLOCK(EXPCUTS_DEF, expcuts)
00074 COMMON_BLOCK_DEF(EXPCUTS_DEF,EXPCUTS);
00075
00076
00077
00078 typedef struct{
00079 int on,ill;
00080 } SWITCH_DEF;
00081 #define SWITCH COMMON_BLOCK(SWITCH_DEF,switchfsr)
00082 COMMON_BLOCK_DEF(SWITCH_DEF,SWITCH);
00083
00084
00085 typedef struct{
00086 int iarun;
00087 } SWITCHARUN_DEF;
00088 #define SWITCHARUN COMMON_BLOCK(SWITCHARUN_DEF,switcharun)
00089 COMMON_BLOCK_DEF(SWITCHARUN_DEF,SWITCHARUN);
00090
00091
00092 typedef struct{
00093 int ich;
00094 } CHANNEL_DEF;
00095 #define CHANNEL COMMON_BLOCK(CHANNEL_DEF,channel)
00096 COMMON_BLOCK_DEF(CHANNEL_DEF,CHANNEL);
00097
00098
00099 typedef struct{
00100 int ires;
00101 } RESONANCES_DEF;
00102 #define RESONANCES COMMON_BLOCK(RESONANCES_DEF,resonances)
00103 COMMON_BLOCK_DEF(RESONANCES_DEF,RESONANCES);
00104
00105
00106 typedef struct{
00107 int seed;
00108 } DECLAREINT_DEF;
00109 #define DECLAREINT COMMON_BLOCK(DECLAREINT_DEF,declareint)
00110 COMMON_BLOCK_DEF(DECLAREINT_DEF, DECLAREINT );
00111
00112
00113 typedef struct{
00114 char tuple,phcut,cutg;
00115 } DECLARESTR_DEF;
00116 #define DECLARESTR COMMON_BLOCK(DECLARESTR_DEF,declarestr)
00117 COMMON_BLOCK_DEF(DECLARESTR_DEF,DECLARESTR);
00118
00119
00120
00121
00122 extern "C" {
00123 extern float flat_();
00124 }
00125
00126 float flat_(){
00127 float rr;
00128 double rd=BabayagaRandom::random();
00129
00130 rr=(float)rd;
00131 return (float)BabayagaRandom::random();
00132
00133 }
00134
00135
00136
00137 PROTOCCALLSFSUB1(BABAYAGA,babayaga,INT)
00138 #define BABAYAGA(NEVENTS) CCALLSFSUB1(BABAYAGA,babayaga,INT,NEVENTS)
00139
00140
00141 Babayaga::Babayaga(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
00142 {
00143
00144 declareProperty("Channel", m_Ich = 1);
00145 declareProperty("Ebeam", m_Ebeam = 1.548);
00146 declareProperty("MinThetaAngle", m_Thmin = 70);
00147 declareProperty("MaxThetaAngle", m_Thmax = 178);
00148 declareProperty("MinimumEnergy", m_Emin = 0.4);
00149 declareProperty("MaximumAcollinearity", m_Zmax = 10);
00150 declareProperty("RunningAlpha", m_Iarun = 1);
00151 declareProperty("HadronicResonance", m_Ires = 0);
00152 declareProperty("FSR_swich", m_on = 0);
00153 declareProperty("MinEnerCutG", m_Egmin = 0.01);
00154 declareProperty("MinAngCutG", m_Thgmin = 5);
00155 declareProperty("MaxAngCutG", m_Thgmax = 21);
00156 declareProperty("HBOOK", HN= 1);
00157 declareProperty("PHCUT", m_PHCUT = 0);
00158 declareProperty("CUTG", m_CUTG = 0);
00159 declareProperty("CutNgam", m_CutNgam = 0);
00160 declareProperty("CutEgam", m_CutEgam = 0);
00161 }
00162
00163 StatusCode Babayaga::initialize(){
00164
00165 MsgStream log(messageService(), name());
00166
00167 log << MSG::INFO << "Babayaga initialize" << endreq;
00168
00169
00170 static const bool CREATEIFNOTTHERE(true);
00171 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00172 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00173 {
00174 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00175 return RndmStatus;
00176 }
00177 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("Babayaga");
00178 std::cout<<"==============================="<<engine<<endl;
00179 BabayagaRandom::setRandomEngine(engine);
00180
00181
00182 if(HN==1) {DECLARESTR.tuple='Y';} else DECLARESTR.tuple='N';
00183 if(m_PHCUT==1){ DECLARESTR.phcut='Y';} else DECLARESTR.cutg='N';
00184 if(m_CUTG ==1){ DECLARESTR.cutg='Y';} else DECLARESTR.cutg='N';
00185 CHANNEL.ich = m_Ich;
00186 BEAMENERGY.ebeam=m_Ebeam;
00187 EXPCUTS.thmin=m_Thmin;
00188 EXPCUTS.thmax=m_Thmax;
00189 EXPCUTS.emin =m_Emin;
00190 EXPCUTS.zmax=m_Zmax;
00191 SWITCHARUN.iarun=m_Iarun;
00192 RESONANCES.ires =m_Ires;
00193 SWITCH.on =m_on;
00194 EXPCUTS.egmin=m_Egmin;
00195 EXPCUTS.thgmin=m_Thgmin;
00196 EXPCUTS.thgmax=m_Thgmax;
00197
00198
00199
00200 getMaxEvent();
00201 std::cout<<"m_evtMax = "<<m_evtMax<<std::endl;
00202 DECLAREINT.seed=m_Int;
00203 BABAYAGA(m_evtMax);
00204
00205 return StatusCode::SUCCESS;
00206 }
00207
00208 static int evtgen=1;
00209 StatusCode Babayaga::execute()
00210 {
00211 MsgStream log(messageService(), name());
00212
00213
00214
00215 int npart = 0;
00216 int pid1,pid2,pst1,pst2;
00217 if(m_Ich==1) {pid1=11; pid2=-11; pst1=1;pst2=1;}
00218 if(m_Ich==2) {pid1=13; pid2=-13; pst1=1;pst2=1;}
00219 if(m_Ich==3) {pid1=22; pid2= 22; pst1=1;pst2=1;}
00220 if(m_Ich==4) {pid1=-211; pid2= 211; pst1=1;pst2=1;}
00221
00222
00223 GenEvent* evt = new GenEvent(1,1);
00224
00225 GenVertex* prod_vtx = new GenVertex();
00226
00227 evt->add_vertex( prod_vtx );
00228
00229
00230 GenParticle* p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p1[1][evtgen], MOMSET.p1[2][evtgen],
00231 MOMSET.p1[3][evtgen], MOMSET.p1[0][evtgen]),
00232 11, 3);
00233 p->suggest_barcode( ++npart );
00234 prod_vtx->add_particle_in(p);
00235
00236
00237
00238 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q1[1][evtgen], MOMSET.q1[2][evtgen],
00239 MOMSET.q1[3][evtgen], MOMSET.q1[4][evtgen]),
00240 -11, 3);
00241
00242 p->suggest_barcode( ++npart );
00243 prod_vtx->add_particle_in(p);
00244
00245
00246 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p2[1][evtgen], MOMSET.p2[2][evtgen],
00247 MOMSET.p2[3][evtgen], MOMSET.p2[0][evtgen]),
00248 pid1,pst1);
00249 p->suggest_barcode( ++npart );
00250 prod_vtx->add_particle_out(p);
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q2[1][evtgen], MOMSET.q2[2][evtgen],
00261 MOMSET.q2[3][evtgen], MOMSET.q2[0][evtgen]),
00262 pid2, pst2);
00263 p->suggest_barcode( ++npart );
00264 prod_vtx->add_particle_out(p);
00265
00266
00267 int iphot=0;
00268
00269 for (iphot=0; iphot<ISRPHOTONS.ncqph[evtgen-1]; iphot++)
00270 {
00271
00272
00273
00274
00275 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.phot[1][iphot][evtgen], MOMSET.phot[2][iphot][evtgen],
00276 MOMSET.phot[3][iphot][evtgen], MOMSET.phot[0][iphot][evtgen]),
00277 22, 1);
00278 p->suggest_barcode( ++npart );
00279 prod_vtx->add_particle_out(p);
00280
00281 }
00282
00283
00284 evtgen++;
00285
00286 if( log.level() < MSG::INFO )
00287 {
00288 evt->print();
00289 }
00290
00291
00292 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00293 if (anMcCol!=0)
00294 {
00295
00296 MsgStream log(messageService(), name());
00297 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00298 McGenEvent* mcEvent = new McGenEvent(evt);
00299 anMcCol->push_back(mcEvent);
00300 }
00301 else
00302 {
00303
00304 McGenEventCol *mcColl = new McGenEventCol;
00305 McGenEvent* mcEvent = new McGenEvent(evt);
00306 mcColl->push_back(mcEvent);
00307 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00308 if (sc != StatusCode::SUCCESS)
00309 {
00310 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00311 delete mcColl;
00312 delete evt;
00313 delete mcEvent;
00314 return StatusCode::FAILURE;
00315 }
00316 else
00317 {
00318
00319 }
00320 }
00321
00322 return StatusCode::SUCCESS;
00323 }
00324
00325 StatusCode Babayaga::finalize()
00326 {
00327 MsgStream log(messageService(), name());
00328 char delcmd[300];
00329 strcpy(delcmd,"cat ");
00330 strcat(delcmd,"CrossSection.txt");
00331 system(delcmd);
00332
00333 std::cout<< "BABAYAGA finalized" << endl;
00334 return StatusCode::SUCCESS;
00335 }
00336
00337 StatusCode Babayaga::getMaxEvent() {
00338 IProperty* appPropMgr=0;
00339 StatusCode status =
00340 serviceLocator()->getService("ApplicationMgr", IProperty::interfaceID(),
00341 reinterpret_cast<IInterface*&>( appPropMgr ));
00342 if( status.isFailure() ) return status;
00343
00344 IntegerProperty evtMax("EvtMax",0);
00345 status = appPropMgr->getProperty( &evtMax );
00346 if (status.isFailure()) return status;
00347
00348 m_evtMax = evtMax.value();
00349 return status;
00350 }
00351