UtilityApps/src/graphicalScan.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 //==========================================================================
0011 //
0012 //  Simple program to make a "CT-scan" of a detector materials
0013 //  makes series of slices perpendicular to X, Y, or Z axis
0014 //  produces a TFile with a TH2F per slice
0015 //    scans the material along a few (2*nSamples, in the 2 slice dimensions) paths across each bin
0016 //    fills histogram bins with:
0017 //       average X0/length across these paths
0018 //       average lambda/length across these paths
0019 //       for each material, fraction of path length which crosses that material
0020 //   (inspired by material scan)
0021 //
0022 //  Author     : D.Jeans, UTokyo
0023 //
0024 //==========================================================================
0025
0026 #include <TError.h>
0027 #include <TFile.h>
0028 #include <TH2F.h>
0029
0030 // Framework include files
0031 #include <DD4hep/Detector.h>
0032 #include <DD4hep/Printout.h>
0033 #include <DDRec/MaterialManager.h>
0034
0035 #include <iostream>
0036 #include <climits>
0037 #include <cerrno>
0038 #include <string>
0039 #include <map>
0040
0041 #undef NDEBUG
0042 #include <cassert>
0043
0044 using namespace dd4hep;
0045 using namespace dd4hep::rec;
0046
0047 using std::cout;
0048 using std::endl;
0049
0050 int main_wrapper(int argc, char** argv)   {
0051   struct Handler  {
0052     Handler() { SetErrorHandler(Handler::print); }
0053     static void print(int level, Bool_t abort, const char *location, const char *msg)  {
0054       if ( level > kInfo || abort ) ::printf("%s: %s\n", location, msg);
0055     }
0056     static void usage()  {
0057       std::cout << " usage: graphicalScan compact.xml axis xMin xMax yMin yMax zMin zMax nSlices nBins nSamples FieldOrMaterial OutfileName" << std::endl
0058                 << " axis (X, Y, or Z)             : perpendicular to the slices" << std::endl
0059                 << " xMin xMax yMin yMax zMin zMax : range of scans " << std::endl
0060                 << " nSlices                       : number of slices (equally spaced along chose axis)" << std::endl
0061                 << " nBins                         : number of bins along each axis of histograms" << std::endl
0062                 << " nSamples                      : the number of times each bin is sampled " << std::endl
0063                 << " FieldOrMaterial               : scan field or material? F = field, M = material, FM or MF = both" << std::endl
0064                 << " OutfileName                   : output root file name" << std::endl
0065                 << "        -> produces graphical scans of material and/or fields defined in a compact xml description"
0066                 << std::endl;
0067       exit(1);
0068     }
0069   } _handler;
0070
0071   // "axis" is the normal to the slices, which are equally spaced between the corresponding max and min
0072   // each slice has nBins x nBins in the specified range
0073   // the material in each bin is sampled along 2*nTests paths
0074
0075   if( argc != 14 ) Handler::usage();
0076
0077   std::string inFile = argv[1]; // input geometry description compact xml file
0078
0079   std::string XYZ    = argv[2]; // the axis
0080
0081   TString labx, laby;
0082   unsigned int index[3] = {99,99,99};
0083   if ( XYZ=="x" || XYZ=="X" ) {
0084     index[0] = 0; // this is the one perpendicular to slices
0085     index[1] = 2;
0086     index[2] = 1;
0087     labx="Z [cm]";
0088     laby="Y [cm]";
0089   } else if ( XYZ=="y" || XYZ=="Y" ) {
0090     index[0] = 1;
0091     index[1] = 2;
0092     index[2] = 0;
0093     labx="Z [cm]";
0094     laby="X [cm]";
0095   } else if ( XYZ=="z" || XYZ=="Z" ) {
0096     index[0] = 2;
0097     index[1] = 0;
0098     index[2] = 1;
0099     labx="X [cm]";
0100     laby="Y [cm]";
0101   } else {
0102     cout << "invalid XYZ" << endl;
0103     return -1;
0104   }
0105
0106   double x0,y0,z0,x1,y1,z1;
0107   unsigned int nslice, nbins, mm_count;
0108   std::stringstream sstr;
0109   sstr <<
0110     argv[3] << " " << argv[4] << " " << argv[5] << " " <<
0111     argv[6] << " " << argv[7] << " " << argv[8] << " " <<
0112     argv[9] << " " << argv[10] << " " << argv[11] << "NONE";
0113   sstr >> x0 >> x1 >> y0 >> y1 >> z0 >> z1 >> nslice >> nbins >> mm_count ;
0114   if ( !sstr.good() )   {
0115     Handler::usage();
0116     ::exit(EINVAL);
0117   }
0118
0119   std::string FM = argv[12];
0120   std::string outFileName = argv[13];
0121
0122   if ( x0>x1 ) { double temp=x0; x0=x1; x1=temp; }
0123   if ( y0>y1 ) { double temp=y0; y0=y1; y1=temp; }
0124   if ( z0>z1 ) { double temp=z0; z0=z1; z1=temp; }
0125
0126   if ( ! (nbins>0 && nbins<USHRT_MAX && nslice>0 && nslice<USHRT_MAX) ) {
0127     cout << "funny # bins/slices " << endl;
0128     ::exit(EINVAL);
0129   }
0130
0131   bool scanField(false);
0132   bool scanMaterial(false);
0133   if ( FM=="f" || FM=="F" ) {
0134     scanField=true;
0135   } else if ( FM=="m" || FM=="M" ) {
0136     scanMaterial=true;
0137   } else if ( FM=="fm" || FM=="FM" || FM=="mf" || FM=="MF" ) {
0138     scanField=true;
0139     scanMaterial=true;
0140   } else {
0141     cout << "invalid field/material flag: use one of f/F/m/M/fm/FM/mf/MF" << endl;
0142     return 1;
0143   }
0144
0145   double mmin[3]={x0,y0,z0};
0146   double mmax[3]={x1,y1,z1};
0147
0148
0149   //------
0150
0151   Detector& description = Detector::getInstance();
0152   description.fromCompact(inFile);
0153
0154   //-----
0155
0156   bool found_tessellated = false;
0157   TFile* f = new TFile(outFileName.c_str(),"recreate");
0158   Vector3D p0, p1; // the two points between which material is calculated
0159   MaterialManager matMgr( description.world().volume() ) ;
0160
0161   for (unsigned int isl=0; isl<nslice; isl++) { // loop over slices
0162
0163     double sz = nslice > 1 ?
0164       mmin[index[0]] + isl*( mmax[index[0]] - mmin[index[0]] )/( nslice - 1 ) :
0165       (mmin[index[0]] + mmax[index[0]])/2. ;
0166
0167     p0.array()[ index[0] ] = sz;
0168     p1.array()[ index[0] ] = sz;
0169
0170     cout << "scanning slice " << isl << " at "+XYZ+" = " << sz << endl;
0171
0172     TString dirn = "Slice"; dirn+=isl;
0173     f->mkdir(dirn);
0174     f->cd(dirn);
0175
0176     std::map < std::string , TH2F* > scanmap;
0177
0178
0179     TString hn = "slice"; hn+=isl; hn+="_X0";
0180     TString hnn = "X0 "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0181
0182     for (int j=1; j<=2; j++) {
0183       if ( mmax[index[j]] - mmin[index[j]] < 1e-4 ) {
0184         cout << "ERROR: max and min of axis are the same!" << endl;
0185         assert(0);
0186       }
0187     }
0188
0189
0190     TH2F* h2slice = 0;
0191
0192     if (scanMaterial) {
0193
0194       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0195       scanmap["x0"] = h2slice;
0196
0197       hn = "slice"; hn+=isl; hn+="_lambda";
0198       hnn = "lambda "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0199       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0200       scanmap["lambda"] = h2slice;
0201     }
0202
0203     if (scanField) {
0204       hn = "slice"; hn+=isl; hn+="_Bx";
0205       hnn = "Bx[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0206       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0207       scanmap["Bx"] = h2slice;
0208
0209       hn = "slice"; hn+=isl; hn+="_By";
0210       hnn = "By[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0211       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0212       scanmap["By"] = h2slice;
0213
0214       hn = "slice"; hn+=isl; hn+="_Bz";
0215       hnn = "Bz[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0216       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0217       scanmap["Bz"] = h2slice;
0218
0219       hn = "slice"; hn+=isl; hn+="_Ex";
0220       hnn = "Ex[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0221       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0222       scanmap["Ex"] = h2slice;
0223
0224       hn = "slice"; hn+=isl; hn+="_Ey";
0225       hnn = "Ey[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0226       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0227       scanmap["Ey"] = h2slice;
0228
0229       hn = "slice"; hn+=isl; hn+="_Ez";
0230       hnn = "Ez[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0231       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0232       scanmap["Ez"] = h2slice;
0233     }
0234
0235
0236     for (int ix=1; ix<=h2slice->GetNbinsX(); ix++) {  // loop over one axis of slice
0237
0238       double xmin = h2slice->GetXaxis()->GetBinLowEdge(ix);
0239       double xmax = h2slice->GetXaxis()->GetBinUpEdge(ix);
0240
0241       for (int iy=1; iy<=h2slice->GetNbinsY(); iy++) { // and the other axis
0242
0243         double ymin = h2slice->GetYaxis()->GetBinLowEdge(iy);
0244         double ymax = h2slice->GetYaxis()->GetBinUpEdge(iy);
0245
0246         if (scanField) {
0247           // first get b field components in centre of bin
0248           double posV[3];
0249           posV[index[0]] = sz;
0250           posV[index[1]] = 0.5*(xmin+xmax);
0251           posV[index[2]] = 0.5*(ymin+ymax);
0252
0253           double fieldV[3] ;
0254           description.field().combinedMagnetic( posV  , fieldV  ) ;
0255           scanmap["Bx"]->SetBinContent(ix, iy, fieldV[0] / dd4hep::tesla );
0256           scanmap["By"]->SetBinContent(ix, iy, fieldV[1] / dd4hep::tesla );
0257           scanmap["Bz"]->SetBinContent(ix, iy, fieldV[2] / dd4hep::tesla );
0258
0259           description.field().combinedElectric( posV  , fieldV  ) ;
0260           scanmap["Ex"]->SetBinContent(ix, iy, fieldV[0] / ( dd4hep::volt/dd4hep::meter ) );
0261           scanmap["Ey"]->SetBinContent(ix, iy, fieldV[1] / ( dd4hep::volt/dd4hep::meter ) );
0262           scanmap["Ez"]->SetBinContent(ix, iy, fieldV[2] / ( dd4hep::volt/dd4hep::meter ) );
0263         }
0264
0265         if (scanMaterial) {
0266
0267           // for this bin, estimate the material
0268           double sum_lambda(0);
0269           double sum_x0(0);
0270           double sum_length(0);
0271           std::map < std::string , float > materialmap;
0272
0273           for (unsigned int jx=0; jx<2*mm_count; jx++) {
0274             if ( jx<mm_count ) {
0275               double xcom = xmin + (1+jx)*( xmax - xmin )/(mm_count+1.);
0276               p0.array()[index[1]] = xcom;  p0.array()[index[2]] = ymin;
0277               p1.array()[index[1]] = xcom;  p1.array()[index[2]] = ymax;
0278             } else {
0279               double ycom =  ymin + (jx-mm_count+1)*( ymax - ymin )/(mm_count+1.);
0280               p0.array()[index[1]] = xmin;  p0.array()[index[2]] = ycom;
0281               p1.array()[index[1]] = xmax;  p1.array()[index[2]] = ycom;
0282             }
0283             const MaterialManager::ScanData scan = matMgr.entriesBetween(p0, p1);
0284             const auto& materials = scan.materials;
0285             const auto& places    = scan.places;
0286
0287             for( unsigned i=0,n=materials.size();i<n;++i){
0288               TGeoMaterial* mat =  materials[i].first->GetMaterial();
0289               double length = materials[i].second;
0290               sum_length += length;
0291               double nx0 = length / mat->GetRadLen();
0292               sum_x0 += nx0;
0293               double nLambda = length / mat->GetIntLen();
0294               sum_lambda += nLambda;
0295
0296               std::string mname = mat->GetName();
0297               if ( materialmap.find( mname )!=materialmap.end() ) {
0298                 materialmap[mname]+=length;
0299               } else {
0300                 materialmap[mname]=length;
0301               }
0302
0303             }
0304             if( !found_tessellated )  {
0305               for( const auto& p : places )  {
0306                 Volume volume(p.first.volume());
0307                 Solid  shape(volume.solid());
0308                 if ( shape->IsA() == TGeoTessellated::Class() )  {
0309                   found_tessellated = true;
0310                 }
0311               }
0312             }
0313           }
0314
0315           scanmap["x0"]->SetBinContent(ix, iy, sum_x0/sum_length); // normalise to cm (ie x0/cm density: indep of bin size)
0316           scanmap["lambda"]->SetBinContent(ix, iy, sum_lambda/sum_length);
0317
0318           for (  std::map < std::string , float >::iterator jj = materialmap.begin(); jj!=materialmap.end(); jj++) {
0319             if ( scanmap.find( jj->first )==scanmap.end() ) {
0320               hn = "slice"; hn+=isl; hn+="_"+jj->first;
0321               hnn = jj->first; hnn += " "+XYZ; hnn+="=";
0322               // hnn+=sz;
0323               hnn += Form("%7.3f",sz);
0324               hnn+=" [cm]";
0325               scanmap[jj->first] = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0326             }
0327             scanmap[jj->first]->SetBinContent(ix, iy, jj->second / sum_length );
0328           }
0329         } // if (scanMaterial)
0330       }
0331     }
0332
0333     for (  std::map < std::string , TH2F* >::iterator jj = scanmap.begin(); jj!=scanmap.end(); jj++) {
0334       jj->second->SetOption("zcol");
0335       jj->second->GetXaxis()->SetTitle(labx);
0336       jj->second->GetYaxis()->SetTitle(laby);
0337     }
0338   }
0339
0340   if ( scanMaterial && found_tessellated )  {
0341     const char* line = " +------------------------------------------------------------"
0342       "--------------------------------------------------------------------------------------\n";
0343     ::printf("%s",line);
0344     ::printf(" |  WARNING: Tessellated shape were encountered during the volume traversal.\n");
0345     ::printf(" |  WARNING: The results of the material scan(s) are unreliable!\n");
0346     ::printf("%s",line);
0347   }
0348   f->Write();
0349   f->Close();
0350   return 0;
0351 }
0352
0353 /// Main entry point as a program
0354 int main(int argc, char** argv)   {
0355   try  {
0356     return main_wrapper(argc, argv);
0357   }
0358   catch(const std::exception& e)  {
0359     std::cout << "Got uncaught exception: " << e.what() << std::endl;
0360   }
0361   catch (...)  {
0362     std::cout << "Got UNKNOWN uncaught exception." << std::endl;
0363   }
0364   return EINVAL;
0365 }