packages/TrackerMillepedeAlignment/MakeMilleFiles.cc

0001 #include "MakeMilleFiles.h"
0002
0003 #include "Mille.h"
0004
0005 /// Tracking includes
0006
0007 #include <trackbase/MvtxDefs.h>
0008 #include <trackbase/TpcDefs.h>  // for side
0009 #include <trackbase/TrackFitUtils.h>
0010 #include <trackbase/TrkrCluster.h>
0011 #include <trackbase/TrkrClusterContainer.h>
0012 #include <trackbase/TrkrDefs.h>
0013
0014 #include <trackbase_historic/ActsTransformations.h>
0015 #include <trackbase_historic/SvtxAlignmentState.h>
0016 #include <trackbase_historic/SvtxTrack.h>
0017 #include <trackbase_historic/SvtxTrackMap.h>
0018 #include <trackbase_historic/SvtxTrackState.h>
0019 #include <trackbase_historic/TrackAnalysisUtils.h>
0020
0021 #include <trackreco/ActsPropagator.h>
0022
0023 #include <g4detectors/PHG4TpcCylinderGeom.h>
0024 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0025
0026 #include <fun4all/Fun4AllReturnCodes.h>
0027
0028 #include <Acts/Definitions/Algebra.hpp>
0029
0030 #include <phool/PHCompositeNode.h>
0031 #include <phool/getClass.h>
0032 #include <phool/phool.h>
0033
0034 #include <TFile.h>
0035 #include <TNtuple.h>
0036 #include <TF1.h>
0037
0038 #include <cmath>
0039 #include <fstream>
0040 #include <iostream>
0041 #include <map>
0042 #include <set>
0043 #include <utility>
0044
0045 namespace
0046 {
0047   /// square
0048   template <class T>
0049   inline constexpr T square(const T& x)
0050   {
0051     return x * x;
0052   }
0053 }  // namespace
0054
0055 //____________________________________________________________________________..
0056 MakeMilleFiles::MakeMilleFiles(const std::string& name)
0057   : SubsysReco(name)
0058   , _mille(nullptr)
0059 {
0060 }
0061
0062 //____________________________________________________________________________..
0063 int MakeMilleFiles::InitRun(PHCompositeNode* topNode)
0064 {
0065   int ret = GetNodes(topNode);
0066   if (ret != Fun4AllReturnCodes::EVENT_OK)
0067   {
0068     return ret;
0069   }
0070
0071   // Instantiate Mille and open output data file
0072   //  _mille = new Mille(data_outfilename.c_str(), false);   // write text in data files, rather than binary, for debugging only
0073   _mille = new Mille(data_outfilename.c_str(), _binary);
0074
0075   // Write the steering file here, and add the data file path to it
0076   std::ofstream steering_file(steering_outfilename);
0077   steering_file << data_outfilename << std::endl;
0078   steering_file << m_constraintFileName << std::endl;
0079   steering_file.close();
0080
0081   m_constraintFile.open(m_constraintFileName);
0082   if (m_useEventVertex)
0083   {
0084     for (int i : AlignmentDefs::glbl_vtx_label)
0085     {
0086       m_constraintFile << " Constraint   0.0" << std::endl;
0087       m_constraintFile << "       " << i << "   1" << std::endl;
0088     }
0089   }
0090   // print grouping setup to log file:
0091   std::cout << "MakeMilleFiles::InitRun: Surface groupings are mvtx " << mvtx_group << " intt " << intt_group << " tpc " << tpc_group << " mms " << mms_group << std::endl;
0092
0093   // Make TTree for cross check
0094   if (!m_tfile_name.empty())
0095   {
0096     m_file = TFile::Open(m_tfile_name.c_str(), "RECREATE");
0097     m_ntuple = new TNtuple (
0098       "ntp", "ntp",
0099       "dXdR:dXdX0:dXdY0:dXdZs:dXdZ0:"
0100       "dXdalpha:dXdbeta:dXdgamma:dXdx:dXdy:dXdz:"
0101       "dYdR:dYdX0:dYdY0:dYdZs:dYdZ0:"
0102       "dYdalpha:dYdbeta:dYdgamma:dYdx:dYdy:dYdz"
0103     );
0104     m_ntuple->SetDirectory(m_file);
0105   }
0106
0107   return ret;
0108 }
0109
0110 //____________________________________________________________________________..
0111 int MakeMilleFiles::process_event(PHCompositeNode* /*topNode*/)
0112 {
0113   // Outline:
0114   //
0115   // loop over track alignment states
0116   //   Make any track cuts here to skip undesirable tracks (maybe low pT?)
0117   //   loop over track states+measurements for each track
0118   //      for each measurement, performed in trackreco/ActsAlignmentStates.cc
0119   //         Get measurement value and error
0120   //         Calculate derivatives and residuals from Acts jacobians
0121   //         These are stored in a map and unpacked for mille
0122   //   Call _mille->mille() with arguments obtained from previous iteration:
0123   //     local pars
0124   //     array of local derivatives
0125   //     global pars
0126   //     array of global derivatives
0127   //     array of integer global par labels
0128   //     residual value (float) z = measurement - track state
0129   //     sigma of measurement
0130   //   After processing all measurements for this track, call _mille->end() to add buffer to file and reset buffer
0131   // After all tracks are processed, file is closed when Mille destructor is called
0132   // Note: all units are in the Acts units of mm and GeV to avoid converting matrices
0133
0134   if (Verbosity() > 0)
0135   {
0136     std::cout << PHWHERE << " track map size " << _track_map->size() << std::endl;
0137     std::cout << "state map size " << _state_map->size() << std::endl;
0138   }
0139
0140   Acts::Vector3 eventVertex = Acts::Vector3::Zero();
0141   if (m_useEventVertex)
0142   {
0143     eventVertex = getEventVertex();
0144   }
0145
0146   ActsPropagator propagator(_tGeometry);
0147
0148   for (auto [key, statevec] : *_state_map)
0149   {
0150     // Check if track was removed from cleaner
0151     auto iter = _track_map->find(key);
0152     if (iter == _track_map->end())
0153     {
0154       continue;
0155     }
0156
0157     SvtxTrack* track = iter->second;
0158
0159     if (Verbosity() > 0)
0160     {
0161       std::cout << std::endl
0162                 << __LINE__ << ": Processing track itrack: " << key << ": nhits: " << track->size_cluster_keys()
0163                 << ": Total tracks: " << _track_map->size() << ": phi: " << track->get_phi() << std::endl;
0164     }
0165
0166     //! Make any desired track cuts here
0167     //! Maybe set a lower pT limit - low pT tracks are not very sensitive to alignment
0168     addTrackToMilleFile(statevec);
0169
0170     //! Only take tracks that have 2 mm within event vertex
0171     if (m_useEventVertex &&
0172         fabs(track->get_z() - eventVertex.z()) < 0.2 &&
0173         fabs(track->get_x()) < 0.2 &&
0174         fabs(track->get_y()) < 0.2)
0175     {
0176       //! set x and y to 0 since we are constraining to the x-y origin
0177       //! and add constraints to pede later
0178       eventVertex(0) = 0;
0179       eventVertex(1) = 0;
0180
0181       auto dcapair = TrackAnalysisUtils::get_dca(track, eventVertex);
0182       Acts::Vector2 vtx_residual(-dcapair.first.first, -dcapair.second.first);
0183       vtx_residual *= Acts::UnitConstants::cm;
0184
0185       float lclvtx_derivative[SvtxAlignmentState::NRES][SvtxAlignmentState::NLOC];
0186       bool success = getLocalVtxDerivativesXY(track, propagator,
0187                                               eventVertex, lclvtx_derivative);
0188
0189       // The global derivs dimensions are [alpha/beta/gamma](x/y/z)
0190       float glblvtx_derivative[SvtxAlignmentState::NRES][3];
0191       getGlobalVtxDerivativesXY(track, eventVertex, glblvtx_derivative);
0192
0193       if (Verbosity() > 2)
0194       {
0195         std::cout << "vertex info for trakc " << track->get_id() << " with charge " << track->get_charge() << std::endl;
0196         std::cout << "vertex is " << eventVertex.transpose() << std::endl;
0197         std::cout << "vertex residuals " << vtx_residual.transpose()
0198                   << std::endl;
0199         std::cout << "global vtx derivatives " << std::endl;
0200         for (auto& i : glblvtx_derivative)
0201         {
0202           for (float j : i)
0203           {
0204             std::cout << j << ", ";
0205           }
0206           std::cout << std::endl;
0207         }
0208       }
0209       if (success)
0210       {
0211         for (int i = 0; i < 2; i++)
0212         {
0213           if (!std::isnan(vtx_residual(i)))
0214           {
0215             _mille->mille(SvtxAlignmentState::NLOC, lclvtx_derivative[i],
0216                           AlignmentDefs::NGLVTX, glblvtx_derivative[i],
0217                           AlignmentDefs::glbl_vtx_label, vtx_residual(i),
0218                           m_vtxSigma(i));
0219
0220           }
0221         }
0222       }
0223     }
0224
0225     //! Finish this track
0226     _mille->end();
0227   }
0228
0229   if (Verbosity() > 0)
0230   {
0231     std::cout << "Finished processing mille file " << std::endl;
0232   }
0233
0234   return Fun4AllReturnCodes::EVENT_OK;
0235 }
0236
0237 int MakeMilleFiles::End(PHCompositeNode* /*unused*/)
0238 {
0239   delete _mille;
0240   m_constraintFile.close();
0241
0242   if (m_file && m_ntuple)
0243   {
0244     m_ntuple->Write();
0245     m_file->Write();
0246     m_file->Close();
0247   }
0248
0249   return Fun4AllReturnCodes::EVENT_OK;
0250 }
0251
0252 int MakeMilleFiles::GetNodes(PHCompositeNode* topNode)
0253 {
0254   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0255   if (!_cluster_map)
0256   {
0257     std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0258     return Fun4AllReturnCodes::ABORTEVENT;
0259   }
0260
0261   _track_map = findNode::getClass<SvtxTrackMap>(topNode, m_track_map_name);
0262   if (!_track_map)
0263   {
0264     std::cout
0265         << PHWHERE << "\n"
0266         << "\tCan't find track map:\n"
0267         << "\t" << m_track_map_name << "\n"
0268         << "\tAborting\n"
0269         << std::endl;
0270     return Fun4AllReturnCodes::ABORTEVENT;
0271   }
0272
0273   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0274   if (!_tGeometry)
0275   {
0276     std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
0277     return Fun4AllReturnCodes::ABORTEVENT;
0278   }
0279
0280   _state_map = findNode::getClass<SvtxAlignmentStateMap>(topNode, m_state_map_name);
0281   if (!_state_map)
0282   {
0283     std::cout
0284         << PHWHERE << "\n"
0285         << "\tCan't find state map:\n"
0286         << "\t" << m_state_map_name << "\n"
0287         << "\tAborting\n"
0288         << std::endl;
0289     return Fun4AllReturnCodes::ABORTEVENT;
0290   }
0291
0292   return Fun4AllReturnCodes::EVENT_OK;
0293 }
0294
0295 bool MakeMilleFiles::getLocalVtxDerivativesXY(SvtxTrack* track,
0296                                               ActsPropagator& propagator,
0297                                               const Acts::Vector3& vertex,
0298                                               float lclvtx_derivative[SvtxAlignmentState::NRES][SvtxAlignmentState::NLOC])
0299 {
0300   //! Get the first track state beyond the vertex, which will be the
0301   //! innermost track state and propagate it to the vertex surface to
0302   //! get the jacobian at the vertex
0303   SvtxTrackState* firststate = (*std::next(track->begin_states(), 1)).second;
0304
0305   TrkrDefs::cluskey ckey = firststate->get_cluskey();
0306   auto cluster = _cluster_map->findCluster(ckey);
0307   auto surf = _tGeometry->maps().getSurface(ckey, cluster);
0308
0309   auto param = propagator.makeTrackParams(firststate, track->get_charge(), surf).value();
0310   auto perigee = propagator.makeVertexSurface(vertex);
0311   auto actspropagator = propagator.makePropagator();
0312
0313   Acts::PropagatorOptions<> options(_tGeometry->geometry().getGeoContext(),
0314                                     _tGeometry->geometry().magFieldContext);
0315
0316   auto result = actspropagator.propagate(param, *perigee, options);
0317
0318   if (result.ok())
0319   {
0320     auto jacobian = *result.value().transportJacobian;
0321
0322     Eigen::Matrix<double, 2, 6> projector = Eigen::Matrix<double, 2, 6>::Zero();
0323     projector(0, 0) = 1;
0324     projector(1, 1) = 1;
0325     auto deriv = projector * jacobian;
0326     if (Verbosity() > 2)
0327     {
0328       std::cout << "local vtxderiv " << std::endl
0329                 << deriv << std::endl;
0330     }
0331
0332     for (int i = 0; i < deriv.rows(); i++)
0333     {
0334       for (int j = 0; j < deriv.cols(); j++)
0335       {
0336         lclvtx_derivative[i][j] = deriv(i, j);
0337       }
0338     }
0339
0340     return true;
0341   }
0342
0343   return false;
0344 }
0345
0346 void MakeMilleFiles::getGlobalVtxDerivativesXY(SvtxTrack* track,
0347                                                const Acts::Vector3& vertex,
0348                                                float glblvtx_derivative[SvtxAlignmentState::NRES][3])
0349 {
0350   Acts::SquareMatrix3 identity = Acts::SquareMatrix3::Identity();
0351
0352   Acts::Vector3 projx = Acts::Vector3::Zero();
0353   Acts::Vector3 projy = Acts::Vector3::Zero();
0354   getProjectionVtxXY(track, vertex, projx, projy);
0355
0356   glblvtx_derivative[0][0] = identity.col(0).dot(projx);
0357   glblvtx_derivative[0][1] = identity.col(1).dot(projx);
0358   glblvtx_derivative[0][2] = identity.col(2).dot(projx);
0359   glblvtx_derivative[1][0] = identity.col(0).dot(projy);
0360   glblvtx_derivative[1][1] = identity.col(1).dot(projy);
0361   glblvtx_derivative[1][2] = identity.col(2).dot(projy);
0362
0363   /// swap sign to fit expectation of pede to have derivative of fit rather than
0364   /// residual
0365   for (int i = 0; i < 2; i++)
0366   {
0367     for (int j = 0; j < 3; j++)
0368     {
0369       glblvtx_derivative[i][j] *= -1;
0370     }
0371   }
0372 }
0373 void MakeMilleFiles::getProjectionVtxXY(SvtxTrack* track,
0374                                         const Acts::Vector3& vertex,
0375                                         Acts::Vector3& projx,
0376                                         Acts::Vector3& projy)
0377 {
0378   Acts::Vector3 tangent(track->get_px(), track->get_py(), track->get_pz());
0379   Acts::Vector3 normal(track->get_px(), track->get_py(), 0);
0380
0381   tangent /= tangent.norm();
0382   normal /= normal.norm();
0383
0384   Acts::Vector3 localx(1, 0, 0);
0385   Acts::Vector3 localz(0, 0, 1);
0386
0387   Acts::Vector3 xglob = localToGlobalVertex(track, vertex, localx);
0388   Acts::Vector3 yglob = localz + vertex;
0389
0390   Acts::Vector3 X = (xglob - vertex) / (xglob - vertex).norm();
0391   Acts::Vector3 Y = (yglob - vertex) / (yglob - vertex).norm();
0392
0393   // see equation 31 of the ATLAS paper (and discussion) for this
0394   projx = X - (tangent.dot(X) / tangent.dot(normal)) * normal;
0395   projy = Y - (tangent.dot(Y) / tangent.dot(normal)) * normal;
0396
0397   return;
0398 }
0399 Acts::Vector3 MakeMilleFiles::localToGlobalVertex(SvtxTrack* track,
0400                                                   const Acts::Vector3& vertex,
0401                                                   const Acts::Vector3& localx) const
0402 {
0403   Acts::Vector3 mom(track->get_px(),
0404                     track->get_py(),
0405                     track->get_pz());
0406
0407   Acts::Vector3 r = mom.cross(Acts::Vector3(0., 0., 1.));
0408   float phi = atan2(r(1), r(0));
0409   Acts::RotationMatrix3 rot;
0410   Acts::RotationMatrix3 rot_T;
0411
0412   rot(0, 0) = cos(phi);
0413   rot(0, 1) = -sin(phi);
0414   rot(0, 2) = 0;
0415   rot(1, 0) = sin(phi);
0416   rot(1, 1) = cos(phi);
0417   rot(1, 2) = 0;
0418   rot(2, 0) = 0;
0419   rot(2, 1) = 0;
0420   rot(2, 2) = 1;
0421
0422   rot_T = rot.transpose();
0423
0424   Acts::Vector3 pos_R = rot * localx;
0425
0426   pos_R += vertex;
0427
0428   return pos_R;
0429 }
0430
0431 Acts::Vector3 MakeMilleFiles::getEventVertex()
0432 {
0433   /**
0434    * Returns event vertex in cm as averaged track positions
0435    */
0436   float xsum = 0;
0437   float ysum = 0;
0438   float zsum = 0;
0439   int nacceptedtracks = 0;
0440
0441   for (auto [key, statevec] : *_state_map)
0442   {
0443     // Check if track was removed from cleaner
0444     auto iter = _track_map->find(key);
0445     if (iter == _track_map->end())
0446     {
0447       continue;
0448     }
0449
0450     SvtxTrack* track = iter->second;
0451
0452     /// The track vertex is given by the fit as the PCA to the beamline
0453     xsum += track->get_x();
0454     ysum += track->get_y();
0455     zsum += track->get_z();
0456
0457     nacceptedtracks++;
0458   }
0459
0460   return Acts::Vector3(xsum / nacceptedtracks,
0461                        ysum / nacceptedtracks,
0462                        zsum / nacceptedtracks);
0463 }
0464
0465 void MakeMilleFiles::addTrackToMilleFile(SvtxAlignmentStateMap::StateVec& statevec)
0466 {
0467   for (auto state : statevec)
0468   {
0469     TrkrDefs::cluskey ckey = state->get_cluster_key();
0470
0471     if (Verbosity() > 2)
0472     {
0473       std::cout << "adding state for ckey " << ckey << " with hitsetkey "
0474                 << (int) TrkrDefs::getHitSetKeyFromClusKey(ckey) << std::endl;
0475     }
0476     // The global alignment parameters are given initial values of zero by default, we do not specify them
0477     // We identify the global alignment parameters for this surface
0478
0479     TrkrCluster* cluster = _cluster_map->findCluster(ckey);
0480     const unsigned int layer = TrkrDefs::getLayer(ckey);
0481     const unsigned int trkrid = TrkrDefs::getTrkrId(ckey);
0482     const SvtxAlignmentState::ResidualVector residual = state->get_residual();
0483     const Acts::Vector3 global = _tGeometry->getGlobalPosition(ckey, cluster);
0484
0485     // need standard deviation of measurements
0486     SvtxAlignmentState::ResidualVector clus_sigma = SvtxAlignmentState::ResidualVector::Zero();
0487
0488     double clusRadius = sqrt(global[0] * global[0] + global[1] * global[1]);
0489     auto para_errors = _ClusErrPara.get_clusterv5_modified_error(cluster, clusRadius, ckey);
0490     double phierror = sqrt(para_errors.first);
0491     double zerror = sqrt(para_errors.second);
0492     clus_sigma(1) = zerror * Acts::UnitConstants::cm;
0493     clus_sigma(0) = phierror * Acts::UnitConstants::cm;
0494
0495     if (std::isnan(clus_sigma(0)) ||
0496         std::isnan(clus_sigma(1)))
0497     {
0498       continue;
0499     }
0500
0501     auto surf = _tGeometry->maps().getSurface(ckey, cluster);
0502
0503     int glbl_label[SvtxAlignmentState::NGL];
0504     if (layer < 3)
0505     {
0506       AlignmentDefs::getMvtxGlobalLabels(surf, glbl_label, mvtx_group);
0507     }
0508     else if (layer > 2 && layer < 7)
0509     {
0510       AlignmentDefs::getInttGlobalLabels(surf, glbl_label, intt_group);
0511     }
0512     else if (layer < 55)
0513     {
0514       AlignmentDefs::getTpcGlobalLabels(surf, ckey, glbl_label, tpc_group);
0515     }
0516     else if (layer < 57)
0517     {
0518       AlignmentDefs::getMMGlobalLabels(surf, glbl_label, mms_group);
0519     }
0520
0521     if (Verbosity() > 1)
0522     {
0523       std::cout << std::endl;
0524     }
0525
0526     float glbl_derivative[SvtxAlignmentState::NRES][SvtxAlignmentState::NGL]{};
0527     float lcl_derivative[SvtxAlignmentState::NRES][SvtxAlignmentState::NLOC]{};
0528
0529     /// For N residual local coordinates x, z
0530     for (int i = 0; i < SvtxAlignmentState::NRES; ++i)
0531     {
0532       // Add the measurement separately for each coordinate direction to Mille
0533       for (int j = 0; j < SvtxAlignmentState::NGL; ++j)
0534       {
0535         glbl_derivative[i][j] = state->get_global_derivative_matrix()(i, j);
0536
0537         if (is_layer_fixed(layer) ||
0538             is_layer_param_fixed(layer, j, fixed_layer_gparams))
0539         {
0540           glbl_derivative[i][j] = 0.;
0541         }
0542
0543         if (trkrid == TrkrDefs::mvtxId)
0544         {
0545           // need stave to get clamshell
0546           auto stave = MvtxDefs::getStaveId(ckey);
0547           auto clamshell = AlignmentDefs::getMvtxClamshell(layer, stave);
0548           if (is_layer_param_fixed(layer, j, fixed_layer_gparams) ||
0549               is_mvtx_layer_fixed(layer, clamshell))
0550           {
0551             glbl_derivative[i][j] = 0;
0552           }
0553         }
0554         else if (trkrid == TrkrDefs::inttId)
0555         {
0556           if (is_layer_param_fixed(layer, j, fixed_layer_gparams) ||
0557               is_layer_fixed(layer))
0558           {
0559             glbl_derivative[i][j] = 0;
0560           }
0561         }
0562         if (trkrid == TrkrDefs::tpcId)
0563         {
0564           auto sector = TpcDefs::getSectorId(ckey);
0565           auto side = TpcDefs::getSide(ckey);
0566           if (is_layer_param_fixed(layer, j, fixed_layer_gparams) ||
0567               is_tpc_sector_fixed(layer, sector, side) ||
0568               is_layer_fixed(layer))
0569           {
0570             glbl_derivative[i][j] = 0.0;
0571           }
0572         }
0573       }
0574
0575       for (int j = 0; j < SvtxAlignmentState::NLOC; ++j)
0576       {
0577         lcl_derivative[i][j] = state->get_local_derivative_matrix()(i, j);
0578
0579         if (is_layer_param_fixed(layer, j, fixed_layer_lparams))
0580         {
0581           lcl_derivative[i][j] = 0.;
0582         }
0583       }
0584       if (Verbosity() > 2)
0585       {
0586         std::cout << "coordinate " << i << " has residual " << residual(i) << " and clus_sigma " << clus_sigma(i) << std::endl
0587                   << "global deriv " << std::endl;
0588
0589         for (float k : glbl_derivative[i])
0590         {
0591           if (k > 0 || k < 0)
0592           {
0593             std::cout << "NONZERO GLOBAL DERIVATIVE" << std::endl;
0594           }
0595           std::cout << k << ", ";
0596         }
0597         std::cout << std::endl
0598                   << "local deriv " << std::endl;
0599         for (float k : lcl_derivative[i])
0600         {
0601           std::cout << k << ", ";
0602         }
0603         std::cout << std::endl;
0604       }
0605
0606       if (clus_sigma(i) < 1.0)  // discards crazy clusters
0607       {
0608         if (Verbosity() > 3)
0609         {
0610           std::cout << "ckey " << ckey << " and layer " << layer << " buffers:" << std::endl;
0611           AlignmentDefs::printBuffers(i, residual, clus_sigma, lcl_derivative[i], glbl_derivative[i], glbl_label);
0612         }
0613         float errinf = 1.0;
0614         if (m_layerMisalignment.find(layer) != m_layerMisalignment.end())
0615         {
0616           errinf = m_layerMisalignment.find(layer)->second;
0617         }
0618
0619         _mille->mille(SvtxAlignmentState::NLOC, lcl_derivative[i], SvtxAlignmentState::NGL, glbl_derivative[i], glbl_label, residual(i), errinf * clus_sigma(i));
0620       }
0621     }
0622
0623     float ntp_data[] = {
0624       lcl_derivative[0][0], lcl_derivative[0][1], lcl_derivative[0][2], lcl_derivative[0][3], lcl_derivative[0][4],
0625       glbl_derivative[0][0], glbl_derivative[0][1], glbl_derivative[0][2], glbl_derivative[0][3], glbl_derivative[0][4], glbl_derivative[0][5],
0626       lcl_derivative[1][0], lcl_derivative[1][1], lcl_derivative[1][2], lcl_derivative[1][3], lcl_derivative[1][4],
0627       glbl_derivative[1][0], glbl_derivative[1][1], glbl_derivative[1][2], glbl_derivative[1][3], glbl_derivative[1][4], glbl_derivative[1][5],
0628     };
0629
0630     if (m_ntuple) m_ntuple->Fill(ntp_data);
0631   }
0632
0633   return;
0634 }
0635
0636 bool MakeMilleFiles::is_layer_fixed(unsigned int layer)
0637 {
0638   bool ret = false;
0639   auto it = fixed_layers.find(layer);
0640   if (it != fixed_layers.end())
0641   {
0642     ret = true;
0643   }
0644
0645   return ret;
0646 }
0647 void MakeMilleFiles::set_layers_fixed(unsigned int minlayer,
0648                                       unsigned int maxlayer)
0649 {
0650   for (unsigned int i = minlayer; i < maxlayer; i++)
0651   {
0652     fixed_layers.insert(i);
0653   }
0654 }
0655 void MakeMilleFiles::set_layer_fixed(unsigned int layer)
0656 {
0657   fixed_layers.insert(layer);
0658 }
0659
0660 bool MakeMilleFiles::is_mvtx_layer_fixed(unsigned int layer, unsigned int clamshell)
0661 {
0662   bool ret = false;
0663
0664   std::pair pair = std::make_pair(layer, clamshell);
0665   auto it = fixed_mvtx_layers.find(pair);
0666   if (it != fixed_mvtx_layers.end())
0667   {
0668     ret = true;
0669   }
0670
0671   return ret;
0672 }
0673 void MakeMilleFiles::set_mvtx_layer_fixed(unsigned int layer, unsigned int clamshell)
0674 {
0675   fixed_mvtx_layers.insert(std::make_pair(layer, clamshell));
0676 }
0677 bool MakeMilleFiles::is_layer_param_fixed(unsigned int layer, unsigned int param, std::set<std::pair<unsigned int, unsigned int>>& param_fixed)
0678 {
0679   bool ret = false;
0680   std::pair<unsigned int, unsigned int> pair = std::make_pair(layer, param);
0681   auto it = param_fixed.find(pair);
0682   if (it != param_fixed.end())
0683   {
0684     ret = true;
0685   }
0686
0687   return ret;
0688 }
0689
0690 void MakeMilleFiles::set_layer_gparam_fixed(unsigned int layer, unsigned int param)
0691 {
0692   fixed_layer_gparams.insert(std::make_pair(layer, param));
0693 }
0694 void MakeMilleFiles::set_layer_lparam_fixed(unsigned int layer, unsigned int param)
0695 {
0696   fixed_layer_lparams.insert(std::make_pair(layer, param));
0697 }
0698 bool MakeMilleFiles::is_tpc_sector_fixed(unsigned int layer, unsigned int sector, unsigned int side)
0699 {
0700   bool ret = false;
0701   unsigned int region = AlignmentDefs::getTpcRegion(layer);
0702   unsigned int subsector = region * 24 + side * 12 + sector;
0703   auto it = fixed_sectors.find(subsector);
0704   if (it != fixed_sectors.end())
0705   {
0706     ret = true;
0707   }
0708
0709   return ret;
0710 }