DD4hep  1.32.1
Detector Description Toolkit for High Energy Physics
graphicalScan.cpp
Go to the documentation of this file.
1 //==========================================================================
2 // AIDA Detector description implementation
3 //--------------------------------------------------------------------------
4 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
5 // All rights reserved.
6 //
7 // For the licensing terms see $DD4hepINSTALL/LICENSE.
8 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
9 //
10 //==========================================================================
11 //
12 // Simple program to make a "CT-scan" of a detector materials
13 // makes series of slices perpendicular to X, Y, or Z axis
14 // produces a TFile with a TH2F per slice
15 // scans the material along a few (2*nSamples, in the 2 slice dimensions) paths across each bin
16 // fills histogram bins with:
17 // average X0/length across these paths
18 // average lambda/length across these paths
19 // for each material, fraction of path length which crosses that material
20 // (inspired by material scan)
21 //
22 // Author : D.Jeans, UTokyo
23 //
24 //==========================================================================
25 
26 #include <TError.h>
27 #include <TFile.h>
28 #include <TH2F.h>
29 
30 // Framework include files
31 #include <DD4hep/Detector.h>
32 #include <DD4hep/Printout.h>
33 #include <DDRec/MaterialManager.h>
34 
35 #include <iostream>
36 #include <climits>
37 #include <cerrno>
38 #include <string>
39 #include <map>
40 
41 #undef NDEBUG
42 #include <cassert>
43 
44 using namespace dd4hep;
45 using namespace dd4hep::rec;
46 
47 using std::cout;
48 using std::endl;
49 
50 int main_wrapper(int argc, char** argv) {
51  struct Handler {
52  Handler() { SetErrorHandler(Handler::print); }
53  static void print(int level, Bool_t abort, const char *location, const char *msg) {
54  if ( level > kInfo || abort ) ::printf("%s: %s\n", location, msg);
55  }
56  static void usage() {
57  std::cout << " usage: graphicalScan compact.xml axis xMin xMax yMin yMax zMin zMax nSlices nBins nSamples FieldOrMaterial OutfileName" << std::endl
58  << " axis (X, Y, or Z) : perpendicular to the slices" << std::endl
59  << " xMin xMax yMin yMax zMin zMax : range of scans " << std::endl
60  << " nSlices : number of slices (equally spaced along chose axis)" << std::endl
61  << " nBins : number of bins along each axis of histograms" << std::endl
62  << " nSamples : the number of times each bin is sampled " << std::endl
63  << " FieldOrMaterial : scan field or material? F = field, M = material, FM or MF = both" << std::endl
64  << " OutfileName : output root file name" << std::endl
65  << " -> produces graphical scans of material and/or fields defined in a compact xml description"
66  << std::endl;
67  exit(1);
68  }
69  } _handler;
70 
71  // "axis" is the normal to the slices, which are equally spaced between the corresponding max and min
72  // each slice has nBins x nBins in the specified range
73  // the material in each bin is sampled along 2*nTests paths
74 
75  if( argc != 14 ) Handler::usage();
76 
77  std::string inFile = argv[1]; // input geometry description compact xml file
78 
79  std::string XYZ = argv[2]; // the axis
80 
81  TString labx, laby;
82  unsigned int index[3] = {99,99,99};
83  if ( XYZ=="x" || XYZ=="X" ) {
84  index[0] = 0; // this is the one perpendicular to slices
85  index[1] = 2;
86  index[2] = 1;
87  labx="Z [cm]";
88  laby="Y [cm]";
89  } else if ( XYZ=="y" || XYZ=="Y" ) {
90  index[0] = 1;
91  index[1] = 2;
92  index[2] = 0;
93  labx="Z [cm]";
94  laby="X [cm]";
95  } else if ( XYZ=="z" || XYZ=="Z" ) {
96  index[0] = 2;
97  index[1] = 0;
98  index[2] = 1;
99  labx="X [cm]";
100  laby="Y [cm]";
101  } else {
102  cout << "invalid XYZ" << endl;
103  return -1;
104  }
105 
106  double x0,y0,z0,x1,y1,z1;
107  unsigned int nslice, nbins, mm_count;
108  std::stringstream sstr;
109  sstr <<
110  argv[3] << " " << argv[4] << " " << argv[5] << " " <<
111  argv[6] << " " << argv[7] << " " << argv[8] << " " <<
112  argv[9] << " " << argv[10] << " " << argv[11] << "NONE";
113  sstr >> x0 >> x1 >> y0 >> y1 >> z0 >> z1 >> nslice >> nbins >> mm_count ;
114  if ( !sstr.good() ) {
115  Handler::usage();
116  ::exit(EINVAL);
117  }
118 
119  std::string FM = argv[12];
120  std::string outFileName = argv[13];
121 
122  if ( x0>x1 ) { double temp=x0; x0=x1; x1=temp; }
123  if ( y0>y1 ) { double temp=y0; y0=y1; y1=temp; }
124  if ( z0>z1 ) { double temp=z0; z0=z1; z1=temp; }
125 
126  if ( ! (nbins>0 && nbins<USHRT_MAX && nslice>0 && nslice<USHRT_MAX) ) {
127  cout << "funny # bins/slices " << endl;
128  ::exit(EINVAL);
129  }
130 
131  bool scanField(false);
132  bool scanMaterial(false);
133  if ( FM=="f" || FM=="F" ) {
134  scanField=true;
135  } else if ( FM=="m" || FM=="M" ) {
136  scanMaterial=true;
137  } else if ( FM=="fm" || FM=="FM" || FM=="mf" || FM=="MF" ) {
138  scanField=true;
139  scanMaterial=true;
140  } else {
141  cout << "invalid field/material flag: use one of f/F/m/M/fm/FM/mf/MF" << endl;
142  return 1;
143  }
144 
145  double mmin[3]={x0,y0,z0};
146  double mmax[3]={x1,y1,z1};
147 
148 
149  //------
150 
151  Detector& description = Detector::getInstance();
152  description.fromCompact(inFile);
153 
154  //-----
155 
156  bool found_tessellated = false;
157  TFile* f = new TFile(outFileName.c_str(),"recreate");
158  Vector3D p0, p1; // the two points between which material is calculated
159  MaterialManager matMgr( description.world().volume() ) ;
160 
161  for (unsigned int isl=0; isl<nslice; isl++) { // loop over slices
162 
163  double sz = nslice > 1 ?
164  mmin[index[0]] + isl*( mmax[index[0]] - mmin[index[0]] )/( nslice - 1 ) :
165  (mmin[index[0]] + mmax[index[0]])/2. ;
166 
167  p0.array()[ index[0] ] = sz;
168  p1.array()[ index[0] ] = sz;
169 
170  cout << "scanning slice " << isl << " at "+XYZ+" = " << sz << endl;
171 
172  TString dirn = "Slice"; dirn+=isl;
173  f->mkdir(dirn);
174  f->cd(dirn);
175 
176  std::map < std::string , TH2F* > scanmap;
177 
178 
179  TString hn = "slice"; hn+=isl; hn+="_X0";
180  TString hnn = "X0 "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
181 
182  for (int j=1; j<=2; j++) {
183  if ( mmax[index[j]] - mmin[index[j]] < 1e-4 ) {
184  cout << "ERROR: max and min of axis are the same!" << endl;
185  assert(0);
186  }
187  }
188 
189 
190  TH2F* h2slice = 0;
191 
192  if (scanMaterial) {
193 
194  h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
195  scanmap["x0"] = h2slice;
196 
197  hn = "slice"; hn+=isl; hn+="_lambda";
198  hnn = "lambda "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
199  h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
200  scanmap["lambda"] = h2slice;
201  }
202 
203  if (scanField) {
204  hn = "slice"; hn+=isl; hn+="_Bx";
205  hnn = "Bx[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
206  h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
207  scanmap["Bx"] = h2slice;
208 
209  hn = "slice"; hn+=isl; hn+="_By";
210  hnn = "By[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
211  h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
212  scanmap["By"] = h2slice;
213 
214  hn = "slice"; hn+=isl; hn+="_Bz";
215  hnn = "Bz[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
216  h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
217  scanmap["Bz"] = h2slice;
218 
219  hn = "slice"; hn+=isl; hn+="_Ex";
220  hnn = "Ex[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
221  h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
222  scanmap["Ex"] = h2slice;
223 
224  hn = "slice"; hn+=isl; hn+="_Ey";
225  hnn = "Ey[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
226  h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
227  scanmap["Ey"] = h2slice;
228 
229  hn = "slice"; hn+=isl; hn+="_Ez";
230  hnn = "Ez[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
231  h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
232  scanmap["Ez"] = h2slice;
233  }
234 
235 
236  for (int ix=1; ix<=h2slice->GetNbinsX(); ix++) { // loop over one axis of slice
237 
238  double xmin = h2slice->GetXaxis()->GetBinLowEdge(ix);
239  double xmax = h2slice->GetXaxis()->GetBinUpEdge(ix);
240 
241  for (int iy=1; iy<=h2slice->GetNbinsY(); iy++) { // and the other axis
242 
243  double ymin = h2slice->GetYaxis()->GetBinLowEdge(iy);
244  double ymax = h2slice->GetYaxis()->GetBinUpEdge(iy);
245 
246  if (scanField) {
247  // first get b field components in centre of bin
248  double posV[3];
249  posV[index[0]] = sz;
250  posV[index[1]] = 0.5*(xmin+xmax);
251  posV[index[2]] = 0.5*(ymin+ymax);
252 
253  double fieldV[3] ;
254  description.field().combinedMagnetic( posV , fieldV ) ;
255  scanmap["Bx"]->SetBinContent(ix, iy, fieldV[0] / dd4hep::tesla );
256  scanmap["By"]->SetBinContent(ix, iy, fieldV[1] / dd4hep::tesla );
257  scanmap["Bz"]->SetBinContent(ix, iy, fieldV[2] / dd4hep::tesla );
258 
259  description.field().combinedElectric( posV , fieldV ) ;
260  scanmap["Ex"]->SetBinContent(ix, iy, fieldV[0] / ( dd4hep::volt/dd4hep::meter ) );
261  scanmap["Ey"]->SetBinContent(ix, iy, fieldV[1] / ( dd4hep::volt/dd4hep::meter ) );
262  scanmap["Ez"]->SetBinContent(ix, iy, fieldV[2] / ( dd4hep::volt/dd4hep::meter ) );
263  }
264 
265  if (scanMaterial) {
266 
267  // for this bin, estimate the material
268  double sum_lambda(0);
269  double sum_x0(0);
270  double sum_length(0);
271  std::map < std::string , float > materialmap;
272 
273  for (unsigned int jx=0; jx<2*mm_count; jx++) {
274  if ( jx<mm_count ) {
275  double xcom = xmin + (1+jx)*( xmax - xmin )/(mm_count+1.);
276  p0.array()[index[1]] = xcom; p0.array()[index[2]] = ymin;
277  p1.array()[index[1]] = xcom; p1.array()[index[2]] = ymax;
278  } else {
279  double ycom = ymin + (jx-mm_count+1)*( ymax - ymin )/(mm_count+1.);
280  p0.array()[index[1]] = xmin; p0.array()[index[2]] = ycom;
281  p1.array()[index[1]] = xmax; p1.array()[index[2]] = ycom;
282  }
283  const MaterialManager::ScanData scan = matMgr.entriesBetween(p0, p1);
284  const auto& materials = scan.materials;
285  const auto& places = scan.places;
286 
287  for( unsigned i=0,n=materials.size();i<n;++i){
288  TGeoMaterial* mat = materials[i].first->GetMaterial();
289  double length = materials[i].second;
290  sum_length += length;
291  double nx0 = length / mat->GetRadLen();
292  sum_x0 += nx0;
293  double nLambda = length / mat->GetIntLen();
294  sum_lambda += nLambda;
295 
296  std::string mname = mat->GetName();
297  if ( materialmap.find( mname )!=materialmap.end() ) {
298  materialmap[mname]+=length;
299  } else {
300  materialmap[mname]=length;
301  }
302 
303  }
304  if( !found_tessellated ) {
305  for( const auto& p : places ) {
306  Volume volume(p.first.volume());
307  Solid shape(volume.solid());
308  if ( shape->IsA() == TGeoTessellated::Class() ) {
309  found_tessellated = true;
310  }
311  }
312  }
313  }
314 
315  scanmap["x0"]->SetBinContent(ix, iy, sum_x0/sum_length); // normalise to cm (ie x0/cm density: indep of bin size)
316  scanmap["lambda"]->SetBinContent(ix, iy, sum_lambda/sum_length);
317 
318  for ( std::map < std::string , float >::iterator jj = materialmap.begin(); jj!=materialmap.end(); jj++) {
319  if ( scanmap.find( jj->first )==scanmap.end() ) {
320  hn = "slice"; hn+=isl; hn+="_"+jj->first;
321  hnn = jj->first; hnn += " "+XYZ; hnn+="=";
322  // hnn+=sz;
323  hnn += Form("%7.3f",sz);
324  hnn+=" [cm]";
325  scanmap[jj->first] = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
326  }
327  scanmap[jj->first]->SetBinContent(ix, iy, jj->second / sum_length );
328  }
329  } // if (scanMaterial)
330  }
331  }
332 
333  for ( std::map < std::string , TH2F* >::iterator jj = scanmap.begin(); jj!=scanmap.end(); jj++) {
334  jj->second->SetOption("zcol");
335  jj->second->GetXaxis()->SetTitle(labx);
336  jj->second->GetYaxis()->SetTitle(laby);
337  }
338  }
339 
340  if ( scanMaterial && found_tessellated ) {
341  const char* line = " +------------------------------------------------------------"
342  "--------------------------------------------------------------------------------------\n";
343  ::printf("%s",line);
344  ::printf(" | WARNING: Tessellated shape were encountered during the volume traversal.\n");
345  ::printf(" | WARNING: The results of the material scan(s) are unreliable!\n");
346  ::printf("%s",line);
347  }
348  f->Write();
349  f->Close();
350  return 0;
351 }
352 
354 int main(int argc, char** argv) {
355  try {
356  return main_wrapper(argc, argv);
357  }
358  catch(const std::exception& e) {
359  std::cout << "Got uncaught exception: " << e.what() << std::endl;
360  }
361  catch (...) {
362  std::cout << "Got UNKNOWN uncaught exception." << std::endl;
363  }
364  return EINVAL;
365 }
dd4hep::rec::MaterialManager::ScanData
Definition: MaterialManager.h:59
dd4hep::OverlayedField::combinedElectric
void combinedElectric(const Position &pos, double *field) const
Returns the 3 electric field components (x, y, z) if many components are present.
Definition: Fields.cpp:155
dd4hep::Detector::world
virtual DetElement world() const =0
Return reference to the top-most (world) detector element.
dd4hep::rec
Namespace for the reconstruction part of the AIDA detector description toolkit.
Definition: SurfaceInstaller.h:28
dd4hep::rec::Vector3D::array
double * array()
direct access to data as double* - allows modification
Definition: Vector3D.h:216
dd4hep::rec::Vector3D
Definition: Vector3D.h:32
Detector.h
dd4hep::exception
void exception(const std::string &src, const std::string &msg)
Definition: RootDictionary.h:69
dd4hep::Volume::solid
Solid solid() const
Access to Solid (Shape)
Definition: Volumes.cpp:1252
main
int main(int argc, char **argv)
Main entry point as a program.
Definition: graphicalScan.cpp:354
dd4hep::Solid_type< TGeoShape >
dd4hep::detail::printf
std::size_t printf(const char *fmt,...)
Definition: Printout.cpp:97
dd4hep::OverlayedField::combinedMagnetic
void combinedMagnetic(const Position &pos, double *field) const
Returns the 3 magnetic field components (x, y, z) if many components are present.
Definition: Fields.cpp:161
dd4hep::Detector::getInstance
static Detector & getInstance(const std::string &name="default")
—Factory method----—
Definition: DetectorImp.cpp:132
dd4hep::Volume
Handle class holding a placed volume (also called physical volume)
Definition: Volumes.h:371
dd4hep::DetElement::volume
Volume volume() const
Access to the logical volume of the detector element's placement.
Definition: DetElement.cpp:352
dd4hep::Detector::field
virtual OverlayedField field() const =0
Return handle to the combined electromagentic field description.
dd4hep::rec::MaterialManager
Definition: MaterialManager.h:41
dd4hep::Detector::fromCompact
virtual void fromCompact(const std::string &fname, DetectorBuildType type=BUILD_DEFAULT)=0
Deprecated call (use fromXML): Read compact geometry description or alignment file.
main_wrapper
int main_wrapper(int argc, char **argv)
Definition: graphicalScan.cpp:50
dd4hep::rec::MaterialManager::entriesBetween
const ScanData entriesBetween(const Vector3D &p0, const Vector3D &p1, double eps=MaterialManager::epsilon)
As above, but optionally allow access to traversed placements.
Definition: MaterialManager.cpp:40
MaterialManager.h
dd4hep
Namespace for the AIDA detector description toolkit.
Definition: AlignmentsCalib.h:28
dd4hep::rec::MaterialManager::ScanData::places
const PlacementVec & places
Definition: MaterialManager.h:62
dd4hep::Detector
The main interface to the dd4hep detector description package.
Definition: Detector.h:90
usage
void usage(std::string argv0)
Definition: listcomponents.cpp:41
dd4hep::rec::MaterialManager::ScanData::materials
const MaterialVec & materials
Definition: MaterialManager.h:61
Printout.h
Generated on Sun Sep 7 2025 08:57:16 for DD4hep by   doxygen 1.8.18