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

Run2-sim08 Enabling depth studies in Ecal crystal #19099

Merged
merged 1 commit into from Jun 7, 2017
Merged

Run2-sim08 Enabling depth studies in Ecal crystal #19099

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
2 changes: 1 addition & 1 deletion SimG4CMS/Calo/interface/ECalSD.h
View file
Open in desktop
Expand Up @@ -53,7 +53,7 @@ class ECalSD : public CaloSD {
bool useWeight, storeTrack, storeRL, storeLayerTimeSim;
bool useBirk, useBirkL3;
double birk1, birk2, birk3, birkSlope, birkCut;
double slopeLY;
double slopeLY, scaleRL;
std::string crystalMat, depth1Name, depth2Name;
std::map<G4LogicalVolume*,double> xtalLMap;
std::vector<G4LogicalVolume*> useDepth1, useDepth2, noWeight;
Expand Down
7 changes: 4 additions & 3 deletions SimG4CMS/Calo/python/CaloSimHitStudy_cfi.py
View file
Open in desktop
Expand Up @@ -7,7 +7,8 @@
EECollection = cms.untracked.string('EcalHitsEE'),
ESCollection = cms.untracked.string('EcalHitsES'),
HCCollection = cms.untracked.string('HcalHits'),
MaxEnergy = cms.untracked.double(200.0),
TimeCut = cms.untracked.double(100.0),
MIPCut = cms.untracked.double(0.70)
MaxEnergy = cms.untracked.double(200.0),
TimeCut = cms.untracked.double(100.0),
MIPCut = cms.untracked.double(0.70),
StoreRL = cms.untracked.bool(False)
)
90 changes: 55 additions & 35 deletions SimG4CMS/Calo/src/ECalSD.cc
View file
Open in desktop
Expand Up @@ -28,7 +28,7 @@

#include<algorithm>

//#define DebugLog
//#define EDM_ML_DEBUG

