Skip to content

Commit

Permalink
Also save cumulative energy deposit for pasive volumes in the list
Browse files Browse the repository at this point in the history
  • Loading branch information
Sunanda committed Nov 9, 2017
1 parent 5ebfbed commit 962adfc
Show file tree
Hide file tree
Showing 6 changed files with 101 additions and 72 deletions.
11 changes: 7 additions & 4 deletions SimDataFormats/CaloHit/interface/PassiveHit.h
View file
Expand Up @@ -10,9 +10,10 @@ class PassiveHit {

public:

PassiveHit(std::string vname, unsigned int id, float e=0., float t=0.,
int it=0) : vname_(vname), id_(id), energy_(e), time_(t), it_(it) {}
PassiveHit() : vname_(""), id_(0), energy_(0), time_(0), it_(0) {}
PassiveHit(std::string vname, unsigned int id, float e=0, float etot=0,
float t=0, int it=0) : vname_(vname), id_(id), energy_(e),
etotal_(etot), time_(t), it_(it) {}
PassiveHit() : vname_(""), id_(0), energy_(0), etotal_(0), time_(0), it_(0) {}

//Names
static const char *name() { return "PassiveHit"; }
Expand All @@ -22,7 +23,8 @@ class PassiveHit {
//Energy deposit of the Hit
double energy() const { return energy_; }
void setEnergy(double e) { energy_ = e; }

double energyTotal() const { return etotal_; }
void setEnergyTotal(double e) { etotal_ = e; }

//Time of the deposit
double time() const { return time_; }
Expand All @@ -49,6 +51,7 @@ class PassiveHit {
std::string vname_;
unsigned int id_;
float energy_;
float etotal_;
float time_;
int it_;
};
Expand Down
5 changes: 2 additions & 3 deletions SimDataFormats/CaloHit/src/PassiveHit.cc
View file
Expand Up @@ -3,9 +3,8 @@

std::ostream & operator<<(std::ostream& o,const PassiveHit& hit) {
o << hit.vname() << " 0x" <<std::hex << hit.id() << std::dec
<< ": Energy " << hit.energy() << " GeV "
<< " Tof " << hit.time() << " ns "
<< " Track #" << hit.trackId();
<< ": Energy " << hit.energy() << " GeV: " << hit.energyTotal()
<< " GeV: Tof " << hit.time() << " ns: " << " Track # " << hit.trackId();

return o;
}
3 changes: 2 additions & 1 deletion SimDataFormats/CaloHit/src/classes_def.xml
View file
Expand Up @@ -28,7 +28,8 @@
<class name="edm::Wrapper<HFShowerLibraryEventInfo>" />
<class name="edm::Wrapper<std::vector<HFShowerLibraryEventInfo> >"/>
<typedef name="CastorShowerEvent::Point" />
<class name="PassiveHit" ClassVersion="11">
<class name="PassiveHit" ClassVersion="12">
<version ClassVersion="12" checksum="1168914655"/>
<version ClassVersion="11" checksum="2866279991"/>
</class>
<class name="std::vector<PassiveHit>"/>
Expand Down
140 changes: 82 additions & 58 deletions SimG4CMS/HGCalTestBeam/plugins/HGCPassive.cc
View file
Expand Up @@ -20,27 +20,26 @@
//#define EDM_ML_DEBUG

HGCPassive::HGCPassive(const edm::ParameterSet &p) : topPV_(nullptr), topLV_(nullptr),
count_(0), init_(false) {
count_(0), init_(false) {

edm::ParameterSet m_Passive = p.getParameter<edm::ParameterSet>("HGCPassive");
LVNames_ = m_Passive.getParameter<std::vector<std::string> >("LVNames");
motherName_= m_Passive.getParameter<std::string>("MotherName");

#ifdef EDM_ML_DEBUG
std::cout << "Name of the mother volume " << motherName_ << std::endl;
edm::LogVerbatim("ValidHGCal") << "Name of the mother volume " <<motherName_;
unsigned k(0);
#endif
for (auto name : LVNames_) {
produces<edm::PassiveHitContainer>(Form("%sPassiveHits",name.c_str()));
#ifdef EDM_ML_DEBUG
std::cout << "Collection name[" << k << "] " << name << std::endl;
edm::LogVerbatim("ValidHGCal") << "Collection name[" << k << "] " << name;
++k;
#endif
}
}

HGCPassive::~HGCPassive() {
}
HGCPassive::~HGCPassive() { }

void HGCPassive::produce(edm::Event& e, const edm::EventSetup&) {

Expand All @@ -64,11 +63,13 @@ void HGCPassive::update(const BeginOfRun * run) {
}

#ifdef EDM_ML_DEBUG
std::cout << "HGCPassive::Finds " << mapLV_.size() << " logical volumes\n";
edm::LogVerbatim("ValidHGCal") << "HGCPassive::Finds " << mapLV_.size()
<< " logical volumes";
unsigned int k(0);
for (const auto& lvs : mapLV_) {
std::cout << "Entry[" << k << "] " << lvs.first << ": ("
<< (lvs.second).first << ", " << (lvs.second).second << ")\n";
edm::LogVerbatim("ValidHGCal") << "Entry[" << k << "] " << lvs.first
<< ": (" << (lvs.second).first << ", "
<< (lvs.second).second << ")";
++k;
}
#endif
Expand All @@ -79,67 +80,84 @@ void HGCPassive::update(const BeginOfRun * run) {
void HGCPassive::update(const BeginOfEvent * evt) {

int iev = (*evt)()->GetEventID();
edm::LogInfo("ValidHGCal") << "HGCPassive: =====> Begin event = "
<< iev << std::endl;
edm::LogVerbatim("ValidHGCal") << "HGCPassive: =====> Begin event = "
<< iev << std::endl;

++count_;
store_.clear();
}

//=================================================================== each STEP
// //=================================================================== each STEP
void HGCPassive::update(const G4Step * aStep) {

if (aStep != nullptr) {

G4VSensitiveDetector* curSD = aStep->GetPreStepPoint()->GetSensitiveDetector();
const G4VTouchable* touchable = aStep->GetPreStepPoint()->GetTouchable();

if (curSD==nullptr) {

G4TouchableHistory* touchable = (G4TouchableHistory*)aStep->GetPreStepPoint()->GetTouchable();
G4LogicalVolume* plv = (G4LogicalVolume*)touchable->GetVolume()->GetLogicalVolume();
G4LogicalVolume* plv = touchable->GetVolume()->GetLogicalVolume();
auto it = (init_) ? mapLV_.find(plv) : findLV(plv);
double time = aStep->GetTrack()->GetGlobalTime();
double energy = (aStep->GetTotalEnergyDeposit())/CLHEP::GeV;

unsigned int copy(0);
if (((aStep->GetPostStepPoint() == nullptr) ||
(aStep->GetTrack()->GetNextVolume() == nullptr)) &&
(aStep->IsLastStepInVolume())) {
#ifdef EDM_ML_DEBUG
std::cout << plv->GetName() << " F|L Step "
<< aStep->IsFirstStepInVolume() << ":"
<< aStep->IsLastStepInVolume() << " Position"
<< aStep->GetPreStepPoint()->GetPosition() << " Track "
<< aStep->GetTrack()->GetDefinition()->GetParticleName()
<< " at" << aStep->GetTrack()->GetPosition() << " Volume "
<< aStep->GetTrack()->GetVolume() << ":"
<< aStep->GetTrack()->GetNextVolume() << " Status "
<< aStep->GetTrack()->GetTrackStatus() << " KE "
<< aStep->GetTrack()->GetKineticEnergy() << " Deposit "
<< aStep->GetTotalEnergyDeposit() << " Map "
<< (it != mapLV_.end()) << std::endl;
edm::LogVerbatim("ValidHGCal") << plv->GetName() << " F|L Step "
<< aStep->IsFirstStepInVolume() << ":"
<< aStep->IsLastStepInVolume()
<< " Position" << aStep->GetPreStepPoint()->GetPosition()
<< " Track " << aStep->GetTrack()->GetDefinition()->GetParticleName()
<< " at" << aStep->GetTrack()->GetPosition()
<< " Volume " << aStep->GetTrack()->GetVolume()
<< ":" << aStep->GetTrack()->GetNextVolume()
<< " Status " << aStep->GetTrack()->GetTrackStatus()
<< " KE " << aStep->GetTrack()->GetKineticEnergy()
<< " Deposit " << aStep->GetTotalEnergyDeposit()
<< " Map " << (it != mapLV_.end());
#endif
double time = aStep->GetTrack()->GetGlobalTime();
double energy = (aStep->GetPreStepPoint()->GetKineticEnergy() +
aStep->GetTotalEnergyDeposit())/CLHEP::GeV;
if (it != mapLV_.end()) {
storeInfo(it, plv, 0, time, energy);
} else if (topLV_ != nullptr) {
auto itr = (init_) ? mapLV_.find(topLV_) : findLV(topLV_);
if (itr != mapLV_.end()) {
storeInfo(itr, topLV_, 0, time, energy);
}
}
} else if (it != mapLV_.end()) {
unsigned int copy = (unsigned int)(touchable->GetReplicaNumber(0) +
1000*touchable->GetReplicaNumber(1));
double time = (aStep->GetPostStepPoint()->GetGlobalTime());
double edeposit = (aStep->GetTotalEnergyDeposit())/CLHEP::GeV;
storeInfo(it, plv, copy, time, edeposit);
energy += (aStep->GetPreStepPoint()->GetKineticEnergy()/CLHEP::GeV);
} else {
time = (aStep->GetPostStepPoint()->GetGlobalTime());
copy = (unsigned int)(touchable->GetReplicaNumber(0) +
1000*touchable->GetReplicaNumber(1));
}
if (it != mapLV_.end()) {
storeInfo(it, plv, copy, time, energy, true);
} else if (topLV_ != nullptr) {
auto itr = findLV(topLV_);
auto itr = (init_) ? mapLV_.find(topLV_) : findLV(topLV_);
if (itr != mapLV_.end()) {
double time = (aStep->GetPostStepPoint()->GetGlobalTime());
double edeposit = (aStep->GetTotalEnergyDeposit())/CLHEP::GeV;
storeInfo(itr, topLV_, 100, time, edeposit);
storeInfo(itr, topLV_, copy, time, energy, true);
}
}
}//if (curSD==NULL)

//Now for the mother volumes
int level = (touchable->GetHistoryDepth());
if (level > 0) {
double energy = (aStep->GetTotalEnergyDeposit())/CLHEP::GeV;
double time = (aStep->GetTrack()->GetGlobalTime());

for (int ii=0; ii<level; ++ii) {
int i = level - ii;
G4LogicalVolume* plv = touchable->GetVolume(i)->GetLogicalVolume();
auto it = (init_) ? mapLV_.find(plv) : findLV(plv);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("ValidHGCal") << "Level: " << ii << ":" << i << " "
<< plv->GetName() <<" flag in the List "
<< (it != mapLV_.end());
#endif
if (it != mapLV_.end()) {
unsigned int copy = (ii == level) ? 0 :
(unsigned int)(touchable->GetReplicaNumber(i) +
1000*touchable->GetReplicaNumber(i+1));
storeInfo(it, plv, copy, time, energy, false);
}
}
}
}//if (aStep != NULL)


Expand All @@ -157,11 +175,12 @@ void HGCPassive::endOfEvent(edm::PassiveHitContainer& hgcPH, unsigned int k) {
if (it != mapLV_.end()) {
if ((it->second).first == k) {
PassiveHit hit((it->second).second,(element.first).second,
(element.second).second,(element.second).first);
(element.second)[1],(element.second)[2],
(element.second)[0]);
hgcPH.push_back(hit);
#ifdef EDM_ML_DEBUG
std::cout << "HGCPassive[" << k << "] Hit[" << kount << "] " << hit
<< std::endl;
edm::LogVerbatim("ValidHGCal") << "HGCPassive[" << k << "] Hit["
<< kount << "] " << hit;
++kount;
#endif
}
Expand Down Expand Up @@ -193,20 +212,25 @@ HGCPassive::volumeIterator HGCPassive::findLV(G4LogicalVolume * plv) {

void HGCPassive::storeInfo(const HGCPassive::volumeIterator it,
G4LogicalVolume* plv, unsigned int copy,
double time, double energy) {
double time, double energy, bool flag) {

std::pair<G4LogicalVolume*,unsigned int> key(plv,copy);
auto itr = store_.find(key);
auto itr = store_.find(key);
double ee = (flag) ? energy : 0;
if (itr == store_.end()) {
store_[key] = std::pair<double,double>(time,energy);
store_[key] = { {time, energy, energy} };
} else {
(itr->second).second += energy;
(itr->second)[1] += ee;
(itr->second)[2] += energy;
}
#ifdef EDM_ML_DEBUG
std::cout << "HGCPassive: Element " << (it->second).first << ":"
<< (it->second).second << ":" << copy << " T "
<< (itr->second).first << " E " << (itr->second).second
<< std::endl;
itr = store_.find(key);
edm::LogVerbatim("ValidHGCal") << "HGCPassive: Element "
<< (it->second).first << ":"
<< (it->second).second << ":" << copy << " T "
<< (itr->second)[0] << " E "
<< (itr->second)[1] << ":"
<< (itr->second)[2];
#endif
}

Expand Down
9 changes: 5 additions & 4 deletions SimG4CMS/HGCalTestBeam/plugins/HGCPassive.h
View file
Expand Up @@ -29,9 +29,10 @@
#include "G4Track.hh"
#include "G4TouchableHistory.hh"

#include <vector>
#include <string>
#include <array>
#include <map>
#include <string>
#include <vector>

class HGCPassive : public SimProducer,
public Observer<const BeginOfRun *>,
Expand Down Expand Up @@ -61,7 +62,7 @@ class HGCPassive : public SimProducer,
G4VPhysicalVolume * getTopPV();
volumeIterator findLV(G4LogicalVolume * plv);
void storeInfo(const volumeIterator itr, G4LogicalVolume* plv,
unsigned int copy, double time, double energy);
unsigned int copy, double time, double energy, bool flag);

private:

Expand All @@ -74,7 +75,7 @@ class HGCPassive : public SimProducer,
// some private members for ananlysis
unsigned int count_;
bool init_;
std::map<std::pair<G4LogicalVolume*,unsigned int>,std::pair<double,double>> store_;
std::map<std::pair<G4LogicalVolume*,unsigned int>,std::array<double,3>> store_;
};


5 changes: 3 additions & 2 deletions SimG4CMS/HGCalTestBeam/test/HGCalTBCERN170_cfg.py
View file
Expand Up @@ -28,6 +28,7 @@
if 'MessageLogger' in process.__dict__:
process.MessageLogger.categories.append('HGCSim')
process.MessageLogger.categories.append('HcalSim')
process.MessageLogger.categories.append('ValidHGCal')

# Input source
process.source = cms.Source("EmptySource")
Expand Down Expand Up @@ -93,8 +94,8 @@
process.g4SimHits.HGCSD.RotatedWafer = True
process.g4SimHits.Watchers = cms.VPSet(cms.PSet(
HGCPassive = cms.PSet(
LVNames = cms.vstring('HGCalEE','HGCalHE','HGCalAH', 'CMSE'),
MotherName = cms.string('CMSE'),
LVNames = cms.vstring('HGCalEE','HGCalHE','HGCalAH', 'HGCalBeam', 'CMSE'),
MotherName = cms.string('OCMS'),
),
type = cms.string('HGCPassive'),
)
Expand Down

0 comments on commit 962adfc

Please sign in to comment.