#include <Ekhara.h>
Public Member Functions | |
Ekhara (const std::string &name, ISvcLocator *pSvcLocator) | |
StatusCode | initialize () |
StatusCode | execute () |
StatusCode | finalize () |
Private Attributes | |
IBesRndmGenSvc * | p_BesRndmGenSvc |
int | m_ngUpperBound |
int | m_initSeed |
int | m_channel_id |
int | m_sw_2pi |
int | m_sw_1pi |
int | m_sw_silent |
int | m_sw |
int | m_piggFFsw |
double | m_E |
double | m_Es |
double | cut_th1min |
double | cut_th1max |
double | cut_E1min |
double | cut_E1max |
double | cut_th2min |
double | cut_th2max |
double | cut_E2min |
double | cut_E2max |
double | cut_Epionmin |
double | cut_Epionmax |
double | cut_thpionmin |
double | cut_thpionmax |
IHistogram1D * | hMCPosiMom |
IHistogram1D * | hMCPosiThe |
IHistogram1D * | hMCPosiPhi |
IHistogram1D * | hMCElecMom |
IHistogram1D * | hMCElecThe |
IHistogram1D * | hMCElecPhi |
IHistogram1D * | hMCEtaPMom |
IHistogram1D * | hMCEtaPThe |
IHistogram1D * | hMCEtaPPhi |
Definition at line 19 of file Ekhara.h.
Ekhara::Ekhara | ( | const std::string & | name, | |
ISvcLocator * | pSvcLocator | |||
) |
Definition at line 83 of file Ekhara.cxx.
References cut_E1max, cut_E1min, cut_E2max, cut_E2min, cut_Epionmax, cut_Epionmin, cut_th1max, cut_th1min, cut_th2max, cut_th2min, cut_thpionmax, cut_thpionmin, m_channel_id, m_E, m_initSeed, m_ngUpperBound, m_piggFFsw, m_sw, m_sw_1pi, m_sw_2pi, and m_sw_silent.
00083 : Algorithm(name, pSvcLocator) 00084 { 00085 declareProperty( "Ecm", m_E = 1.02 ); // CMS-energy (GeV) 00086 declareProperty( "InitialSeed", m_initSeed = 100); 00087 declareProperty( "InitialEvents", m_ngUpperBound = 50000 ); // # of events to determine the maximum 00088 declareProperty( "Channel", m_channel_id = 2 ); // pi+pi-(1),pi0(2),, eta(3), eta' (4) 00089 declareProperty( "Thetaminpositron", cut_th1min = 0.00 ); //minimum angle e+ 00090 declareProperty( "Thetamaxpositron", cut_th1max = 180.00 ); //maximum angle e+ 00091 declareProperty( "Energyminpositron", cut_E1min = 0.00 ); //minimum Energy e+ 00092 declareProperty( "Energymaxpositron", cut_E1max = 110.00 ); //maximum Energy e+ 00093 declareProperty( "Thetaminelectron", cut_th2min = 0.00 ); //minimum angle e- 00094 declareProperty( "Thetamaxelectron", cut_th2max = 180.00 ); //maximum angle e- 00095 declareProperty( "Energyminelectron", cut_E2min = 0.00 ); //minimum Energy e- 00096 declareProperty( "Energymaxelectron", cut_E2max = 110.00 ); //maximum Energy e- 00097 declareProperty( "Thetaminpion", cut_thpionmin = 0.00 ); //minimum angle pi0 00098 declareProperty( "Thetamaxpion", cut_thpionmax = 180.00 ); //maximum angle pi0 00099 declareProperty( "Eminpion", cut_Epionmin = 0.00 ); //minimum Energy pi0 00100 declareProperty( "Emaxpion", cut_Epionmax = 100.00 ); //maximum Energy pi0 00101 declareProperty( "sw_silent", m_sw_silent = 1 ); //suppress output to stdout 00102 declareProperty( "sw_2pi", m_sw_2pi = 2 ); // s (1); t (2); s+t (3); s = signal 00103 declareProperty( "sw_1pi", m_sw_1pi = 2 ); // s (1); t (2); s+t (3); s = signal 00104 declareProperty( "sw", m_sw = 2 ); 00105 declareProperty( "Formfactor", m_piggFFsw = 5 ); 00106 // 00107 }
StatusCode Ekhara::execute | ( | ) |
Definition at line 199 of file Ekhara.cxx.
References EKHARA, calibUtil::ERROR, GET_MOMENTUM, hMCElecMom, hMCElecPhi, hMCElecThe, hMCEtaPMom, hMCEtaPPhi, hMCEtaPThe, hMCPosiMom, hMCPosiPhi, hMCPosiThe, genRecEmupikp::i, Bes_Common::INFO, m_channel_id, and msgSvc().
00199 { 00200 MsgStream log(msgSvc(), name()); 00201 log << MSG::INFO << "Ekhara in execute()" << endreq; 00202 00203 EKHARA(0); 00204 00205 //std::ofstream file1; 00206 //file1.open("file1.txt"); 00207 00208 00209 double p1[4],p2[4],q1[4],q2[4]; 00210 double pi1[4],pi2[4],qpion[4]; 00211 00212 for (int i=0;i<4;i++) { 00213 p1[i]=0.0; 00214 p2[i]=0.0; 00215 q1[i]=0.0; 00216 q2[i]=0.0; 00217 pi1[i]=0.0; 00218 pi2[i]=0.0; 00219 qpion[i]=0.0; 00220 } 00221 00222 GET_MOMENTUM(p1,p2,q1,q2,pi1,pi2,qpion); 00223 00224 double PosiMom = sqrt(q1[1]*q1[1]+q1[2]*q1[2]+q1[3]*q1[3]); 00225 double PosiThe = acos(q1[3]/PosiMom); 00226 double PosiPhi = atan2(q1[2],q1[1]); 00227 00228 double ElecMom = sqrt(q2[1]*q2[1]+q2[2]*q2[2]+q2[3]*q2[3]); 00229 double ElecThe = acos(q2[3]/ElecMom); 00230 double ElecPhi = atan2(q2[2],q2[1]); 00231 00232 double EtaPMom = sqrt(qpion[1]*qpion[1]+qpion[2]*qpion[2]+qpion[3]*qpion[3]); 00233 double EtaPThe = acos(qpion[3]/EtaPMom); 00234 double EtaPPhi = atan2(qpion[2],qpion[1]); 00235 00236 hMCPosiMom->fill(PosiMom); 00237 hMCElecMom->fill(ElecMom); 00238 hMCEtaPMom->fill(EtaPMom); 00239 00240 hMCPosiPhi->fill(PosiPhi*TMath::RadToDeg()); 00241 hMCElecPhi->fill(ElecPhi*TMath::RadToDeg()); 00242 hMCEtaPPhi->fill(EtaPPhi*TMath::RadToDeg()); 00243 00244 hMCPosiThe->fill(PosiThe*TMath::RadToDeg()); 00245 hMCElecThe->fill(ElecThe*TMath::RadToDeg()); 00246 hMCEtaPThe->fill(EtaPThe*TMath::RadToDeg()); 00247 00248 00249 00250 // Fill event information 00251 GenEvent* evt = new GenEvent(1,1); 00252 00253 GenVertex* prod_vtx = new GenVertex(); 00254 evt->add_vertex( prod_vtx ); 00255 00256 // incoming beam e+ 00257 GenParticle* p = new GenParticle( HepLorentzVector(p1[1],p1[2],p1[3],p1[0]), -11, 3); 00258 p->suggest_barcode(1 ); 00259 prod_vtx->add_particle_in(p); 00260 00261 // incoming beam e- 00262 p = new GenParticle(HepLorentzVector(p2[1],p2[2],p2[3],p2[0]), 11, 3); 00263 p->suggest_barcode( 2 ); 00264 prod_vtx->add_particle_in(p); 00265 00266 00267 // outcoming beam e+ 00268 p = new GenParticle(HepLorentzVector(q1[1],q1[2],q1[3],q1[0]), -11, 1); 00269 p->suggest_barcode(3 ); 00270 prod_vtx->add_particle_out(p); 00271 00272 // outcoming beam e- 00273 p = new GenParticle( HepLorentzVector(q2[1],q2[2],q2[3],q2[0] ), 11, 1); 00274 p->suggest_barcode( 4 ); 00275 prod_vtx->add_particle_out(p); 00276 00277 if (m_channel_id == 2){ //e+e- pi0 00278 p = new GenParticle( HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 111, 1); 00279 p->suggest_barcode( 5 ); 00280 prod_vtx->add_particle_out(p); 00281 } 00282 else if (m_channel_id == 3){ //e+e- eta 00283 p = new GenParticle( HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 221, 1); 00284 p->suggest_barcode( 5 ); 00285 prod_vtx->add_particle_out(p); 00286 } 00287 else if (m_channel_id == 4){ //e+e- eta' 00288 p = new GenParticle( HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 331, 1); 00289 p->suggest_barcode( 5 ); 00290 prod_vtx->add_particle_out(p); 00291 } 00292 else if (m_channel_id == 1) { // e+e-pi+ pi- 00293 // pi+ 00294 p = new GenParticle( HepLorentzVector( pi1[1],pi1[2],pi1[3],pi1[0]),211, 1); 00295 p->suggest_barcode( 5 ); 00296 prod_vtx->add_particle_out(p); 00297 // pi- 00298 p = new GenParticle( HepLorentzVector(pi2[1],pi2[2],pi2[3],pi2[0]), -211, 1); 00299 p->suggest_barcode( 6 ); 00300 prod_vtx->add_particle_out(p); 00301 } 00302 00303 if( log.level() <= MSG::INFO ) 00304 { 00305 00306 evt->print(); 00307 00308 // file1<<qpion[1]<<endl; 00309 00310 } 00311 00312 // Check if the McCollection already exists 00313 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen"); 00314 if (anMcCol!=0) 00315 { 00316 // Add event to existing collection 00317 MsgStream log(messageService(), name()); 00318 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq; 00319 McGenEvent* mcEvent = new McGenEvent(evt); 00320 anMcCol->push_back(mcEvent); 00321 } 00322 else 00323 { 00324 // Create Collection and add to the transient store 00325 McGenEventCol *mcColl = new McGenEventCol; 00326 McGenEvent* mcEvent = new McGenEvent(evt); 00327 mcColl->push_back(mcEvent); 00328 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl); 00329 if (sc != StatusCode::SUCCESS) 00330 { 00331 log << MSG::ERROR << "Could not register McGenEvent" << endreq; 00332 delete mcColl; 00333 delete evt; 00334 delete mcEvent; 00335 return StatusCode::FAILURE; 00336 } 00337 else 00338 { 00339 log << MSG::INFO << "McGenEventCol created " << endreq; 00340 } 00341 //file1.close(); 00342 } 00343 00344 log<<MSG::INFO<< "before execute() return"<<endreq; 00345 return StatusCode::SUCCESS; 00346 00347 }
StatusCode Ekhara::finalize | ( | ) |
Definition at line 349 of file Ekhara.cxx.
References DIAGNOSE, EKHARA, Bes_Common::INFO, and msgSvc().
00349 { 00350 MsgStream log(msgSvc(), name()); 00351 log << MSG::INFO << "Ekhara in finalize()" << endreq; 00352 log << MSG::INFO << "Ekhara is terminated now successfully" << endreq; 00353 DIAGNOSE(); 00354 EKHARA(1); 00355 return StatusCode::SUCCESS; 00356 }
StatusCode Ekhara::initialize | ( | ) |
Definition at line 109 of file Ekhara.cxx.
References BOSS_INIT_READ, CHANNELSEL, cut_E1max, cut_E1min, cut_E2max, cut_E2min, cut_Epionmax, cut_Epionmin, cut_th1max, cut_th1min, cut_th2max, cut_th2min, cut_thpionmax, cut_thpionmin, Bes_Common::DEBUG, DIAGNOSE, EKHARA, calibUtil::ERROR, IBesRndmGenSvc::GetEngine(), histoSvc(), hMCElecMom, hMCElecPhi, hMCElecThe, hMCEtaPMom, hMCEtaPPhi, hMCEtaPThe, hMCPosiMom, hMCPosiPhi, hMCPosiThe, Bes_Common::INFO, m_channel_id, m_E, m_ngUpperBound, m_piggFFsw, m_sw_1pi, m_sw_2pi, msgSvc(), p_BesRndmGenSvc, PIONFFSW, EkharaRandom::setRandomEngine(), and SWDIAG.
00109 { 00110 MsgStream log(msgSvc(), name()); 00111 log << MSG::INFO << "Ekhara in initialize()" << endreq; 00112 00113 00114 hMCPosiMom = histoSvc()->book("MonteCarlo",1,"PosiMom",250,0.,5.); 00115 hMCPosiThe = histoSvc()->book("MonteCarlo",2,"PosiThe",45,0.,180.); 00116 hMCPosiPhi = histoSvc()->book("MonteCarlo",3,"PosiPhi",90,-180.,180.); 00117 hMCElecMom = histoSvc()->book("MonteCarlo",4,"ElecMom",250,0.,5.); 00118 hMCElecThe = histoSvc()->book("MonteCarlo",5,"ElecThe",45,0.,180.); 00119 hMCElecPhi = histoSvc()->book("MonteCarlo",6,"ElecPhi",90,-180.,180.); 00120 hMCEtaPMom = histoSvc()->book("MonteCarlo",7,"EtaPMom",250,0.,5.); 00121 hMCEtaPThe = histoSvc()->book("MonteCarlo",8,"EtaPThe",45,0.,180.); 00122 hMCEtaPPhi = histoSvc()->book("MonteCarlo",9,"EtaPPhi",90,-180.,180.); 00123 00124 00125 //set Bes unified random engine 00126 static const bool CREATEIFNOTTHERE(true); 00127 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE); 00128 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc) 00129 { 00130 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq; 00131 return RndmStatus; 00132 } 00133 HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("Ekhara"); 00134 EkharaRandom::setRandomEngine(engine); 00135 00136 CHANNELSEL.channel_id = m_channel_id; 00137 SWDIAG.sw_2pi=m_sw_2pi; 00138 SWDIAG.sw_1pi=m_sw_1pi; 00139 PIONFFSW.piggFFsw=m_piggFFsw; 00140 /* 00141 if(m_channel_id==1){ 00142 m_sw=m_sw_2pi;} 00143 else if(m_channel_id==2||m_channel_id==3||m_channel_id==4){ 00144 m_sw=m_sw_1pi;} 00145 */ 00146 00147 // --- print run data ------------------------------ 00148 cout << "-------------------------------------------------------------" << endl; 00149 if (CHANNELSEL.channel_id == 2) 00150 cout << " EKHARA 2.1 : e^+ e^- -> e^+ e^- pi^0" << endl; 00151 else if (CHANNELSEL.channel_id == 1) 00152 cout << " EKHARA 2.1: e^+ e^- -> e^+ e^- pi^+ pi^-" << endl; 00153 else if (CHANNELSEL.channel_id == 3) 00154 cout << " EKHARA 2.1: e^+ e^- -> e^+ e^- eta" << endl; 00155 else if (CHANNELSEL.channel_id == 4) 00156 cout << " EKHARA 2.1: e^+ e^- -> e^+ e^- eta'" << endl; 00157 00158 00159 if(m_channel_id ==2||m_channel_id ==3||m_channel_id ==4) 00160 printf(" %s %f,%f\n","angular cuts on e+ = ",cut_th1min,cut_th1max); 00161 printf(" %s %f,%f\n","angular cuts on e- = ",cut_th2min,cut_th2max); 00162 printf(" %s %f,%f\n","Energy cuts on e+ = ",cut_E1min,cut_E1max); 00163 printf(" %s %f,%f\n","Energy cuts on e- = ",cut_E2min,cut_E2max); 00164 //--------------------------------- 00165 00166 double xpar[100]; 00167 xpar[0]=m_E; 00168 00169 xpar[1]=m_ngUpperBound; 00170 xpar[2]=m_channel_id; 00171 xpar[3]=cut_th1min; 00172 xpar[4]=cut_th1max; 00173 xpar[5]=cut_E1min; 00174 xpar[6]=cut_E1max; 00175 xpar[7]=cut_th2min; 00176 xpar[8]=cut_th2max; 00177 xpar[9]=cut_E2min; 00178 xpar[10]=cut_E2max; 00179 xpar[11]=cut_thpionmin; 00180 xpar[12]=cut_thpionmax; 00181 xpar[13]=cut_Epionmin; 00182 xpar[14]=cut_Epionmax; 00183 00184 log << MSG::DEBUG << "Options filled in array" << endreq; 00185 BOSS_INIT_READ(xpar); 00186 log << MSG::DEBUG << "Options read by FORTRAN routine" << endreq; 00187 DIAGNOSE(); 00188 log << MSG::DEBUG << "FORTRAN diagnose passed" << endreq; 00189 EKHARA(-1); 00190 log << MSG::DEBUG << "Ekhara Fortran routines initialized" << endreq; 00191 //EKHARA_SET_ONE_EVENT_MODE(); 00192 //log << MSG::DEBUG << "Ekhara set to single event mode" << endreq; 00193 // EKHARA_SET_SILENT(); 00194 //log << MSG::DEBUG << "Ekhara set to silent operation" << endreq; 00195 00196 return StatusCode::SUCCESS; 00197 }
double Ekhara::cut_E1max [private] |
double Ekhara::cut_E1min [private] |
double Ekhara::cut_E2max [private] |
double Ekhara::cut_E2min [private] |
double Ekhara::cut_Epionmax [private] |
double Ekhara::cut_Epionmin [private] |
double Ekhara::cut_th1max [private] |
double Ekhara::cut_th1min [private] |
double Ekhara::cut_th2max [private] |
double Ekhara::cut_th2min [private] |
double Ekhara::cut_thpionmax [private] |
double Ekhara::cut_thpionmin [private] |
IHistogram1D* Ekhara::hMCElecMom [private] |
IHistogram1D* Ekhara::hMCElecPhi [private] |
IHistogram1D* Ekhara::hMCElecThe [private] |
IHistogram1D* Ekhara::hMCEtaPMom [private] |
IHistogram1D* Ekhara::hMCEtaPPhi [private] |
IHistogram1D* Ekhara::hMCEtaPThe [private] |
IHistogram1D* Ekhara::hMCPosiMom [private] |
IHistogram1D* Ekhara::hMCPosiPhi [private] |
IHistogram1D* Ekhara::hMCPosiThe [private] |
int Ekhara::m_channel_id [private] |
double Ekhara::m_E [private] |
double Ekhara::m_Es [private] |
int Ekhara::m_initSeed [private] |
int Ekhara::m_ngUpperBound [private] |
int Ekhara::m_piggFFsw [private] |
int Ekhara::m_sw [private] |
int Ekhara::m_sw_1pi [private] |
int Ekhara::m_sw_2pi [private] |
int Ekhara::m_sw_silent [private] |
IBesRndmGenSvc* Ekhara::p_BesRndmGenSvc [private] |