Skip to content

Commit

Permalink
Merge pull request #616 from michaelmackenzie/Mu2eII_MergeMain
Browse files Browse the repository at this point in the history 
Mu2eII_SM21: Merge in main at tag v10_07_00
  • Loading branch information
@brownd1978
brownd1978 authored Oct 12, 2021
2 parents 5bc9d27 + 398fcb3 commit 38b8236
Show file tree
Hide file tree
Showing 129 changed files with 8,394 additions and 7,116 deletions.
2 changes: 1 addition & 1 deletion .muse
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# prefer to build with this environment
ENVSET p012
ENVSET p013
# add these to python path
PYTHONPATH TrackerAlignment/scripts
PYTHONPATH Trigger/python
Expand Down
43 changes: 43 additions & 0 deletions Analyses/fcl/DSBFieldDump.fcl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
#
# Graph the DS field on its axis
#
#include "Offline/fcl/standardServices.fcl"

process_name : DSField

#services : @local::Services.Reco
services : {
message : @local::default_message
GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" }
GeometryService : {
inputFile : "Offline/Mu2eG4/geom/geom_common.txt"
simulatedDetector : { tool_type: "Mu2e" }
}
TFileService: {
fileName : "DSField.root"
}
}
source : {
module_type : EmptyEvent
}

physics : {
analyzers: {
BFieldPlotter : {
module_type : BFieldPlotter
plane : "x"
planeValue : -3904.0
axisOneMin : 0.0
axisOneMax : 0.0
axisTwoMin : 4100
axisTwoMax : 13500
mapBinSize : 10.0
dump : true
dumpName : "Test"
detector : true
}
}
e1 : [BFieldPlotter]
end_paths : [e1]
}
services.GeometryService.inputFile : "Offline/Mu2eG4/geom/geom_common_reco.txt"
31 changes: 31 additions & 0 deletions Analyses/fcl/DSField.fcl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
#
# Graph the DS field on its axis
#
#include "Offline/fcl/standardServices.fcl"

process_name : DSField

services : @local::Services.Reco
source : {
module_type : EmptyEvent
}

physics : {
analyzers: {
BFieldPlotter : {
module_type : BFieldPlotter
plane : "x"
planeValue : -3904.0
axisOneMin : 0.0
axisOneMax : 0.0
axisTwoMin : 4100
axisTwoMax : 13500
mapBinSize : 10.0
dumpFile : true
detector : true
}
}
e1 : [BFieldPlotter]
end_paths : [e1]
}
services.GeometryService.inputFile : "Offline/Mu2eG4/geom/geom_common_no_tsu_ps.txt"
192 changes: 114 additions & 78 deletions Analyses/src/BFieldPlotter_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,17 @@

// Mu2e includes
#include "Offline/BFieldGeom/inc/BFieldManager.hh"
#include "Offline/GeometryService/inc/GeomHandle.hh"

// art includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "fhiclcpp/types/Atom.h"
#include "art_root_io/TFileService.h"
#include "art/Framework/Services/Registry/ServiceDefinitionMacros.h"
#include "Offline/GeometryService/inc/GeomHandle.hh"
#include "Offline/GeometryService/inc/DetectorSystem.hh"

// ROOT includes
#include "TFile.h"
Expand All @@ -26,46 +28,53 @@
#include <iomanip>
#include <vector>
#include <string>
#include <fstream>
#include <cmath>

using namespace std;

