00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 #include "CosmicGenerator/CosmicGenerator.h"
00074 #include "CosmicGenerator/CosmicGun.h"
00075 #include "CosmicGenerator/CosmicEventParser.h"
00076
00077
00078 #include "GaudiKernel/MsgStream.h"
00080
00081
00082 #include "CLHEP/Vector/TwoVector.h"
00083 #include "CLHEP/Vector/ThreeVector.h"
00084 #include "CLHEP/Geometry/Normal3D.h"
00085 #include "CLHEP/Geometry/Point3D.h"
00086
00087
00089 #include "CLHEP/Units/PhysicalConstants.h"
00090
00091 #include "HepMC/GenEvent.h"
00092 #include "HepMC/GenVertex.h"
00093 #include "HepMC/GenParticle.h"
00094 #include "HepMC/Polarization.h"
00095
00096 #include "GaudiKernel/Bootstrap.h"
00097 #include "GaudiKernel/ISvcLocator.h"
00098 #include "GaudiKernel/IMessageSvc.h"
00099 #include "GaudiKernel/GaudiException.h"
00100 #include "GaudiKernel/AlgFactory.h"
00101 #include "GaudiKernel/DataSvc.h"
00102 #include "GaudiKernel/SmartDataPtr.h"
00103 #include "GaudiKernel/MsgStream.h"
00104
00105
00106 #include "GaudiKernel/IDataProviderSvc.h"
00107 #include "GaudiKernel/PropertyMgr.h"
00108
00109 #include "GaudiKernel/INTupleSvc.h"
00110 #include "GaudiKernel/NTuple.h"
00111 #include "GaudiKernel/Bootstrap.h"
00112 #include "GaudiKernel/IHistogramSvc.h"
00113
00114 #include "GeneratorObject/McGenEvent.h"
00115
00116
00120
00121 #include "BesKernel/IBesRndmGenSvc.h"
00122
00123 #include "CLHEP/Random/Ranlux64Engine.h"
00124 #include "CLHEP/Random/RandFlat.h"
00125
00126 #include <limits.h>
00127
00128 #include <cmath>
00129 #include <vector>
00130 #include <string>
00131 #include <fstream>
00132 using namespace CLHEP;
00133
00134 static inline int sign(double x) { return (x>0 ? 1: -1);}
00135
00136
00137 float PI = 3.1415927;
00138 double mass =0.1055658;
00139
00141
00142 IBesRndmGenSvc* CosmicGenerator::p_BesRndmGenSvc = 0;
00143
00144 extern "C" float cosmicrndm_(int* )
00145 {
00146
00147 HepRandomEngine* engine = CosmicGenerator::p_BesRndmGenSvc->GetEngine("PYTHIA");
00148
00149 return RandFlat::shoot(engine);
00150 }
00151
00152
00153 CosmicGenerator::CosmicGenerator(const std::string& name,
00154 ISvcLocator* pSvcLocator): Algorithm(name,pSvcLocator)
00155
00156 {
00157
00158
00159
00160
00161
00162
00163 m_GeV = 1000;
00164 m_mm = 10;
00165 m_readfile = false;
00166 m_exzCut=false ;
00167 m_tuple=0;
00168 m_tuple1=0;
00169 declareProperty("eventfile", m_infile = "NONE" );
00170 declareProperty("emin", m_emin =0.1 );
00171 declareProperty("emax", m_emax =10. );
00172 declareProperty("xvert_low", m_xlow =-110.7);
00173 declareProperty("xvert_hig", m_xhig =110.7 );
00174 declareProperty("zvert_low", m_zlow =-171.2);
00175 declareProperty("zvert_hig", m_zhig = 171.2 );
00176 declareProperty("yvert_val", m_yval = 0 );
00177 declareProperty("tmin", m_tmin =0. );
00178 declareProperty("tmax", m_tmax =24. );
00179 declareProperty("tcor", m_tcor =0. );
00180 declareProperty("IPx", m_IPx =0. );
00181 declareProperty("IPy", m_IPy =0. );
00182 declareProperty("IPz", m_IPz =0. );
00183 declareProperty("Radius", m_radius =0. );
00184 declareProperty("ExzCut", m_exzCut );
00185 declareProperty("CubeProjection", m_cubeProj = true );
00186 declareProperty("OptimizeSphere", m_sphereOpt = false );
00187
00188
00189 declareProperty("SwapYZAxis", m_swapYZAxis = false);
00190 declareProperty("ctcut", m_ctcut =0.0 );
00191
00192
00193 declareProperty("PrintEvent", m_printEvent=10);
00194 declareProperty("PrintMod", m_printMod=100);
00195 declareProperty("RMax", m_rmax = 10000000. );
00196 declareProperty("ThetaMin", m_thetamin = 0.);
00197 declareProperty("ThetaMax", m_thetamax = 1.);
00198 declareProperty("PhiMin", m_phimin = -1*PI);
00199 declareProperty("PhiMax", m_phimax = PI);
00200 declareProperty("Xposition", m_xpos = 263.5-0.0001);
00201 declareProperty("Yposition", m_ypos = 263.5-0.0001);
00202 declareProperty("Zposition", m_zpos = 287.5-0.0001);
00203
00204 declareProperty("DumpMC", m_dumpMC = 0);
00205 }
00206
00207
00208 CosmicGenerator::~CosmicGenerator()
00209
00210 {}
00211
00212
00213 StatusCode CosmicGenerator::initialize() {
00214
00215
00216
00217
00218 m_events = 0;
00219 m_tried=0;
00220 m_accepted=0;
00221 m_rejected=0;
00222 if(m_emin<0.1) {m_emin=0.1;std::cout<<"WARNING: set emin to 0.1 GeV"<<std::endl;}
00223 m_emin =m_emin *m_GeV;
00224 m_emax =m_emax *m_GeV;
00225 m_xlow =m_xlow *m_mm;
00226 m_xhig =m_xhig *m_mm;
00227 m_zlow =m_zlow *m_mm;
00228 m_zhig =m_zhig *m_mm;
00229 m_yval =m_yval *m_mm;
00230 m_xpos =m_xpos *m_mm;
00231 m_ypos =m_ypos *m_mm;
00232 m_zpos =m_zpos *m_mm;
00233 m_radius=m_radius*m_mm;
00234 m_tcor=m_tcor
00235 +sqrt((m_xpos-m_xlow)*(m_xpos-m_xlow)+(m_zpos-m_zlow)*(m_zpos-m_zlow)+(m_ypos-m_yval)*(m_ypos-m_yval))
00236 /(m_emin/sqrt(m_emin*m_emin+mass*mass*m_GeV*m_GeV));
00237
00238 ISvcLocator* svcLocator = Gaudi::svcLocator();
00239 p_msgSvc = 0;
00240 StatusCode result = svcLocator->service("MessageSvc", p_msgSvc);
00241
00242 if ( !result.isSuccess() || p_msgSvc == 0 )
00243 {
00244 std::cerr << "VGenCommand(): could not initialize MessageSvc!" << std::endl;
00245 throw GaudiException("Could not initalize MessageSvc","CosmicGenerator",StatusCode::FAILURE);
00246 }
00247
00248 MsgStream log (p_msgSvc,"CosmicGenerator");
00249 if(m_dumpMC==1) {
00250 StatusCode status;
00251 NTuplePtr nt(ntupleSvc(), "FILE1/comp");
00252 if(nt) m_tuple = nt;
00253 else {
00254 m_tuple = ntupleSvc()->book ("FILE1/comp", CLID_ColumnWiseTuple, "ks N-Tuple example");
00255 if ( m_tuple ) {
00256 status = m_tuple->addItem ("cosmic_e", m_cosmicE);
00257 status = m_tuple->addItem ("cosmic_theta", m_cosmicTheta);
00258 status = m_tuple->addItem ("cosmic_phi", m_cosmicPhi);
00259 status = m_tuple->addItem ("cosmic_charge", m_cosmicCharge);
00260 }
00261 else {
00262 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple) << endmsg;
00263 return StatusCode::FAILURE;
00264 }
00265 }
00266 NTuplePtr nt1(ntupleSvc(), "FILE1/compp");
00267 if(nt1) m_tuple1 = nt1;
00268 else {
00269 m_tuple1 = ntupleSvc()->book ("FILE1/compp", CLID_ColumnWiseTuple, "ks N-Tuple example");
00270 if ( m_tuple1 ) {
00271 status = m_tuple1->addItem ("mc_theta", mc_theta);
00272 status = m_tuple1->addItem ("mc_phi", mc_phi);
00273 status = m_tuple1->addItem ("mc_px", mc_px);
00274 status = m_tuple1->addItem ("mc_py", mc_py);
00275 status = m_tuple1->addItem ("mc_pz", mc_pz);
00276 }
00277 else {
00278 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00279 return StatusCode::FAILURE;
00280 }
00281 }
00282 }
00283 if(m_infile=="NONE")
00284
00285 {
00286
00287
00288 static const bool CREATEIFNOTTHERE(true);
00289 StatusCode RndmStatus = svcLocator->service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00290
00291 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00292 {
00293 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00294 return RndmStatus;
00295 }
00296
00297 CosmicGun* gun = CosmicGun::GetCosmicGun();
00298
00299 gun->SetEnergyRange(m_emin/m_GeV,m_emax/m_GeV);
00300 gun->SetCosCut(m_ctcut);
00301 gun->PrintLevel(m_printEvent, m_printMod);
00302 float flux_withCT = gun->InitializeGenerator();
00303
00304 log << MSG::INFO << "Initialisation cosmic gun done." << endreq;
00305 log << MSG::INFO << "Accepted diff flux after E and cos(theta) cuts = " <<
00306 flux_withCT << " /cm^2/s" << endreq;
00307 log << MSG::INFO << "Accepted total flux after E and cos(theta) cuts = " <<
00308 flux_withCT*(m_xhig-m_xlow)/m_mm*(m_zhig-m_zlow)/m_mm << " /s" << endreq;
00309
00310 std::cout<< "Accepted diff flux after E and cos(theta) cuts = " <<
00311 flux_withCT << " /cm^2/s" << std::endl;
00312 std::cout << "Accepted total flux after E and cos(theta) cuts = " <<
00313 flux_withCT*(m_xhig-m_xlow)/m_mm*(m_zhig-m_zlow)/m_mm << " /s" << std::endl;
00314
00315 }
00316 else
00317 {
00318 log << MSG::INFO << "Cosmics are read from file " << m_infile << endreq;
00319 m_ffile.open(m_infile.c_str());
00320 if(!m_ffile)
00321 {
00322 log << MSG::FATAL << "Could not open input file - stop! " << endreq;
00323 return StatusCode::FAILURE;
00324 }
00325 m_readfile = true;
00326 }
00327
00328 m_center=Hep3Vector(m_IPx, m_IPy, m_IPz);
00329 log << MSG::INFO << "Initialisation done" << endreq;
00330 return StatusCode::SUCCESS;
00331
00332 }
00333
00334 HepLorentzVector CosmicGenerator::generateVertex(void) {
00335
00336 MsgStream log(messageService(), name());
00337
00338
00339
00340
00341 HepRandomEngine* engine = CosmicGenerator::p_BesRndmGenSvc->GetEngine("PYTHIA");
00342
00343
00344 float x_val = RandFlat::shoot(engine, m_xlow, m_xhig);
00345 float z_val = RandFlat::shoot(engine, m_zlow, m_zhig);
00346
00347
00348 float t_val=0;
00349 if(m_tmin < m_tmax){
00350 t_val = RandFlat::shoot(engine, m_tmin, m_tmax);
00351 }
00352 else if(m_tmin == m_tmax){
00353 t_val = m_tmin;
00354 }
00355 else log << MSG::FATAL << " You specified m_tmin = " << m_tmin << " and m_tmax " << m_tmax << endreq;
00356 HepLorentzVector p(x_val,m_yval,z_val, t_val*c_light);
00357
00358 return p;
00359
00360 }
00361
00362
00363 StatusCode CosmicGenerator::execute() {
00364
00365 MsgStream log(messageService(), name());
00366 log << MSG::INFO << "CosmicGenerator executing" << endreq;
00367
00368
00369 ++m_events;
00370
00371 log << MSG::DEBUG << "Event #" << m_events << endreq;
00372
00373
00374
00375 m_fourPos.clear();
00376 m_fourMom.clear();
00377 m_polarization.clear();
00378 m_pdgCode.clear();
00379
00380
00381 if(m_readfile)
00382 {
00383 if(!m_ffile.eof())
00384 {
00385 CosmicEventParser evt;
00386 m_ffile >> evt;
00387
00388 std::cout << evt << std::endl;
00389
00390 double polx = 0;
00391 double poly = 0;
00392 double polz = 0;
00393 Polarization thePolarization;
00394 thePolarization.set_normal3d(HepNormal3D(polx,poly,polz));
00395 m_polarization.push_back(thePolarization);
00396
00397
00398
00399
00400 m_fourPos.push_back(evt.Vertex());
00401 m_fourMom.push_back(evt.Momentum());
00402 m_pdgCode.push_back(evt.pdgID());
00403
00404 }
00405 else
00406 {
00407 log << MSG::FATAL << "End of file reached - stop " << endreq;
00408 exit(1);
00409 return StatusCode::FAILURE;
00410 }
00411 }
00412
00413 else
00414 {
00415
00416 bool accepted=false;
00417 bool projected=false;
00418 HepLorentzVector pp;
00419 CosmicGun* gun = CosmicGun::GetCosmicGun();
00420 HepLorentzVector vert;
00421 HepLorentzVector vert_proj;
00422 while(!accepted){
00423 m_tried++;
00424 vert = generateVertex();
00425 Hep3Vector vert3(vert.x(),vert.y(),vert.z());
00426
00427 pp = gun->GenerateEvent();
00428 if(m_dumpMC==1) {
00429 m_cosmicE=pp.e()*m_GeV;
00430 m_cosmicTheta=pp.theta();
00431 m_cosmicPhi=pp.phi();
00432 m_cosmicCharge=gun->GetMuonCharge();
00433 m_tuple->write();
00434 }
00435 double theta1=pp.theta();
00436 double phi1=pp.phi();
00437 double mag1=pp.vect().mag();
00438
00439 Hep3Vector pp_corr(mag1*sin(theta1)*cos(phi1),
00440 -mag1*cos(theta1),
00441 mag1*sin(theta1)*sin(phi1));
00442 Hep3Vector direction(pp_corr.x(),pp_corr.y(), pp_corr.z());
00443 Hep3Vector center_dir=m_center-vert3;
00444
00445
00446
00447 if (m_cubeProj) {
00448 double alpha,beta;
00449 if(m_sphereOpt){
00450
00451 beta=direction.angle(center_dir);
00452 alpha=asin(m_radius/center_dir.mag());
00453 if(fabs(beta)>alpha){
00454 log<<MSG::DEBUG<<alpha<<", "<<beta<<" rejected muon due to sphere cut! "<<endreq;
00455 m_rejected++;
00456 continue;
00457 }
00458 }
00459
00460 double l_fight0,l_fight1, l_fight2;
00461 double coor_x0, coor_y0, coor_z0;
00462 double coor_x1, coor_y1, coor_z1;
00463 double coor_x2, coor_y2, coor_z2;
00464 bool isIn0=false;
00465 bool isIn1=false;
00466 bool isIn2=false;
00467 HepPoint3D vert_pro0;
00468 HepPoint3D vert_pro1;
00469 HepPoint3D vert_pro2;
00470 HepPoint3D vert_org(vert3);
00471
00472 coor_y0 = m_ypos;
00473 coor_x0 = direction.x()*(coor_y0 - vert.y())/direction.y() +vert.x();
00474 coor_z0 = direction.z()*(coor_y0 - vert.y())/direction.y() +vert.z();
00475
00476
00477
00478 if( fabs(coor_x0) <= m_xpos && fabs(coor_z0) <= m_zpos){
00479 vert_pro0=HepPoint3D (coor_x0, coor_y0, coor_z0);
00480 l_fight0=vert_org.distance(vert_pro0);
00481 isIn0=true;
00482 }
00483 else{
00484
00485 coor_z1 = sign(coor_z0-vert.z())*m_zpos;
00486 coor_x1 = direction.x()*(coor_z1 - vert.z())/direction.z() +vert.x();
00487 coor_y1 = direction.y()*(coor_z1 - vert.z())/direction.z() +vert.y();
00488
00489 vert_pro1=HepPoint3D (coor_x1, coor_y1, coor_z1);
00490 l_fight1=vert_org.distance(vert_pro1);
00491
00492
00493 coor_x2 = sign(coor_x0-vert.x())*m_xpos;
00494 coor_z2 = direction.z()*(coor_x2 - vert.x())/direction.x() +vert.z();
00495 coor_y2 = direction.y()*(coor_x2 - vert.x())/direction.x() +vert.y();
00496 vert_pro2=HepPoint3D (coor_x2, coor_y2, coor_z2);
00497 l_fight2=vert_org.distance(vert_pro2);
00498 if(l_fight1<l_fight2)
00499 isIn1=true;
00500 else
00501 isIn2=true;
00502 }
00503
00504 if(isIn0){
00505 vert_proj=HepLorentzVector (vert_pro0.x(),vert_pro0.y(),vert_pro0.z() , l_fight0);
00506
00507 projected =true;
00508 accepted = true;
00509 m_accepted++;
00510 }
00511 else if(isIn1){
00512 vert_proj=HepLorentzVector (vert_pro1.x(),vert_pro1.y(),vert_pro1.z() , l_fight1);
00513
00514 projected =true;
00515 accepted = true;
00516 m_accepted++;
00517 }
00518 else if(isIn2){
00519 vert_proj=HepLorentzVector (vert_pro2.x(),vert_pro2.y(),vert_pro2.z() , l_fight2);
00520
00521 projected =true;
00522 accepted = true;
00523 m_accepted++;
00524 } else {
00525 log<<MSG::DEBUG<<" rejected muon due to cube projection request! "<<endreq;
00526 m_rejected++;
00527 }
00528
00529
00530
00531
00532 }
00533 else accepted=true;
00534 }
00535
00536
00537
00538
00539
00540
00541 pp.setX(pp.x());
00542 pp.setY(pp.y());
00543 pp.setZ(pp.z());
00544 pp.setT(pp.t());
00545
00546 int charge = gun->GetMuonCharge();
00547
00548 m_pdgCode.push_back(charge*-13);
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559 double polx = 0;
00560 double poly = 0;
00561 double polz = 0;
00562
00563 Polarization thePolarization;
00564
00565
00566
00567
00568
00569
00570
00571 if(!m_swapYZAxis){
00572
00573 thePolarization.set_normal3d(HepNormal3D(polx,poly,polz));
00574 }
00575 else
00576 thePolarization.set_normal3d(HepNormal3D(polx,polz,-poly));
00577
00578 m_polarization.push_back(thePolarization);
00579
00580
00581
00582
00583 double e = pp.e();
00584 double theta = pp.theta();
00585 double phi = pp.phi();
00586
00587
00588
00589
00590 double p2 = e*e - mass*mass;
00591 if ( p2 < 0 )
00592 {
00593 log << MSG::ERROR
00594 << "Event #" << m_events
00595 << " E=" << e << ", mass=" << mass
00596 << " -- you have generated a tachyon! Increase energy or change particle ID."
00597 << endmsg;
00598 return StatusCode::FAILURE;
00599 }
00600
00601 double p = sqrt(p2);
00602 double px = p*sin(theta)*cos(phi);
00603 double pz = p*sin(theta)*sin(phi);
00604 double py = -p*cos(theta);
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614 if(!m_swapYZAxis) {
00615
00616
00617 m_fourMom.push_back(HepLorentzVector(px,py,pz,pp.e()));
00618 }
00619 else
00620 m_fourMom.push_back(HepLorentzVector(px,pz,-py,pp.e()));
00621 double x = vert.x();
00622 double y = vert.y();
00623 double z = vert.z();
00624 double t = vert.t();
00625 if(projected){
00626
00627 x = vert_proj.x();
00628 y = vert_proj.y();
00629 z = vert_proj.z();
00630 t = vert.t()-vert_proj.t()/HepLorentzVector(px,py,pz,pp.e()).beta()
00631 +m_tcor;
00632
00633
00634
00635 }
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647 if(!m_swapYZAxis)
00648 m_fourPos.push_back(HepLorentzVector(x,y,z,t));
00649 else
00650 m_fourPos.push_back(HepLorentzVector(x,z,y,t));
00651
00652 log << MSG::DEBUG
00653 << " (x,y,z,t) = ("
00654 << m_fourPos.back().x() << ","
00655 << m_fourPos.back().y() << ","
00656 << m_fourPos.back().z() << ","
00657 << m_fourPos.back().t() << "), (Px,Py,Pz,E) = ("
00658 << m_fourMom.back().px() << ","
00659 << m_fourMom.back().py() << ","
00660 << m_fourMom.back().pz() << ","
00661 << m_fourMom.back().e() << ")"
00662 << endreq;
00663 log << MSG::DEBUG
00664 << " (theta,phi) = (" << theta << "," << phi << "), "
00665 << "polarization(x,y,z) = ("
00666 << m_polarization.back().normal3d().x() << ","
00667 << m_polarization.back().normal3d().y() << ","
00668 << m_polarization.back().normal3d().z() << ")"
00669 << endreq;
00670 if(m_dumpMC==1) {
00671
00672
00673
00674
00675 mc_theta=m_fourMom.back().theta();
00676 mc_phi=m_fourMom.back().phi();
00677 mc_px=m_fourMom.back().px();
00678 mc_py=m_fourMom.back().py();
00679 mc_pz=m_fourMom.back().pz();
00680
00681 m_tuple1->write();
00682 }
00683 }
00684
00685
00686
00687 GenEvent* event = new GenEvent(1,1);
00688
00689 if(m_fourMom.size()==m_fourPos.size()&&m_fourMom.size()==m_polarization.size()){
00690
00691 for(std::size_t v=0;v<m_fourMom.size();++v){
00692
00693
00694
00695
00696
00697
00698 GenParticle* particle = new GenParticle( m_fourMom[v], m_pdgCode[v], 1);
00699 particle->set_polarization( m_polarization[v] );
00700
00701
00702 GenVertex* vertex = new GenVertex(m_fourPos[v]);
00703 vertex->add_particle_out( particle );
00704
00705
00706 event->add_vertex( vertex );
00707
00708 }
00709
00710 event->set_event_number(m_events);
00711
00712
00713
00714 } else {
00715
00716 log<<MSG::ERROR<<"Wrong different number of vertexes/momenta/polaritazions!"<<endreq;
00717 return StatusCode::FAILURE;
00718 }
00719
00720 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00721 if (anMcCol!=0)
00722 {
00723
00724
00725 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00726 McGenEvent* mcEvent = new McGenEvent(event);
00727 anMcCol->push_back(mcEvent);
00728 }
00729 else
00730 {
00731
00732 McGenEventCol *mcColl = new McGenEventCol;
00733 McGenEvent* mcEvent = new McGenEvent(event);
00734 mcColl->push_back(mcEvent);
00735 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00736 if (sc != StatusCode::SUCCESS)
00737 {
00738 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00739 delete mcColl;
00740 delete event;
00741 delete mcEvent;
00742 return StatusCode::FAILURE;
00743 }
00744
00745 }
00746
00747
00748
00749
00750
00751 return StatusCode::SUCCESS;
00752
00753 }
00754
00755
00756
00757 StatusCode CosmicGenerator::finalize() {
00758
00759
00760
00761 if(m_cubeProj) {
00762
00763 MsgStream log(messageService(), name());
00764
00765 log << MSG::INFO<<"***************************************************"<<endreq;
00766 log << MSG::INFO <<"** you have run CosmicGenerator with some " <<endreq;
00767 log << MSG::INFO <<"** filters for cosmic muon simulation" <<endreq;
00768 log << MSG::INFO <<"** "<<m_tried<<" muons were tried" <<endreq;
00769 log << MSG::INFO <<"** "<<m_accepted<<" muons were accepted" <<endreq;
00770 log << MSG::INFO <<"** "<<m_rejected<<" muons were rejected" <<endreq;
00771 log << MSG::INFO<<"***************************************************"<<endreq;
00772 std::cout<<"***************************************************"<<std::endl;
00773 std::cout <<"** you have run CosmicGenerator with some " <<std::endl;
00774 std::cout <<"** filters for cosmic muon simulation" <<std::endl;
00775 std::cout <<"** "<<m_tried<<" muons were tried" <<std::endl;
00776 std::cout <<"** "<<m_accepted<<" muons were accepted" <<std::endl;
00777 std::cout <<"** "<<m_rejected<<" muons were rejected" <<std::endl;
00778 std::cout<<"***************************************************"<<std::endl;
00779
00780 }
00781
00782 return StatusCode::SUCCESS;
00783 }
00784
00785
00786
00787
00788
00789 bool CosmicGenerator::exzCut(const Hep3Vector& pos,const HepLorentzVector& p)
00790 {
00791
00792
00793 double r =0;
00794 bool cut = false;
00795 if(pos.z()<0){
00796 r = sqrt((pow(pos.x(),2)+pow(pos.z()+28000,2))) ;
00797 double e = 0.45238*r+5000 ;
00798 cut = p.e()*m_GeV>e;
00799 }
00800 else
00801 {
00802 r = sqrt((pow(pos.x(),2)+pow(pos.z()-20000,2))) ;
00803 if(r<15000) {
00804 cut = true;
00805 } else
00806 {
00807 double e = 0.461538*(r-15000)+10000 ;
00808 cut = p.e()*m_GeV>e;
00809
00810 }
00811 }
00812
00813
00814 cut = cut && r < m_rmax ;
00815
00816 return cut;
00817
00818
00819 }