24 #include <Acts/Geometry/GeometryIdentifier.hpp>
31 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
32 #include <ActsExamples/EventData/TrkrClusterSourceLink.hpp>
48 template<
class T>
inline constexpr
T square(
const T&
x ) {
return x*
x; }
53 template<
class T>
inline constexpr
T deltaPhi(
const T&
phi)
56 return phi - 2. *
M_PI;
57 else if (phi <= -
M_PI)
58 return phi + 2.*
M_PI;
89 std::cout <<
"Starting PHTpcResiduals event "
93 std::cout <<
"PHTpcResiduals processed " <<
m_event
94 <<
" events" << std::endl;
99 std::cout <<
"Finished PHTpcResiduals event "
109 std::cout <<
"PHTpcResiduals::End - writing matrices to " <<
m_outputfile << std::endl;
111 { std::cout <<
"PHTpcResiduals::End - Number of bad SL propagations " <<
m_nBadProps << std::endl; }
116 std::unique_ptr<TFile> outputfile( TFile::Open(
m_outputfile.c_str(),
"RECREATE" ) );
125 for(
const auto o:std::initializer_list<TObject*>({
h_rphiResid,
h_zResid,
h_etaResidLayer,
h_zResidLayer,
h_etaResid,
h_index,
h_alpha,
h_beta,
h_deltarphi_layer,
h_deltaz_layer,
residTup }) )
126 {
if( o ) o->Write(); }
137 { std::cout <<
"PHTpcResiduals::processTracks - proto track size " <<
m_trackMap->
size() <<std::endl; }
142 std::cout <<
"Processing track key " << trackKey
158 std::cout <<
"Track has pt " << track->
get_pt() << std::endl;
168 auto key = *clusIter;
179 std::cout <<
"Number of mvtx/intt/MM hits "
180 << nMvtxHits <<
"/" << nInttHits <<
"/"
181 << nMMHits << std::endl;
186 if(nMvtxHits<2 || nInttHits<2 || nMMHits<2)
189 if(nMvtxHits<2 || nInttHits<2 )
253 auto iter = surfMap.find(hitsetkey);
254 if(iter != surfMap.end())
271 auto surfvec = iter->second;
272 return surfvec.at(surfkey);
293 double p = track->
get_p();
296 for(
int i =0; i<6; i++)
298 for(
int j =0; j<6; j++)
307 const auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(
position);
321 std::cout <<
"Propagating silicon+MM fit params momentum: "
322 << track->
get_p() <<
" and position "
324 <<
", " << track->
get_z() <<
" cm "
348 auto pathLength = (*result).first;
349 auto trackStateParams = std::move(*(*result).second);
352 std::cout <<
"PHTpcResiduals::processTrack -"
353 <<
" path length: " << pathLength
354 <<
" track momentum : "
355 << trackParams.momentum()
356 <<
" propagator momentum : "
357 << trackStateParams.momentum()
373 std::cout <<
"Starting track params position/momentum: "
375 << std::endl << trackParams.momentum().transpose()
377 <<
"Track params phi/eta "
378 << std::atan2(trackParams.momentum().y(),
379 trackParams.momentum().x())
381 << std::atanh(trackParams.momentum().z() /
382 trackParams.momentum().norm())
396 return std::visit([params, sl,
this]
398 using InputMagneticField =
399 typename std::decay_t<decltype(inputField)>::element_type;
404 MagneticField field(inputField);
416 Acts::LoggerWrapper{*logger});
418 auto result = propagator.propagate(params, sl.referenceSurface(),
427 return result.error();
430 std::move(m_tGeometry->magField));
449 const auto momentum = params.momentum();
457 for (
int i = 0; i < 6; ++i)
458 for (
int j = 0; j < 6; ++j)
459 { state.
set_error(i, j, globalCov(i,j)); }
474 clusPhi = std::atan2(globClusPos(1), globClusPos(0));
475 clusZ = globClusPos(2);
482 std::cout <<
"cluster key is " <<
cluskey <<std::endl;
483 std::cout <<
"Cluster r phi and z " <<
clusR <<
" "
494 const auto globalStateMom = params.momentum();
495 const auto globalStateCov = *params.covariance();
508 const auto globStateZ =
stateZ;
513 const auto trackDrDt = (globStateX * globalStateMom(0) +
514 globStateY * globalStateMom(1)) / stateR;
515 const auto trackDxDr = globalStateMom(0) / trackDrDt;
516 const auto trackDyDr = globalStateMom(1) / trackDrDt;
517 const auto trackDzDr = globalStateMom(2) / trackDrDt;
519 const auto trackX = globStateX + dr * trackDxDr;
520 const auto trackY = globStateY + dr * trackDyDr;
521 const auto trackZ = globStateZ + dr * trackDzDr;
524 std::cout <<
"State Calculations: " << stateR <<
", "
525 << dr <<
", " << trackDrDt <<
", " << trackDxDr
526 <<
", " << trackDyDr <<
", " << trackDzDr
527 <<
" , " << trackX <<
", " << trackY <<
", "
528 << trackZ << std::endl;
530 statePhi = std::atan2(trackY, trackX);
534 std::cout <<
"State r phi and z "
549 std::cout <<
"TPC residuals " <<
drphi <<
" " <<
dz << std::endl;
552 = std::atanh(params.momentum().z() / params.absoluteMomentum());
553 const auto clusEta = std::atanh(
clusZ / std::sqrt(
558 const auto trackPPhi = -params.momentum()(0) * std::sin(statePhi) +
559 params.momentum()(1) * std::cos(statePhi);
560 const auto trackPR = params.momentum()(0) * std::cos(statePhi) +
561 params.momentum()(1) * std::sin(statePhi);
562 const auto trackPZ = params.momentum()(2);
563 const auto trackAlpha = -trackPPhi / trackPR;
564 const auto trackBeta = -trackPZ / trackPR;
567 std::cout <<
"Track angles " << trackPPhi <<
", " << trackPR
568 <<
", " << trackPZ <<
", " << trackAlpha <<
", " << trackBeta
575 const auto index =
getCell(globClusPos);
577 { std::cout <<
"Bin index found is " << index << std::endl; }
579 if(index < 0 )
return;
588 h_etaResid->Fill(trackEta, clusEta - trackEta);
641 float phi = std::atan2(loc(1), loc(0));
647 const auto r =
get_r( loc(0), loc(1) );
648 if( r < m_rMin || r >=
m_rMax )
return -1;
652 const auto z = loc(2);
653 if( z < m_zMin || z >=
m_zMax )
return -1;
670 m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
674 std::cout <<
PHWHERE <<
"No Acts surface maps, exiting"
683 std::cout <<
PHWHERE <<
"No TRKR_CLUSTER node on node tree. Exiting."
688 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
691 std::cout <<
"ActsTrackingGeometry not on node tree. Exiting."
697 m_trackMap = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxSiliconMMTrackMap");
701 std::cout <<
PHWHERE <<
"SvtxSiliconMMTrackMap not on node tree. Exiting."
717 h_beta =
new TH2F(
"betadz",
";tan#beta; #Deltaz [cm]",100,-0.5,0.5,100,-0.5,0.5);
718 h_alpha =
new TH2F(
"alphardphi",
";tan#alpha; r#Delta#phi [cm]", 100,-0.5,0.5,100,-0.5,0.5);
719 h_index =
new TH1F(
"index",
";index",total_bins, 0, total_bins);
720 h_rphiResid =
new TH2F(
"rphiResid",
";r [cm]; #Deltar#phi [cm]",
721 60, 20, 80, 500, -2, 2);
722 h_zResid =
new TH2F(
"zResid",
";z [cm]; #Deltaz [cm]",
723 200, -100, 100, 1000, -2, 2);
724 h_etaResid =
new TH2F(
"etaResid",
";#eta;#Delta#eta",
725 20, -1, 1, 500, -0.2, 0.2);
727 60, 20, 80, 500, -0.2, 0.2);
728 h_zResidLayer =
new TH2F(
"zResidLayer",
";r [cm]; #Deltaz [cm]",
729 60, 20, 80, 1000, -2, 2);
731 h_deltarphi_layer =
new TH2F(
"deltarphi_layer",
";layer; r.#Delta#phi_{track-cluster} (cm)", 57, 0, 57, 500, -2, 2 );
732 h_deltaz_layer =
new TH2F(
"deltaz_layer",
";layer; #Deltaz_{track-cluster} (cm)", 57, 0, 57, 100, -2, 2 );
734 residTup =
new TTree(
"residTree",
"tpc residual info");