11 #include "GaudiKernel/MsgStream.h"
12 #include "GaudiKernel/ITHistSvc.h"
37 #include <TDecompLU.h>
39 #define MIN_TRACK_SEC 50 // min number of tracks per sector
44 template<
typename ...Ptr>
46 anyNullPtr(Ptr&&...
p){
47 return ((
p!=
nullptr) && ...) ;
104 std::string tree_title =
"Ambank " +
std::to_string(0) +
" matrices";
105 m_good_tree =
new TTree(tree_name.c_str(), tree_title.c_str());
111 return StatusCode::SUCCESS;
119 auto h_vc = std::make_unique<TH1F>(
"h_vc",
"h_vc",500,-1
e-4,1
e-4);
120 auto h_vd = std::make_unique<TH1F>(
"h_vd",
"h_vd",500,-1
e-2,1
e-2);
121 auto h_vf = std::make_unique<TH1F>(
"h_vf",
"h_vf",500,-1
e-2,1
e-2);
122 auto h_vz = std::make_unique<TH1F>(
"h_vz",
"h_vz",500,-1
e-2,1
e-2);
123 auto h_veta = std::make_unique<TH1F>(
"h_veta",
"h_veta",500,-1
e-2,1
e-2);
130 return StatusCode::SUCCESS;
137 TTree *old_tree = (TTree*)
file->Get(
"slice");
158 old_tree->GetEntry(0);
161 TTree *new_tree =
new TTree(
"slice",
"Slice boundaries");
186 return StatusCode::SUCCESS;
198 std::string tree_title =
"Ambank " +
std::to_string(0) +
" constants";
199 m_ctree =
new TTree(tree_name.c_str(), tree_title.c_str());
207 m_ctree->Branch(
"sectorID", &anInt,
"sectorID[nplane]/I");
208 m_ctree->Branch(
"hashID", &anInt,
"hashID[nplane]/I");
209 m_ctree->Branch(
"nhit", &aFloat,
"nhit/F");
211 m_ctree->Branch(
"Cd", &aDouble,
"Cd/D");
212 m_ctree->Branch(
"Cc", &aDouble,
"Cc/D");
213 m_ctree->Branch(
"Cf", &aDouble,
"Cf/D");
214 m_ctree->Branch(
"Cz0", &aDouble,
"Cz0/D");
215 m_ctree->Branch(
"Co", &aDouble,
"Co/D");
217 m_ctree->Branch(
"Vc", &aDouble,
"Vc[ndim]/D");
218 m_ctree->Branch(
"Vd", &aDouble,
"Vd[ndim]/D");
219 m_ctree->Branch(
"Vf", &aDouble,
"Vf[ndim]/D");
220 m_ctree->Branch(
"Vz0", &aDouble,
"Vz0[ndim]/D");
221 m_ctree->Branch(
"Vo", &aDouble,
"Vo[ndim]/D");
223 m_ctree->Branch(
"kernel", &aDouble,
"kernel[nkernel]/D");
224 m_ctree->Branch(
"kaverages", &aDouble,
"kaverages[nkaverages]/D");
227 return StatusCode::SUCCESS;
236 if (!
file.is_open())
return;
240 while (
file >> sector)
242 if (sector >= nSectors)
266 std::vector<module_t> & modules =
reader.getModules();
282 if (!success)
continue;
306 unsigned int nusable = acc_norm.
nusable() - 1;
308 coordsToUse[missing] =
false;
312 coordsToUse[missing+1] =
false;
346 std::vector<double> inv_covariance =
invert(
m_nCoords, mtx_reduced, coordsToUse);
349 std::vector<double> eigvals;
352 eigen(nusable,
m_nCoords, mtx_reduced, coordsToUse, eigvals, eigvecs);
355 geo =
makeConsts(acc_norm, coordsToUse, inv_covariance, eigvals, eigvecs);
368 float coverage =
static_cast<float>(
acc.track_bins.size());
369 m_ctree->SetBranchAddress(
"sectorID",
acc.FTK_modules.data());
370 m_ctree->SetBranchAddress(
"hashID", modules.data());
371 m_ctree->SetBranchAddress(
"nhit", &coverage);
373 m_ctree->SetBranchAddress(
"Cc", &
geo.pars.qOverPt);
374 m_ctree->SetBranchAddress(
"Cd", &
geo.pars.d0);
375 m_ctree->SetBranchAddress(
"Cf", &
geo.pars.phi);
376 m_ctree->SetBranchAddress(
"Cz0", &
geo.pars.z0);
377 m_ctree->SetBranchAddress(
"Co", &
geo.pars.eta);
379 m_ctree->SetBranchAddress(
"Vc",
geo.Vcurvature.data());
380 m_ctree->SetBranchAddress(
"Vd",
geo.Vd0.data());
381 m_ctree->SetBranchAddress(
"Vf",
geo.Vphi.data());
382 m_ctree->SetBranchAddress(
"Vz0",
geo.Vz0.data());
383 m_ctree->SetBranchAddress(
"Vo",
geo.Veta.data());
385 m_ctree->SetBranchAddress(
"kaverages",
geo.kaverages.data());
386 m_ctree->SetBranchAddress(
"kernel",
geo.kernel.data());
392 const std::string
prefix{
"/TRIGFPGATrackSimTREEGOODOUT/"};
393 auto getHistogram = [&](
const std::string &
suffix)->TH1*{
396 return (
sc == StatusCode::SUCCESS) ?
ptr:
nullptr;
398 auto h_vc = getHistogram(
"h_vc");
399 auto h_vd = getHistogram(
"h_vd");
400 auto h_vf = getHistogram(
"h_vf");
401 auto h_vz = getHistogram(
"h_vz");
402 auto h_veta = getHistogram(
"h_veta");
403 if (anyNullPtr(h_vc, h_vd, h_vf, h_vz, h_veta)){
404 ATH_MSG_ERROR(
"FPGATrackSimConstGenAlgo::fillConstTree; nullptr");
409 h_vc->Fill(
geo.Vcurvature[
i]);
410 h_vd->Fill(
geo.Vd0[
i]);
411 h_vf->Fill(
geo.Vphi[
i]);
412 h_vz->Fill(
geo.Vz0[
i]);
413 h_veta->Fill(
geo.Veta[
i]);
421 if (TMath::IsNaN(
value))
429 #define CHECK_NAN(var) (isNAN((var), #var))
486 double coverage =
static_cast<double>(
acc.track_bins.size());
487 size_t nCoords =
acc.hit_coords.size();
490 acc.pars[
i] /= coverage;
492 for (
unsigned i = 0;
i < nCoords;
i++)
494 acc.hit_coords[
i] /= coverage;
496 else acc.coords_usable[
i] =
true;
500 for (
unsigned i = 0;
i < nCoords;
i++)
502 acc.hit_x_QoP[
i] = (
acc.hit_x_QoP[
i] -
acc.hit_coords[
i] *
acc.pars.qOverPt * coverage) / (coverage-1);
503 acc.hit_x_d0[
i] = (
acc.hit_x_d0[
i] -
acc.hit_coords[
i] *
acc.pars.d0 * coverage) / (coverage-1);
504 acc.hit_x_z0[
i] = (
acc.hit_x_z0[
i] -
acc.hit_coords[
i] *
acc.pars.z0 * coverage) / (coverage-1);
505 acc.hit_x_phi[
i] = (
acc.hit_x_phi[
i] -
acc.hit_coords[
i] *
acc.pars.phi * coverage) / (coverage-1);
506 acc.hit_x_eta[
i] = (
acc.hit_x_eta[
i] -
acc.hit_coords[
i] *
acc.pars.eta * coverage) / (coverage-1);
509 for (
unsigned j =
i ; j < nCoords; ++j)
511 acc.covariance[j * nCoords +
i] =
acc.covariance[
i * nCoords + j] =
512 (
acc.covariance[
i * nCoords + j] -
acc.hit_coords[
i] *
acc.hit_coords[j] * coverage) / (coverage-1);
521 std::vector<double>
const & inv_covariance,
524 size_t nCoords =
acc.hit_coords.size();
535 if (!usable[
i])
continue;
536 for (
size_t j = 0; j < nCoords; j++)
538 if (!usable[j])
continue;
540 geo.kernel(
i, j) = eigvecs(
i, j) / sqrt(eigvals[
i]);
542 geo.kaverages[
i] += -
geo.kernel(
i, j) *
acc.hit_coords[j];
547 for (
size_t i = 0;
i < nCoords;
i++)
550 geo.pars.qOverPt += -
geo.Vcurvature[
i] *
acc.hit_coords[
i];
556 geo.real =
static_cast<int>(
acc.track_bins.size());
571 std::vector<double>
x(
n);
573 for (
size_t i = 0;
i <
n;
i++)
574 for (
size_t j = 0; j <
n; j++)
575 x[
i] +=
A[
i *
n + j] *
b[j];
599 TMatrixD mtx(
n,
n, mtx_v.data());
600 TMatrixD newmtx(nDimToUse, nDimToUse);
603 for (
size_t i = 0;
i <
n;
i++)
605 if (!usable[
i])
continue;
608 for (
size_t j = 0; j <
n; j++)
610 if (!usable[j])
continue;
611 newmtx[counteri][counterj] = mtx[
i][j];
614 assert(counterj == nDimToUse);
617 assert(counteri == nDimToUse);
633 std::vector<double> inv(n_full * n_full);
638 for (
size_t i = 0;
i < n_full;
i++)
640 if (!usable[
i])
continue;
643 for (
size_t j = 0; j < n_full; j++)
645 if (!usable[j])
continue;
646 inv[
i*n_full + j] = mtx[counteri][counterj];
673 TVectorD eigvals_redu;
674 TMatrixD eigvecs_redu = mtx.EigenVectors(eigvals_redu);
677 eigvals_full.resize(n_full, 0);
678 eigvecs_full.
resize(n_full, n_full, 0);
682 size_t j_redu = n_redu - 1;
683 for (
size_t i_full = 0; i_full < n_full; i_full++)
685 if (!usable[i_full])
continue;
688 for (
size_t j_full = 0; j_full < n_full; j_full++)
690 if (!usable[j_full])
continue;
691 eigvecs_full(i_full, j_full) = eigvecs_redu[i_redu][j_redu];
694 eigvals_full[i_full] = eigvals_redu[j_redu];
703 for (
size_t i = 0;
i <
size;
i++)
704 total += vec1[
i] *
vec2[
i];
713 if (!usable[
i])
continue;
716 auto project = [&](std::vector<double> & hit_x_par,
double &
par)
718 double pr =
dot(hit_x_par.data(),
geo.kernel[
i], nCoords) /
norm;
719 for (
int j = 0; j < nCoords; j++) hit_x_par[j] += -
geo.kernel(
i,j) *
pr;
738 return StatusCode::SUCCESS;
756 for (
int missing = 0; missing <
m_nLayers; missing++) {
768 return StatusCode::SUCCESS;
782 FILE *sector_file = fopen(sector_filename.c_str(),
"w");
783 FILE *sectorHW_file = fopen(sectorHW_filename.c_str(),
"w");
790 while (
reader.nextEntry())
793 size_t sector =
reader.getEntry();
795 fprintf(sector_file,
"%zu ", sector);
796 fprintf(sectorHW_file,
"%zu ", sector);
799 fprintf(sector_file,
"%d ",
acc.FTK_modules[
i]);
800 fprintf(sectorHW_file,
"%d ",
reader.getModules()[
i]);
802 fprintf(sector_file,
"0 %zu",
acc.track_bins.size());
803 fprintf(sector_file,
"\n");
804 fprintf(sectorHW_file,
"0 %zu",
acc.track_bins.size());
805 fprintf(sectorHW_file,
"\n");
812 fclose(sectorHW_file);
814 slice.saveSlices(slice_filename);
821 FILE *const_file = fopen(
filename.c_str(),
"w");
823 fprintf(const_file,
"! *** RECONSTRUCTION GEOMETRY CONSTANTS ***\n");
824 fprintf(const_file,
"\n");
825 fprintf(const_file,
"Version 2 ! File format version number\n");
826 fprintf(const_file,
"\n");
827 fprintf(const_file,
"! *** PHI is now in GLOBAL coordinate system ***\n");
828 fprintf(const_file,
" NPLANES\n");
830 fprintf(const_file,
" NSECTORS\n");
831 fprintf(const_file,
"%zu\n",geo_consts.size());
832 fprintf(const_file,
" NDIM\n");
833 fprintf(const_file,
" 2\n");
837 for (
size_t sector = 0; sector < geo_consts.size(); sector++)
840 fprintf(const_file,
"sector\n");
841 fprintf(const_file,
"%zu\n", sector);
843 fprintf(const_file,
" Vc \n");
845 fprintf(const_file,
"%e\n",geo_consts[sector].Vcurvature[
i]);
848 fprintf(const_file,
" Vd \n");
850 fprintf(const_file,
"%e\n",geo_consts[sector].Vd0[
i]);
853 fprintf(const_file,
" Vf \n");
855 fprintf(const_file,
"%e\n",geo_consts[sector].Vphi[
i]);
858 fprintf(const_file,
" Vz0 \n");
860 fprintf(const_file,
"%e\n",geo_consts[sector].Vz0[
i]);
863 fprintf(const_file,
" Vo \n");
865 fprintf(const_file,
"%e\n",geo_consts[sector].Veta[
i]);
868 fprintf(const_file,
"kaverages\n");
870 fprintf(const_file,
"%e\n",
m_geo_consts[sector].kaverages[
i]);
873 fprintf(const_file,
"kernel\n");
876 fprintf(const_file,
"%e\n",geo_consts[sector].kernel(
i, j));
880 fprintf(const_file,
"Cc\n");
881 fprintf(const_file,
"%e\n",geo_consts[sector].
pars.qOverPt);
883 fprintf(const_file,
"Cd\n");
884 fprintf(const_file,
"%e\n",geo_consts[sector].
pars.d0);
886 fprintf(const_file,
"Cf\n");
887 fprintf(const_file,
"%e\n",geo_consts[sector].
pars.phi);
889 fprintf(const_file,
"Cz0\n");
890 fprintf(const_file,
"%e\n",geo_consts[sector].
pars.z0);
892 fprintf(const_file,
"Co\n");
893 fprintf(const_file,
"%e\n",geo_consts[sector].
pars.eta);