tpc/jobB/DistortionCorrectionMatrixInversion.C

0001 #include <g4detectors/PHG4CylinderGeomContainer.h>
0002 #include <g4detectors/PHG4CylinderGeom.h>
0003
0004 #include <tpccalib/TpcSpaceChargeMatrixInversion.h>
0005 #include <tpccalib/TpcSpaceChargeReconstructionHelper.h>
0006
0007 #include <micromegas/CylinderGeomMicromegas.h>
0008 #include <micromegas/MicromegasDefs.h>
0009
0010 #include <trackreco/MakeActsGeometry.h>
0011
0012 #include <ffamodules/CDBInterface.h>
0013
0014 #include <fun4all/Fun4AllRunNodeInputManager.h>
0015 #include <fun4all/Fun4AllServer.h>
0016
0017 #include <phool/getClass.h>
0018 #include <phool/recoConsts.h>
0019
0020 #include <TString.h>
0021
0022 #include <cstdio>
0023 #include <sstream>
0024 #include <string>
0025
0026 R__LOAD_LIBRARY(libffamodules.so)
0027 R__LOAD_LIBRARY(libfun4all.so)
0028 R__LOAD_LIBRARY(libtpccalib.so)
0029 R__LOAD_LIBRARY(libtrack_reco.so)
0030
0031 namespace
0032 {
0033
0034   // range adaptor to be able to use range-based for loop
0035   template<class T> class range_adaptor
0036   {
0037     public:
0038     explicit range_adaptor( const T& range ):m_range(range){}
0039     const typename T::first_type& begin() {return m_range.first;}
0040     const typename T::second_type& end() {return m_range.second;}
0041     private:
0042     T m_range;
0043   };
0044
0045   // generic pair printout
0046   template<class T, class U>
0047   std::ostream& operator << (std::ostream& out, const std::pair<T, U> pair)
0048   {
0049     out << "{" << pair.first << "," << pair.second << "}";
0050     return out;
0051   }
0052
0053   // generic container printout
0054   template<template<typename,typename> class Container, class T, class A>
0055     std::ostream& operator << (std::ostream& out, const Container<T,A>& container)
0056   {
0057     out << "{";
0058     bool first = true;
0059     for( const auto& value:container )
0060     {
0061       if( !first ) { out << ", "; };
0062       first = false;
0063       out << value;
0064     }
0065     out << "}";
0066     return out;
0067   }
0068 }
0069
0070 //_______________________________________________
0071 // load tpot geometry ranges from current calibrations
0072 void load_tpot_geometry( int runnumber )
0073 {
0074
0075   // initialization
0076   auto *rc = recoConsts::instance();
0077   auto *se = Fun4AllServer::instance();
0078   auto *topNode = se->topNode();
0079
0080   rc->set_IntFlag("RUNNUMBER", runnumber);
0081   rc->set_IntFlag("RUNSEGMENT", 0);
0082   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0083   rc->set_uint64Flag("TIMESTAMP", runnumber);
0084
0085   // load geometry file
0086   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0087   std::cout << "Geometry - geofile: " <<  geofile.c_str() << std::endl;
0088   auto *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0089   ingeo->AddFile(geofile);
0090   ingeo->run(0);
0091
0092   // acts geometry initialization
0093   auto *geom = new MakeActsGeometry;
0094   geom->set_mvtx_applymisalign(true);
0095   geom->InitRun( topNode );
0096
0097   // get relevant nodes
0098   // micromegas geometry
0099   auto *m_micromegas_geomcontainer = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0100
0101   // ACTS geometry
0102   auto *m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0103
0104   using range_list_t = std::vector<CylinderGeomMicromegas::range_t>;
0105   range_list_t theta_range_list;
0106   range_list_t phi_range_list;
0107
0108   // loop over layers
0109   for( const auto& [layer, layergeom] : range_adaptor( m_micromegas_geomcontainer->get_begin_end() ) )
0110   {
0111     // sanity
0112     assert(layer == layergeom->get_layer());
0113
0114     auto *micromegas_layergeom = static_cast<CylinderGeomMicromegas*>(layergeom);
0115
0116     // tiles
0117     const unsigned int tile_count = micromegas_layergeom->get_tiles_count();
0118
0119     // segmentation
0120     const auto segmentation = micromegas_layergeom->get_segmentation_type();
0121
0122     for( unsigned int itile = 0; itile < tile_count; ++itile )
0123     {
0124       switch(segmentation)
0125       {
0126         case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0127         phi_range_list.emplace_back(micromegas_layergeom->get_phi_range(itile, m_tGeometry));
0128         break;
0129
0130         case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0131         theta_range_list.emplace_back(micromegas_layergeom->get_theta_range(itile, m_tGeometry));
0132         break;
0133       }
0134     }
0135   }
0136
0137   // calculate inner phi range for each sector
0138   auto get_phi_range = [phi_range_list]( const std::vector<int>& indexes )
0139   {
0140     CylinderGeomMicromegas::range_t out{0,0};
0141     for(const auto& i:indexes)
0142     {
0143       const auto& window = phi_range_list[i];
0144       if( out.first == 0 || window.first > out.first ) out.first = window.first;
0145       if( out.second == 0 || window.second < out.second ) out.second = window.second;
0146     }
0147     return out;
0148   };
0149
0150   // calculate phi ranges
0151   // central (= bottommost) sector corresponds to tiles 0 to 3
0152   const auto phi_range_central = get_phi_range( {0,1,2,3} );
0153
0154   // east sector corresponds to tiles 4 and 5
0155   const auto phi_range_east = get_phi_range( {4,5} );
0156
0157   // west sector corresponds to tiles 6 and 7
0158   const auto phi_range_west = get_phi_range( {6,7} );
0159
0160   // print
0161   std::cout << "phi_range_central=" << phi_range_central << std::endl;
0162   std::cout << "phi_range_east=" << phi_range_east << std::endl;
0163   std::cout << "phi_range_west=" << phi_range_west << std::endl;
0164   std::cout << std::endl;
0165
0166   // assign
0167   TpcSpaceChargeReconstructionHelper::set_phi_range_central(phi_range_central);
0168   TpcSpaceChargeReconstructionHelper::set_phi_range_east(phi_range_east);
0169   TpcSpaceChargeReconstructionHelper::set_phi_range_west(phi_range_west);
0170
0171   // generate theta ranges
0172   const range_list_t theta_range_central{theta_range_list[0], theta_range_list[1], theta_range_list[2], theta_range_list[3]};
0173   const range_list_t theta_range_east{theta_range_list[4], theta_range_list[5]};
0174   const range_list_t theta_range_west{theta_range_list[6], theta_range_list[7]};
0175
0176   // print
0177   std::cout << "theta_range_central=" << theta_range_central << std::endl;
0178   std::cout << "theta_range_east=" << theta_range_east << std::endl;
0179   std::cout << "theta_range_west=" << theta_range_west << std::endl;
0180
0181   // assign
0182   TpcSpaceChargeReconstructionHelper::set_theta_range_central(theta_range_central);
0183   TpcSpaceChargeReconstructionHelper::set_theta_range_east(theta_range_east);
0184   TpcSpaceChargeReconstructionHelper::set_theta_range_west(theta_range_west);
0185 }
0186
0187 //_______________________________________________
0188 // get list of files matching selection
0189 std::vector<std::string> list_files( const std::string& selection )
0190 {
0191   std::vector<std::string> out;
0192
0193   std::cout << "list_files - selection: " << selection.c_str() << std::endl;
0194   if( selection.empty() ) return out;
0195
0196   const std::string command = std::string("ls -1 ") + selection;
0197   auto *tmp = popen( command.c_str(), "r" );
0198   char line[512];
0199   while( fgets( line, 512, tmp ) )
0200   {
0201
0202     std::istringstream istr( line );
0203     std::string filename;
0204     istr >> filename;
0205
0206     if( filename.empty() ) continue;
0207     if( access( filename.c_str(), R_OK ) ) continue;
0208
0209     out.push_back( filename );
0210   }
0211   pclose( tmp );
0212   return out;
0213 }
0214
0215 //_______________________________________________
0216 void DistortionCorrectionMatrixInversion()
0217 {
0218
0219   // input files
0220   /*
0221    * this is the list of distortion correction matrices files coming out from Job A
0222    * that needs to be inverted, to get track-based, beam-induced distortions inside
0223    * TPOT acceptance
0224    */
0225   std::string tag = "_flat_genfit_truth_notpc_distorted-new";
0226   std::string inputFile = "DST/CONDOR" + tag + "/TpcSpaceChargeMatrices" + tag + std::string("_*.root");
0227
0228   // Central membrane distortion corrections
0229   /*
0230    * this is the 2D distortion corrections measured at the central membrane using diffuse lasers
0231    * it is used to extrapolate the distortions measured in the TPOT acceptance to the rest of the TPC acceptance
0232    * see: https://indico.bnl.gov/event/22887/contributions/90413/attachments/54020/92443/distortion_extrapolation_hp.pdf
0233    */
0234   const std::string inputfile_cm = "distortion_maps/average_minus_static_distortion_cm.root";
0235
0236   // output file
0237   std::string outputFile = "Rootfiles/Distortions_full" + tag + "_mm.root";
0238
0239   std::cout << "DistortionCorrectionMatrixInversion - inputFile: " << inputFile.c_str() << std::endl;
0240   std::cout << "DistortionCorrectionMatrixInversion - inputfile_cm: " << inputfile_cm.c_str() << std::endl;
0241   std::cout << "DistortionCorrectionMatrixInversion - outputFile: " << outputFile .c_str()<< std::endl;
0242
0243   auto filenames = list_files( inputFile );
0244   std::cout << "SpaceChargeMatrixInversion - loaded " << filenames.size() << " files" << std::endl;
0245
0246   // update TPOT phi range, needed for the interpolation
0247   load_tpot_geometry(52077);
0248   std::cout << "SpaceChargeMatrixInversion - done adjusting TPOT phi and theta ranges" << std::endl;
0249
0250   // perform matrix inversion
0251   TpcSpaceChargeMatrixInversion spaceChargeMatrixInversion;
0252
0253   // load input files
0254   for( const auto& file:filenames )
0255   { spaceChargeMatrixInversion.add_from_file( file ); }
0256
0257   // calculate the distortions in TPOT acceptance
0258   spaceChargeMatrixInversion.calculate_distortion_corrections();
0259
0260   // load central membrane corrections
0261   spaceChargeMatrixInversion.load_cm_distortion_corrections( inputfile_cm );
0262   spaceChargeMatrixInversion.extrapolate_distortion_corrections();
0263
0264   // write to output
0265   spaceChargeMatrixInversion.save_distortion_corrections( outputFile );
0266
0267   std::cout << "DistortionCorrectionMatrixInversion - all done." << std::endl;
0268 }