17 #include <Eigen/Dense>
24 static constexpr
float m_phimin = 0;
25 static constexpr
float m_phimax = 2.*M_PI;
29 static constexpr
float m_rmin = 20;
30 static constexpr
float m_rmax = 78;
33 static constexpr
float m_zmin = -105.5;
34 static constexpr
float m_zmax = 105.5;
57 std::cout <<
"TpcSpaceChargeMatrixInversion::add_from_file - could not open file " <<
filename << std::endl;
62 std::unique_ptr<TpcSpaceChargeMatrixContainer>
source( dynamic_cast<TpcSpaceChargeMatrixContainer*>( inputfile->Get( objectname.c_str() ) ) );
65 std::cout <<
"TpcSpaceChargeMatrixInversion::add_from_file - could not find object name " << objectname <<
" in file " <<
filename << std::endl;
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 ) );
112 for(
const auto&
h:{ hentries, hphi,
hz, hr } )
114 h->GetXaxis()->SetTitle(
"#phi (rad)" );
115 h->GetYaxis()->SetTitle(
"r (cm)" );
116 h->GetZaxis()->SetTitle(
"z (cm)" );
121 static constexpr
int ncoord = 3;
122 using matrix_t = Eigen::Matrix<float, ncoord, ncoord >;
123 using column_t = Eigen::Matrix<float, ncoord, 1 >;
126 for(
int iphi = 0; iphi <
phibins; ++iphi )
127 for(
int ir = 0; ir < rbins; ++ir )
128 for(
int iz = 0; iz < zbins; ++iz )
135 static constexpr
int min_cluster_count = 2;
137 if( cell_entries < min_cluster_count )
continue;
141 for(
int i = 0;
i < ncoord; ++
i )
142 for(
int j = 0;
j < ncoord; ++
j )
146 for(
int i = 0;
i < ncoord; ++
i )
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;
159 const auto cov = lhs.inverse();
160 auto partialLu = lhs.partialPivLu();
161 const auto result = partialLu.solve( rhs );
164 hentries->SetBinContent( iphi+1, ir+1, iz+1, cell_entries );
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) ) );
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) ) );
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) ) );
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;
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" ) );
192 for(
const auto&
h: {hentries, hphi, hr,
hz} )
202 for(
const auto&
h: { hentries, hphi, hr,
hz } ) {
h->Write(); }
206 auto process_histogram = []( TH3*
h,
const TString&
name )
218 process_histogram( hentries,
"hentries" );
219 process_histogram( hphi,
"hIntDistortionP" );
220 process_histogram( hr,
"hIntDistortionR" );
221 process_histogram(
hz,
"hIntDistortionZ" );