Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Run3-sim46 Add step length information for passive hits #27964

Merged
merged 3 commits into from Sep 12, 2019
Merged

Run3-sim46 Add step length information for passive hits #27964

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ba7c2cc
Add step length information for passive hits
Sep 10, 2019
ada32c2
Code check
Sep 10, 2019
05f8f92
Put the macro as well
Sep 10, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
18 changes: 15 additions & 3 deletions SimDataFormats/CaloHit/interface/PassiveHit.h
View file
Open in desktop
Expand Up @@ -8,9 +8,16 @@

class PassiveHit {
public:
PassiveHit(std::string vname, unsigned int id, float e = 0, float etot = 0, float t = 0, int it = 0, int ip = 0)
: vname_(vname), id_(id), energy_(e), etotal_(etot), time_(t), it_(it), ip_(ip) {}
PassiveHit() : vname_(""), id_(0), energy_(0), etotal_(0), time_(0), it_(0), ip_(0) {}
PassiveHit(std::string vname,
unsigned int id,
float e = 0,
float etot = 0,
float t = 0,
int it = 0,
int ip = 0,
float stepl = 0)
: vname_(vname), id_(id), energy_(e), etotal_(etot), time_(t), it_(it), ip_(ip), stepl_(stepl) {}
PassiveHit() : vname_(""), id_(0), energy_(0), etotal_(0), time_(0), it_(0), ip_(0), stepl_(0) {}

//Names
static const char *name() { return "PassiveHit"; }
Expand Down Expand Up @@ -43,6 +50,10 @@ class PassiveHit {
int pdgId() const { return ip_; }
void setPDGId(int ip) { ip_ = ip; }

//Step length for the current Hit
float stepLength() const { return stepl_; }
void setStepLength(float stepl) { stepl_ = stepl; }

//Comparisons
bool operator<(const PassiveHit &d) const { return energy_ < d.energy_; }

Expand All @@ -57,6 +68,7 @@ class PassiveHit {
float time_;
int it_;
int ip_;
float stepl_;
};

namespace edm {
Expand Down
3 changes: 2 additions & 1 deletion SimDataFormats/CaloHit/src/classes_def.xml
View file
Open in desktop
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="13">
<class name="PassiveHit" ClassVersion="14">
<version ClassVersion="14" checksum="3344618767"/>
<version ClassVersion="13" checksum="1735297430"/>
<version ClassVersion="12" checksum="1168914655"/>
<version ClassVersion="11" checksum="2866279991"/>
Expand Down
5 changes: 2 additions & 3 deletions SimG4CMS/Calo/interface/CaloSteppingAction.h
View file
Open in desktop
Expand Up @@ -100,9 +100,8 @@ class CaloSteppingAction : public SimProducer,
double birkCutEC_, birkC1HC_, birkC2HC_;
double birkC3HC_, timeSliceUnit_;
std::map<std::pair<int, CaloHitID>, CaloGVHit> hitMap_[nSD_];
typedef std::tuple<const G4LogicalVolume *, uint32_t, int, int> PassiveKey;
typedef std::tuple<int, double, double, double> PassiveData;
std::map<PassiveKey, PassiveData> store_;
typedef std::tuple<const G4LogicalVolume *, uint32_t, int, int, double, double, double, double> PassiveData;
std::vector<PassiveData> store_;
};

#endif
61 changes: 43 additions & 18 deletions SimG4CMS/Calo/macros/MakeGVPlots.C
View file
Open in desktop
Expand Up @@ -9,7 +9,7 @@
// To plot on the same Canvas plots from Geant4 and GeantV done using
// CaloSimHitAnalysis
// makeGV2Plots(fnmG4, fnmGV, todomin, todomax, normalize, tag, text,
// save, dirnm)
// save, mode, dirnm)
//
// To plot from passive hit collection produced using CaloSimHitAnalysis
// makeGVSPlots(fnmG4, fnmGV, todomin, todomax, tag, text, save, dirnm)
Expand All @@ -23,10 +23,12 @@
// ifGV bool If the file created using Geant4 (true)
// todomin int The first one in the series to be created (0)
// todomax int The last one in the series to be created
// (3:2:5 for GVPlots:GV2Plots:GVSPlots)
// (3:2:7 for GVPlots:GV2Plots:GVSPlots)
// tag std::string To be added to the name of the canvas ("")
// text std::string To be added to the title of the histogram ("")
// save bool If the canvas is to be saved as jpg file (false)
// mode int Flag if one file is G4 and other GV (0), or
// both files are GV (1) or G4 (2) (defualt 0)
// dirnm std::string Name of the directory ("caloSimHitAnalysis")
//
// The histogram series have different meanings for each function
Expand All @@ -46,9 +48,10 @@
// "Energy deposit for EM particles", "Energy deposit for non-EM particles",
// "R", "Z", "#eta", "phi"
//
// GVSPlots (6 in total):
// "Total Hits", "Tracks", "Energy Deposit in all components",
// "Time of all hits", "Energy Deposit in tracker subdetectors",
// GVSPlots (8 in total):
// "Total Hits", "Tracks", "Step Length',
// "Energy Deposit in all components", "Time of all hits",
// "Energy Deposit in tracker subdetectors",
// "Time of hits in tracker subdetectors"
//
// All plots in GVPlots or GV2Plots happen for EB, EE, HB and HE
Expand Down Expand Up @@ -172,7 +175,7 @@ void makeGVPlots(std::string fname="analG4.root", bool ifG4=true,
void makeGV2Plots(std::string fnmG4="analG4.root",
std::string fnmGV="analGV.root", int todomin=0,
int todomax=2, bool normalize=true, std::string tag="",
std::string text="", bool save=false,
std::string text="", bool save=false, int mode=0,
std::string dirnm="caloSimHitAnalysis") {

std::string names[13] = {"Edep", "Time", "Etot", "Edep15", "EdepT", "EdepT15",
Expand All @@ -197,7 +200,7 @@ void makeGV2Plots(std::string fnmG4="analG4.root",
gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite);

if (normalize) gStyle->SetOptStat(0);
else gStyle->SetOptStat(111100);
else gStyle->SetOptStat(111110);
TFile *file1 = new TFile(fnmG4.c_str());
TFile *file2 = new TFile(fnmGV.c_str());
if (file1 && file2) {
Expand All @@ -208,7 +211,12 @@ void makeGV2Plots(std::string fnmG4="analG4.root",
for (int i2=todomin; i2<=todomax; ++i2) {
sprintf (name, "%s%d", names[i2].c_str(), dets[i1]);
sprintf (cname, "%s%d%s", names[i2].c_str(), dets[i1], tag.c_str());
sprintf (title, "%s %s (Geant4 vs GeantV)", text.c_str(), detName[i1].c_str());
if (mode == 1)
sprintf (title, "%s %s (2 runs of GeantV)", text.c_str(), detName[i1].c_str());
else if (mode == 2)
sprintf (title, "%s %s (2 runs of Geant4)", text.c_str(), detName[i1].c_str());
else
sprintf (title, "%s %s (Geant4 vs GeantV)", text.c_str(), detName[i1].c_str());
TH1D *hist[2];
hist[0] = (TH1D*)dir1->FindObjectAny(name);
hist[1] = (TH1D*)dir2->FindObjectAny(name);
Expand All @@ -235,7 +243,13 @@ void makeGV2Plots(std::string fnmG4="analG4.root",
hist[i]->SetLineStyle(isty[i]);
hist[i]->SetNdivisions(505,"X");
total[i] = hist[i]->GetEntries();
legend->AddEntry(hist[i],type[i].c_str(),"lp");
if (mode == 1)
sprintf (name, "%s (run %d)", type[1].c_str(), i);
else if (mode == 2)
sprintf (name, "%s (run %d)", type[0].c_str(), i);
else
sprintf (name, "%s", type[i].c_str());
legend->AddEntry(hist[i], name, "lp");
if (i == 0) {
if (normalize) hist[i]->DrawNormalized("hist");
else hist[i]->Draw();
Expand Down Expand Up @@ -335,18 +349,23 @@ void makeGV2Plots(std::string fnmG4="analG4.root",

void makeGVSPlots(std::string fnmG4="analG4.root",
std::string fnmGV="analGV.root",
int todomin=0, int todomax=5,
int todomin=0, int todomax=7,
std::string tag="", std::string text="", bool save=false,
std::string dirnm="caloSimHitAnalysis") {

std::string names[6] = {"hitp", "trackp", "edepp", "timep", "edept", "timet"};
std::string xtitle[6] = {"Hits", "Tracks", "Energy Deposit (MeV)",
"Time (ns)", "Energy Deposit (MeV)", "Time (ns)"};
std::string ytitle[6] = {"Events", "Events", "Hits", "Hits", "Hits", "Hits"};
std::string names[8] = {"hitp", "trackp", "stepp", "edepp", "timep", "volp",
"edept", "timet"};
std::string xtitle[8] = {"Hits", "Tracks", "Step Length (cm)",
"Energy Deposit (MeV)", "Time (ns)",
"Volume elements", "Energy Deposit (MeV)",
"Time (ns)"};
std::string ytitle[8] = {"Events", "Events", "Hits", "Hits", "Hits",
"Events", "Hits", "Hits"};
std::string detName[6] = {"Barrel Pixel", "Forward Pixel", "TIB", "TID",
"TOB", "TEC"};
int boxp[6] = {0, 1, 0, 0, 0, 0};
int nhis[6] = {1, 1, 1, 1, 6, 6};
std::string particles[3] = {"Electrons", "Positrons", "Photons"};
int boxp[8] = {0, 1, 0, 0, 0, 0, 0, 0};
int nhis[8] = {-4, -4, -4, 1, 1, 1, 6, 6};

gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite);
Expand All @@ -358,11 +377,17 @@ void makeGVSPlots(std::string fnmG4="analG4.root",
TDirectory *dir2 = (TDirectory*)file2->FindObjectAny(dirnm.c_str());
char name[100], cname[100], title[100];
for (int i2=todomin; i2<=todomax; ++i2) {
for (int i1=0; i1<nhis[i2]; ++ i1) {
if (nhis[i2] == 1) {
int nhism = (nhis[i2] >= 0) ? nhis[i2] : -nhis[i2];
for (int i1=0; i1<nhism; ++ i1) {
if (nhis[i2] <= 1 && i1 == 0) {
sprintf (name, "%s", names[i2].c_str());
sprintf (cname, "%s%s", names[i2].c_str(), tag.c_str());
sprintf (title, "%s (Geant4 vs GeantV)", text.c_str());
} else if (nhis[i2] < 0) {
sprintf (name, "%s%d", names[i2].c_str(), i1-1);
sprintf (cname, "%s%d%s", names[i2].c_str(), i1-1, tag.c_str());
sprintf (title, "%s in %s (Geant4 vs GeantV)",
particles[i1-1].c_str(), text.c_str());
} else {
sprintf (name, "%s%d", names[i2].c_str(), i1);
sprintf (cname, "%s%d%s", names[i2].c_str(), i1, tag.c_str());
Expand Down
10 changes: 7 additions & 3 deletions SimG4CMS/Calo/plugins/CaloSimHitAnalysis.cc
View file
Open in desktop
Expand Up @@ -76,7 +76,7 @@ class CaloSimHitAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm:
TH1F *h_edepEM_[nCalo_], *h_edepHad_[nCalo_], *h_rr_[nCalo_], *h_zz_[nCalo_];
TH1F *h_eta_[nCalo_], *h_phi_[nCalo_], *h_etot_[nCalo_], *h_etotg_[nCalo_];
TH2F *h_rz_, *h_rz1_, *h_etaphi_;
TH1F *h_hitp_, *h_trackp_, *h_edepp_, *h_timep_;
TH1F *h_hitp_, *h_trackp_, *h_edepp_, *h_timep_, *h_stepp_;
std::vector<TH1F*> h_edepTk_, h_timeTk_;
};

Expand Down Expand Up @@ -213,8 +213,11 @@ CaloSimHitAnalysis::CaloSimHitAnalysis(const edm::ParameterSet& ps)
h_edepp_->GetXaxis()->SetTitle("Energy Deposit (MeV)");
h_edepp_->GetYaxis()->SetTitle("Hits");
h_timep_ = tfile->make<TH1F>("timep", "All Steps", 100, 0.0, 100.0);
h_hitp_->GetXaxis()->SetTitle("Hits");
h_hitp_->GetYaxis()->SetTitle("Hit Time (ns)");
h_timep_->GetXaxis()->SetTitle("Hits");
h_timep_->GetYaxis()->SetTitle("Hit Time (ns)");
h_stepp_ = tfile->make<TH1F>("stepp", "All Steps", 1000, 0.0, 100.0);
h_stepp_->GetXaxis()->SetTitle("Hits");
h_stepp_->GetYaxis()->SetTitle("Step length (cm)");
for (unsigned int k = 0; k < detNames_.size(); ++k) {
sprintf(name, "edept%d", k);
sprintf(title, "Energy Deposit (MeV) in %s", detNames_[k].c_str());
Expand Down Expand Up @@ -449,6 +452,7 @@ void CaloSimHitAnalysis::analyzePassiveHits(std::vector<PassiveHit>& hits) {
++(ktr->second);
h_edepp_->Fill(hit.energy());
h_timep_->Fill(hit.time());
h_stepp_->Fill(hit.stepLength());
if ((name.find(active) != std::string::npos) || (name.find(sensor) != std::string::npos)) {
unsigned idet = detNames_.size();
for (unsigned int k = 0; k < detNames_.size(); ++k) {
Expand Down
52 changes: 13 additions & 39 deletions SimG4CMS/Calo/src/CaloSteppingAction.cc
View file
Open in desktop
Expand Up @@ -109,16 +109,17 @@ void CaloSteppingAction::fillHits(edm::PCaloHitContainer& cc, int type) {
void CaloSteppingAction::fillPassiveHits(edm::PassiveHitContainer& cc) {
edm::LogVerbatim("Step") << "CaloSteppingAction::fillPassiveHits with " << store_.size() << " hits";
for (const auto& element : store_) {
auto lv = std::get<0>(element.first);
auto lv = std::get<0>(element);
auto it = mapLV_.find(lv);
if (it != mapLV_.end()) {
PassiveHit hit(it->second,
std::get<1>(element.first),
std::get<2>(element.second),
std::get<3>(element.second),
std::get<1>(element.second),
std::get<2>(element.first),
std::get<0>(element.second));
std::get<1>(element),
std::get<5>(element),
std::get<6>(element),
std::get<4>(element),
std::get<2>(element),
std::get<3>(element),
std::get<7>(element));
cc.emplace_back(hit);
}
}
Expand Down Expand Up @@ -206,7 +207,7 @@ void CaloSteppingAction::update(const BeginOfEvent* evt) {
slave_[k].get()->Initialize();
}
if (allSteps_ > 0)
store_.erase(store_.begin(), store_.end());
store_.clear();
}

//=================================================================== each STEP
Expand Down Expand Up @@ -274,14 +275,14 @@ void CaloSteppingAction::update(const G4Step* aStep) {
double energy = aStep->GetTotalEnergyDeposit() / CLHEP::MeV;
auto const touch = aStep->GetPreStepPoint()->GetTouchable();
double time = aStep->GetTrack()->GetGlobalTime() / CLHEP::nanosecond;
int itime = static_cast<int>(1000.0 * time);
int trackId = aStep->GetTrack()->GetTrackID();
int pdg = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
double stepl = (aStep->GetStepLength() / CLHEP::cm);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("Step") << "CaloSteppingAction: Volume " << lv->GetName() << " History "
<< touch->GetHistoryDepth() << " Pointers " << aStep->GetPostStepPoint() << ":"
<< aStep->GetTrack()->GetNextVolume() << ":" << aStep->IsLastStepInVolume() << " E "
<< energy << " T " << time << ":" << itime << " PDG " << pdg;
<< energy << " T " << time << " PDG " << pdg << " step " << stepl;
#endif
uint32_t copy = 0;
if (((aStep->GetPostStepPoint() == nullptr) || (aStep->GetTrack()->GetNextVolume() == nullptr)) &&
Expand All @@ -294,35 +295,8 @@ void CaloSteppingAction::update(const G4Step* aStep) {
: static_cast<uint32_t>(touch->GetReplicaNumber(0) + 1000 * touch->GetReplicaNumber(1));
}
if (it != mapLV_.end()) {
PassiveKey key(std::make_tuple(lv, copy, trackId, itime));
if (allSteps_ > 1) {
bool flag(true);
while (store_.find(key) != store_.end()) {
copy += 1000000;
key = std::make_tuple(lv, copy, trackId, itime);
if (copy > 100000000) {
flag = false;
break;
}
}
if (flag) {
store_[key] = std::make_tuple(pdg, time, energy, energy);
} else {
auto itr = store_.find(key);
double e1 = std::get<2>(itr->second) + energy;
double e2 = std::get<3>(itr->second) + energy;
store_[key] = std::make_tuple(pdg, time, e1, e2);
}
} else {
auto itr = store_.find(key);
if (itr == store_.end()) {
store_[key] = std::make_tuple(pdg, time, energy, energy);
} else {
double e1 = std::get<2>(itr->second) + energy;
double e2 = std::get<3>(itr->second) + energy;
store_[key] = std::make_tuple(pdg, time, e1, e2);
}
}
PassiveData key(std::make_tuple(lv, copy, trackId, pdg, time, energy, energy, stepl));
store_.push_back(key);
}
}
}
Expand Down