drich-dev/src/generate_material_tables.cpp

0001 #include <cstdlib>
0002 #include <iostream>
0003 #include <fmt/format.h>
0004
0005 #include "g4dRIChOptics.hh"
0006 #include "surfaceEnums.h"
0007
0008 #include <TCanvas.h>
0009 #include <TGraphErrors.h>
0010 #include <TAxis.h>
0011 #include <TFile.h>
0012
0013 ///////////////////////////////////
0014 // SETTINGS
0015 const int    aerOptModel  = 3;      // aerogel optical model used to estimate the refractive Index
0016 const double filter_thr   = 300*nm; // wavelength filter cutoff
0017 const bool   vsWavelength = true;   // if true, make plots vs. wavelength
0018 const std::string xmlOutput      = "out/optical_materials_drich.xml";
0019 const TString     root_file_name = "out/optical_materials_drich.root";
0020 ///////////////////////////////////
0021
0022 // other global vars
0023 std::FILE *xmlFile;
0024
0025 // ===========================================================================
0026 // structures for preferred units
0027 class UnitDef {
0028   public:
0029     UnitDef(double divisor_, std::string name_) : divisor(divisor_), name(name_), xml(""), title("") {
0030       if(name!="") {
0031         xml   = "*" + name;
0032         title = " [" + name + "]";
0033       }
0034     }
0035     double divisor;
0036     std::string name, xml, title;
0037 };
0038
0039
0040 // ===========================================================================
0041 // extends g4dRIChOptics class with plots and XML printing
0042 template<class MAT> class MaterialTable {
0043   public:
0044     MAT *mpt;
0045     std::string name;
0046     MaterialTable(MAT *mpt_, std::string name_) : mpt(mpt_), name(name_) {
0047       if(name=="C2F6" or name=="C4F10") PreferredUnits = &PreferredUnits2;
0048       else PreferredUnits = &PreferredUnits1;
0049     }
0050     ~MaterialTable() { if(mpt!=nullptr) delete mpt; };
0051
0052     // print XML matrices, for `optical_materials.xml`
0053     void PrintXML(bool isSurface=false, G4String detectorName="DRICH") {
0054       fmt::print(xmlFile,"\n<!-- {:_^60} -->\n\n",name);
0055       // function to print a row
0056       auto PrintRow = [] (int indentation, UnitDef units) {
0057         return [indentation,&units] (G4double energy, G4double value) {
0058           fmt::print(xmlFile,"{:{}}{:<#.6g}{}   {:>#.7g}{}\n",
0059               "",                  indentation,
0060               energy/eV,           "*eV",
0061               value/units.divisor, units.xml
0062               );
0063         };
0064       };
0065       if(isSurface) {
0066         auto surf = mpt->getSurface();
0067         fmt::print(xmlFile,"{:4}<opticalsurface name=\"{}_{}\" model=\"{}\" finish=\"{}\" type=\"{}\">\n",
0068             "",
0069             mpt->getLogicalVName(),
0070             detectorName,
0071             surfaceEnum::GetModel(surf),
0072             surfaceEnum::GetFinish(surf),
0073             surfaceEnum::GetType(surf)
0074             );
0075         for(const auto& propName : mpt->getMaterialPropertyNames()) {
0076           fmt::print(xmlFile,"{:6}<property name=\"{}\" coldim=\"2\" values=\"\n", "", propName);
0077           mpt->loopMaterialPropertyTable(propName,PrintRow(8,PreferredUnits->at(propName)));
0078           fmt::print(xmlFile,"{:8}\"/>\n","");
0079         }
0080         fmt::print(xmlFile,"{:4}</opticalsurface>\n","");
0081       } else {
0082         for(const auto& propName : mpt->getMaterialPropertyNames()) {
0083           fmt::print(xmlFile,"{:4}<matrix name=\"{}__{}_{}\" coldim=\"2\" values=\"\n",
0084               "",
0085               propName,
0086               mpt->getMaterialName(),
0087               detectorName
0088               );
0089           mpt->loopMaterialPropertyTable(propName,PrintRow(6,PreferredUnits->at(propName)));
0090           fmt::print(xmlFile,"{:6}\"/>\n","");
0091         }
0092       }
0093     } // PrintXML
0094
0095     // draw plots of material property tables
0096     void DrawPlots() {
0097       TString canvName("materialProperties_"+name);
0098       TString pngName("out/"+canvName+".png");
0099       int ncols = 2;
0100       int nrows = 1+(mpt->getMaterialPropertyTableSize()-1)/ncols;
0101       if(nrows==1) ncols=mpt->getMaterialPropertyTableSize();
0102       auto canv = new TCanvas(canvName,canvName,ncols*800,nrows*600);
0103       canv->Divide(ncols,nrows);
0104       int pad=0;
0105       for(const auto& propName : mpt->getMaterialPropertyNames()) {
0106         fmt::print("Plotting property {}\n", propName);
0107         auto units = PreferredUnits->at(propName);
0108         canv->cd(++pad);
0109         canv->GetPad(pad)->SetGrid(1,1);
0110         canv->GetPad(pad)->SetLeftMargin(0.2);
0111         canv->GetPad(pad)->SetRightMargin(0.15);
0112         auto graph = new TGraphErrors();
0113         graph->SetName(TString("graph_"+name+"_"+propName));
0114         std::string xTitle = vsWavelength ? "wavelength [nm]" : "energy [eV]";
0115         graph->SetTitle(TString(name+" "+propName+units.title+";"+xTitle));
0116         auto makeGraph = [&graph,&units] (G4double energy, G4double value) {
0117           auto xval = vsWavelength ? g4dRIChOptics::e2wl(energy)/nm : energy/eV;
0118           graph->AddPoint(xval, value / units.divisor);
0119         };
0120         mpt->loopMaterialPropertyTable(propName,makeGraph, vsWavelength);
0121         graph->SetMarkerStyle(kFullCircle);
0122         graph->SetMarkerColor(kAzure);
0123         graph->GetXaxis()->SetLabelSize(0.05);
0124         graph->GetYaxis()->SetLabelSize(0.05);
0125         graph->GetXaxis()->SetTitleOffset(1.3);
0126         graph->Draw("AP");
0127         graph->Write();
0128       }
0129       canv->SaveAs(pngName);
0130     }
0131
0132     // sets of preferred units
0133     const std::map<G4String,UnitDef> *PreferredUnits;
0134     const std::map<G4String,UnitDef> PreferredUnits1 = {
0135       { "RINDEX",          UnitDef( 1.,    ""      )},
0136       { "GROUPVEL",        UnitDef( mm/ns, "mm/ns" )},
0137       { "RAYLEIGH",        UnitDef( mm,    "mm"    )},
0138       { "ABSLENGTH",       UnitDef( mm,    "mm"    )},
0139       { "REFLECTIVITY",    UnitDef( 1.,    ""      )},
0140       { "REALRINDEX",      UnitDef( 1.,    ""      )},
0141       { "IMAGINARYRINDEX", UnitDef( 1.,    ""      )},
0142       { "EFFICIENCY",      UnitDef( 1.,    ""      )}
0143     };
0144     const std::map<G4String,UnitDef> PreferredUnits2 = {
0145       { "RINDEX",          UnitDef( 1.,    ""      )},
0146       { "GROUPVEL",        UnitDef( mm/ns, "mm/ns" )},
0147       { "RAYLEIGH",        UnitDef( m,     "m"     )},
0148       { "ABSLENGTH",       UnitDef( m,     "m"     )},
0149       { "REFLECTIVITY",    UnitDef( 1.,    ""      )},
0150       { "REALRINDEX",      UnitDef( 1.,    ""      )},
0151       { "IMAGINARYRINDEX", UnitDef( 1.,    ""      )},
0152       { "EFFICIENCY",      UnitDef( 1.,    ""      )}
0153     };
0154
0155 }; // class MaterialTable
0156
0157 // ===========================================================================
0158
0159 int main(int argc, char** argv) {
0160
0161   // start XML file and ROOT file
0162   xmlFile = std::fopen(xmlOutput.c_str(),"w");
0163   auto root_file = new TFile(root_file_name,"RECREATE");
0164
0165   // build detector by text file
0166   fmt::print("[+] read model text file\n");
0167   G4tgbVolumeMgr *volmgr = G4tgbVolumeMgr::GetInstance();
0168   auto model_file = G4String("text/drich-materials.txt");
0169   volmgr->AddTextFile(model_file);
0170   fmt::print("[+] construct detector from text file\n");
0171   G4VPhysicalVolume *vesselPhysVol = volmgr->ReadAndConstructDetector();
0172   fmt::print("[+] done construction\n");
0173
0174   // produce material property tables ///////////////////////
0175
0176   // aerogel
0177   MaterialTable Aerogel(new g4dRIChAerogel("Aerogel"),"Aerogel");
0178   Aerogel.mpt->setOpticalParams(aerOptModel);
0179   Aerogel.PrintXML();
0180   Aerogel.DrawPlots();
0181
0182   // acrylic filter
0183   MaterialTable Acrylic(new g4dRIChFilter("Acrylic"),"Acrylic");
0184   Acrylic.mpt->setOpticalParams(filter_thr);
0185   Acrylic.PrintXML();
0186   Acrylic.DrawPlots();
0187
0188   // gas
0189   // - C2F6
0190   MaterialTable C2F6(new g4dRIChGas("C2F6"),"C2F6");
0191   C2F6.mpt->setOpticalParams();
0192   C2F6.PrintXML();
0193   C2F6.DrawPlots();
0194   // - C4F10
0195   MaterialTable C4F10(new g4dRIChGas("C4F10"),"C4F10");
0196   C4F10.mpt->setOpticalParams();
0197   C4F10.PrintXML(false,"PFRICH");
0198   C4F10.DrawPlots();
0199
0200   // mirror surface
0201   MaterialTable MirrorSurface(new g4dRIChMirror("MirrorSurface"),"MirrorSurface");
0202   MirrorSurface.mpt->setOpticalParams("ciDRICH");
0203   MirrorSurface.PrintXML(true);
0204   MirrorSurface.DrawPlots();
0205
0206   // photo sensor surface
0207   MaterialTable SensorSurface(new g4dRIChPhotosensor("SensorSurface"),"SensorSurface");
0208   SensorSurface.mpt->setOpticalParams("ciDRICH");
0209   SensorSurface.PrintXML(true);
0210   SensorSurface.DrawPlots();
0211
0212   // cleanup
0213   std::fclose(xmlFile);
0214   root_file->Close();
0215   fmt::print("\nwrote XML nodes to {}\n\n",xmlOutput);
0216
0217 } // main