namespace mu2e {

class BFieldPlotter : public art::EDAnalyzer {
public:
struct Config {
using Name=fhicl::Name;
using Comment=fhicl::Comment;
fhicl::Atom<std::string> plane {Name("plane" ), Comment("Axis plane is defined by (x, y, or z)")};
fhicl::Atom<double> planeValue{Name("planeValue"), Comment("Value on the axis the plane intersects (mm)")};
fhicl::Atom<double> axisOneMin{Name("axisOneMin"), Comment("Axis one lower edge (bin centered) for plotting (mm) (plane = x/y/z --> axis one = y/x/x)")};
fhicl::Atom<double> axisOneMax{Name("axisOneMax"), Comment("Axis one upper edge (bin centered) for plotting (mm) (plane = x/y/z --> axis one = y/x/x)")};
fhicl::Atom<double> axisTwoMin{Name("axisTwoMin"), Comment("Axis two lower edge (bin centered) for plotting (mm) (plane = x/y/z --> axis two = z/z/y)")};
fhicl::Atom<double> axisTwoMax{Name("axisTwoMax"), Comment("Axis two upper edge (bin centered) for plotting (mm) (plane = x/y/z --> axis two = z/z/y)")};
fhicl::Atom<double> mapBinSize{Name("mapBinSize"), Comment("Map bin size (mm) (must be a divisor of both axis lengths)"), 10.};
};
typedef art::EDAnalyzer::Table<Config> Parameters;

explicit BFieldPlotter(const Parameters& pset);
virtual ~BFieldPlotter() { }

void analyze(const art::Event& e);

void fillHistogram(BFMap const *map, art::ServiceHandle<art::TFileService> tfs);

private:

std::string _plane ; // x, y, or z map plane
double _planeValue; // value to set the plane at (mm)
double _axisOneMin; // sampling axis values (mm)
double _axisOneMax; // sampling axis values (mm)
double _axisTwoMin; // sampling axis values (mm)
double _axisTwoMax; // sampling axis values (mm)
double _mapBinSize; // histogram bin size (mm)

std::map<std::string,TH2F*> _hMap; //histogram of the map

public:
struct Config {
using Name=fhicl::Name;
using Comment=fhicl::Comment;
fhicl::Atom<std::string> plane {Name("plane" ), Comment("Axis plane is defined by (x, y, or z)")};
fhicl::Atom<double> planeValue{Name("planeValue"), Comment("Value on the axis the plane intersects (mm)")};
fhicl::Atom<double> axisOneMin{Name("axisOneMin"), Comment("Axis one lower edge (bin centered) for plotting (mm) (plane = x/y/z --> axis one = y/x/x)")};
fhicl::Atom<double> axisOneMax{Name("axisOneMax"), Comment("Axis one upper edge (bin centered) for plotting (mm) (plane = x/y/z --> axis one = y/x/x)")};
fhicl::Atom<double> axisTwoMin{Name("axisTwoMin"), Comment("Axis two lower edge (bin centered) for plotting (mm) (plane = x/y/z --> axis two = z/z/y)")};
fhicl::Atom<double> axisTwoMax{Name("axisTwoMax"), Comment("Axis two upper edge (bin centered) for plotting (mm) (plane = x/y/z --> axis two = z/z/y)")};
fhicl::Atom<double> mapBinSize{Name("mapBinSize"), Comment("Map bin size (mm) (must be a divisor of both axis lengths)"), 10.};
fhicl::Atom<bool> dump {Name("dump"), Comment("Dump map samples to a CSV file"), false};
fhicl::Atom<std::string> dumpName {Name("dumpName"), Comment("Suffix to the dump file name"),"Dump"};
fhicl::Atom<bool> detector {Name("detector"), Comment("Use Detector coordinate system"), false};
};
typedef art::EDAnalyzer::Table<Config> Parameters;

explicit BFieldPlotter(const Parameters& pset);
virtual ~BFieldPlotter() { }

void analyze(const art::Event& e);

void fillHistogram(BFMap const *map, art::ServiceHandle<art::TFileService> tfs);

private:

std::string _plane ; // x, y, or z map plane
double _planeValue; // value to set the plane at (mm)
double _axisOneMin; // sampling axis values (mm)
double _axisOneMax; // sampling axis values (mm)
double _axisTwoMin; // sampling axis values (mm)
double _axisTwoMax; // sampling axis values (mm)
double _mapBinSize; // histogram bin size (mm)
bool _dump; // dump to CSV?
std::string _dumpname;
bool _detector; // use detector system?

std::map<std::string,TH2F*> _hMap; //histogram of the map

};

BFieldPlotter::BFieldPlotter(const Parameters& pset) :
Expand All @@ -77,19 +86,22 @@ namespace mu2e {
, _axisTwoMin(pset().axisTwoMin())
, _axisTwoMax(pset().axisTwoMax())
, _mapBinSize(pset().mapBinSize())
, _dump(pset().dump())
, _dumpname(pset().dumpName())
, _detector(pset().detector())
{
if(_axisOneMin >= _axisOneMax || _axisTwoMin >= _axisTwoMax) {
if(_axisOneMin > _axisOneMax || _axisTwoMin > _axisTwoMax) {
throw cet::exception("BADCONFIG") << "BField mapping plane ill defined: "
<< _axisOneMin << " < axisOneValues < " << _axisOneMax << ", "
<< _axisTwoMin << " < axisTwoValues < " << _axisTwoMax;
<< _axisOneMin << " < axisOneValues < " << _axisOneMax << ", "
<< _axisTwoMin << " < axisTwoValues < " << _axisTwoMax;
}
if(_mapBinSize <= 0.) {
throw cet::exception("BADCONFIG") << "BField map binning should be >= 0 but given as " << _mapBinSize;
}

if(!(_plane == "x" || _plane == "y" || _plane == "z")) {
throw cet::exception("BADCONFIG") << "BField map plane not recognized! Options are x, y, or z but given "
<< _plane.c_str();
throw cet::exception("BADCONFIG") << "BField map plane not recognized! Options are x, y, or z but given "
<< _plane.c_str();
}

} //end constructor
Expand All @@ -100,16 +112,16 @@ namespace mu2e {
GeomHandle<BFieldManager> bf;
const BFieldManager::MapContainerType& innerMaps = bf->getInnerMaps();
const BFieldManager::MapContainerType& outerMaps = bf->getOuterMaps();

//plot inner maps
for(std::shared_ptr<BFMap> const & map : innerMaps) {
for(std::shared_ptr<BFMap> const & map : innerMaps) {
fillHistogram(map.get(), tfs);
}
}

//plot outer maps
for(std::shared_ptr<BFMap> const & map : outerMaps) {
for(std::shared_ptr<BFMap> const & map : outerMaps) {
fillHistogram(map.get(), tfs);
}
}

//plot the default field returned from the manager
fillHistogram(NULL, tfs);
Expand All @@ -118,64 +130,88 @@ namespace mu2e {

//fill a histogram with the magnetic field
void BFieldPlotter::fillHistogram(BFMap const *map, art::ServiceHandle<art::TFileService> tfs) {
GeomHandle<DetectorSystem> det;
//define histogram binning such that map edges are bin centers and inside histogram
long nbinsOne = std::lround((_axisOneMax - _axisOneMin)/_mapBinSize) + 1;
long nbinsTwo = std::lround((_axisTwoMax - _axisTwoMin)/_mapBinSize) + 1;
//check that the map bin step size works with the edges given
if(std::abs(_axisOneMin + (nbinsOne-1)*_mapBinSize - _axisOneMax > (_axisOneMax-_axisOneMin)/(100.*(nbinsOne-1))))
if(nbinsOne > 1 && std::abs(_axisOneMin + (nbinsOne-1)*_mapBinSize - _axisOneMax > (_axisOneMax-_axisOneMin)/(100.*(nbinsOne-1))))
throw cet::exception("BADCONFIG") << "BField mapping axis values not steppable with step size given: "
<< _axisOneMin << " < axisOneValues < " << _axisOneMax << ", "
<< "step size = " << _mapBinSize;
if(std::abs(_axisTwoMin + (nbinsTwo-1)*_mapBinSize - _axisTwoMax > (_axisTwoMax-_axisTwoMin)/(100.*(nbinsTwo-1))))
<< _axisOneMin << " < axisOneValues < " << _axisOneMax << ", "
<< "step size = " << _mapBinSize;

if(nbinsTwo > 1 && std::abs(_axisTwoMin + (nbinsTwo-1)*_mapBinSize - _axisTwoMax > (_axisTwoMax-_axisTwoMin)/(100.*(nbinsTwo-1))))
throw cet::exception("BADCONFIG") << "BField mapping axis values not steppable with step size given: "
<< _axisTwoMin << " < axisTwoValues < " << _axisTwoMax << ", "
<< "step size = " << _mapBinSize;
<< _axisTwoMin << " < axisTwoValues < " << _axisTwoMax << ", "
<< "step size = " << _mapBinSize;

//if a null map, plot the default field from the manager
const std::string name = (map) ? map->getKey() : "default";

if(_hMap[name]) return; //already filled the maps in a previous event

// if dumping, create the file
std::ofstream fs;
if(_dump){
string dumpfile= name+_dumpname+".csv";
fs.open (dumpfile.c_str(), fstream::out);
if(fs.is_open()){
fs << "# Dump of " << name << endl;
fs << "# x , y , z , Bx , By , Bz" << endl;
} else
throw cet::exception("BADCONFIG") << "Cannot open file " << dumpfile << endl;
}

//make a directory for the map
art::TFileDirectory tfdir = tfs->mkdir( ("BFieldMapper_"+name).c_str() );

//define a histogram
_hMap[name] = tfdir.make<TH2F>(("hMap"+name).c_str() , (name + " Magnetic field map").c_str(),
// add offsets so all values are bin centers and fit edge into the map
nbinsOne, _axisOneMin - _mapBinSize/2.,_axisOneMax + _mapBinSize/2.,
nbinsTwo, _axisTwoMin - _mapBinSize/2.,_axisTwoMax + _mapBinSize/2.);
_hMap[name] = tfdir.make<TH2F>(("hMap"+name).c_str() , (name + " Magnetic field map").c_str(),
// add offsets so all values are bin centers and fit edge into the map
nbinsOne, _axisOneMin - _mapBinSize/2.,_axisOneMax + _mapBinSize/2.,
nbinsTwo, _axisTwoMin - _mapBinSize/2.,_axisTwoMax + _mapBinSize/2.);

GeomHandle<BFieldManager> bf; //only needed in null map case
//Loop through the points in the magnetic field
double axisOne, axisTwo;
for(int binOne = 0; binOne < nbinsOne; ++binOne) {
axisOne = _axisOneMin + binOne*_mapBinSize;
if(binOne == 0 || binOne == nbinsOne-1) //if first or last bin, at slight step into map to avoid edge issues
axisOne += (binOne) ? -_mapBinSize/100. : _mapBinSize/100.;
axisOne += (binOne) ? -_mapBinSize/100. : _mapBinSize/100.;
for(int binTwo = 0; binTwo < nbinsTwo; ++binTwo) {
axisTwo = _axisTwoMin + binTwo*_mapBinSize;
if(binTwo == 0 || binTwo == nbinsTwo-1) //if first or last bin, at slight step into map to avoid edge issues
axisTwo += (binTwo) ? -_mapBinSize/100. : _mapBinSize/100.;
CLHEP::Hep3Vector point;
CLHEP::Hep3Vector field;
if(_plane == "x")
point = CLHEP::Hep3Vector(_planeValue, axisOne, axisTwo);
else if(_plane == "y")
point = CLHEP::Hep3Vector(axisOne, _planeValue, axisTwo);
else if(_plane == "z")
point = CLHEP::Hep3Vector(axisOne, axisTwo, _planeValue);
if(map)
map->getBFieldWithStatus(point,field);
else //if null, get the manager and plot the field returned by it
bf->getBFieldWithStatus(point,field);

_hMap[name]->Fill(axisOne, axisTwo, field.mag()); //fill with weight of the field magnitude
axisTwo = _axisTwoMin + binTwo*_mapBinSize;
if(binTwo == 0 || binTwo == nbinsTwo-1) //if first or last bin, at slight step into map to avoid edge issues
axisTwo += (binTwo) ? -_mapBinSize/100. : _mapBinSize/100.;
CLHEP::Hep3Vector point;
CLHEP::Hep3Vector field;
if(_plane == "x")
point = CLHEP::Hep3Vector(_planeValue, axisOne, axisTwo);
else if(_plane == "y")
point = CLHEP::Hep3Vector(axisOne, _planeValue, axisTwo);
else if(_plane == "z")
point = CLHEP::Hep3Vector(axisOne, axisTwo, _planeValue);
if(map)
map->getBFieldWithStatus(point,field);
else //if null, get the manager and plot the field returned by it
bf->getBFieldWithStatus(point,field);

_hMap[name]->Fill(axisOne, axisTwo, field.mag()); //fill with weight of the field magnitude

// if there's a dump file, fill it.
if(_dump){
CLHEP::Hep3Vector dpoint = _detector ? det->toDetector(point) : point;
fs << std::setw(10) << std::setprecision(5)
<< dpoint.x() << ", " << std::setprecision(5)
<< dpoint.y() << ", " << std::setprecision(5)
<< dpoint.z() << ", " << std::setprecision(5)
<< field.x() << ", " << std::setprecision(5)
<< field.y() << ", " << std::setprecision(5)
<< field.z() << endl;
}

} //end axis two loop
} //end axis one loop
if(_dump && fs.is_open())fs.close();
} //end fillHistorgram


} // end namespace mu2e

DEFINE_ART_MODULE(mu2e::BFieldPlotter);
1 change: 1 addition & 0 deletions Analyses/src/SConscript
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,7 @@ helper.make_plugins( [ mainlib,
'tbb',
'cetlib',
'cetlib_except',
'General',
'CLHEP',
'HepPDT',
'HepPID',
Expand Down
Loading

0 comments on commit 38b8236

Please sign in to comment.