From 140bde0a829363f631f4ef038438531883ae6c2e Mon Sep 17 00:00:00 2001 From: Markus Frank <markus.frank@cern.ch> Date: Mon, 3 Nov 2014 20:15:27 +0000 Subject: [PATCH] First implementation of dipole field, HepMC reader for DDG4 and fix to import constants into DDG4 python --- DDCore/src/Evaluator/setStdMath.cpp | 2 ++ DDCore/src/FieldTypes.cpp | 1 + DDCore/src/plugins/Compact2Objects.cpp | 7 ++-- DDG4/examples/CLICSidSimuMarkus.py | 20 ++++------- DDG4/examples/readHEPMC.py | 17 +++++----- DDG4/include/DDG4/Geant4Particle.h | 1 + DDG4/plugins/Geant4EventReaderHepMC.cpp | 25 ++++++++++---- DDG4/python/DDG4.py | 45 +++++++++++++++++-------- DDG4/src/Geant4InputAction.cpp | 2 +- DDG4/src/Geant4InputHandling.cpp | 2 +- DDG4/src/Geant4Particle.cpp | 27 +++++++++++++-- examples/CLICSiD/compact/compact.xml | 20 +++++++++-- 12 files changed, 118 insertions(+), 51 deletions(-) diff --git a/DDCore/src/Evaluator/setStdMath.cpp b/DDCore/src/Evaluator/setStdMath.cpp index d72ead330..db1857f60 100644 --- a/DDCore/src/Evaluator/setStdMath.cpp +++ b/DDCore/src/Evaluator/setStdMath.cpp @@ -27,6 +27,7 @@ static double eval_tanh (double a) { return std::tanh(a); } static double eval_exp (double a) { return std::exp(a); } static double eval_log (double a) { return std::log(a); } static double eval_log10(double a) { return std::log10(a); } +static double eval_double(double a) { return a; } static double eval_int (double a) { return (double)int(a); } static double eval_nint (double a) { return std::floor(a); } static double eval_floor(double a) { return std::floor(a); } @@ -54,6 +55,7 @@ void Evaluator::setStdMath() { // S E T S T A N D A R D F U N C T I O N S setFunction("floor", eval_floor); + setFunction("double",eval_double); setFunction("int", eval_int); setFunction("nint", eval_nint); setFunction("abs", eval_abs); diff --git a/DDCore/src/FieldTypes.cpp b/DDCore/src/FieldTypes.cpp index df631bc28..8413b0789 100644 --- a/DDCore/src/FieldTypes.cpp +++ b/DDCore/src/FieldTypes.cpp @@ -51,6 +51,7 @@ void SolenoidField::fieldComponents(const double* pos, double* field) { /// Initializing constructor DipoleField::DipoleField() : zmax(INFINITY), zmin(-INFINITY), rmax(INFINITY) { + type = CartesianField::MAGNETIC; } /// Call to access the field components at a given location diff --git a/DDCore/src/plugins/Compact2Objects.cpp b/DDCore/src/plugins/Compact2Objects.cpp index b16dc6460..d64a7fd9b 100644 --- a/DDCore/src/plugins/Compact2Objects.cpp +++ b/DDCore/src/plugins/Compact2Objects.cpp @@ -10,6 +10,7 @@ // Framework includes #include "DD4hep/DetFactoryHelper.h" #include "DD4hep/IDDescriptor.h" +#include "DD4hep/DD4hepUnits.h" #include "DD4hep/FieldTypes.h" #include "DD4hep/Printout.h" #include "DD4hep/Plugins.h" @@ -138,7 +139,9 @@ static Ref_t create_DipoleField(lcdd_t& /* lcdd */, xml_h e) { xml_comp_t c(e); CartesianField obj; DipoleField* ptr = new DipoleField(); - double val, lunit = c.attr<double>(_U(lunit)), funit = c.attr<double>(_U(funit)); + double val; + double lunit = c.hasAttr(_U(lunit)) ? c.attr<double>(_U(lunit)) : 1.0; + double funit = c.hasAttr(_U(funit)) ? c.attr<double>(_U(funit)) : 1.0; if (c.hasAttr(_U(zmin))) ptr->zmin = _multiply<double>(c.attr < string > (_U(zmin)), lunit); @@ -148,7 +151,7 @@ static Ref_t create_DipoleField(lcdd_t& /* lcdd */, xml_h e) { ptr->rmax = _multiply<double>(c.attr < string > (_U(rmax)), lunit); for (xml_coll_t coll(c, _U(dipole_coeff)); coll; ++coll) { val = funit / pow(lunit, (int) ptr->coefficents.size()); - val = _multiply<double>(coll.value(), val); + val = _multiply<double>(coll.text(), val); ptr->coefficents.push_back(val); } obj.assign(ptr, c.nameStr(), c.typeStr()); diff --git a/DDG4/examples/CLICSidSimuMarkus.py b/DDG4/examples/CLICSidSimuMarkus.py index 4642ca2be..28e5eeebf 100644 --- a/DDG4/examples/CLICSidSimuMarkus.py +++ b/DDG4/examples/CLICSidSimuMarkus.py @@ -112,9 +112,9 @@ def run(): #gen.Input = "LCIOFileReader|/home/frankm/SW/data/mcparticles_pi-_5GeV.slcio" #gen.Input = "LCIOFileReader|/home/frankm/SW/data/mcparticles_mu-_5GeV.slcio" #gen.Input = "LCIOFileReader|/home/frankm/SW/data/bbbb_3TeV.slcio" - #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/FCC-eh.stdhep" - gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt" - gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/sherpa-2.1.1_zjets.hepmc2g" + gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/FCC-eh.stdhep" + #gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt" + #gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/sherpa-2.1.1_zjets.hepmc2g" gen.OutputLevel = 4 # generator_output_level gen.MomentumScale = 1.0 gen.Mask = 1 @@ -201,12 +201,12 @@ def run(): # First the tracking detectors seq,act = simple.setupTracker('SiVertexBarrel') - seq.adopt(f1) + #seq.adopt(f1) #seq.adopt(f4) - act.adopt(f1) + #act.adopt(f1) seq,act = simple.setupTracker('SiVertexEndcap') - seq.adopt(f1) + #seq.adopt(f1) #seq.adopt(f4) seq,act = simple.setupTracker('SiTrackerBarrel') @@ -225,14 +225,6 @@ def run(): # Now build the physics list: phys = simple.setupPhysics('QGSP_BERT') - ph = DDG4.PhysicsList(kernel,'Geant4PhysicsList/Myphysics') - ph.addParticleConstructor('G4BosonConstructor') - ph.addParticleConstructor('G4LeptonConstructor') - ph.addParticleProcess('e[+-]','G4eMultipleScattering',-1,1,1) - ph.addPhysicsConstructor('G4OpticalPhysics') - ph.enableUI() - phys.adopt(ph) - phys.dump() kernel.configure() diff --git a/DDG4/examples/readHEPMC.py b/DDG4/examples/readHEPMC.py index 1098db629..e40ef0789 100644 --- a/DDG4/examples/readHEPMC.py +++ b/DDG4/examples/readHEPMC.py @@ -1,6 +1,6 @@ # # -import os, time, DDG4 +import os, time, getopt, DDG4 from DDG4 import OutputLevel as Output from SystemOfUnits import * # @@ -16,17 +16,18 @@ DD4hep simulation example setup using the python configuration def run(): kernel = DDG4.Kernel() lcdd = kernel.lcdd() - simple = DDG4.Simple(kernel) gen = DDG4.GeneratorAction(kernel,"Geant4InputAction/Input") kernel.generatorAction().adopt(gen) - gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt"; + gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt" + gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/Atlas_Pythia8.hepmc" + gen.OutputLevel = Output.DEBUG parts = gen.new_particles() - gen.readParticles(0,parts); - parts.clear() - print 132*'*' - print 132*'*' - gen.readParticles(1,parts); + ret = 1 + while ret: + ret = gen.readParticles(0,parts) + parts.clear() + print 132*'*',ret if __name__ == "__main__": run() diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h index ed3614f86..252cd5f8d 100644 --- a/DDG4/include/DDG4/Geant4Particle.h +++ b/DDG4/include/DDG4/Geant4Particle.h @@ -198,6 +198,7 @@ namespace DD4hep { /// Output type 3:+++ "tag" ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] \#Par: 0 \#Dau: 4 void dumpWithVertex(int level, const std::string& src, const char* tag) const; void dumpWithMomentum(int level, const std::string& src, const char* tag) const; + void dumpWithMomentumAndVertex(int level, const std::string& src, const char* tag) const; void dump4(int level, const std::string& src, const char* tag) const; /// Handlers diff --git a/DDG4/plugins/Geant4EventReaderHepMC.cpp b/DDG4/plugins/Geant4EventReaderHepMC.cpp index 48a500911..815f54e44 100644 --- a/DDG4/plugins/Geant4EventReaderHepMC.cpp +++ b/DDG4/plugins/Geant4EventReaderHepMC.cpp @@ -239,7 +239,7 @@ void HepMC::fix_particles(EventStream& info) { } } for(vector<Geant4Particle*>::iterator i=beam.begin(); i!=beam.end();++i) { - cout << "Clear parents of " << (*i)->id << endl; + //cout << "Clear parents of " << (*i)->id << endl; (*i)->parents.clear(); (*i)->status = G4PARTICLE_GEN_DECAYED; } @@ -264,6 +264,7 @@ char HepMC::get_input(istream& is, istringstream& iline) { is.clear(ios::badbit); return -1; } + iline.clear(); iline.str(line); if ( line.length()>0 && line[1] == ' ' ) iline >> firstc; return iline ? value : -1; @@ -490,8 +491,8 @@ int HepMC::read_heavy_ion(EventStream &, istringstream & input) { float impact = 0., plane = 0., xcen = 0., inel = 0.; input >> nh >> np >> nt >> nc >> neut >> prot >> nw >> nwn >> nwnw; input >> impact >> plane >> xcen >> inel; - cerr << "Reading heavy ion, but igoring data!" << endl; /* + cerr << "Reading heavy ion, but igoring data!" << endl; ion->set_Ncoll_hard(nh); ion->set_Npart_proj(np); ion->set_Npart_targ(nt); @@ -527,8 +528,8 @@ int HepMC::read_pdf(EventStream &, istringstream & input) { if( !input.eof() ) { input >> pdf_id1 >> pdf_id2; } - cerr << "Reading pdf, but igoring data!" << endl; /* + cerr << "Reading pdf, but igoring data!" << endl; pdf->set_id1( id1 ); pdf->set_id2( id2 ); pdf->set_pdf_id1( pdf_id1 ); @@ -560,16 +561,26 @@ void HepMC::EventStream::clear() { bool HepMC::EventStream::read() { EventStream& info = *this; bool event_read = false; + static int num_evt = 0; + int num_line = 0, num_line_accepted = 0; releaseObjects(vertices())(); releaseObjects(particles())(); + ++num_evt; + if ( num_evt == 998 ) { + cout << "Hallo"; + } while( instream.good() ) { char value = instream.peek(); istringstream input_line; - - if ( value == 'E' && event_read ) break; - else if ( instream.eof() ) return false; + ++num_line; + if ( value == 'E' && event_read ) + break; + else if ( instream.eof() && event_read ) + goto Done; + else if ( instream.eof() ) + return false; else if ( value=='#' || ::isspace(value) ) { get_input(instream,input_line); continue; @@ -580,6 +591,7 @@ bool HepMC::EventStream::read() { if( !input_line || value < 0 ) goto Skip; + ++num_line_accepted; switch( value ) { case 'H': { int iotype = 0; @@ -670,6 +682,7 @@ bool HepMC::EventStream::read() { event_read = false; if ( instream.eof() ) return false; } + Done: fix_particles(info); releaseObjects(vertices())(); return true; diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py index 9002881b3..7a331166d 100644 --- a/DDG4/python/DDG4.py +++ b/DDG4/python/DDG4.py @@ -41,12 +41,9 @@ LCDD.globalVal = _constant #--------------------------------------------------------------------------- """ -import DDG4 -l=DDG4.LCDD.getInstance() -l.addConstant(DDG4.Ref_t(DDG4.NamedObject('AA','10*cm'))) -DDG4.importConstants(l,'sub-name-space') + Import the LCDD constants into the DDG4 namespace """ -def importConstants(lcdd,namespace=None): +def importConstants(lcdd,namespace=None,debug=False): scope = current ns = current if namespace is not None and not hasattr(current,namespace): @@ -54,15 +51,35 @@ def importConstants(lcdd,namespace=None): m = imp.new_module('DDG4.'+namespace) setattr(current,namespace,m) ns = m - values = {} - for c in lcdd.constants(): - values[c.first] = c.second.GetTitle() - evaluator = DD4hep.evaluator() - for k in values.keys(): - setattr(ns,k,values[k]) - #print 'Imported global value:',c.first,'=',c.second.GetTitle(),'into namespace',ns.__name__ -#--------------------------------------------------------------------------- - + evaluator = DD4hep.g4Evaluator() + cnt = 0 + num = 0 + todo = {} + for c in lcdd.constants(): todo[c.first] = c.second.GetTitle().replace('(int)','') + while len(todo) and cnt<100: + cnt = cnt + 1 + if cnt == 100: + print '%s %d out of %d %s "%s": [%s]\n+++ %s'\ + %('+++ FAILED to import', + len(todo),len(todo)+num, + 'global values into namespace', + ns.__name__,'Try to cintinue anyway',100*'=',) + for k,v in todo.items(): + if not hasattr(ns,k): + val = evaluator.evaluate(v) + status = evaluator.status() + if status == 0: + evaluator.setVariable(k,val) + setattr(ns,k,val) + if debug: print 'Imported global value:',k,'=',val,'into namespace',ns.__name__ + del todo[k] + num = num + 1 + elif cnt == 100: + print 'FAILED to import:',k,'=',v + if cnt<100: + print '+++ Imported %d global values to namespace:%s'%(num,ns.__name__,) + +#--------------------------------------------------------------------------- def _registerGlobalAction(self,action): self.get().registerGlobalAction(Interface.toAction(action)) def _registerGlobalFilter(self,filter): diff --git a/DDG4/src/Geant4InputAction.cpp b/DDG4/src/Geant4InputAction.cpp index bfbf58ef9..b649d35f2 100644 --- a/DDG4/src/Geant4InputAction.cpp +++ b/DDG4/src/Geant4InputAction.cpp @@ -170,6 +170,6 @@ void Geant4InputAction::operator()(G4Event* event) { vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters } inter->particles.insert(make_pair(p->id,p)); - p.dumpWithVertex(outputLevel()-1,name(),"+->"); + p.dumpWithMomentumAndVertex(outputLevel()-1,name(),"->"); } } diff --git a/DDG4/src/Geant4InputHandling.cpp b/DDG4/src/Geant4InputHandling.cpp index 2e09c3a0d..53fa7e617 100644 --- a/DDG4/src/Geant4InputHandling.cpp +++ b/DDG4/src/Geant4InputHandling.cpp @@ -325,7 +325,7 @@ int DD4hep::Simulation::generatePrimaries(const Geant4Action* caller, PropertyMask reason(r->reason); reason.set(G4PARTICLE_PRIMARY); v4->SetPrimary(p4); - ::snprintf(text,sizeof(text),"+-> G4Primary[%3d]",num_part); + ::snprintf(text,sizeof(text),"-> G4Primary[%3d]",num_part); r.dumpWithMomentum(caller->outputLevel()-1,caller->name(),text); ++num_part; } diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp index e9f58bcb5..996dedda3 100644 --- a/DDG4/src/Geant4Particle.cpp +++ b/DDG4/src/Geant4Particle.cpp @@ -269,7 +269,7 @@ void Geant4ParticleHandle::dumpWithVertex(int level, const std::string& src, con ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i); } printout((DD4hep::PrintLevel)level,src, - "+++ %s ID:%3d %-12s status:%08X typ:%9d Vtx:(%+.2e,%+.2e,%+.2e)[mm] " + "+++ %s ID:%3d %-12s status:%08X PDG:%6d Vtx:(%+.2e,%+.2e,%+.2e)[mm] " "time: %+.2e [ns] #Dau:%3d #Par:%1d%-6s", tag,p->id,p.particleName().c_str(),p->status,p->pdgID, p->vsx/mm,p->vsy/mm,p->vsz/mm,p->time/ns, @@ -292,7 +292,7 @@ void Geant4ParticleHandle::dumpWithMomentum(int level, const std::string& src, c ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i); } printout((DD4hep::PrintLevel)level,src, - "+++ %s ID:%3d %-12s status:%08X typ:%9d Mom:(%+.2e,%+.2e,%+.2e)[mm] " + "+++%s ID:%3d %-12s stat:%08X PDG:%6d Mom:(%+.2e,%+.2e,%+.2e)[MeV] " "time: %+.2e [ns] #Dau:%3d #Par:%1d%-6s", tag,p->id,p.particleName().c_str(),p->status,p->pdgID, p->psx/MeV,p->psy/MeV,p->psz/MeV,p->time/ns, @@ -301,6 +301,29 @@ void Geant4ParticleHandle::dumpWithMomentum(int level, const std::string& src, c text); } +/// Output type 3:+++ <tag> ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] \#Par: 0 \#Dau: 4 +void Geant4ParticleHandle::dumpWithMomentumAndVertex(int level, const std::string& src, const char* tag) const { + char text[256]; + Geant4ParticleHandle p(*this); + text[0]=0; + if ( p->parents.size() == 1 ) + ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin())); + else if ( p->parents.size() > 1 ) { + text[0]='/';text[1]=0; + for(std::set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i) + ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i); + } + printout((DD4hep::PrintLevel)level,src, + "+++%s %3d %-12s stat:%08X PDG:%6d Mom:(%+.2e,%+.2e,%+.2e)[MeV] " + "Vtx:(%+.2e,%+.2e,%+.2e)[mm] #Dau:%3d #Par:%1d%-6s", + tag,p->id,p.particleName().c_str(),p->status,p->pdgID, + p->psx/MeV,p->psy/MeV,p->psz/MeV, + p->vsx/mm,p->vsy/mm,p->vsz/mm, + int(p->daughters.size()), + int(p->parents.size()), + text); +} + void Geant4ParticleHandle::dump4(int level, const std::string& src, const char* tag) const { Geant4ParticleHandle p(*this); char equiv[32]; diff --git a/examples/CLICSiD/compact/compact.xml b/examples/CLICSiD/compact/compact.xml index 6b84eab79..5856ed685 100644 --- a/examples/CLICSiD/compact/compact.xml +++ b/examples/CLICSiD/compact/compact.xml @@ -28,10 +28,12 @@ <constant name="CaloSides" value="12"/> <constant name="MuonSides" value="8"/> + <constant name="EcalBarrel_rmin" value="126.50*cm + world_side - world_side"/> <constant name="EcalBarrel_rmin" value="126.50*cm"/> <constant name="EcalBarrel_zmax" value="176.50*cm"/> <constant name="EcalEndcap_rmin" value="21.0*cm"/> - <constant name="EcalEndcap_rmax" value="(EcalBarrel_rmin - 1.5*cm) / (cos(pi/CaloSides))"/> <!-- Correction from going from inner circle to outer circle --> + <!-- Correction from going from inner circle to outer circle --> + <constant name="EcalEndcap_rmax" value="(EcalBarrel_rmin - 1.5*cm) / (cos(pi/CaloSides))"/> <constant name="EcalEndcap_zmin" value="165.70*cm"/> <constant name="HcalBarrel_rmin" value="141.90*cm"/> @@ -39,7 +41,8 @@ <constant name="HcalBarrel_layer_thickness" value="1.0*cm + 0.65*cm"/> <constant name="HcalEndcap_zmin" value="EcalBarrel_zmax + 4.0*cm"/> <!-- Gap for cables --> <constant name="HcalEndcap_rmin" value="50.0*cm"/> - <constant name="HcalEndcap_rmax" value="(HcalBarrel_rmin + HcalBarrel_layers * HcalBarrel_layer_thickness) / (cos(pi/CaloSides))"/> <!-- Correction from going from inner circle to outer circle --> + <!-- Correction from going from inner circle to outer circle --> + <constant name="HcalEndcap_rmax" value="(HcalBarrel_rmin + HcalBarrel_layers * HcalBarrel_layer_thickness) / (cos(pi/CaloSides))"/> <constant name="HcalEndcap_layers" value="60"/> <constant name="HcalEndcap_layer_thickness" value="2.0*cm + 0.65*cm"/> <constant name="HcalEndcap_zmax" value="HcalEndcap_zmin + HcalEndcap_layers * HcalEndcap_layer_thickness"/> @@ -71,7 +74,8 @@ <constant name="MuonBarrel_layer_thickness" value="10.0*cm + 4.0*cm"/> <constant name="MuonEndcap_zmin" value="MuonBarrel_zmax + 10.0*cm"/> <!-- Space for cables etc. --> <constant name="MuonEndcap_rmin" value="69.0*cm"/> <!-- Space for QD0 and anti-solenoid--> - <constant name="MuonEndcap_rmax" value="(MuonBarrel_rmin + 57.0*cm + MuonBarrel_layers * MuonBarrel_layer_thickness) / (cos(pi/MuonSides))"/> <!-- Correction from going from inner circle to outer circle --> + <!-- Correction from going from inner circle to outer circle --> + <constant name="MuonEndcap_rmax" value="(MuonBarrel_rmin + 57.0*cm + MuonBarrel_layers * MuonBarrel_layer_thickness) / (cos(pi/MuonSides))"/> <constant name="MuonEndcap_layers" value="18"/> <constant name="MuonEndcap_layer_thickness" value="10.0*cm + 4.0*cm"/> <constant name="MuonEndcap_zmax" value="MuonEndcap_zmin + MuonEndcap_layers * MuonEndcap_layer_thickness"/> @@ -1534,6 +1538,7 @@ <id>system:8,layer:8,barrel:3,layer:8,slice:5,x:32:-16,y:-16</id> </readout> </readouts> + <fields> <field name="GlobalSolenoid" type="solenoid" inner_field="5.0*tesla" @@ -1541,5 +1546,14 @@ zmax="SolenoidCoilOuterZ" outer_radius="SolenoidalFieldRadius"> </field> + + <field name="GlobalDipole" type="DipoleMagnet" + rmax="50*cm" + zmin="0*cm" + zmax="50*cm"> + <dipole_coeff>1.0*tesla</dipole_coeff> + <dipole_coeff>0.1*tesla/pow(cm,1)</dipole_coeff> + <dipole_coeff>0.01*tesla/pow(cm,2)</dipole_coeff> + </field> </fields> </lccdd> -- GitLab