ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcSpaceChargeMatrixInversion.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcSpaceChargeMatrixInversion.cc
1 
10 
11 #include <frog/FROG.h>
12 
13 #include <TFile.h>
14 #include <TH3.h>
15 
16 #include <Eigen/Core>
17 #include <Eigen/Dense>
18 
19 #include <memory>
20 
21 namespace
22 {
23  // phi range
24  static constexpr float m_phimin = 0;
25  static constexpr float m_phimax = 2.*M_PI;
26 
27  // TODO: could try to get the r and z range from TPC geometry
28  // r range
29  static constexpr float m_rmin = 20;
30  static constexpr float m_rmax = 78;
31 
32  // z range
33  static constexpr float m_zmin = -105.5;
34  static constexpr float m_zmax = 105.5;
35 }
36 
37 //_____________________________________________________________________
39  Fun4AllBase( name)
40 {}
41 
42 //_____________________________________________________________________
45 
46 //_____________________________________________________________________
47 bool TpcSpaceChargeMatrixInversion::add_from_file( const std::string& shortfilename, const std::string& objectname )
48 {
49  // get filename from frog
50  FROG frog;
51  const auto filename = frog.location( shortfilename );
52 
53  // open TFile
54  std::unique_ptr<TFile> inputfile( TFile::Open( filename ) );
55  if( !inputfile )
56  {
57  std::cout << "TpcSpaceChargeMatrixInversion::add_from_file - could not open file " << filename << std::endl;
58  return false;
59  }
60 
61  // load object from input file
62  std::unique_ptr<TpcSpaceChargeMatrixContainer> source( dynamic_cast<TpcSpaceChargeMatrixContainer*>( inputfile->Get( objectname.c_str() ) ) );
63  if( !source )
64  {
65  std::cout << "TpcSpaceChargeMatrixInversion::add_from_file - could not find object name " << objectname << " in file " << filename << std::endl;
66  return false;
67  }
68 
69  // add object
70  return add( *source.get() );
71 }
72 
73 //_____________________________________________________________________
75 {
76  // check internal container, create if necessary
77  if( !m_matrix_container )
78  {
80 
81  // get grid dimensions from source
82  int phibins = 0;
83  int rbins = 0;
84  int zbins = 0;
85  source.get_grid_dimensions( phibins, rbins, zbins );
86 
87  // assign
88  m_matrix_container->set_grid_dimensions( phibins, rbins, zbins );
89  }
90 
91  // add content
92  return m_matrix_container->add( source );
93 }
94 
95 //_____________________________________________________________________
97 {
98 
99  // get grid dimensions from matrix container
100  int phibins = 0;
101  int rbins = 0;
102  int zbins = 0;
103  m_matrix_container->get_grid_dimensions( phibins, rbins, zbins );
104 
105  // create output histograms
106  auto hentries( new TH3F( "hentries_rec", "hentries_rec", phibins, m_phimin, m_phimax, rbins, m_rmin, m_rmax, zbins, m_zmin, m_zmax ) );
107  auto hphi( new TH3F( "hDistortionP_rec", "hDistortionP_rec", phibins, m_phimin, m_phimax, rbins, m_rmin, m_rmax, zbins, m_zmin, m_zmax ) );
108  auto hz( new TH3F( "hDistortionZ_rec", "hDistortionZ_rec", phibins, m_phimin, m_phimax, rbins, m_rmin, m_rmax, zbins, m_zmin, m_zmax ) );
109  auto hr( new TH3F( "hDistortionR_rec", "hDistortionR_rec", phibins, m_phimin, m_phimax, rbins, m_rmin, m_rmax, zbins, m_zmin, m_zmax ) );
110 
111  // set axis labels
112  for( const auto& h:{ hentries, hphi, hz, hr } )
113  {
114  h->GetXaxis()->SetTitle( "#phi (rad)" );
115  h->GetYaxis()->SetTitle( "r (cm)" );
116  h->GetZaxis()->SetTitle( "z (cm)" );
117  }
118 
119  // matrix convenience definition
120  /* number of coordinates must match that of the matrix container */
121  static constexpr int ncoord = 3;
122  using matrix_t = Eigen::Matrix<float, ncoord, ncoord >;
123  using column_t = Eigen::Matrix<float, ncoord, 1 >;
124 
125  // loop over bins
126  for( int iphi = 0; iphi < phibins; ++iphi )
127  for( int ir = 0; ir < rbins; ++ir )
128  for( int iz = 0; iz < zbins; ++iz )
129  {
130 
131  // get cell index
132  const auto icell = m_matrix_container->get_cell_index( iphi, ir, iz );
133 
134  // minimum number of entries per bin
135  static constexpr int min_cluster_count = 2;
136  const auto cell_entries = m_matrix_container->get_entries(icell);
137  if( cell_entries < min_cluster_count ) continue;
138 
139  // build eigen matrices from container
140  matrix_t lhs;
141  for( int i = 0; i < ncoord; ++i )
142  for( int j = 0; j < ncoord; ++j )
143  { lhs(i,j) = m_matrix_container->get_lhs( icell, i, j ); }
144 
145  column_t rhs;
146  for( int i = 0; i < ncoord; ++i )
147  { rhs(i) = m_matrix_container->get_rhs( icell, i ); }
148 
149  if (Verbosity())
150  {
151  // print matrices and entries
152  std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortions - inverting bin " << iz << ", " << ir << ", " << iphi << std::endl;
153  std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortions - entries: " << cell_entries << std::endl;
154  std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortions - lhs: \n" << lhs << std::endl;
155  std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortions - rhs: \n" << rhs << std::endl;
156  }
157 
158  // calculate result using linear solving
159  const auto cov = lhs.inverse();
160  auto partialLu = lhs.partialPivLu();
161  const auto result = partialLu.solve( rhs );
162 
163  // fill histograms
164  hentries->SetBinContent( iphi+1, ir+1, iz+1, cell_entries );
165 
166  hphi->SetBinContent( iphi+1, ir+1, iz+1, result(0) );
167  hphi->SetBinError( iphi+1, ir+1, iz+1, std::sqrt( cov(0,0) ) );
168 
169  hz->SetBinContent( iphi+1, ir+1, iz+1, result(1) );
170  hz->SetBinError( iphi+1, ir+1, iz+1, std::sqrt( cov(1,1) ) );
171 
172  hr->SetBinContent( iphi+1, ir+1, iz+1, result(2) );
173  hr->SetBinError( iphi+1, ir+1, iz+1, std::sqrt( cov(2,2) ) );
174 
175  if (Verbosity())
176  {
177  std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortions - drphi: " << result(0) << " +/- " << std::sqrt( cov(0,0) ) << std::endl;
178  std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortions - dz: " << result(1) << " +/- " << std::sqrt( cov(1,1) ) << std::endl;
179  std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortions - dr: " << result(2) << " +/- " << std::sqrt( cov(2,2) ) << std::endl;
180  std::cout << std::endl;
181  }
182  }
183 
184  // save everything to root file
185  std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortions - writing histograms to " << m_outputfile << std::endl;
186  std::unique_ptr<TFile> outputfile( TFile::Open( m_outputfile.c_str(), "RECREATE" ) );
187  outputfile->cd();
188 
189  // when using migromegas, one needs to extrapolate to the rest of the acceptance
190  if( m_use_micromegas )
191  {
192  for( const auto& h: {hentries, hphi, hr, hz} )
193  {
194  if( !h ) continue;
198  }
199  }
200 
201  // write source histograms
202  for( const auto& h: { hentries, hphi, hr, hz } ) { h->Write(); }
203 
204  // split histograms in two along z axis and write
205  // also write histograms suitable for space charge reconstruction
206  auto process_histogram = []( TH3* h, const TString& name )
207  {
208  TH3* hneg;
209  TH3* hpos;
210  std::tie( hneg, hpos ) = TpcSpaceChargeReconstructionHelper::split( h );
211  hneg->Write();
212  hpos->Write();
214  TpcSpaceChargeReconstructionHelper::copy_histogram( hneg, Form( "%s_negz", name.Data() ) )->Write();
215  TpcSpaceChargeReconstructionHelper::copy_histogram( hpos, Form( "%s_posz", name.Data() ) )->Write();
216  };
217 
218  process_histogram( hentries, "hentries" );
219  process_histogram( hphi, "hIntDistortionP" );
220  process_histogram( hr, "hIntDistortionR" );
221  process_histogram( hz, "hIntDistortionZ" );
222 
223  // close output file
224  outputfile->Close();
225 
226 }