template <class T>
bool any(const std::vector<T> & v, const T &what) {
Expand Down Expand Up @@ -66,6 +66,7 @@ ECalSD::ECalSD(G4String name, const DDCompactView & cpv,
bool isItTB = m_EC.getUntrackedParameter<bool>("TestBeam", false);
bool nullNS = m_EC.getUntrackedParameter<bool>("NullNumbering", false);
storeRL = m_EC.getUntrackedParameter<bool>("StoreRadLength", false);
scaleRL = m_EC.getUntrackedParameter<double>("ScaleRadLength",1.0);

//Changes for improved timing simulation
storeLayerTimeSim = m_EC.getUntrackedParameter<bool>("StoreLayerTimeSim", false);
Expand All @@ -84,13 +85,20 @@ ECalSD::ECalSD(G4String name, const DDCompactView & cpv,
useWeight= true;
std::vector<double> tempD = getDDDArray("EnergyWeight",sv);
if (tempD.size() > 0) { if (tempD[0] < 0.1) useWeight = false; }
#ifdef EDM_ML_DEBUG
edm::LogInfo("EcalSim") << "ECalSD:: useWeight " << tempD.size() << ":"
<< useWeight << std::endl;
#endif
std::vector<std::string> tempS = getStringArray("Depth1Name",sv);
if (tempS.size() > 0) depth1Name = tempS[0];
else depth1Name = " ";
tempS = getStringArray("Depth2Name",sv);
if (tempS.size() > 0) depth2Name = tempS[0];
else depth2Name = " ";

#ifdef EDM_ML_DEBUG
edm::LogInfo("EcalSim") << "Names (Depth 1):" << depth1Name << " (Depth 2):"
<< depth2Name << std::endl;
#endif
EcalNumberingScheme* scheme=0;
if (nullNS) scheme = 0;
else if (name == "EcalHitsEB")
Expand All @@ -113,7 +121,7 @@ ECalSD::ECalSD(G4String name, const DDCompactView & cpv,


if (scheme) setNumberingScheme(scheme);
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "Constructing a ECalSD with name " << GetName();
#endif
if (useWeight) {
Expand All @@ -136,7 +144,7 @@ ECalSD::ECalSD(G4String name, const DDCompactView & cpv,
<< "\tions below " << kmaxIon << " MeV"
<< "\n\tDepth1 Name = " << depth1Name
<< "\tDepth2 Name = " << depth2Name
<< "\n\tstoreRL" << storeRL
<< "\n\tstoreRL" << storeRL << ":" << scaleRL
<< "\tstoreLayerTimeSim " << storeLayerTimeSim
<< "\n\ttime Granularity " << p.getParameter<edm::ParameterSet>("ECalSD").getParameter<double>("TimeSliceUnit") << " ns";
if (useWeight) initMap(name,cpv);
Expand Down Expand Up @@ -168,7 +176,7 @@ double ECalSD::getEnergyDeposit(G4Step * aStep) {
((pdg/10)%100) > 0)) && (ke<kmaxIon)) weight = 0;
if ((pdg == 2212) && (ke < kmaxProton)) weight = 0;
if ((pdg == 2112) && (ke < kmaxNeutron)) weight = 0;
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
if (weight == 0)
LogDebug("EcalSim") << "Ignore Track " << theTrack->GetTrackID()
<< " Type " << theTrack->GetDefinition()->GetParticleName()
Expand Down Expand Up @@ -204,7 +212,7 @@ double ECalSD::getEnergyDeposit(G4Step * aStep) {
}
*/
if(wt2 > 0.0) { edep *= wt2; }
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD:: " << nameVolume
<<" Light Collection Efficiency " <<weight << ":" <<wt1
<< " Weighted Energy Deposit " << edep/MeV << " MeV";
Expand Down Expand Up @@ -236,11 +244,15 @@ int ECalSD::getTrackID(G4Track* aTrack) {

uint16_t ECalSD::getDepth(G4Step * aStep) {
G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetTouchable()->GetVolume(0)->GetLogicalVolume();
if (any(useDepth1,lv)) return 1;
else if (any(useDepth2,lv)) return 2;
else if (storeRL) return getRadiationLength(aStep);
else if (storeLayerTimeSim) return getLayerIDForTimeSim(aStep);
return 0;
if (storeRL) {
return getRadiationLength(aStep);
} else if (storeLayerTimeSim) {
return getLayerIDForTimeSim(aStep);
} else {
if (any(useDepth1,lv)) return 1;
else if (any(useDepth2,lv)) return 2;
else return 0;
}
}

uint16_t ECalSD::getRadiationLength(G4Step * aStep) {
Expand All @@ -249,14 +261,20 @@ uint16_t ECalSD::getRadiationLength(G4Step * aStep) {
if (aStep != NULL) {
G4StepPoint* hitPoint = aStep->GetPreStepPoint();
G4LogicalVolume* lv = hitPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume();

#ifdef EDM_ML_DEBUG
edm::LogInfo("EcalSim") << lv->GetName() << " useWight " << useWeight;
#endif
if (useWeight) {
G4ThreeVector localPoint = setToLocal(hitPoint->GetPosition(),
hitPoint->GetTouchable());
double crlength = crystalLength(lv);
double radl = hitPoint->GetMaterial()->GetRadlen();
double detz = (float)(0.5*crlength + localPoint.z());
thisX0 = (uint16_t)floor(detz/radl);
thisX0 = (uint16_t)floor(scaleRL*detz/radl);
#ifdef EDM_ML_DEBUG
edm::LogInfo("EcalSim") << "Crystal Length " << crlength << " Radl "
<< radl << " DetZ " << detz << " " << thisX0;
#endif
}
}
return thisX0;
Expand Down Expand Up @@ -342,17 +360,18 @@ void ECalSD::initMap(G4String sd, const DDCompactView & cpv) {
if (strncmp(lvname.c_str(), depth1Name.c_str(), 4) == 0) {
if (!any(useDepth1, lv)) {
useDepth1.push_back(lv);
#ifdef DebugLog
LogDebug("EcalSim") << "ECalSD::initMap Logical Volume " << lvname
<<" in Depth 1 volume list";
#ifdef EDM_ML_DEBUG
edm::LogInfo("EcalSim") << "ECalSD::initMap Logical Volume "
<< lvname <<" in Depth 1 volume list\n";
#endif
}
G4LogicalVolume* lvr = nameMap[lvname + "_refl"];
if (lvr != 0 && !any(useDepth1, lvr)) {
useDepth1.push_back(lvr);
#ifdef DebugLog
LogDebug("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << "_refl"
<<" in Depth 1 volume list";
#ifdef EDM_ML_DEBUG
edm::LogInfo("EcalSim") << "ECalSD::initMap Logical Volume "
<< lvname << "_refl"
<<" in Depth 1 volume list\n";
#endif
}
}
Expand All @@ -361,17 +380,18 @@ void ECalSD::initMap(G4String sd, const DDCompactView & cpv) {
if (strncmp(lvname.c_str(), depth2Name.c_str(), 4) == 0) {
if (!any(useDepth2, lv)) {
useDepth2.push_back(lv);
#ifdef DebugLog
LogDebug("EcalSim") << "ECalSD::initMap Logical Volume " << lvname
<<" in Depth 2 volume list";
#ifdef EDM_ML_DEBUG
edm::LogInfo("EcalSim") << "ECalSD::initMap Logical Volume "
<< lvname <<" in Depth 2 volume list\n";
#endif
}
G4LogicalVolume* lvr = nameMap[lvname + "_refl"];
if (lvr != 0 && !any(useDepth2,lvr)) {
useDepth2.push_back(lvr);
#ifdef DebugLog
LogDebug("EcalSim") << "ECalSD::initMap Logical Volume " << lvname << "_refl"
<<" in Depth 2 volume list";
#ifdef EDM_ML_DEBUG
edm::LogInfo("EcalSim") << "ECalSD::initMap Logical Volume "
<< lvname << "_refl"
<<" in Depth 2 volume list\n";
#endif
}
}
Expand All @@ -382,7 +402,7 @@ void ECalSD::initMap(G4String sd, const DDCompactView & cpv) {
lvused.push_back(lv);
const DDSolid & sol = fv.logicalPart().solid();
const std::vector<double> & paras = sol.parameters();
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD::initMap (for " << sd << "): Solid "
<< lvname << " Shape " << sol.shape()
<< " Parameter 0 = "<< paras[0]
Expand All @@ -399,15 +419,15 @@ void ECalSD::initMap(G4String sd, const DDCompactView & cpv) {
} else {
if (!any(noWeight,lv)) {
noWeight.push_back(lv);
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD::initMap Logical Volume " << lvname
<< " Material " << matname <<" in noWeight list";
#endif
}
lv = nameMap[lvname];
if (lv != 0 && !any(noWeight,lv)) {
noWeight.push_back(lv);
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD::initMap Logical Volume " << lvname
<< " Material " << matname <<" in noWeight list";
#endif
Expand All @@ -416,7 +436,7 @@ void ECalSD::initMap(G4String sd, const DDCompactView & cpv) {
}
dodet = fv.next();
}
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD: Length Table for " << attribute << " = "
<< sd << ":";
std::map<G4LogicalVolume*,double>::const_iterator ite = xtalLMap.begin();
Expand Down Expand Up @@ -482,7 +502,7 @@ void ECalSD::getBaseNumber(const G4Step* aStep) {
if ( theSize > 1 ) {
for (int ii = 0; ii < theSize ; ii++) {
theBaseNumber.addLevel(touch->GetVolume(ii)->GetName(),touch->GetReplicaNumber(ii));
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD::getBaseNumber(): Adding level " << ii
<< ": " << touch->GetVolume(ii)->GetName() << "["
<< touch->GetReplicaNumber(ii) << "]";
Expand All @@ -507,7 +527,7 @@ double ECalSD::getBirkL3(G4Step* aStep) {
if (weight < birkCut) weight = birkCut;
else if (weight > 1.) weight = 1.;
}
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD::getBirkL3 in " << mat->GetName()
<< " Charge " << charge << " dE/dx " << dedx
<< " Birk Const " << rkb << " Weight = " << weight
Expand All @@ -521,12 +541,12 @@ double ECalSD::getBirkL3(G4Step* aStep) {
std::vector<double> ECalSD::getDDDArray(const std::string & str,
const DDsvalues_type & sv) {

#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD:getDDDArray called for " << str;
#endif
DDValue value(str);
if (DDfetch(&sv,value)) {
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << value;
#endif
const std::vector<double> & fvec = value.doubles();
Expand All @@ -540,12 +560,12 @@ std::vector<double> ECalSD::getDDDArray(const std::string & str,
std::vector<std::string> ECalSD::getStringArray(const std::string & str,
const DDsvalues_type & sv) {

#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << "ECalSD:getStringArray called for " << str;
#endif
DDValue value(str);
if (DDfetch(&sv,value)) {
#ifdef DebugLog
#ifdef EDM_ML_DEBUG
LogDebug("EcalSim") << value;
#endif
const std::vector<std::string> & fvec = value.strings();
Expand Down
7 changes: 5 additions & 2 deletions SimG4CMS/Calo/test/CaloSimHitStudy.cc
View file
Open in desktop
Expand Up @@ -17,6 +17,7 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
double maxEnergy_= ps.getUntrackedParameter<double>("MaxEnergy", 200.0);
tmax_ = ps.getUntrackedParameter<double>("TimeCut", 100.0);
eMIP_ = ps.getUntrackedParameter<double>("MIPCut", 0.70);
storeRL_ = ps.getUntrackedParameter<bool>("StoreRL", false);

muonLab[0] = "MuonRPCHits";
muonLab[1] = "MuonCSCHits";
Expand Down Expand Up @@ -267,8 +268,10 @@ void CaloSimHitStudy::analyzeHits (std::vector<PCaloHit>& hits, int indx) {
}
int idx = -1;
if (indx != 3) {
if (indx == 0) idx = hits[i].depth();
else idx = indx+2;
if (indx == 0) {
if (storeRL_) idx = 0;
else idx = hits[i].depth();
} else idx = indx+2;
time_[idx]->Fill(time);
edep_[idx]->Fill(edep);
edepEM_[idx]->Fill(edepEM);
Expand Down
1 change: 1 addition & 0 deletions SimG4CMS/Calo/test/CaloSimHitStudy.h
View file
Open in desktop
Expand Up @@ -58,6 +58,7 @@ class CaloSimHitStudy: public edm::EDAnalyzer{
edm::EDGetTokenT<edm::PSimHitContainer> toks_tkLow_[6];
std::string muonLab[3], tkHighLab[6], tkLowLab[6];
double tmax_, eMIP_;
bool storeRL_;

TH1F *hit_[9], *time_[9], *edepEM_[9], *edepHad_[9], *edep_[9];
TH1F *etot_[9], *etotg_[9], *timeAll_[9], *hitMu, *hitHigh;
Expand Down
43 changes: 23 additions & 20 deletions SimG4CMS/Calo/test/python/runWithGun_cfg.py
View file
Open in desktop
Expand Up @@ -27,7 +27,7 @@
cout = cms.untracked.PSet(
# threshold = cms.untracked.string('DEBUG'),
INFO = cms.untracked.PSet(
limit = cms.untracked.int32(-1)
limit = cms.untracked.int32(0)
),
DEBUG = cms.untracked.PSet(
limit = cms.untracked.int32(0)
Expand All @@ -39,7 +39,7 @@
limit = cms.untracked.int32(-1)
),
SimTrackManager = cms.untracked.PSet(
limit = cms.untracked.int32(-1)
limit = cms.untracked.int32(0)
),
SimG4CoreApplication = cms.untracked.PSet(
limit = cms.untracked.int32(0)
Expand All @@ -63,7 +63,7 @@
limit = cms.untracked.int32(0)
),
HcalSim = cms.untracked.PSet(
limit = cms.untracked.int32(0)
limit = cms.untracked.int32(-1)
)
)
)
Expand Down Expand Up @@ -110,10 +110,10 @@
ignoreTotal = cms.untracked.int32(1)
)

process.Tracer = cms.Service("Tracer")
#process.Tracer = cms.Service("Tracer")

process.TFileService = cms.Service("TFileService",
fileName = cms.string('runWithGun_QGSP_FTFP_BERT_EML.root')
fileName = cms.string('runWithGun_FTFP_BERT_EMM.root')
)

process.generation_step = cms.Path(process.pgen)
Expand All @@ -122,17 +122,20 @@
process.out_step = cms.EndPath(process.output)

process.caloSimHitStudy.MaxEnergy = 1000.0
#process.g4SimHits.Physics.type = 'SimG4Core/Physics/QGSP_FTFP_BERT_EML'
process.g4SimHits.Physics.type = 'SimG4Core/Physics/FTFP_BERT_EMM'
process.g4SimHits.Physics.MonopoleCharge = 1
process.g4SimHits.Physics.Verbosity = 0
process.g4SimHits.CaloSD.UseResponseTables = [1,1,0,1]
process.g4SimHits.CaloSD.EminHits[0] = 0
process.g4SimHits.ECalSD.StoreSecondary = True
process.g4SimHits.ECalSD.StoreRadLength = True
process.g4SimHits.ECalSD.ScaleRadLength = 100.0
process.g4SimHits.CaloTrkProcessing.PutHistory = True
process.g4SimHits.CaloResponse.UseResponseTable = True
process.g4SimHits.CaloResponse.ResponseScale = 1.0
process.g4SimHits.CaloResponse.ResponseFile = 'SimG4CMS/Calo/data/responsTBpim50.dat'
process.g4SimHits.G4Commands = ['/run/verbose 2']
process.caloSimHitStudy.StoreRL = True
process.common_maximum_timex = cms.PSet(
MaxTrackTime = cms.double(1000.0),
MaxTimeNames = cms.vstring(),
Expand Down Expand Up @@ -182,20 +185,20 @@
EkinParticles = cms.vstring(),
Verbosity = cms.untracked.int32(2)
)
process.g4SimHits.Watchers = cms.VPSet(cms.PSet(
CheckForHighEtPhotons = cms.untracked.bool(False),
TrackMin = cms.untracked.int32(0),
TrackMax = cms.untracked.int32(0),
TrackStep = cms.untracked.int32(1),
EventMin = cms.untracked.int32(0),
EventMax = cms.untracked.int32(0),
EventStep = cms.untracked.int32(1),
PDGids = cms.untracked.vint32(),
VerboseLevel = cms.untracked.int32(0),
G4Verbose = cms.untracked.bool(True),
DEBUG = cms.untracked.bool(False),
type = cms.string('TrackingVerboseAction')
))
#process.g4SimHits.Watchers = cms.VPSet(cms.PSet(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hi @bsunanda - do we need the commented out lines?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That watcher is not working right now. I shall work on it later - make it work and then uncomment these lines

# CheckForHighEtPhotons = cms.untracked.bool(False),
# TrackMin = cms.untracked.int32(0),
# TrackMax = cms.untracked.int32(0),
# TrackStep = cms.untracked.int32(1),
# EventMin = cms.untracked.int32(0),
# EventMax = cms.untracked.int32(0),
# EventStep = cms.untracked.int32(1),
# PDGids = cms.untracked.vint32(),
# VerboseLevel = cms.untracked.int32(0),
# G4Verbose = cms.untracked.bool(True),
# DEBUG = cms.untracked.bool(False),
# type = cms.string('TrackingVerboseAction')
#))

# Schedule definition
process.schedule = cms.Schedule(process.generation_step,
Expand Down