BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
BabayagaNLO.cxx
Go to the documentation of this file.
1//#include "HepMC/HEPEVT_Wrapper.h"
2#include "HepMC/IO_HEPEVT.h"
3//#include "HepMC/IO_Ascii.h"
4#include "HepMC/GenEvent.h"
5
6#include "GaudiKernel/MsgStream.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/AlgFactory.h"
9#include "GaudiKernel/DataSvc.h"
10#include "GaudiKernel/SmartDataPtr.h"
11
16#include "CLHEP/Vector/LorentzVector.h"
17
18#include <iostream>
19#include "boost/filesystem/path.hpp"
20
21#define fs boost::filesystem
22using std::cout;
23using std::endl;
24using namespace HepMC;
25
26
27extern "C" {
28
29 extern void bossinterface_ (int *xpari, double *xpard);
30 extern void init_babayaga_ ();
31 extern void generate_event_(bool &use_unweighted);
32 extern void print_statistics_ ();
33 extern void set_ecms_spread_ (double &ecms, double &bspread);
34
35 // common/momentainitial/pin1,pin2
36 extern struct {
37 double pin1[4];
38 double pin2[4];
40
41// common/event_mom/p1,p2,qph
42 extern struct {
43 double p1[4],p2[4],qph[4][40];
45
46
47 //common/weights/use_weighted,sdif
48 extern struct {
49 double sdif;
52
53
54 // common/babayagainitint/in_conf_spin,nphmax,naccepted,nwhenmax,nover,istopsearch,nneg,ng
55 extern struct {
60
61 // common/vpcnskf/nskfilepath
62 extern struct {
63 char nskfilepath[400];
65
66}
67
68
69DECLARE_COMPONENT(BabayagaNLO)
70BabayagaNLO::BabayagaNLO(const std::string& name, ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator)
71{
72// enum int_pars {kChannel,kVerbose,kNsearch,kAlphaRun,kPhMode,kModelVP,kPrecision,kUseWeighted,kDarkMode};
73// enum dbl_pars {kEcm,kEspread,kCutTh0,kCutTh1,kCutAColl,kCutEmin,kCutInvMass0,kCutInvMass1,kDarkMass,kDarkWidth,kDarkVec,kDarkAx,kMaxSdif,kScaleVPerr};
74
75 declareProperty("Channel", m_ch = 0); // 0=ee, 1=mumu, 2 = gg
76 declareProperty("Verbose", m_iverbose = 0);
77 declareProperty("Nsearch", m_nsearch = 4000000); // event used for maximum search
78 declareProperty("RunningAlpha", m_arun = 0); // alpha running
79 declareProperty("VPparam", m_iteubn = 0); // 0= Jegerlehner, 1 = Teubner, 2 = Novosibirsk
80 declareProperty("PhotonNumber", m_photmode = -1); // maximum number of photon (-1 = unlimited (40))
81 declareProperty("Precision", m_precision = 0); // 0=best, 1=alpha2, 2=born
82
83 declareProperty("CMSEnergy", m_ecmsinput = 3.686);
84 declareProperty("BeamEnergySpread", m_beamspread = 0.0013);
85
86 declareProperty("ThetaMin", m_thmin = 10); //| Range of charged particles (ee,mumu)
87 declareProperty("ThetaMax", m_thmax = 170);//| or one photon pair (gg)
88 declareProperty("AcollMax", m_zmax = 180.0);// Maximum acollinearity
89 declareProperty("Emin", m_emin = 0.001); // minimum energy of charged tracks/most energetic photons
90 declareProperty("MinInvMass", m_Minv_min = 0.0); //| Range of invariant mass to be generated
91 declareProperty("MaxInvMass", m_Minv_max = -1.0);//| charged (ee,mumu) or at least one photon pair (gg)
92
93 declareProperty("MaximumVal", m_sdif_max = 1.e-18); //starting value for maximum search
94 declareProperty("ScaleVPerr", m_scale_vperr = 0.0); // add hadronic vp error scaled by this value ( vp_hadr = vp_hadr + scale*d_hvp, default: no addition)
95
96
97 declareProperty("UseWeightedEv", m_weighted = 0); //enable generation of weighted events (default: unweighted)
98
99 declareProperty("DarkMode", m_darkmode = 0); // use dark photon (A') model
100 declareProperty("DarkMass", m_DarkMass = 0.4); // A' mass
101 declareProperty("DarkWidth", m_DarkWidth = -1.0); // A' width (default: point-like)
102 declareProperty("DarkVectorCoupling", m_DarkVec = 1.e-3); //| A'
103 declareProperty("DarkAxialCoupling" , m_DarkAxial = 0); //| couplings
104
105 std::string path = getenv("BABAYAGANLOROOT");
106 path.append("/share/vpol_novosibirsk.dat");
107 declareProperty("NskFilePath" , m_nsk_filepath = path); // path to nsk VP input file
108
109
110
111
112
113}
114
116 MsgStream log(msgSvc(), name());
117 log << MSG::DEBUG << "BabayagaNLO in initialize()" << endreq;
118
119 //set Bes unified random engine
120 static const bool CREATEIFNOTTHERE(true);
121 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
122 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
123 {
124 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
125 return RndmStatus;
126 }
127 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("BabayagaNLO");
129
130
131 //start initialize
132 int xpari[10];
133 double xpard[14];
134
135 xpari[kChannel]=m_ch;
136 xpari[kVerbose]=m_iverbose;
137 xpari[kNsearch]=m_nsearch;
138 xpari[kAlphaRun]=m_arun;
139 xpari[kPhMode]=m_photmode;
140 xpari[kModelVP]=m_iteubn;
141 xpari[kPrecision]=m_precision;
142 xpari[kUseWeighted]=m_weighted;
143 xpari[kDarkMode]=m_darkmode;
144
145 xpard[kEcm]=m_ecmsinput;
146 xpard[kEspread]=m_beamspread;
147 xpard[kCutTh0]=m_thmin;
148 xpard[kCutTh1]=m_thmax;
149 xpard[kCutAColl]=m_zmax;
150 xpard[kCutEmin]=m_emin;
151 xpard[kCutInvMass0]=m_Minv_min;
152 xpard[kCutInvMass1]=m_Minv_max;
153 xpard[kDarkMass]=m_DarkMass;
154 xpard[kDarkWidth]=m_DarkWidth;
155 xpard[kDarkVec]=m_DarkVec;
156 xpard[kDarkAx]=m_DarkAxial;
157 xpard[kMaxSdif]=m_sdif_max;
158 xpard[kScaleVPerr]=m_scale_vperr;
159
160 using namespace boost::filesystem;
161 path p_nskfile (m_nsk_filepath.data());
162 //if (!exists(p_nskfile) && m_arun>0 && m_iteubn==2){
163 // log <<MSG::FATAL<<"Cannot file Novosibirsk VP file"<<endreq;
164 // return StatusCode::ERROR;
165 //}
166
167 m_nsk_filepath = p_nskfile.native();
168 strcpy(vpcnskf_.nskfilepath,m_nsk_filepath.data());
169
170 bossinterface_(xpari,xpard);
171
173
174 // if (!m_weighted) {
175 // log <<MSG::INFO<< "Looking for the maximum" <<endreq;
176 // bool gen_unw = 0;
177 // while (babayagainitint_.istopsearch==0){
178 // generate_event_(gen_unw);
179 // }
180 // log <<MSG::INFO<< "Search for maximum ended"<<endreq;
181 // }
182
183 log <<MSG::DEBUG<< "Finish BabayagaNLO initialize()" <<endreq;
184 return StatusCode::SUCCESS;
185}
186
188 MsgStream log(msgSvc(), name());
189 log << MSG::DEBUG << "BabayagaNLO in execute()" << endreq;
190
191 // Fill event information
192 GenEvent* evt = new GenEvent(1,1);
193 GenVertex* prod_vtx = new GenVertex();
194 evt->add_vertex( prod_vtx );
195
196 log << MSG::DEBUG << "check point 1" << endreq;
197 bool unw = !(m_weighted);
198
199 log << MSG::DEBUG << "check point 2" << endreq;
200 generate_event_(unw);
201
202 log << MSG::DEBUG << "check point 3" << endreq;
203 int finalpidm,finalpidp;
204 if (m_ch == 0) {
205 finalpidm=11;
206 finalpidp=-11;
207 } else if (m_ch == 1) {
208 finalpidm=13;
209 finalpidp=-13;
210 } else if (m_ch == 2) {
211 finalpidm=22;
212 finalpidp=22;
213 } else {
214 finalpidm=11;
215 finalpidp=-11;
216 }
217
218 double *p1 = momentainitial_.pin1;
219 double *p2 = momentainitial_.pin2;
220 double *p3 = event_mom_.p1;
221 double *p4 = event_mom_.p2;
222
223
224 //incoming particle -
225 GenParticle *p = new GenParticle(CLHEP::HepLorentzVector(p1[1],p1[2],-p1[3],p1[0]), 11, 1);
226 p->suggest_barcode(3);
227 prod_vtx->add_particle_in(p);
228
229 //incoming particle +
230 p = new GenParticle(CLHEP::HepLorentzVector(p2[1],p2[2],-p2[3],p2[0]), -11, 1);
231 p->suggest_barcode(4);
232 prod_vtx->add_particle_in(p);
233
234 log << MSG::DEBUG << "check point 4" << endreq;
235 //outgoing particle -
236 p = new GenParticle(CLHEP::HepLorentzVector(p3[1],p3[2],-p3[3],p3[0]), finalpidm, 1);
237 p->suggest_barcode(3);
238 prod_vtx->add_particle_out(p);
239
240 //outgoing particle +
241 p = new GenParticle(CLHEP::HepLorentzVector(p4[1],p4[2],-p4[3],p4[0]), finalpidp, 1);
242 p->suggest_barcode(4);
243 prod_vtx->add_particle_out(p);
244
245 log << MSG::DEBUG << "check point 5" << endreq;
246 //photons
247 int id_cntr = 5;
248 for (int i=0;i<babayagainitint_.ng;i++) {
249 double px = event_mom_.qph[1][i];
250 double py = event_mom_.qph[2][i];
251 double pz = event_mom_.qph[3][i];
252 double eph = event_mom_.qph[0][i];
253 p = new GenParticle(CLHEP::HepLorentzVector(px,py,pz,eph), 22, 1);
254 p->suggest_barcode(id_cntr);
255 prod_vtx->add_particle_out(p);
256 id_cntr++;
257 }
258
259 log << MSG::DEBUG << "check point 6" << endreq;
260 if (!unw){ //Add dummy particle with diffential cross section as mass
261 p = new GenParticle(CLHEP::HepLorentzVector(0,0,0,weights_.sdif), 99, 1);
262 p->suggest_barcode(id_cntr);
263 prod_vtx->add_particle_out(p);
264 id_cntr++;
265 }
266
267 if( log.level() <= MSG::DEBUG ){
268 evt->print();
269 }
270
271 // Check if the McCollection already exists
272 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
273 if (anMcCol!=0){
274 // Add event to existing collection
275 log<<MSG::WARNING<<"add event"<<endreq;
276 MsgStream log(messageService(), name());
277 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
278 McGenEvent* mcEvent = new McGenEvent(evt);
279 anMcCol->push_back(mcEvent);
280 } else {
281 // Create Collection and add to the transient store
282 log<<MSG::WARNING<<"create collection"<<endreq;
283 McGenEventCol *mcColl = new McGenEventCol;
284 McGenEvent* mcEvent = new McGenEvent(evt);
285 mcColl->push_back(mcEvent);
286 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
287 if (sc != StatusCode::SUCCESS) {
288 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
289 delete mcColl;
290 delete evt;
291 delete mcEvent;
292 return StatusCode::FAILURE;
293 } else {
294 log << MSG::INFO << "McGenEventCol created!" << endreq;
295 }
296 }
297
298
299 log<<MSG::DEBUG<< "before execute() return"<<endreq;
300 return StatusCode::SUCCESS;
301
302}
303
305 MsgStream log(msgSvc(), name());
306 log << MSG::INFO << "BabayagaNLO in finalize()" << endreq;
307
309 return StatusCode::SUCCESS;
310}
311
312
int ng
int in_conf_spin
void init_babayaga_()
double pin1[4]
int nover
void set_ecms_spread_(double &ecms, double &bspread)
double sdif
double pin2[4]
bool use_unweighted
struct @7 momentainitial_
long naccepted
double qph[4][40]
void print_statistics_()
double p1[4]
int nwhenmax
struct @11 vpcnskf_
struct @9 weights_
void generate_event_(bool &use_unweighted)
char nskfilepath[400]
struct @8 event_mom_
int istopsearch
double p2[4]
int nphmax
void bossinterface_(int *xpari, double *xpard)
struct @10 babayagainitint_
int nneg
*********Class see also m_nmax DOUBLE PRECISION m_emin
Definition KarFin.h:12
ObjectVector< McGenEvent > McGenEventCol
Definition McGenEvent.h:39
IMessageSvc * msgSvc()
static void setRandomEngine(CLHEP::HepRandomEngine *randomEngine)
StatusCode finalize()
StatusCode execute()
StatusCode initialize()
virtual CLHEP::HepRandomEngine * GetEngine(const std::string &StreamName)=0
Interface to the CLHEP engine.
const double ecms
Definition inclkstar.cxx:42
double double double * p4
Definition qcdloop1.h:77
double double * p3
Definition qcdloop1.h:76