Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 23 additions & 3 deletions Analyses/src/StepPointMC1stHitDumper_module.cc
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,11 @@ namespace mu2e {
float x;
float y;
float z;

float InitX;
float InitY;
float InitZ;

float time;

float px;
Expand All @@ -83,11 +88,17 @@ namespace mu2e {
int pdgId;
unsigned particleId;
unsigned volumeCopyNumber;
int eventId;
int subrunId;

VDHit() : x(std::numeric_limits<double>::quiet_NaN())
, y(std::numeric_limits<double>::quiet_NaN())
, z(std::numeric_limits<double>::quiet_NaN())

, InitX(std::numeric_limits<double>::quiet_NaN())
, InitY(std::numeric_limits<double>::quiet_NaN())
, InitZ(std::numeric_limits<double>::quiet_NaN())

, time(std::numeric_limits<double>::quiet_NaN())

, px(std::numeric_limits<double>::quiet_NaN())
Expand All @@ -100,14 +111,20 @@ namespace mu2e {
, pdgId(0)
, particleId(-1U)
, volumeCopyNumber(-1U)
, eventId(0)
, subrunId(0)
{}

//----------------------------------------------------------------
VDHit(const SimParticleTimeOffset& toff, const StepPointMC& hit)
VDHit(const SimParticleTimeOffset& toff, const art::Event& event, const StepPointMC& hit)
: x(hit.position().x())
, y(hit.position().y())
, z(hit.position().z())

, InitX(hit.simParticle()->startPosition().x())
, InitY(hit.simParticle()->startPosition().y())
, InitZ(hit.simParticle()->startPosition().z())

, time(toff.timeWithOffsetsApplied(hit))

, px(hit.momentum().x())
Expand All @@ -122,6 +139,9 @@ namespace mu2e {
, pdgId(hit.simParticle()->pdgId())
, particleId(hit.simParticle()->id().asUint())
, volumeCopyNumber(hit.volumeId())

, eventId(event.id().event())
, subrunId(event.id().subRun())
{}

}; // struct VDHit
Expand Down Expand Up @@ -169,7 +189,7 @@ namespace mu2e {
//================================================================
void StepPointMC1stHitDumper::beginJob() {
art::ServiceHandle<art::TFileService> tfs;
static const char branchDesc[] = "x/F:y/F:z/F:time/F:px/F:py/F:pz/F:pmag/F:ek/F:charge/F:pdgId/I:particleId/i:volumeCopy/i";
static const char branchDesc[] = "x/F:y/F:z/F:InitX/F:InitY/F:InitZ/F:time/F:px/F:py/F:pz/F:pmag/F:ek/F:charge/F:pdgId/I:particleId/i:volumeCopy/i:eventId/I:subrunId/I";
nt_ = tfs->make<TTree>( "nt", "StepPointMC1stHitDumper ntuple");
nt_->Branch("hits", &hit_, branchDesc);
if(writeProperTime_) {
Expand All @@ -195,7 +215,7 @@ namespace mu2e {
// std::cout << event.id() << ": " << stepPoints.at(i) << std::endl;
// }
if (stepPoints.size()) {
hit_ = VDHit(toff_, stepPoints.at(0));
hit_ = VDHit(toff_, event, stepPoints.at(0));

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

The event info seems unrelated to the momentum cutoff: what is it used for? Just diagnostics?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

Hi Dave,

The StepPointMC1stHitDumper module is just an analyzer module that makes a simple ntuple from the first step point MC. I found it handy to make simple plots for diagnostics purposes.
The upper momentum cutoff is done in the filter module: FilterStepPointMomentum.
The upper momentum cut of 0.5 GeV removes all the straight tracks and it reduces the size of the files by about a factor 10. It will be useful to have this functionality for a large scale cosmics MC production.

I've submitted these changes in one commit, and they are completely unrelated.
Sorry, if it adds confusion.
Would it help to separate commits in the future?

// std::cout << event.id() << ", " << hit_.pdgId
// << ", " << hit_.time << ", " << hit_.charge
// << std::endl;
Expand Down
44 changes: 36 additions & 8 deletions Filters/src/FilterStepPointMomentum_module.cc
Original file line number Diff line number Diff line change
Expand Up @@ -23,28 +23,56 @@ namespace mu2e {
typedef std::vector<art::InputTag> InputTags;
InputTags inputs_;
double cutMomentumMin_;
double cutMomentumMax_;

// statistics counters
unsigned numInputEvents_;
unsigned numPassedEvents_;
public:
explicit FilterStepPointMomentum(const fhicl::ParameterSet& pset);

struct Config {
using Name=fhicl::Name;
using Comment=fhicl::Comment;

fhicl::Sequence<art::InputTag> inputs {
Name("inputs"),
Comment("Particles and StepPointMCs mentioned in thise collections will be preserved.")
};

fhicl::Atom<double> cutMomentumMin {
Name("cutMomentumMin"),
Comment("The filter passes events if any of the step points satisties pmag>cutMomentumMin\n"
"By default cutMomentumMin=-inf\n"),
-std::numeric_limits<double>::max()
};

fhicl::Atom<double> cutMomentumMax {
Name("cutMomentumMax"),
Comment("The filter passes events if any of the step points satisties pmag<cutMomentumMax\n"
"By default cutMomentumMax=inf\n"),
std::numeric_limits<double>::max()
};

};

using Parameters = art::EDFilter::Table<Config>;
explicit FilterStepPointMomentum(const Parameters& conf);
virtual bool filter(art::Event& event) override;
virtual void endJob() override;
};

//================================================================
FilterStepPointMomentum::FilterStepPointMomentum(const fhicl::ParameterSet& pset)
: art::EDFilter{pset}
, cutMomentumMin_(pset.get<double>("cutMomentumMin"))
FilterStepPointMomentum::FilterStepPointMomentum(const Parameters& conf)
: art::EDFilter{conf}
, cutMomentumMin_(conf().cutMomentumMin())
, cutMomentumMax_(conf().cutMomentumMax())
, numInputEvents_(0)
, numPassedEvents_(0)
{
typedef std::vector<std::string> VS;
const VS in(pset.get<VS>("inputs"));
for(const auto& i : in) {
for(const auto& i : conf().inputs()) {
inputs_.emplace_back(i);
}

}

//================================================================
Expand All @@ -53,7 +81,7 @@ namespace mu2e {
for(const auto& cn : inputs_) {
auto ih = event.getValidHandle<StepPointMCCollection>(cn);
for(const auto& hit : *ih) {
if(hit.momentum().mag() > cutMomentumMin_) {
if(hit.momentum().mag() > cutMomentumMin_ && hit.momentum().mag() < cutMomentumMax_) {
passed = true;
break;
}
Expand Down
19 changes: 9 additions & 10 deletions JobConfig/cosmic/dsstops_resampler.fcl
Original file line number Diff line number Diff line change
Expand Up @@ -109,13 +109,6 @@ physics.filters.dsFilter: {
vetoDaughters : []
}

physics.filters.stepPointMomentumFilter: {
module_type: FilterStepPointMomentum
inputs : [ "dsResample:crvStage1" ]
cutMomentumMin: 45. // Drop low momentum tracks
cutMomentumMax: 500. // Drop straight tracks
}

physics.filters.trackerStepPointFilter: {
module_type : TrackerStepPointFilter
}
Expand All @@ -128,14 +121,21 @@ physics.filters.detectorFilter: {
vetoDaughters: []
}

physics.filters.stepPointMomentumFilter: {
module_type: FilterStepPointMomentum
inputs : [ "detectorFilter:tracker" ]
cutMomentumMin: 45. // Drop low momentum tracks
cutMomentumMax: 500. // Drop straight tracks
}

physics.filters.g4status: {
module_type: FilterStatusG4
input: "g4run"
maxAcceptedStatus: 1 # status 10 and above means StepPointMCCollection may have non-dereferencable pointers
}


physics.trig: [genCounter, dsResample, g4run, g4consistent, stepPointMomentumFilter, dsFilter, trackerStepPointFilter, detectorFilter]
physics.trig: [genCounter, dsResample, g4run, g4consistent, dsFilter, trackerStepPointFilter, detectorFilter, stepPointMomentumFilter]
physics.g4StatusFilter : [dsResample, g4run, "!g4status"]
physics.outputs: [truncatedEvtsOutput, DSVacuumOut]
physics.an: [genCountLogger]
Expand Down Expand Up @@ -189,5 +189,4 @@ physics.producers.g4run.Mu2eG4CommonCut:{
physics.producers.g4run.ResourceLimits.maxSimParticleCollectionSize: 100000
// Test:
//physics.filters.dsResample.fileNames: ["/pnfs/mu2e/persistent/users/oksuzian/workflow/cry_rs1_1019_g4_10_5/outstage/363672.fcllist_191019181515/00/00000/sim.oksuzian.minedep-filter-cat.cryconcat_rs1_1019.002701_00000000.art"]
//physics.filters.dsResample.fileNames: ["/pnfs/mu2e/scratch/users/oksuzian/workflow/cry_rs1_1019_g4_10_5/good/682690.fcllist_191024104002/00/00047/sim.oksuzian.minedep-filter.cry-minedep_filter-10-5.002701_00020094.art"]
//physics.filters.dsResample.fileNames: ["/pnfs/mu2e/scratch/users/oksuzian/workflow/cry_rs1_1019_g4_10_5/good/683306.fcllist_191024104850/00/00000/sim.oksuzian.minedep-filter.cry-minedep_filter-10-5.002701_00000385.art"]
//physics.filters.dsResample.fileNames: ["/pnfs/mu2e/persistent/users/oksuzian/workflow/cry_rs1_1019_g4_10_5/outstage/683306.fcllist_191024104850/00/00000/sim.oksuzian.minedep-filter.cry-minedep_filter-10-5.002701_00000385.art"]
1 change: 1 addition & 0 deletions JobConfig/mixing/DS-flash-TrkCal-cut.fcl
Original file line number Diff line number Diff line change
Expand Up @@ -51,3 +51,4 @@ services.TFileService.fileName : "nts.owner.DS-flash-TrkCal-cut.version.sequence
services.message.destinations.log.categories.ArtReport.reportEvery : 1
services.message.destinations.log.categories.ArtReport.limit : 1
services.message.destinations.log.categories.ArtReport.timespan : 300
services.GeometryService.simulatedDetector.tool_type : "Mu2e"
163 changes: 4 additions & 159 deletions Mu2eG4/geom/geom_common_phase1.txt
Original file line number Diff line number Diff line change
@@ -1,164 +1,9 @@
//
// Top level geometry file first used in pass 2 of (stage 4 and after) of the CD3c simulation campagin.
//
// Warning: do not write 10000. as 10,000.; it will be read as two numbers (10., 0.).
// A geometry file for using all regular concrete shielding design.
// Most of the geometry is inherited from the existing shielding geometry, except as below.

string detector.name = "g4geom_v00";

bool hasHall = true;
bool hasTarget = true;
bool hasProtonAbsorber = true;
bool hasTSdA = true;
bool hasExternalShielding = true;
bool hasDiskCalorimeter = true;
bool hasBFieldManager = true;
bool hasBeamline = true;
bool hasVirtualDetector = true; // some components, e.g. ProtonAbsorber assume vd presence now;
bool hasCosmicRayShield = true;
bool hasSTM = true;
bool hasMBS = true; // note the two subcomponents, see mbs section below;
// no MBS implies no downstream hole in Cosmic Ray Passive Shield
// Magnetic field may be affected as well

#include "Mu2eG4/geom/g4_visOptions.txt"

//-------------------------------------------
// Mu2e geometry includes
//-------------------------------------------

// X-offset of the PS(+x) and DS(-x) from the Mu2e origin.
// The origin of the detector coordinate system is on the DS axis at the specified z.
double mu2e.solenoidOffset = 3904.; // mm
double mu2e.detectorSystemZ0 = 10171.; // mm G4BL: (17730-7292=9801 mm)

#include "Mu2eG4/geom/mu2eWorld.txt"
// mu2eHall.txt should be used with protonBeamDump_v02.txt, below
//#include "Mu2eG4/geom/mu2eHall.txt"
// whereas mu2eHall_v02.txt should be used with protonBeamDump_v03.txt, below
#include "Mu2eG4/geom/mu2eHall_v02.txt"

// Solenoids
#include "Mu2eG4/geom/DetectorSolenoid_v03.txt"
#include "Mu2eG4/geom/DSShielding_v01.txt"
#include "Mu2eG4/geom/ProductionSolenoid_v02.txt"
#include "Mu2eG4/geom/psEnclosure_v04.txt"
#include "Mu2eG4/geom/PSShield_v06.txt"
#include "Mu2eG4/geom/PSExternalShielding_v01.txt"
#include "Mu2eG4/geom/TransportSolenoid_v03.txt"
// First, the base geometry to start from...
#include "Mu2eG4/geom/geom_common_current.txt"

// External Shielding
#include "Mu2eG4/geom/ExtShieldUpstream_v06.txt"
#include "Mu2eG4/geom/ExtShieldDownstream_v06.txt"
#include "Mu2eG4/geom/Saddle_v02.txt"
#include "Mu2eG4/geom/Pipe_v03.txt"
#include "Mu2eG4/geom/ElectronicRack_v01.txt"

//#include "Mu2eG4/geom/stoppingTarget_TDR.txt" // 17 foil tapered muon stopping target used for the TDR simulations
//#include "Mu2eG4/geom/stoppingTarget_CD3C_34foils.txt" // 34 foil muon stopping target to be used for the CDC3 simulations
#include "Mu2eG4/geom/stoppingTargetHoles_DOE_review_2017.txt" // 37 foil muon stopping target with holes

#include "Mu2eG4/geom/TSdA_v01.txt"
#include "Mu2eG4/geom/muonBeamStop_v05.txt"

//#include "Mu2eG4/geom/MSTM_v01.txt" // muon stopping target monitor (deprecated)
//#include "Mu2eG4/geom/STM_v01.txt" // (muon) stopping target monitor
//#include "Mu2eG4/geom/STM_v02.txt" // (muon) stopping target monitor
#include "Mu2eG4/geom/STM_v03.txt" // (muon) stopping target monitor

// Proton Absorber
#include "Mu2eG4/geom/protonAbsorber_cylindrical_v02.txt"
#include "Mu2eG4/geom/degrader_v01.txt"

#include "Mu2eG4/geom/ProductionTargetInPS.txt"
//#include "Mu2eG4/geom/protonBeamDump_v02.txt"
#include "Mu2eG4/geom/protonBeamDump_v03.txt"
#include "Mu2eG4/geom/extmon_fnal_v02.txt"

#include "Mu2eG4/geom/tracker_v5.txt"

// Crystal calorimeter
#include "Mu2eG4/geom/calorimeter_CsI.txt"

//CRV counters
#include "Mu2eG4/geom/crv_counters_v07.txt"


//-------------------------------------------
// Magnetic field
//-------------------------------------------

// Form of DS field: 0 is full field;
// 1 is full upstream, const downstream;
// 2 is const throughout
int detSolFieldForm = 0;
vector<string> bfield.innerMaps = {
"BFieldMaps/Mau13/DSMap.header",
"BFieldMaps/Mau13/PSMap.header",
"BFieldMaps/Mau13/TSuMap.header",
"BFieldMaps/Mau13/TSdMap.header",
"BFieldMaps/Mau13/PStoDumpAreaMap.header",
"BFieldMaps/Mau13/ProtonDumpAreaMap.header",
"BFieldMaps/Mau13/DSExtension.header"
};

// Value of the uniform magnetic field with the DS volume (only for
// detSolFieldForm>0)
double toyDS.bz = 1.0;

// Gradient of field in DS2 volume. Applied only in the case
// of detSolFieldForm=1 or detSolFieldForm=2.
double toyDS.gradient = 0.0; // Tesla/m

// This is recommended field map. See geom_mecofield.txt to use the meco field.
string bfield.format = "G4BL";

// The other option is "meco"
string bfield.interpolationStyle = trilinear;

int bfield.verbosityLevel = 0;
bool bfield.writeG4BLBinaries = false;

vector<string> bfield.outerMaps = {
"BFieldMaps/Mau9/ExtMonUCIInternal1AreaMap.header",
"BFieldMaps/Mau9/ExtMonUCIInternal2AreaMap.header",
"BFieldMaps/Mau9/ExtMonUCIAreaMap.header",
"BFieldMaps/Mau9/PSAreaMap.header"
};

// This scale factor is of limited use.
// It can make approximate sense to scale the PS field to get a rough
// answer; the answer will be wrong in detail.
// It never makes sense to scale the TS field.
// Not sure if it ever makes sense to scale the PS field.
double bfield.scaleFactor = 1.0;

//---------------------------------------
// Virtual detectors
//---------------------------------------
double vd.halfLength = 0.01; //mm
int vd.verbosityLevel = 0;
bool vd.visible = true;
bool vd.solid = false;

// // VD right in front of a hall wall
// double vd.ExtMonCommonPlane.z = -11999.99;


//---------------------------------------
// Region visualization
//---------------------------------------
#include "Mu2eG4/geom/visualization_regions.txt"


//
//
// End notes:
//
// 1) Sources of information:
//
//
//
// This tells emacs to view this file in c++ mode.
// Local Variables:
// mode:c++
// End: