tpc/jobB/DistortionCorrectionMatrixInversion.C

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