BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
EeToeeV/EeToeeV-00-00-01/src/EeToeeV/Eepipi.cxx
Go to the documentation of this file.
1//*****************************************************************************
2//
3// Eepipi.cxx
4//
5// Algorithm runs e+e- ->e+e- rho0, rho0->pi+pi- precess
6//
7// July 2016-4-29 Rong-Gang Ping to create package for BES3
8// The original fortran code is generated with FDC, consult Prof. Wang Jianxiong
9//*****************************************************************************
10
11// Header for this module:-
12#include "../Eepipi/Eepipi.h"
13#include "../Eepipi/EepipiRandom.h"
14#include "BesKernel/IBesRndmGenSvc.h"
15
16// Framework Related Headers:-
17#include "HepMC/GenEvent.h"
18#include "HepMC/GenVertex.h"
19#include "HepMC/GenParticle.h"
20
21#include "GaudiKernel/MsgStream.h"
22#include "GaudiKernel/ISvcLocator.h"
23#include "GaudiKernel/AlgFactory.h"
24#include "GaudiKernel/DataSvc.h"
25#include "GaudiKernel/SmartDataPtr.h"
26
27#include "GeneratorObject/McGenEvent.h"
28#include "CLHEP/Vector/LorentzVector.h"
29#include "cfortran/cfortran.h"
30
31#include <stdlib.h>
32#include <time.h>
33
34 typedef struct {
35 double q1[4][160001]; //e+ beam
36 double p1[4][160001]; //e- beam
37 double q2[4][160001]; // e-
38 double p2[4][160001]; // e+
39 double q3[4][160001]; // pi+
40 double p3[4][160001]; // pi-
41 } MOMSET_DEF;
42//MOMSET_DEF MOMSET;
43#define MOMSET COMMON_BLOCK(MOMSET_DEF, momset)
45
46//common/beamEnergy/ebeam
47typedef struct {
48 double ecms;
50#define BEAMENERGY COMMON_BLOCK(BEAMENERGY_DEF, beamenergy)
52//BEAMENERGY_DEF,BEAMENERGY
53
54
55//--//Unify Eepipi random engine with Bes random servive
56extern "C" {
57 extern float flat_();
58}
59
60float flat_(){
61 float rr;
62 double rd=EepipiRandom::random();
63// std::cout<<"EepipiRandom::random()="<<rd<<endl;
64 rr=(float)rd;
65 return (float)EepipiRandom::random();
66 // return rr;
67}
68extern "C" {
69 extern void intxs_();
70}
71
72
73PROTOCCALLSFSUB1(GEVENT,gevent,INT)
74#define GEVENT(NEVENTS) CCALLSFSUB1(GEVENT,gevent,INT,NEVENTS)
75
76
77Eepipi::Eepipi(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
78{
79 declareProperty("Ecms", m_Ecms = 3.65); // Ecm = sqrt(s) [GeV]
80}
81
82StatusCode Eepipi::initialize(){
83
84 MsgStream log(messageService(), name());
85
86 log << MSG::INFO << "Eepipi initialize" << endreq;
87
88 //set Bes unified random engine
89 static const bool CREATEIFNOTTHERE(true);
90 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
91 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
92 {
93 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
94 return RndmStatus;
95 }
96 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("Eepipi");
97 std::cout<<"==============================="<<engine<<endl;
99 // *****************
100
101 BEAMENERGY.ecms=m_Ecms;
102 // std::cout<<"m_Ires= "<<m_Ires<<endl;
103
104 getMaxEvent();
105 std::cout<<"m_evtMax = "<<m_evtMax<<std::endl;
106 intxs_();
107 GEVENT(m_evtMax);
108
109 return StatusCode::SUCCESS;
110}
111
112static int evtgen=1;
113StatusCode Eepipi::execute()
114{
115 MsgStream log(messageService(), name());
116
117
118 int npart = 0;
119 int pid1,pid2,pid3,pid4,pst1,pst2;
120 pid1 = 11;
121 pid2 =-11;
122 pid3 = 211;
123 pid4 =-211;
124 pst1 = 1;
125 pst2 = 1;
126
127
128
129 // Fill event information
130 GenEvent* evt = new GenEvent(1,1);
131
132 GenVertex* prod_vtx = new GenVertex();
133// prod_vtx->add_particle_out( p );
134 evt->add_vertex( prod_vtx );
135
136 // incoming beam e- HepLorentzVector(px,py,pz,energy) required!
137 GenParticle* p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p1[1][evtgen], MOMSET.p1[2][evtgen],
138 MOMSET.p1[3][evtgen], MOMSET.p1[0][evtgen]),
139 11, 3);
140 p->suggest_barcode( ++npart );
141 prod_vtx->add_particle_in(p);
142
143// std::cout<<"incoming beam e+"<<endl;
144 // incoming beam e+, e+ moving along the z-direction
145 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q1[1][evtgen], MOMSET.q1[2][evtgen],
146 MOMSET.q1[3][evtgen], MOMSET.q1[4][evtgen]),
147 -11, 3);
148
149 p->suggest_barcode( ++npart );
150 prod_vtx->add_particle_in(p);
151
152 // scattered lepton +
153 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p2[1][evtgen], MOMSET.p2[2][evtgen],
154 MOMSET.p2[3][evtgen], MOMSET.p2[0][evtgen]),
155 pid1,pst1);
156 p->suggest_barcode( ++npart );
157 prod_vtx->add_particle_out(p);
158
159 // scattered lepton -
160 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q2[1][evtgen], MOMSET.q2[2][evtgen],
161 MOMSET.q2[3][evtgen], MOMSET.q2[0][evtgen]),
162 pid2, pst2);
163 p->suggest_barcode( ++npart );
164 prod_vtx->add_particle_out(p);
165
166 // outgoing pi+
167 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p3[1][evtgen], MOMSET.p3[2][evtgen],
168 MOMSET.p3[3][evtgen], MOMSET.p3[0][evtgen]),
169 pid3,pst1);
170 p->suggest_barcode( ++npart );
171 prod_vtx->add_particle_out(p);
172
173 //outgoing pi-
174 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q3[1][evtgen], MOMSET.q3[2][evtgen],
175 MOMSET.q3[3][evtgen], MOMSET.q3[0][evtgen]),
176 pid4, pst2);
177 p->suggest_barcode( ++npart );
178 prod_vtx->add_particle_out(p);
179
180
181
182 evtgen++;
183
184 if( log.level() < MSG::INFO )
185 {
186 evt->print();
187 }
188
189 // Check if the McCollection already exists
190 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
191 if (anMcCol!=0)
192 {
193 // Add event to existing collection
194 MsgStream log(messageService(), name());
195 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
196 McGenEvent* mcEvent = new McGenEvent(evt);
197 anMcCol->push_back(mcEvent);
198 }
199 else
200 {
201 // Create Collection and add to the transient store
202 McGenEventCol *mcColl = new McGenEventCol;
203 McGenEvent* mcEvent = new McGenEvent(evt);
204 mcColl->push_back(mcEvent);
205 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
206 if (sc != StatusCode::SUCCESS)
207 {
208 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
209 delete mcColl;
210 delete evt;
211 delete mcEvent;
212 return StatusCode::FAILURE;
213 }
214 else
215 {
216 // log << MSG::INFO << "McGenEventCol created and " << npart <<" particles stored in McGenEvent" << endreq;
217 }
218 }
219
220 return StatusCode::SUCCESS;
221}
222
223StatusCode Eepipi::finalize()
224{
225 MsgStream log(messageService(), name());
226 char delcmd[300];
227 strcpy(delcmd,"cat ");
228 strcat(delcmd,"fresult.dat");
229 system(delcmd);
230
231 std::cout<< "Eepipi finalized" << endl;
232 return StatusCode::SUCCESS;
233}
234
235StatusCode Eepipi::getMaxEvent() {
236 IProperty* appPropMgr=0;
237 StatusCode status =
238 serviceLocator()->getService("ApplicationMgr", IProperty::interfaceID(),
239 reinterpret_cast<IInterface*&>( appPropMgr ));
240 if( status.isFailure() ) return status;
241
242 IntegerProperty evtMax("EvtMax",0);
243 status = appPropMgr->getProperty( &evtMax );
244 if (status.isFailure()) return status;
245
246 m_evtMax = evtMax.value();
247 return status;
248}
249
#define MOMSET
Definition: Babayaga.cxx:45
#define BEAMENERGY
Definition: Babayaga.cxx:65
void intxs_()
#define GEVENT(NEVENTS)
void intxs_()
#define COMMON_BLOCK_DEF(DEFINITION, NAME)
#define PROTOCCALLSFSUB1(UN, LN, T1)
static void setRandomEngine(CLHEP::HepRandomEngine *randomEngine)
static double random()
Eepipi(const string &name, ISvcLocator *pSvcLocator)
virtual CLHEP::HepRandomEngine * GetEngine(const std::string &StreamName)=0
Interface to the CLHEP engine.
const double ecms
Definition: inclkstar.cxx:42