Skip to content
GitLab
Commit bb59cd38 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

new file: examples/scripts/mag_field.cxx 

	new file:   src/ConceptDetectors/topside/compact/elements.xml
parent 142b2485
Expand all Hide whitespace changes
Inline Side-by-side
examples/scripts/mag_field.cxx 0 → 100644
View file @ bb59cd38
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/VectorUtil.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH1D.h"
#include "TTree.h"
#include "TF1.h"
#include "ROOT/TDataFrame.hxx"
#include <vector>
#include <tuple>
#include <algorithm>
#include <iterator>
// DD4hep
// -----
// In .rootlogon.C
// gSystem->Load("libDDDetectors");
// gSystem->Load("libDDG4IO");
// gInterpreter->AddIncludePath("/opt/software/local/include");
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "DD4hep/DD4hepUnits.h"
/** Cell size example.
*
*/
void mag_field(double z0 = 0.0 ){
using namespace ROOT::Math;
using namespace dd4hep;
// -------------------------
// Get the DD4hep instance
// Load the compact XML file
dd4hep::Detector& detector = dd4hep::Detector::getInstance();
detector.fromCompact("beamline_test.xml");
double x_min = -100.0;
double x_max = 100.0;
double y_min = -100.0;
double y_max = 100.0;
double z_min = -100.0;
double z_max = 100.0;
int nx_bins = 200;
int ny_bins = 200;
int nz_bins = 200;
double dx = (x_max -x_min)/nx_bins;
double dy = (y_max -y_min)/ny_bins;
double dz = (z_max -z_min)/nz_bins;
double y0 = 20.0*cm;
TH2F* hBx = new TH2F("hBx","; z [cm]; x [cm]; B_{x} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
TH2F* hBz = new TH2F("hBz","; z [cm]; x [cm]; B_{z} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
TH2F* hBr = new TH2F("hBr","; z [cm]; x [cm]; B_{r} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
TH2F* hBphi = new TH2F("hBphi","; z [cm]; x [cm]; B_{#phi} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
for(int ix = 0; ix< nx_bins; ix++) {
for(int iz = 0; iz< nz_bins; iz++) {
double x = (x_min + ix*dx)*cm;
double z = (z_min + iz*dz)*cm;
XYZVector pos{ x, y0, z };
XYZVector field = detector.field().magneticField(pos);
hBz->Fill(z,x, field.z()/tesla);
hBx->Fill(z,x, field.x()/tesla);
hBr->Fill(z,x, field.rho()/tesla);
hBphi->Fill(z,x, field.phi()/tesla);
}
//for(int iy = 0; iy< ny_bins; iy++) {
// double x = (x_min + ix*dx)*cm;
// double y = (y_min + iy*dy)*cm;
// XYyVector pos{ x, y, z0 };
// XYyVector field = detector.field().magneticField(pos);
// hBz->Fill(z,x, field.z()/tesla);
// hBx->Fill(z,x, field.x()/tesla);
// hBr->Fill(z,x, field.rho()/tesla);
//}
}
auto c = new TCanvas();
c->Divide(2,2);
c->cd(1);
hBz->Draw("colz");
c->cd(2);
hBx->Draw("colz");
c->cd(3);
hBr->Draw("colz");
c->cd(4);
hBphi->Draw("colz");
}
src/ConceptDetectors/topside/compact/elements.xml 0 → 100644
View file @ bb59cd38
This diff is collapsed.
src/ConceptDetectors/topside/compact/materials.xml 0 → 100644
View file @ bb59cd38
<?xml version="1.0" encoding="UTF-8"?>
<materials>
<!--
Air by weight from
http://www.engineeringtoolbox.com/air-composition-24_212.html
-->
<material name="Air">
<D type="density" unit="g/cm3" value="0.0012"/>
<fraction n="0.754" ref="N"/>
<fraction n="0.234" ref="O"/>
<fraction n="0.012" ref="Ar"/>
</material>
<!-- We model vakuum just as very thin air -->
<material name="Vacuum">
<D type="density" unit="g/cm3" value="0.0000000001"/>
<fraction n="0.754" ref="N"/>
<fraction n="0.234" ref="O"/>
<fraction n="0.012" ref="Ar"/>
</material>
<material name="Epoxy">
<D type="density" value="1.3" unit="g/cm3"/>
<composite n="44" ref="H"/>
<composite n="15" ref="C"/>
<composite n="7" ref="O"/>
</material>
<material name="Quartz">
<D type="density" value="2.2" unit="g/cm3"/>
<composite n="1" ref="Si"/>
<composite n="2" ref="O"/>
</material>
<material name="G10">
<D type="density" value="1.7" unit="g/cm3"/>
<fraction n="0.08" ref="Cl"/>
<fraction n="0.773" ref="Quartz"/>
<fraction n="0.147" ref="Epoxy"/>
</material>
<material name="Polystyrene">
<D value="1.032" unit="g/cm3"/>
<composite n="19" ref="C"/>
<composite n="21" ref="H"/>
</material>
<material name="Steel235">
<D value="7.85" unit="g/cm3"/>
<fraction n="0.998" ref="Fe"/>
<fraction n=".002" ref="C"/>
</material>
<material name="SiliconOxide">
<D type="density" value="2.65" unit="g/cm3"/>
<composite n="1" ref="Si"/>
<composite n="2" ref="O"/>
</material>
<material name="BoronOxide">
<D type="density" value="2.46" unit="g/cm3"/>
<composite n="2" ref="B"/>
<composite n="3" ref="O"/>
</material>
<material name="SodiumOxide">
<D type="density" value="2.65" unit="g/cm3"/>
<composite n="2" ref="Na"/>
<composite n="1" ref="O"/>
</material>
<material name="AluminumOxide">
<D type="density" value="3.89" unit="g/cm3"/>
<composite n="2" ref="Al"/>
<composite n="3" ref="O"/>
</material>
<material name="PyrexGlass">
<D type="density" value="2.23" unit="g/cm3"/>
<fraction n="0.806" ref="SiliconOxide"/>
<fraction n="0.130" ref="BoronOxide"/>
<fraction n="0.040" ref="SodiumOxide"/>
<fraction n="0.023" ref="AluminumOxide"/>
</material>
<material name="CarbonFiber">
<D type="density" value="1.5" unit="g/cm3"/>
<fraction n="0.65" ref="C"/>
<fraction n="0.35" ref="Epoxy"/>
</material>
<material name="CarbonFiber_50D">
<D type="density" value="0.75" unit="g/cm3"/>
<fraction n="0.65" ref="C"/>
<fraction n="0.35" ref="Epoxy"/>
</material>
<material name="Rohacell31">
<D type="density" value="0.032" unit="g/cm3"/>
<composite n="9" ref="C"/>
<composite n="13" ref="H"/>
<composite n="2" ref="O"/>
<composite n="1" ref="N"/>
</material>
<material name="Rohacell31_50D">
<D type="density" value="0.016" unit="g/cm3"/>
<composite n="9" ref="C"/>
<composite n="13" ref="H"/>
<composite n="2" ref="O"/>
<composite n="1" ref="N"/>
</material>
<material name="RPCGasDefault" state="gas">
<D type="density" value="0.0037" unit="g/cm3"/>
<composite n="209" ref="C"/>
<composite n="239" ref="H"/>
<composite n="381" ref="F"/>
</material>
<material name="PolystyreneFoam">
<D type="density" value="0.0056" unit="g/cm3"/>
<fraction n="1.0" ref="Polystyrene"/>
</material>
<material name="Kapton">
<D value="1.43" unit="g/cm3"/>
<composite n="22" ref="C"/>
<composite n="10" ref="H"/>
<composite n="2" ref="N"/>
<composite n="5" ref="O"/>
</material>
<material name="PEEK">
<D value="1.37" unit="g/cm3"/>
<composite n="19" ref="C"/>
<composite n="12" ref="H"/>
<composite n="3" ref="O"/>
</material>
<material name="TungstenDens23">
<D value="17.7" unit="g / cm3"/>
<fraction n="0.925" ref="W"/>
<fraction n="0.066" ref="Ni"/>
<fraction n="0.009" ref="Fe"/>
</material>
<material name="TungstenDens24">
<D value="17.8" unit="g / cm3"/>
<fraction n="0.93" ref="W"/>
<fraction n="0.061" ref="Ni"/>
<fraction n="0.009" ref="Fe"/>
</material>
<material name="TungstenDens25">
<D value="18.2" unit="g / cm3"/>
<fraction n="0.950" ref="W"/>
<fraction n="0.044" ref="Ni"/>
<fraction n="0.006" ref="Fe"/>
</material>
<material name="CarbonFiber_25percent">
<D type="density" value="0.375" unit="g / cm3"/>
<fraction n="1.0" ref="CarbonFiber"/>
</material>
<material name="CarbonFiber_15percent">
<D type="density" value="0.225" unit="g / cm3"/>
<fraction n="1.0" ref="CarbonFiber"/>
</material>
<material name="Rohacell31_50percent">
<D type="density" value="0.016" unit="g / cm3"/>
<fraction n="1.0" ref="Rohacell31"/>
</material>
<material name="Rohacell31_15percent">
<D type="density" value="0.0048" unit="g / cm3"/>
<fraction n="1.0" ref="Rohacell31"/>
</material>
<material name="BoratedPolyethylene5">
<D value="0.93" unit="g / cm3"/>
<fraction n="0.612" ref="C"/>
<fraction n="0.222" ref="O"/>
<fraction n="0.116" ref="H"/>
<fraction n="0.050" ref="B"/>
</material>
<material name="SiliconCarbide">
<D value="3.1" unit="g / cm3"/>
<composite n="1" ref="Si"/>
<composite n="1" ref="C"/>
</material>
<material name="SiliconCarbide_6percent">
<D value="0.186" unit="g / cm3"/>
<fraction n="1.0" ref="SiliconCarbide"/>
</material>
</materials>
src/ConceptDetectors/topside/compact/topside.xml 0 → 100644
View file @ bb59cd38
This diff is collapsed.
src/ConceptDetectors/topside/digi/CellInfo.cpp 0 → 100644
View file @ bb59cd38
#include <DD4hep/MultiSegmentation.h>
#include "CellInfo.h"
using namespace dd4hep;
namespace topside {
namespace digi {
static const float sqrt12 = sqrt(12.);
static const float sqrt3 = sqrt(3.);
void cartesianCellInfo(CellInfo &info, uint64_t cellID, VolumeManagerContext *ctxt) {
DetElement element = ctxt->element;
Alignment align = element.survey();
SensitiveDetector sd = ctxt->volumePlacement().volume().sensitiveDetector();
Segmentation seg = sd.readout().segmentation();
Position cellPosLocal = seg.position(cellID);
info.meanPos = align.localToWorld(cellPosLocal);
std::vector<double> cellDims = seg.cellDimensions(cellID);
double origin[3] = {0};
std::string type = seg.type();
while (true) {
if (type.compare("CartesianGridXY") == 0) {
info.noiseContrib.push_back(
align.localToWorld(Position(cellDims[0] / sqrt12, 0, 0) + cellPosLocal) - info.meanPos);
info.noiseContrib.push_back(
align.localToWorld(Position(0, cellDims[1] / sqrt12, 0) + cellPosLocal) - info.meanPos);
double zdir[3] = {0, 0, 1};
info.noiseContrib.push_back(
align.localToWorld(
Position(
0, 0,
ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, zdir) / sqrt3) +
cellPosLocal) -
info.meanPos);
break;
} else if (type.compare("CartesianGridXZ") == 0) {
info.noiseContrib.push_back(
align.localToWorld(Position(cellDims[0] / sqrt12, 0, 0) + cellPosLocal) - info.meanPos);
double ydir[3] = {0, 1, 0};
info.noiseContrib.push_back(
align.localToWorld(
Position(
0, ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, ydir) / sqrt3,
0) +
cellPosLocal) -
info.meanPos);
info.noiseContrib.push_back(
align.localToWorld(Position(0, 0, cellDims[1] / sqrt12) + cellPosLocal) - info.meanPos);
break;
} else if (type.compare("CartesianGridYZ") == 0) {
double xdir[3] = {1, 0, 0};
info.noiseContrib.push_back(
align.localToWorld(
Position(
ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, xdir) / sqrt3, 0,
0) +
cellPosLocal) -
info.meanPos);
info.noiseContrib.push_back(
align.localToWorld(Position(0, cellDims[0] / sqrt12, 0) + cellPosLocal) - info.meanPos);
info.noiseContrib.push_back(
align.localToWorld(Position(0, 0, cellDims[1] / sqrt12) + cellPosLocal) - info.meanPos);
break;
} else if (type.compare("CartesianStripX") == 0) {
info.noiseContrib.push_back(
align.localToWorld(Position(cellDims[0] / sqrt12, 0, 0) + cellPosLocal) - info.meanPos);
double ydir[3] = {0, 1, 0};
info.noiseContrib.push_back(
align.localToWorld(
Position(
0, ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, ydir) / sqrt3,
0) +
cellPosLocal) -
info.meanPos);
double zdir[3] = {0, 0, 1};
info.noiseContrib.push_back(
align.localToWorld(
Position(
0, 0,
ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, zdir) / sqrt3) +
cellPosLocal) -
info.meanPos);
break;
} else if (type.compare("CartesianStripY") == 0) {
double xdir[3] = {1, 0, 0};
info.noiseContrib.push_back(
align.localToWorld(
Position(
ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, xdir) / sqrt3, 0,
0) +
cellPosLocal) -
info.meanPos);
info.noiseContrib.push_back(
align.localToWorld(Position(0, cellDims[0] / sqrt12, 0) + cellPosLocal) - info.meanPos);
double zdir[3] = {0, 0, 1};
info.noiseContrib.push_back(
align.localToWorld(
Position(
0, 0,
ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, zdir) / sqrt3) +
cellPosLocal) -
info.meanPos);
break;
} else if (type.compare("CartesianStripZ") == 0) {
double xdir[3] = {1, 0, 0};
info.noiseContrib.push_back(
align.localToWorld(
Position(
ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, xdir) / sqrt3, 0,
0) +
cellPosLocal) -
info.meanPos);
double ydir[3] = {0, 1, 0};
info.noiseContrib.push_back(
align.localToWorld(
Position(
0, ctxt->volumePlacement().volume()->GetShape()->DistFromInside(origin, ydir) / sqrt3,
0) +
cellPosLocal) -
info.meanPos);
info.noiseContrib.push_back(
align.localToWorld(Position(0, 0, cellDims[0] / sqrt12) + cellPosLocal) - info.meanPos);
break;
} else if (type.compare("MultiSegmentation") == 0) {
type = ((DDSegmentation::MultiSegmentation *)seg.segmentation())->subsegmentation(cellID).type();
} else {
std::cout << "Unknown seg type: " << type << std::endl;
break;
}
}
}
} // namespace digi
} // namespace topside
src/ConceptDetectors/topside/digi/CellInfo.h 0 → 100644
View file @ bb59cd38
#ifndef TOPSIDE_DIGI_CELL_INFO_H
#define TOPSIDE_DIGI_CELL_INFO_H
#include <vector>
#include <DD4hep/VolumeManager.h>
namespace topside {
namespace digi {
struct CellInfo {
dd4hep::Position meanPos;
std::vector<dd4hep::Position> noiseContrib;
};
void cartesianCellInfo(CellInfo &info, uint64_t cellID, dd4hep::VolumeManagerContext *ctxt);
} // namespace digi
} // namespace topside
#endif // TOPSIDE_DIGI_CELL_INFO_H
src/ConceptDetectors/topside/digi/Context.h 0 → 100644
View file @ bb59cd38
#ifndef TOPSIDE_DIGI_CONTEXT_H
#define TOPSIDE_DIGI_CONTEXT_H
#include <TRandom.h>
#include <proio/reader.h>
#include <proio/writer.h>
namespace topside {
namespace digi {
struct Context {
proio::Reader *reader;
proio::Writer *writer;
TRandom *RNG;
};
} // namespace digi
} // namespace topside
#endif // TOPSIDE_DIGI_CONTEXT_H
src/ConceptDetectors/topside/digi/README.md 0 → 100644
View file @ bb59cd38
# digi
source code for the `topside_digi` executable which digitizes TOPSiDE simulation output
## Contents
* ReadoutTypes: digitization algorithms that implement `ReadoutFactory` by calling the `READOUT_TYPE` macro
* topside_digi.cpp: main code that calls necessary readout types
src/ConceptDetectors/topside/digi/Readout.h 0 → 100644
View file @ bb59cd38
#ifndef TOPSIDE_DIGI_READOUT_TYPE_H
#define TOPSIDE_DIGI_READOUT_TYPE_H
namespace topside {
namespace digi {
class Readout {
public:
Readout() {}
virtual ~Readout() {}
virtual std::vector<std::string> Types() { return {}; }
// perform set up at the beginning of each event
virtual void BeginEvent(Context *ctxt, proio::Event *event) {}
// process a simulated hit
virtual void ProcessSimHit(Context *ctxt, proio::model::eic::SimHit *simHit, uint64_t simHitID) {}
// clean up at end of each event
virtual void EndEvent(Context *ctxt, proio::Event *event) {}
};
class ReadoutFactory {
public:
virtual Readout *Create() = 0;
};
class ReadoutManager {
public:
static void RegisterReadout(std::string name, ReadoutFactory *factory) { factories()[name] = factory; }
static class Readout *Readout(std::string name) {
if (factories().count(name))
return factories()[name]->Create();
else
return NULL;
}
static std::vector<class Readout *> AvailableReadouts() {
std::vector<class Readout *> types;
for (auto keyValuePair : factories()) {
types.push_back(keyValuePair.second->Create());
}
return types;
}
private:
static std::map<std::string, ReadoutFactory *> &factories() {
static std::map<std::string, ReadoutFactory *> f;
return f;
}
};
} // namespace digi
} // namespace topside
#define READOUT_TYPE(klass) \
namespace topside { \
namespace digi { \
class klass##Factory : public ReadoutFactory { \
public: \
klass##Factory() { ReadoutManager::RegisterReadout(#klass, this); } \
Readout *Create() { return new klass(); } \
}; \
static klass##Factory global_##klass##Factory; \
} \
}
#endif // TOPSIDE_DIGI_READOUT_TYPE_H
src/ConceptDetectors/topside/digi/ReadoutTypes/RPC.cpp 0 → 100644
View file @ bb59cd38
#include <unistd.h>
#include <iostream>
#include <stdexcept>
#include <CLHEP/Units/SystemOfUnits.h>
#include <proio/event.h>
#include <proio/model/eic/eic.pb.h>
#include "CellInfo.h"
#include "Context.h"
#include "Readout.h"
#include "global.h"
using namespace dd4hep;
using namespace proio::model::eic;
namespace topside {
namespace digi {
class RPCReadout : public Readout {
public:
RPCReadout() {}
// This readout is currently disabled because of the lack of timing. Need
// to add timing to the digitization for later analysis that relies on
// invertible 4D covariance matrices.
// virtual std::vector<std::string> Types() { return {"HcalBarrelHits", "HcalEndcapHits"}; }
virtual std::vector<std::string> Types() { return {}; }
virtual void BeginEvent(Context *, proio::Event *) { depRPC.clear(); }
virtual void ProcessSimHit(Context *, SimHit *simHit, uint64_t simHitID) {
auto globalprepos = simHit->globalprepos();
double globalPrePos[3] = {globalprepos.x() * 0.1, globalprepos.y() * 0.1, globalprepos.z() * 0.1};
auto globalpostpos = simHit->globalpostpos();
double globalPostPos[3] = {globalpostpos.x() * 0.1, globalpostpos.y() * 0.1, globalpostpos.z() * 0.1};
double localPrePos[3];
double localPostPos[3];
Segmentation seg;
{
pthread_mutex_lock(&globalMutex);
VolumeManagerContext *vmanCtxt = vman.lookupContext(simHit->volumeid());
std::string typeName =
((SensitiveDetector)vmanCtxt->volumePlacement().volume().sensitiveDetector())