From 226701c6df848e2edf51a9933a76973ea75879ef Mon Sep 17 00:00:00 2001 From: Markus Frank <markus.frank@cern.ch> Date: Mon, 17 Nov 2014 09:28:41 +0000 Subject: [PATCH] Fix TrackerCombine sensitive detector -- thanks Andre --- DDG4/examples/CLICSidSimu.py | 3 +- DDG4/examples/CLICSidSimuMarkus.py | 125 ++++------------------------- DDG4/plugins/Geant4SDActions.cpp | 42 +++++----- DDG4/src/Geant4Action.cpp | 16 ++-- DDG4/src/Geant4Converter.cpp | 14 ++-- DDG4/src/Geant4ParticleHandler.cpp | 6 +- 6 files changed, 54 insertions(+), 152 deletions(-) diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py index 4557ddea3..f05e9845d 100644 --- a/DDG4/examples/CLICSidSimu.py +++ b/DDG4/examples/CLICSidSimu.py @@ -18,7 +18,7 @@ def run(): lcdd = kernel.lcdd() install_dir = os.environ['DD4hepINSTALL'] example_dir = install_dir+'/examples/DDG4/examples'; - kernel.loadGeometry("file:"+install_dir+"/examples/CLICSiD/compact/compact.xml") + kernel.loadGeometry("file:"+install_dir+"/DDDetectors/compact/SiD.xml") kernel.loadXML("file:"+example_dir+"/DDG4_field.xml") DDG4.importConstants(lcdd) @@ -58,7 +58,6 @@ def run(): """ Generation of isotrope tracks of a given multiplicity with overlay: """ - # First particle generator: pi+ gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropPi+"); gen.Particle = 'pi+' diff --git a/DDG4/examples/CLICSidSimuMarkus.py b/DDG4/examples/CLICSidSimuMarkus.py index 28e5eeebf..e8cfcc140 100644 --- a/DDG4/examples/CLICSidSimuMarkus.py +++ b/DDG4/examples/CLICSidSimuMarkus.py @@ -18,87 +18,37 @@ def run(): lcdd = kernel.lcdd() install_dir = os.environ['DD4hepINSTALL'] example_dir = install_dir+'/examples/DDG4/examples'; - kernel.loadGeometry("file:"+install_dir+"/examples/CLICSiD/compact/compact.xml") + kernel.loadGeometry("file:"+install_dir+"/DDDetectors/compact/SiD.xml") kernel.loadXML("file:"+example_dir+"/DDG4_field.xml") DDG4.importConstants(lcdd) - simple = DDG4.Simple(kernel) + simple = DDG4.Simple(kernel,tracker='Geant4TrackerCombineAction') simple.printDetectors() # Configure UI simple.setupCshUI() # Configure Run actions run1 = DDG4.RunAction(kernel,'Geant4TestRunAction/RunInit') - run1.Property_int = 12345 - run1.Property_double = -5e15*keV - run1.Property_string = 'Startrun: Hello_2' - print run1.Property_string, run1.Property_double, run1.Property_int run1.enableUI() kernel.registerGlobalAction(run1) kernel.runAction().adopt(run1) # Configure Event actions prt = DDG4.EventAction(kernel,'Geant4ParticlePrint/ParticlePrint') - prt.OutputLevel = Output.INFO + prt.OutputLevel = Output.DEBUG prt.OutputType = 3 # Print both: table and tree kernel.eventAction().adopt(prt) + generator_output_level = Output.DEBUG + # Configure I/O evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) - evt_lcio.OutputLevel = Output.ERROR - + evt_lcio.OutputLevel = generator_output_level evt_root = simple.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) - generator_output_level = Output.INFO - gen = DDG4.GeneratorAction(kernel,"Geant4GeneratorActionInit/GenerationInit") kernel.generatorAction().adopt(gen) - #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV - """ - Generation of isotrope tracks using the DDG4 partcle gun: - - # Setup particle gun - gun = simple.setupGun('Gun','pi-',energy=10*GeV,isotrop=True,multiplicity=1) - gun.OutputLevel = 5 # generator_output_level - gun.Mask = 1 - #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - """ - - #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV - """ - Generation of isotrope tracks of a given multiplicity with overlay: - - # First particle generator: pi+ - gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropPi+"); - gen.Particle = 'pi+' - gen.Energy = 100 * GeV - gen.Multiplicity = 2 - gen.Mask = 1 - kernel.generatorAction().adopt(gen) - # Install vertex smearing for this interaction - gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearPi+"); - gen.Mask = 1 - gen.Offset = (20*mm, 10*mm, 10*mm, 0*ns) - gen.Sigma = (4*mm, 1*mm, 1*mm, 0*ns) - kernel.generatorAction().adopt(gen) - - # Second particle generator: e- - gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropE-"); - gen.Particle = 'e-' - gen.Energy = 25 * GeV - gen.Multiplicity = 3 - gen.Mask = 2 - kernel.generatorAction().adopt(gen) - # Install vertex smearing for this interaction - gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearE-"); - gen.Mask = 2 - gen.Offset = (-20*mm, -10*mm, -10*mm, 0*ns) - gen.Sigma = (12*mm, 8*mm, 8*mm, 0*ns) - kernel.generatorAction().adopt(gen) - #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - """ - #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV """ Generation of primary particles from LCIO input files @@ -112,52 +62,27 @@ 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 = "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.Input = "LCIOFileReader|/afs/cern.ch/user/n/nikiforo/public/Markus/muons.slcio" + + gen.OutputLevel = generator_output_level gen.MomentumScale = 1.0 gen.Mask = 1 gen.enableUI() kernel.generatorAction().adopt(gen) - # Install vertex smearing for this interaction - gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/Smear1"); - gen.OutputLevel = 4 #generator_output_level - gen.Mask = 1 - gen.Offset = (-20*mm, -10*mm, -10*mm, 0*ns) - gen.Sigma = (12*mm, 8*mm, 8*mm, 0*ns) - gen.enableUI() - kernel.generatorAction().adopt(gen) - """ - # Second particle file reader - gen = DDG4.GeneratorAction(kernel,"LCIOInputAction/LCIO2"); - gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_2.stdhep" - #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/bbnn_3TeV_01.stdhep" - gen.OutputLevel = 4 # generator_output_level - gen.Mask = 2 - gen.MomentumScale = 0.1 - gen.enableUI() - kernel.generatorAction().adopt(gen) - # Install vertex smearing for this interaction - gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/Smear2"); - gen.OutputLevel = generator_output_level - gen.Mask = 2 - gen.Offset = (20*mm, 10*mm, 10*mm, 0*ns) - gen.Sigma = (2*mm, 1*mm, 1*mm, 0*ns) - gen.enableUI() - kernel.generatorAction().adopt(gen) - """ #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # Merge all existing interaction records gen = DDG4.GeneratorAction(kernel,"Geant4InteractionMerger/InteractionMerger") - gen.OutputLevel = 4 #generator_output_level + gen.OutputLevel = generator_output_level gen.enableUI() kernel.generatorAction().adopt(gen) # Finally generate Geant4 primaries gen = DDG4.GeneratorAction(kernel,"Geant4PrimaryHandler/PrimaryHandler") - gen.OutputLevel = 4 #generator_output_level + gen.OutputLevel = generator_output_level gen.enableUI() kernel.generatorAction().adopt(gen) @@ -167,7 +92,7 @@ def run(): #part.SaveProcesses = ['conv','Decay'] part.SaveProcesses = ['Decay'] part.MinimalKineticEnergy = 100*MeV - part.OutputLevel = 5 # generator_output_level + part.OutputLevel = Output.INFO #generator_output_level part.enableUI() user = DDG4.Action(kernel,"Geant4TCUserParticleHandler/UserParticleHandler") user.TrackingVolume_Zmax = DDG4.EcalEndcap_zmin @@ -185,30 +110,11 @@ def run(): kernel.generatorAction().adopt(rdr) """ - # Setup global filters fur use in sensntive detectors - f1 = DDG4.Filter(kernel,'GeantinoRejectFilter/GeantinoRejector') - f2 = DDG4.Filter(kernel,'ParticleRejectFilter/OpticalPhotonRejector') - f2.particle = 'opticalphoton' - f3 = DDG4.Filter(kernel,'ParticleSelectFilter/OpticalPhotonSelector') - f3.particle = 'opticalphoton' - f4 = DDG4.Filter(kernel,'EnergyDepositMinimumCut') - f4.Cut = 10*MeV - f4.enableUI() - kernel.registerGlobalFilter(f1) - kernel.registerGlobalFilter(f2) - kernel.registerGlobalFilter(f3) - kernel.registerGlobalFilter(f4) - + seq,act = simple.setupTracker('SiTrackerBarrel') + """ # First the tracking detectors seq,act = simple.setupTracker('SiVertexBarrel') - #seq.adopt(f1) - #seq.adopt(f4) - #act.adopt(f1) - seq,act = simple.setupTracker('SiVertexEndcap') - #seq.adopt(f1) - #seq.adopt(f4) - seq,act = simple.setupTracker('SiTrackerBarrel') seq,act = simple.setupTracker('SiTrackerEndcap') seq,act = simple.setupTracker('SiTrackerForward') @@ -222,6 +128,7 @@ def run(): seq,act = simple.setupCalorimeter('MuonEndcap') seq,act = simple.setupCalorimeter('LumiCal') seq,act = simple.setupCalorimeter('BeamCal') + """ # Now build the physics list: phys = simple.setupPhysics('QGSP_BERT') diff --git a/DDG4/plugins/Geant4SDActions.cpp b/DDG4/plugins/Geant4SDActions.cpp index 0450a7533..a3333a8a7 100644 --- a/DDG4/plugins/Geant4SDActions.cpp +++ b/DDG4/plugins/Geant4SDActions.cpp @@ -183,25 +183,22 @@ namespace DD4hep { * \ingroup DD4HEP_SIMULATION */ struct TrackerCombine { - typedef Geant4Tracker::Hit Hit; typedef Geant4HitCollection HitCollection; - Hit pre; - Hit post; + Geant4Tracker::Hit pre, post; Geant4Sensitive* sensitive; - G4Track* track; double e_cut; int current; int combined; - TrackerCombine() : pre(), post(), sensitive(0), track(0), e_cut(0.0), current(-1), combined(0) { + TrackerCombine() : pre(), post(), sensitive(0), e_cut(0.0), current(-1), combined(0) { } /// Start a new hit void start(G4Step* step, G4StepPoint* point) { pre.storePoint(step,point); + pre.truth.deposit = 0.0; current = pre.truth.trackID; - track = step->GetTrack(); - sensitive->mark(track); + sensitive->mark(step->GetTrack()); post = pre; combined = 0; } @@ -219,25 +216,24 @@ namespace DD4hep { pre.clear(); current = -1; combined = 0; - track = 0; } bool mustSaveTrack(const G4Track* tr) const { - return track && current != tr->GetTrackID(); + return current > 0 && current != tr->GetTrackID(); } /// Extract hit information and add the created hit to the collection - Hit* extractHit(HitCollection* collection) { - if ( current == -1 || !track ) { - return 0; + void extractHit(Geant4HitCollection* collection) { + if ( current == -1 ) { + return; } Position pos = 0.5 * (pre.position + post.position); Momentum mom = 0.5 * (pre.momentum + post.momentum); double path_len = (post.position - pre.position).R(); - Hit* hit = new Hit(pre.truth.trackID, - pre.truth.pdgID, - pre.truth.deposit, - pre.truth.time); + Geant4Tracker::Hit* hit = new Geant4Tracker::Hit(pre.truth.trackID, + pre.truth.pdgID, + pre.truth.deposit, + pre.truth.time); hit->position = pos; hit->momentum = mom; hit->length = path_len; @@ -247,7 +243,6 @@ namespace DD4hep { pre.truth.trackID,sensitive->c_name(),combined,pre.truth.deposit/MeV, pos.X()/mm,pos.Y()/mm,pos.Z()/mm); clear(); - return hit; } @@ -265,13 +260,12 @@ namespace DD4hep { if ( h.deposit()/keV <= 0 ) { return false; } - /// We can now start collecting the deposits of the next hit. - if ( 0 == track ) { + /// Initialize the deposits of the next hit. + if ( current < 0 ) { start(step, h.pre); } - - // ....update ..... - update(step); + /// ....update ..... + update(step); if ( prePV != postPV ) { void* postSD = h.sd(h.post); @@ -283,7 +277,7 @@ namespace DD4hep { } } } - else if ( track->GetTrackStatus() == fStopAndKill ) { + else if ( h.track->GetTrackStatus() == fStopAndKill ) { extractHit(coll); } return true; @@ -297,7 +291,7 @@ namespace DD4hep { // // Alternatively the 'update' method would become rather CPU consuming, // beacuse the extract action would have to be recalculated over and over. - if ( track ) { + if ( current > 0 ) { Geant4HitCollection* coll = sensitive->collection(0); extractHit(coll); } diff --git a/DDG4/src/Geant4Action.cpp b/DDG4/src/Geant4Action.cpp index f98389b19..4b5136e96 100644 --- a/DDG4/src/Geant4Action.cpp +++ b/DDG4/src/Geant4Action.cpp @@ -200,7 +200,7 @@ void Geant4Action::printP2(const char* fmt, ...) const { void Geant4Action::debug(const char* fmt, ...) const { va_list args; va_start(args, fmt); - DD4hep::printout(DD4hep::DEBUG, "Geant4Action", fmt, args); + DD4hep::printout(DD4hep::DEBUG, m_name, fmt, args); va_end(args); } @@ -208,7 +208,7 @@ void Geant4Action::debug(const char* fmt, ...) const { void Geant4Action::info(const char* fmt, ...) const { va_list args; va_start(args, fmt); - DD4hep::printout(DD4hep::INFO, "Geant4Action", fmt, args); + DD4hep::printout(DD4hep::INFO, m_name, fmt, args); va_end(args); } @@ -216,7 +216,7 @@ void Geant4Action::info(const char* fmt, ...) const { void Geant4Action::warning(const char* fmt, ...) const { va_list args; va_start(args, fmt); - DD4hep::printout(DD4hep::WARNING, "Geant4Action", fmt, args); + DD4hep::printout(DD4hep::WARNING, m_name, fmt, args); va_end(args); } @@ -224,7 +224,7 @@ void Geant4Action::warning(const char* fmt, ...) const { void Geant4Action::error(const char* fmt, ...) const { va_list args; va_start(args, fmt); - DD4hep::printout(DD4hep::ERROR, "Geant4Action", fmt, args); + DD4hep::printout(DD4hep::ERROR, m_name, fmt, args); va_end(args); } @@ -232,7 +232,7 @@ void Geant4Action::error(const char* fmt, ...) const { bool Geant4Action::error(bool return_value, const char* fmt, ...) const { va_list args; va_start(args, fmt); - DD4hep::printout(DD4hep::ERROR, "Geant4Action", fmt, args); + DD4hep::printout(DD4hep::ERROR, m_name, fmt, args); va_end(args); return return_value; } @@ -242,7 +242,7 @@ void Geant4Action::fatal(const char* fmt, ...) const { string err; va_list args; va_start(args, fmt); - DD4hep::printout(DD4hep::FATAL, "Geant4Action", fmt, args); + DD4hep::printout(DD4hep::FATAL, m_name, fmt, args); va_end(args); } @@ -251,8 +251,8 @@ void Geant4Action::except(const char* fmt, ...) const { string err; va_list args; va_start(args, fmt); - DD4hep::printout(DD4hep::FATAL, "Geant4Action", fmt, args); - err = DD4hep::format("Geant4Action", fmt, args); + DD4hep::printout(DD4hep::FATAL, m_name, fmt, args); + err = DD4hep::format(m_name, fmt, args); va_end(args); throw runtime_error(err); } diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp index a33e70a70..51430a98f 100644 --- a/DDG4/src/Geant4Converter.cpp +++ b/DDG4/src/Geant4Converter.cpp @@ -11,6 +11,7 @@ #include "DD4hep/Plugins.h" #include "DD4hep/Volumes.h" #include "DD4hep/Printout.h" +#include "DD4hep/DD4hepUnits.h" #include "DD4hep/objects/ObjectsInterna.h" #include "DD4hep/objects/DetectorInterna.h" @@ -77,6 +78,7 @@ #include "G4ElectroMagneticField.hh" #include "G4FieldManager.hh" #include "G4ReflectionFactory.hh" +#include "CLHEP/Units/SystemOfUnits.h" #include <iostream> #include <iomanip> @@ -724,13 +726,13 @@ void* Geant4Converter::handleRegion(Region region, const set<const TGeoVolume*>& g4 = new G4Region(r.name()); // set production cut G4ProductionCuts* cuts = new G4ProductionCuts(); - cuts->SetProductionCut(r.cut()); + cuts->SetProductionCut(r.cut()*CLHEP::MeV/dd4hep::MeV); g4->SetProductionCuts(cuts); // create region info with storeSecondaries flag G4UserRegionInformation* info = new G4UserRegionInformation(); info->region = r; - info->threshold = r.threshold(); + info->threshold = r.threshold()*CLHEP::MeV/dd4hep::MeV; info->storeSecondaries = r.storeSecondaries(); g4->SetUserInformation(info); @@ -769,13 +771,13 @@ void* Geant4Converter::handleLimitSet(LimitSet limitset, const set<const TGeoVol for (LimitSet::Object::const_iterator i = limits.begin(); i != limits.end(); ++i) { const Limit& l = *i; if (l.name == "step_length_max") - g4->SetMaxAllowedStep(l.value); + g4->SetMaxAllowedStep(l.value*CLHEP::mm/dd4hep::mm); else if (l.name == "track_length_max") - g4->SetMaxAllowedStep(l.value); + g4->SetMaxAllowedStep(l.value*CLHEP::mm/dd4hep::mm); else if (l.name == "time_max") - g4->SetUserMaxTime(l.value); + g4->SetUserMaxTime(l.value*CLHEP::ns/dd4hep::ns); else if (l.name == "ekin_min") - g4->SetUserMinEkine(l.value); + g4->SetUserMinEkine(l.value*CLHEP::MeV/dd4hep::MeV); else if (l.name == "range_min") g4->SetUserMinRange(l.value); else diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp index c39f4b905..09742b6c0 100644 --- a/DDG4/src/Geant4ParticleHandler.cpp +++ b/DDG4/src/Geant4ParticleHandler.cpp @@ -167,7 +167,7 @@ void Geant4ParticleHandler::mark(const G4Track* track) { /// Event generation action callback void Geant4ParticleHandler::operator()(G4Event* event) { typedef Geant4MonteCarloTruth _MC; - info("+++ Event:%d Add EVENT extension of type Geant4ParticleHandler.....",event->GetEventID()); + debug("+++ Event:%d Add EVENT extension of type Geant4ParticleHandler.....",event->GetEventID()); context()->event().addExtension((_MC*)this, typeid(_MC), 0); clear(); /// Call the user particle handler @@ -374,13 +374,13 @@ void Geant4ParticleHandler::endEvent(const G4Event* event) { int level = outputLevel(); do { if ( level <= VERBOSE ) dumpMap("Particle"); - print("+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size()); + debug("+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size()); } while( recombineParents() > 0 ); if ( level <= VERBOSE ) dumpMap("Recombined"); // Rebase the simulated tracks, so that they fit to the generator particles rebaseSimulatedTracks(0); - if ( level <= DEBUG ) dumpMap("Rebased"); + if ( level <= VERBOSE ) dumpMap("Rebased"); // Consistency check.... checkConsistency(); /// Call the user particle handler -- GitLab