From 6c754d61834d92ce48e545e14acd3f97c135f91d Mon Sep 17 00:00:00 2001 From: Markus Frank <markus.frank@cern.ch> Date: Wed, 9 Dec 2015 15:06:08 +0000 Subject: [PATCH] Fix Geant4TrackerWeightedSD for energy weighted positions and remove CppCheck warnings --- DDCore/src/Conditions.cpp | 7 ++++- DDCore/src/Evaluator/hash_map.src | 4 ++- DDCore/src/Volumes.cpp | 3 ++- DDCore/src/XML/tinyxml_inl.h | 8 ++---- DDEve/src/DDG4EventHandler.cpp | 11 ++++---- DDG4/g4gdmlDisplay.cpp | 1 - DDG4/include/DDG4/Geant4StepHandler.h | 27 ++++--------------- .../Geant4DetectorSensitivesConstruction.cpp | 4 +-- DDG4/plugins/Geant4EventReaderHepMC.cpp | 11 ++++---- DDG4/plugins/Geant4TrackerWeightedSD.cpp | 18 ++++++++----- DDG4/src/Geant4StepHandler.cpp | 22 +++++++++++++++ 11 files changed, 63 insertions(+), 53 deletions(-) diff --git a/DDCore/src/Conditions.cpp b/DDCore/src/Conditions.cpp index b7964ed3a..baa24ca64 100644 --- a/DDCore/src/Conditions.cpp +++ b/DDCore/src/Conditions.cpp @@ -139,7 +139,12 @@ const type_info& Condition::typeInfo() const { /// Access to the grammar type const DD4hep::BasicGrammar& Condition::descriptor() const { const BasicGrammar* g = access()->data.grammar; - if ( !g ) invalidHandleError<Condition>(); + if ( !g ) { + invalidHandleError<Condition>(); + // This code is never reached, since function above throws exception! + // Needed to satisfay CppCheck + throw runtime_error("Null pointer in Grammar object"); + } return *g; } diff --git a/DDCore/src/Evaluator/hash_map.src b/DDCore/src/Evaluator/hash_map.src index ee1dc2abd..6ca44e7a9 100644 --- a/DDCore/src/Evaluator/hash_map.src +++ b/DDCore/src/Evaluator/hash_map.src @@ -29,7 +29,8 @@ template<class K, class T> class hash_map { private: hash_map(const hash_map& c) - : table(0), cur_size(0), max_size(0), max_load(0), default_value(c.default_value) + : table(0), cur_size(0), max_size(0), max_load(0), + grow(0.0), default_value(c.default_value) { } hash_map& operator=(const hash_map& c) { @@ -37,6 +38,7 @@ private: cur_size = 0; max_size = 0; max_load = 0; + grow = 0.0; default_value = c.default_value; return *this; } diff --git a/DDCore/src/Volumes.cpp b/DDCore/src/Volumes.cpp index da1338db7..583abce50 100644 --- a/DDCore/src/Volumes.cpp +++ b/DDCore/src/Volumes.cpp @@ -595,7 +595,8 @@ const Volume& Volume::setVisAttributes(const VisAttr& attr) const { m_element->SetFillColor(bright); #endif m_element->SetFillStyle(1001); // Root: solid - m_element->GetMedium()->GetMaterial()->SetTransparency((1-vis->alpha)*100); + // Suggested by Nikiforos. Not optimal. + //m_element->GetMedium()->GetMaterial()->SetTransparency((1-vis->alpha)*100); } else { printout(DEBUG,"setVisAttributes","Set to wireframe vis:%s",name()); diff --git a/DDCore/src/XML/tinyxml_inl.h b/DDCore/src/XML/tinyxml_inl.h index d6e822250..6d9c6272b 100644 --- a/DDCore/src/XML/tinyxml_inl.h +++ b/DDCore/src/XML/tinyxml_inl.h @@ -137,11 +137,9 @@ TiXmlNode::TiXmlNode( NodeType _type ) : TiXmlBase() TiXmlNode::~TiXmlNode() { TiXmlNode* node = firstChild; - TiXmlNode* temp = 0; - while ( node ) { - temp = node; + TiXmlNode* temp = node; node = node->next; delete temp; } @@ -158,11 +156,9 @@ void TiXmlNode::CopyTo( TiXmlNode* target ) const void TiXmlNode::Clear() { TiXmlNode* node = firstChild; - TiXmlNode* temp = 0; - while ( node ) { - temp = node; + TiXmlNode* temp = node; node = node->next; delete temp; } diff --git a/DDEve/src/DDG4EventHandler.cpp b/DDEve/src/DDG4EventHandler.cpp index 29170f9d6..72d9532de 100644 --- a/DDEve/src/DDG4EventHandler.cpp +++ b/DDEve/src/DDG4EventHandler.cpp @@ -36,11 +36,11 @@ namespace { DDG4EventHandler::ParticleAccessor_t particles; void* ptr; }; -} - -static void* _create(const char*) { - EventHandler* h = new DDG4EventHandler(); - return h; + /// Factory entry point + void* _create(const char*) { + EventHandler* h = new DDG4EventHandler(); + return h; + } } using namespace DD4hep::Geometry; DECLARE_CONSTRUCTOR(DDEve_DDG4EventHandler,_create) @@ -189,7 +189,6 @@ Int_t DDG4EventHandler::ReadEvent(Long64_t event_number) { /// Open new data file bool DDG4EventHandler::Open(const std::string&, const std::string& name) { - string err; if ( m_file.first ) m_file.first->Close(); m_hasFile = false; m_hasEvent = false; diff --git a/DDG4/g4gdmlDisplay.cpp b/DDG4/g4gdmlDisplay.cpp index 9aa3d9157..bfe48f24c 100644 --- a/DDG4/g4gdmlDisplay.cpp +++ b/DDG4/g4gdmlDisplay.cpp @@ -64,7 +64,6 @@ int main(int argc, char** argv) { string setup = argv[2]; const char* args[] = {"cmd"}; for(int i=1; i<argc;++i) { - string nam = get_arg(argc,argv,i); if ( argv[i][0]=='-' ) { string n = argv[i]+1; if ( ::strncmp(n.c_str(),"gdml",4) == 0 ) diff --git a/DDG4/include/DDG4/Geant4StepHandler.h b/DDG4/include/DDG4/Geant4StepHandler.h index be86c7f1e..03ec80437 100644 --- a/DDG4/include/DDG4/Geant4StepHandler.h +++ b/DDG4/include/DDG4/Geant4StepHandler.h @@ -73,7 +73,7 @@ namespace DD4hep { /// Returns total energy deposit double totalEnergy() const { if(applyBirksLaw == true) - return BirkAttenuation(step); + return this->birkAttenuation(); else return step->GetTotalEnergyDeposit(); } @@ -129,6 +129,9 @@ namespace DD4hep { double deposit() const { return step->GetTotalEnergyDeposit(); } + double stepLength() const { + return step->GetStepLength(); + } int trkID() const { return track->GetTrackID(); } @@ -234,27 +237,7 @@ namespace DD4hep { G4ThreeVector globalToLocalG4(const G4ThreeVector& loc) const; /// Apply BirksLaw - double BirkAttenuation(const G4Step* aStep) const - { - double energyDeposition = aStep->GetTotalEnergyDeposit(); - double length = aStep->GetStepLength(); - double niel = aStep->GetNonIonizingEnergyDeposit(); - const G4Track* trk = aStep->GetTrack(); - const G4ParticleDefinition* particle = trk->GetDefinition(); - const G4MaterialCutsCouple* couple = trk->GetMaterialCutsCouple(); -#if G4VERSION_NUMBER >= 1001 - G4EmSaturation* emSaturation = new G4EmSaturation(0); -#else - G4EmSaturation* emSaturation = new G4EmSaturation(); -#endif - double engyVis = emSaturation->VisibleEnergyDeposition(particle, - couple, - length, - energyDeposition, - niel); - delete emSaturation; - return engyVis; - } + double birkAttenuation() const; /// Set applyBirksLaw to ture void doApplyBirksLaw(void) { applyBirksLaw = true; diff --git a/DDG4/plugins/Geant4DetectorSensitivesConstruction.cpp b/DDG4/plugins/Geant4DetectorSensitivesConstruction.cpp index 70620c637..b4bf026b6 100644 --- a/DDG4/plugins/Geant4DetectorSensitivesConstruction.cpp +++ b/DDG4/plugins/Geant4DetectorSensitivesConstruction.cpp @@ -84,13 +84,13 @@ void Geant4DetectorSensitivesConstruction::constructSensitives(Geant4DetectorCon typedef Geometry::GeoHandlerTypes::SensitiveVolumes _SV; typedef Geometry::GeoHandlerTypes::ConstVolumeSet VolSet; Geant4GeometryInfo* p = Geant4Mapping::instance().ptr(); - G4VSensitiveDetector* g4sd = 0; _SV& vols = p->sensitives; for(_SV::const_iterator iv=vols.begin(); iv != vols.end(); ++iv) { Geometry::SensitiveDetector sd = (*iv).first; string typ = sd.type(), nam = sd.name(); - g4sd = PluginService::Create<G4VSensitiveDetector*>(typ, nam, &ctxt->lcdd); + G4VSensitiveDetector* g4sd = + PluginService::Create<G4VSensitiveDetector*>(typ, nam, &ctxt->lcdd); if (!g4sd) { string tmp = typ; tmp[0] = ::toupper(tmp[0]); diff --git a/DDG4/plugins/Geant4EventReaderHepMC.cpp b/DDG4/plugins/Geant4EventReaderHepMC.cpp index 6a2cf82cf..306158256 100644 --- a/DDG4/plugins/Geant4EventReaderHepMC.cpp +++ b/DDG4/plugins/Geant4EventReaderHepMC.cpp @@ -374,19 +374,18 @@ int HepMC::read_particle(EventStream &info, istringstream& input, Geant4Particle } int HepMC::read_vertex(EventStream &info, istream& is, istringstream & input) { - int id=0, dummy=0, num_orphans_in=0, num_particles_out=0, weights_size=0; + int id=0, dummy = 0, num_orphans_in=0, num_particles_out=0, weights_size=0; vector<float> weights; Geant4Vertex* v = new Geant4Vertex(); Geant4Particle* p; - if ( v ) { - input >> id >> dummy >> v->x >> v->y >> v->z >> v->time - >> num_orphans_in >> num_particles_out >> weights_size; - } - if(!input) { + if( !input ) { delete v; return 0; } + input >> id >> dummy >> v->x >> v->y >> v->z >> v->time + >> num_orphans_in >> num_particles_out >> weights_size; + v->x *= info.pos_unit; v->y *= info.pos_unit; v->z *= info.pos_unit; diff --git a/DDG4/plugins/Geant4TrackerWeightedSD.cpp b/DDG4/plugins/Geant4TrackerWeightedSD.cpp index 83f872fde..3d21129cc 100644 --- a/DDG4/plugins/Geant4TrackerWeightedSD.cpp +++ b/DDG4/plugins/Geant4TrackerWeightedSD.cpp @@ -62,6 +62,7 @@ namespace DD4hep { double distance_to_inside; double distance_to_outside; double mean_time; + double step_length; double e_cut; int current, parent; int combined; @@ -73,7 +74,7 @@ namespace DD4hep { bool single_deposit_mode; TrackerWeighted() : pre(), post(), sensitive(0), thisSD(0), distance_to_inside(0.0), distance_to_outside(0.0), mean_time(0.0), - e_cut(0.0), current(-1), parent(0), combined(0), + step_length(0.0), e_cut(0.0), current(-1), parent(0), combined(0), hit_position_type(POSITION_MIDDLE), hit_flag(0), g4ID(0), cell(0), single_deposit_mode(false) { @@ -85,6 +86,7 @@ namespace DD4hep { distance_to_inside = 0; distance_to_outside = 0; mean_time = 0; + step_length = 0; post.clear(); pre.clear(); current = -1; @@ -114,11 +116,14 @@ namespace DD4hep { /// Update energy and track information during hit info accumulation TrackerWeighted& update(const G4Step* step) { post.storePoint(step,step->GetPostStepPoint()); + Position mean = (post.position+pre.position)*0.5; + double mean_tm = (post.truth.time+pre.truth.time)*0.5; pre.truth.deposit += post.truth.deposit; - mean_pos.SetX(mean_pos.x()+post.position.x()*post.truth.deposit); - mean_pos.SetY(mean_pos.y()+post.position.y()*post.truth.deposit); - mean_pos.SetZ(mean_pos.z()+post.position.z()*post.truth.deposit); - mean_time += post.truth.time*post.truth.deposit; + mean_pos.SetX(mean_pos.x()+mean.x()*post.truth.deposit); + mean_pos.SetY(mean_pos.y()+mean.y()*post.truth.deposit); + mean_pos.SetZ(mean_pos.z()+mean.z()*post.truth.deposit); + mean_time += mean_tm*post.truth.deposit; + step_length += step->GetStepLength(); if ( 0 == cell ) { cell = sensitive->cellID(step); if ( 0 == cell ) { @@ -163,7 +168,6 @@ namespace DD4hep { Position pos; Momentum mom = 0.5 * (pre.momentum + post.momentum); double time = deposit != 0 ? mean_time / deposit : mean_time; - double path = (post.position - pre.position).R(); char dist_in[64], dist_out[64]; switch(hit_position_type) { @@ -195,7 +199,7 @@ namespace DD4hep { hit->flag = hit_flag; hit->position = pos; hit->momentum = mom; - hit->length = path; + hit->length = step_length; hit->cellID = cell; hit->g4ID = g4ID; diff --git a/DDG4/src/Geant4StepHandler.cpp b/DDG4/src/Geant4StepHandler.cpp index 2410a0635..bae7724ce 100644 --- a/DDG4/src/Geant4StepHandler.cpp +++ b/DDG4/src/Geant4StepHandler.cpp @@ -110,3 +110,25 @@ G4ThreeVector Geant4StepHandler::globalToLocalG4(const G4ThreeVector& global) c G4TouchableHandle t = step->GetPreStepPoint()->GetTouchableHandle(); return t->GetHistory()->GetTopTransform().TransformPoint(global); } + +/// Apply BirksLaw +double Geant4StepHandler::birkAttenuation() const { +#if G4VERSION_NUMBER >= 1001 + static G4EmSaturation s_emSaturation(0); +#else + static G4EmSaturation s_emSaturation(); +#endif + + double energyDeposition = step->GetTotalEnergyDeposit(); + double length = step->GetStepLength(); + double niel = step->GetNonIonizingEnergyDeposit(); + const G4Track* trk = step->GetTrack(); + const G4ParticleDefinition* particle = trk->GetDefinition(); + const G4MaterialCutsCouple* couple = trk->GetMaterialCutsCouple(); + double engyVis = s_emSaturation.VisibleEnergyDeposition(particle, + couple, + length, + energyDeposition, + niel); + return engyVis; +} -- GitLab