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   TFile* f = new TFile(outFileName.c_str(),"recreate");
0157
0158   Vector3D p0, p1; // the two points between which material is calculated
0159
0160   MaterialManager matMgr(description.world().volume() ) ;
0161
0162   for (unsigned int isl=0; isl<nslice; isl++) { // loop over slices
0163
0164     double sz = nslice > 1 ?
0165       mmin[index[0]] + isl*( mmax[index[0]] - mmin[index[0]] )/( nslice - 1 ) :
0166       (mmin[index[0]] + mmax[index[0]])/2. ;
0167
0168     p0.array()[ index[0] ] = sz;
0169     p1.array()[ index[0] ] = sz;
0170
0171     cout << "scanning slice " << isl << " at "+XYZ+" = " << sz << endl;
0172
0173     TString dirn = "Slice"; dirn+=isl;
0174     f->mkdir(dirn);
0175     f->cd(dirn);
0176
0177     std::map < std::string , TH2F* > scanmap;
0178
0179
0180     TString hn = "slice"; hn+=isl; hn+="_X0";
0181     TString hnn = "X0 "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0182
0183     for (int j=1; j<=2; j++) {
0184       if ( mmax[index[j]] - mmin[index[j]] < 1e-4 ) {
0185         cout << "ERROR: max and min of axis are the same!" << endl;
0186         assert(0);
0187       }
0188     }
0189
0190
0191     TH2F* h2slice = 0;
0192
0193     if (scanMaterial) {
0194
0195       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0196       scanmap["x0"] = h2slice;
0197
0198       hn = "slice"; hn+=isl; hn+="_lambda";
0199       hnn = "lambda "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0200       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0201       scanmap["lambda"] = h2slice;
0202     }
0203
0204     if (scanField) {
0205       hn = "slice"; hn+=isl; hn+="_Bx";
0206       hnn = "Bx[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0207       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0208       scanmap["Bx"] = h2slice;
0209
0210       hn = "slice"; hn+=isl; hn+="_By";
0211       hnn = "By[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0212       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0213       scanmap["By"] = h2slice;
0214
0215       hn = "slice"; hn+=isl; hn+="_Bz";
0216       hnn = "Bz[T] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0217       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0218       scanmap["Bz"] = h2slice;
0219
0220       hn = "slice"; hn+=isl; hn+="_Ex";
0221       hnn = "Ex[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0222       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0223       scanmap["Ex"] = h2slice;
0224
0225       hn = "slice"; hn+=isl; hn+="_Ey";
0226       hnn = "Ey[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0227       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0228       scanmap["Ey"] = h2slice;
0229
0230       hn = "slice"; hn+=isl; hn+="_Ez";
0231       hnn = "Ez[V/m] "; hnn += XYZ; hnn+="="; hnn += Form("%7.3f",sz); hnn+=" [cm]";
0232       h2slice = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0233       scanmap["Ez"] = h2slice;
0234     }
0235
0236
0237     for (int ix=1; ix<=h2slice->GetNbinsX(); ix++) {  // loop over one axis of slice
0238
0239       double xmin = h2slice->GetXaxis()->GetBinLowEdge(ix);
0240       double xmax = h2slice->GetXaxis()->GetBinUpEdge(ix);
0241
0242       for (int iy=1; iy<=h2slice->GetNbinsY(); iy++) { // and the other axis
0243
0244         double ymin = h2slice->GetYaxis()->GetBinLowEdge(iy);
0245         double ymax = h2slice->GetYaxis()->GetBinUpEdge(iy);
0246
0247         if (scanField) {
0248           // first get b field components in centre of bin
0249           double posV[3];
0250           posV[index[0]] = sz;
0251           posV[index[1]] = 0.5*(xmin+xmax);
0252           posV[index[2]] = 0.5*(ymin+ymax);
0253
0254           double fieldV[3] ;
0255           description.field().combinedMagnetic( posV  , fieldV  ) ;
0256           scanmap["Bx"]->SetBinContent(ix, iy, fieldV[0] / dd4hep::tesla );
0257           scanmap["By"]->SetBinContent(ix, iy, fieldV[1] / dd4hep::tesla );
0258           scanmap["Bz"]->SetBinContent(ix, iy, fieldV[2] / dd4hep::tesla );
0259
0260           description.field().combinedElectric( posV  , fieldV  ) ;
0261           scanmap["Ex"]->SetBinContent(ix, iy, fieldV[0] / ( dd4hep::volt/dd4hep::meter ) );
0262           scanmap["Ey"]->SetBinContent(ix, iy, fieldV[1] / ( dd4hep::volt/dd4hep::meter ) );
0263           scanmap["Ez"]->SetBinContent(ix, iy, fieldV[2] / ( dd4hep::volt/dd4hep::meter ) );
0264         }
0265
0266         if (scanMaterial) {
0267
0268           // for this bin, estimate the material
0269           double sum_lambda(0);
0270           double sum_x0(0);
0271           double sum_length(0);
0272
0273           std::map < std::string , float > materialmap;
0274
0275           for (unsigned int jx=0; jx<2*mm_count; jx++) {
0276             if ( jx<mm_count ) {
0277               double xcom = xmin + (1+jx)*( xmax - xmin )/(mm_count+1.);
0278               p0.array()[index[1]] = xcom;  p0.array()[index[2]] = ymin;
0279               p1.array()[index[1]] = xcom;  p1.array()[index[2]] = ymax;
0280             } else {
0281               double ycom =  ymin + (jx-mm_count+1)*( ymax - ymin )/(mm_count+1.);
0282               p0.array()[index[1]] = xmin;  p0.array()[index[2]] = ycom;
0283               p1.array()[index[1]] = xmax;  p1.array()[index[2]] = ycom;
0284             }
0285
0286             const MaterialVec& materials = matMgr.materialsBetween(p0, p1);
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
0305           }
0306
0307           scanmap["x0"]->SetBinContent(ix, iy, sum_x0/sum_length); // normalise to cm (ie x0/cm density: indep of bin size)
0308           scanmap["lambda"]->SetBinContent(ix, iy, sum_lambda/sum_length);
0309
0310           for (  std::map < std::string , float >::iterator jj = materialmap.begin(); jj!=materialmap.end(); jj++) {
0311             if ( scanmap.find( jj->first )==scanmap.end() ) {
0312               hn = "slice"; hn+=isl; hn+="_"+jj->first;
0313               hnn = jj->first; hnn += " "+XYZ; hnn+="=";
0314               // hnn+=sz;
0315               hnn += Form("%7.3f",sz);
0316               hnn+=" [cm]";
0317               scanmap[jj->first] = new TH2F( hn, hnn, nbins, mmin[index[1]], mmax[index[1]], nbins, mmin[index[2]], mmax[index[2]] );
0318             }
0319             scanmap[jj->first]->SetBinContent(ix, iy, jj->second / sum_length );
0320           }
0321
0322         } // if (scanMaterial)
0323
0324       }
0325     }
0326
0327     for (  std::map < std::string , TH2F* >::iterator jj = scanmap.begin(); jj!=scanmap.end(); jj++) {
0328       jj->second->SetOption("zcol");
0329       jj->second->GetXaxis()->SetTitle(labx);
0330       jj->second->GetYaxis()->SetTitle(laby);
0331     }
0332   }
0333   f->Write();
0334   f->Close();
0335   return 0;
0336 }
0337
0338 /// Main entry point as a program
0339 int main(int argc, char** argv)   {
0340   try  {
0341     return main_wrapper(argc, argv);
0342   }
0343   catch(const std::exception& e)  {
0344     std::cout << "Got uncaught exception: " << e.what() << std::endl;
0345   }
0346   catch (...)  {
0347     std::cout << "Got UNKNOWN uncaught exception." << std::endl;
0348   }
0349   return EINVAL;
0350 }