21 #include "TGraphErrors.h"
23 #include "TVirtualFitter.h"
36 (
const std::string&
t,
const std::string&
n,
const IInterface*
p)
39 declareInterface<ITRT_TrackSegmentsMaker>(
this);
86 return StatusCode::FAILURE;
93 return StatusCode::FAILURE;
102 return StatusCode::FAILURE;
116 std::string fit_name1=
"ztanphi";
117 std::string fit_name2=
"zphi";
118 std::string fit_name3=
"zphi_ap";
155 std::unique_ptr<InDet::ITRT_TrackSegmentsMaker::IEventData>
159 std::unique_ptr<TRT_TrackSegmentsMaker_ECcosmics::EventData>
160 event_data = std::make_unique<TRT_TrackSegmentsMaker_ECcosmics::EventData>();
167 it=event_data->m_goodHits.begin();
168 itE=event_data->m_goodHits.end();
180 else if(
phi>15 )
phi = 0 ;
185 int zslice=
int((std::abs(
sc.z())-800.)/100.);
187 event_data->m_sectors[
z][zslice][
phi].push_back(*
it);
191 ATH_MSG_DEBUG(
"Number of good driftcircles: "<<event_data->m_goodHits.size());
192 ATH_MSG_DEBUG(
"Number of noise driftcircles: "<<event_data->m_noiseHits.size());
194 return std::unique_ptr<InDet::ITRT_TrackSegmentsMaker::IEventData>(event_data.release());
201 std::unique_ptr<InDet::ITRT_TrackSegmentsMaker::IEventData>
203 (
const EventContext& ,
const std::vector<IdentifierHash>& vTRT)
const
206 std::unique_ptr<TRT_TrackSegmentsMaker_ECcosmics::EventData>
207 event_data = std::make_unique<TRT_TrackSegmentsMaker_ECcosmics::EventData>();
215 it=event_data->m_goodHits.begin();
216 itE=event_data->m_goodHits.end();
228 else if(
phi>15 )
phi = 0 ;
232 int zslice=
int((std::abs(
sc.z())-800.)/100.);
234 event_data->m_sectors[
z][zslice][
phi].push_back(*
it);
238 ATH_MSG_DEBUG(
"Number of good driftcircles: "<<event_data->m_goodHits.size());
239 ATH_MSG_DEBUG(
"Number of noise driftcircles: "<<event_data->m_noiseHits.size());
241 return std::unique_ptr<InDet::ITRT_TrackSegmentsMaker::IEventData>(event_data.release());
312 for(
int endcap=0;endcap<2;endcap++){
324 for(
int j=0;j<20;j++){
325 for(
int i=0;
i<16;
i++){
339 if(sector<0)
continue;
341 max=event_data.
m_sectors[endcap][zslice][sector].size();
353 event_data.
m_sectors[endcap][zslice][sector].clear();
365 std::vector<std::vector<const InDet::TRT_DriftCircle*> *>
::iterator sit,sitE;
367 sit=event_data.
m_seeds.begin();
370 for(;sit!=sitE;++sit){
395 std::vector<const InDet::TRT_DriftCircle*> dc;
398 dc.reserve(event_data.
m_sectors[endcap][zslice][sector].size());
413 for(;vit!=vitE;++vit){
424 if(dc.size()>=40)
return false;
432 for(;vit!=vitE;++vit){
443 shift=(phimax+phimin)/2.-M_PI_4;
444 ATH_MSG_VERBOSE(
"Phimin="<<phimin<<
" Phimax="<<phimax<<
" -> shift = "<<shift);
447 if(phimin<-M_PI_2 && phimax>M_PI_2){
449 ATH_MSG_VERBOSE(
"Boundary between pi and -pi!!! use the following shift: "<<shift);
454 for(vit=dc.begin();vit!=vitE;++vit){
458 double orig=
sc.phi();
460 double corr=orig-shift;
475 double p_best[3]={0,0,0};
479 for(
int j=
i+1;j<
count;j++){
516 for(
auto &
i : muster)
520 std::vector<const InDet::TRT_DriftCircle*> *seed=
new std::vector<const InDet::TRT_DriftCircle*>;
522 for(
int zsl=zslice-1;zsl<=zslice+1;zsl++){
523 for(
int ps=sector-1;ps<=sector+1;ps++){
526 if(zsl<0 || zsl>=20)
continue;
536 double pred_phi=
sc.z()*p_best[0]+p_best[1];
537 double corr=
sc.phi()-p_best[2];
540 if(
sc.z()==p_best[1]) {
550 if(std::abs(
phidiff(pred_phi,corr))<=4./650.){
551 seed->push_back(*
it);
554 muster[zsl-zslice+1][ps-sector+1]=1;
566 ATH_MSG_VERBOSE(
"\t"<<muster[0][0]<<
" "<<muster[1][0]<<
" "<<muster[2][0]);
567 ATH_MSG_VERBOSE(
"\t"<<muster[0][1]<<
" "<<muster[1][1]<<
" "<<muster[2][1]);
568 ATH_MSG_VERBOSE(
"\t"<<muster[0][2]<<
" "<<muster[1][2]<<
" "<<muster[2][2]);
571 if(muster[2][0] || muster[2][1] || muster[2][2]){
575 for(
int ps=sector-2;ps<=sector+2;ps++){
578 if(zsl<0 || zsl>=20)
continue;
588 double pred_phi=
sc.z()*p_best[0]+p_best[1];
589 double corr=
sc.phi()-p_best[2];
592 if(
sc.z()==p_best[1]) {
602 if(std::abs(
phidiff(pred_phi,corr))<=4./650.){
606 seed->push_back(*
it);
612 if(muster[0][0] || muster[0][1] || muster[0][2]){
616 for(
int ps=sector-2;ps<=sector+2;ps++){
619 if(zsl<0 || zsl>=20)
continue;
629 double pred_phi=
sc.z()*p_best[0]+p_best[1];
630 double corr=
sc.phi()-p_best[2];
633 if(
sc.z()==p_best[1]) {
643 if(std::abs(
phidiff(pred_phi,corr))<=4./650.){
647 seed->push_back(*
it);
665 for(;vit!=vitE;++vit){
671 if ( phi<0 ) phi = 15;
672 else if( phi>15 ) phi = 0 ;
676 int zslic=
int((std::abs(
sc.z())-800.)/100.);
678 event_data.
m_sectors[
z][zslic][phi].remove(*vit);
688 for(;lit!=litE;++lit){
698 if ( phi<0 ) phi = 15;
699 else if( phi>15 ) phi = 0 ;
702 int zsl=
int((std::abs(
sc.z())-800.)/100.);
704 if(std::abs(zsl-zslice)<2){
705 double pred_phi=
sc.z()*p_best[0]+p_best[1];
706 double corr=
sc.phi()-p_best[2];
709 if(
sc.z()==p_best[1]) {
719 if(std::abs(
phidiff(pred_phi,corr))<=4./650.){
721 seed->push_back(*lit);
734 event_data.
m_seeds.push_back(seed);
757 for(
int zsl=zslice-1;zsl<=zslice+1;zsl++){
758 for(
int ps=sector-1;ps<=sector+1;ps++){
761 if(zsl<0 || zsl>=20)
continue;
773 double pred_phi=
sc.z()*
p[0]+
p[1];
774 double corr=
sc.phi()-
p[2];
787 if(std::abs(
phidiff(pred_phi,corr))<=4./650.){
818 for(;vit!=vitE;++vit){
846 shift=(phimax+phimin)/2.-M_PI_4;
847 ATH_MSG_VERBOSE(
"Phimin="<<phimin<<
" Phimax="<<phimax<<
" shift="<<shift);
850 if(phimin<-M_PI_2 && phimax>M_PI_2){
852 ATH_MSG_VERBOSE(
"Boundary between pi and -pi!!! use the following shift: "<<shift);
861 double corr=orig-shift;
868 event_data.
m_a_tan_err[
i]=(1.15/700.)/(cosCorr*cosCorr);
880 std::vector<const InDet::TRT_DriftCircle*> tobeadded;
883 for(;lit!=litE;++lit){
887 double fitted_phi=
fit->Eval(
sc.z(),0.,0.);
890 if (fitted_phi >
M_PI) fitted_phi = std::fmod(fitted_phi+
M_PI,2.*
M_PI)-
M_PI;
891 else if(fitted_phi <-
M_PI) fitted_phi = std::fmod(fitted_phi-
M_PI,2.*
M_PI)+
M_PI;
899 for(;vit!=vitE && new_hit;++vit){
905 ATH_MSG_VERBOSE(
"\t\t good hit matched! fitted_phi="<<fitted_phi<<
" z=" <<
sc.z());
911 seed->push_back(*lit);
914 tobeadded.push_back(*lit);
920 std::sort(tobeadded.begin(),tobeadded.end(),
CompareDCz);
928 for(vit=tobeadded.begin();vit!=tobeadded.end();++vit){
951 seed->push_back(tobeadded[
i]);
965 seed->push_back(tobeadded[
i]);
978 std::sort(seed->begin(),seed->end(),
CompareDCz);
984 for(;vit!=vitE;++vit){
1009 shift=(phimax+phimin)/2.-M_PI_4;
1010 ATH_MSG_VERBOSE(
"Phimin="<<phimin<<
" Phimax="<<phimax<<
" shift="<<shift);
1013 if(phimin<-M_PI_2 && phimax>M_PI_2){
1015 ATH_MSG_VERBOSE(
"Boundary between pi and -pi!!! use the following shift: "<<shift);
1022 double corr=orig-shift;
1029 const auto cosCorr{
std::cos(corr)};
1030 event_data.
m_a_tan_err[
i]=(1.15/700.)/(cosCorr*cosCorr);
1041 for(;lit!=litE;++lit){
1049 double phi=
sc.phi();
1050 double fitted_phi=
fit->Eval(
sc.z(),0.,0.);
1053 if (fitted_phi >
M_PI) fitted_phi = std::fmod(fitted_phi+
M_PI,2.*
M_PI)-
M_PI;
1054 else if(fitted_phi <-
M_PI) fitted_phi = std::fmod(fitted_phi-
M_PI,2.*
M_PI)+
M_PI;
1057 if(
z*event_data.
m_a_z[0]<0)
continue;
1063 for(;vit!=vitE && new_hit;++vit){
1070 ATH_MSG_VERBOSE(
"\t\t noise hit matched! fitted_phi="<<fitted_phi<<
" z=" <<
sc.z());
1071 seed->push_back(*lit);
1073 ATH_MSG_VERBOSE(
"\t\t noise hit matched but too far away! fitted_phi="<<fitted_phi<<
" z=" <<
sc.z());
1080 std::sort(seed->begin(),seed->end(),
CompareDCy);
1084 Amg::Vector3D laststraw=((*seed)[seed->size()-1]->detectorElement()->strawCenter(
m_trtid->
straw((*seed)[seed->size()-1]->identify()) ));
1086 double z2=laststraw.z();
1087 double yfirst=firststraw.y();
1088 double ylast=laststraw.y();
1096 if (std::abs(yfirst-ylast)>50){
1100 std::sort(seed->begin(),seed->end(),
CompareDCz);
1104 if ((yfirst>0 && z1<0) || (yfirst<0 && z1>0)) std::sort(seed->begin(),seed->end(),
CompareDCz);
1107 firststraw=(*seed)[0]->detectorElement()->strawCenter(
m_trtid->
straw((*seed)[0]->identify()) );
1108 laststraw=(*seed)[seed->size()-1]->detectorElement()->strawCenter(
m_trtid->
straw((*seed)[seed->size()-1]->identify()) );
1114 Theta=
M_PI-std::atan2(400.,std::abs(z1-z2));
1117 Theta=std::atan2(400.,std::abs(z1-z2));
1121 double globaly=((*seed)[0]->detectorElement()->strawCenter(
m_trtid->
straw((*seed)[0]->identify()) )).
y();
1122 double firstphi=((*seed)[0]->detectorElement()->strawCenter(
m_trtid->
straw((*seed)[0]->identify()) )).phi();
1125 ATH_MSG_VERBOSE(
"Number of tube hits matching straight line: "<<seed->size());
1131 const Amg::Vector3D& sc_last = ((*seed)[seed->size()-1]->detectorElement()->strawCenter(
m_trtid->
straw((*seed)[seed->size()-1]->identify()) ));
1134 double tx1,tx2,ty1,ty2;
1142 double tphi=std::atan2(ty2-ty1,tx2-tx1);
1167 ATH_MSG_VERBOSE(
"globaly = "<<globaly<<
" phi="<<firstphi<<
" theta="<<Theta<<
" --> state = "<<state);
1174 for(;vit!=vitE;++vit){
1194 for(;vit!=vitE;++vit){
1198 double fitted_phi=
fit->Eval(
sc.z(),0.,0.);
1201 if (fitted_phi >
M_PI) fitted_phi = std::fmod(fitted_phi+
M_PI,2.*
M_PI)-
M_PI;
1202 else if(fitted_phi <-
M_PI) fitted_phi = std::fmod(fitted_phi-
M_PI,2.*
M_PI)+
M_PI;
1205 if(fitted_phi<
sc.phi())
sign=-1;
1209 double drifte = 2*std::sqrt((*vit)->localCovariance()(0,0));
1211 double straw_r=640+400./std::abs(z2-z1)*std::abs(
sc.z()-z1);
1213 double dphi=
sign*driftr/straw_r;
1242 for(;vit!=vitE;++vit){
1245 double fitted_phi=
fit->Eval(
sc.z(),0.,0.);
1249 double drifte = 2*std::sqrt((*vit)->localCovariance()(0,0));
1252 if(fitted_phi<
sc.phi())
sign=-1;
1254 double straw_r=640+400./std::abs(z2-z1)*std::abs(
sc.z()-z1);
1255 double dphi=
sign*driftr/straw_r;
1293 for(;vit!=vitE;++vit){
1297 double fitted_phi=
fit->Eval(
sc.z(),0.,0.);
1300 if (fitted_phi >
M_PI) fitted_phi = std::fmod(fitted_phi+
M_PI,2.*
M_PI)-
M_PI;
1301 else if(fitted_phi <-
M_PI) fitted_phi = std::fmod(fitted_phi-
M_PI,2.*
M_PI)+
M_PI;
1304 double drifte = 2*std::sqrt((*vit)->localCovariance()(0,0));
1307 if(fitted_phi<
sc.phi())
sign=-1;
1309 double straw_r=640+400./std::abs(z2-z1)*std::abs(
sc.z()-z1);
1310 double dphi=
sign*driftr/straw_r;
1331 if(
count>=400)
return;
1344 double initial_r=-10.;
1345 double initial_locz=0.;
1348 for(;vit!=vitE;++vit){
1350 (&((*vit)->detectorElement())->surface( (*vit)->identify()));
1361 if(std::abs(z2-z1)>0.000001){
1363 locz=400./std::abs(z2-z1)*std::abs(
z-z1)-200.;
1365 locz=200.-400./std::abs(z2-z1)*std::abs(
z-z1);
1368 double fitted_phi=
fit->Eval(
z,0.,0.);
1371 if (fitted_phi >
M_PI) fitted_phi = std::fmod(fitted_phi+
M_PI,2.*
M_PI)-
M_PI;
1372 else if(fitted_phi <-
M_PI) fitted_phi = std::fmod(fitted_phi-
M_PI,2.*
M_PI)+
M_PI;
1374 double fitted_r=(
phidiff(fitted_phi,
sc.phi()))*(840+locz);
1375 double temp_phi=fitted_phi;
1377 if(initial_r==-10. && std::abs(fitted_r)<4.0){
1400 bool useDrift=
false;
1404 double drifte = 2*std::sqrt((*vit)->localCovariance()(0,0));
1407 if(fitted_phi<
sc.phi())
sign=-1;
1409 double straw_r=640+400./std::abs(z2-z1)*std::abs(
sc.z()-z1);
1410 double dphi=
sign*driftr/straw_r;
1412 if(std::abs(
phidiff(fitted_phi,
sc.phi()-dphi)) < 3* drifte/straw_r){
1414 chi2+=(std::abs(fitted_r)-std::abs(driftr))*(std::abs(fitted_r)-std::abs(driftr))/(drifte*drifte);
1422 rio.push_back(
m_riomakerD->correct(*(*vit),Tp, Gaudi::Hive::currentContext()));
1425 chi2+=(fitted_r/1.15)*(fitted_r/1.15);
1427 rio.push_back(
m_riomakerN->correct(*(*vit),Tp, Gaudi::Hive::currentContext()));
1442 if(state==2 || state==3){
1458 rio.push_back(pseudo);
1467 const Trk::Surface * sur = &((*seed)[0]->detectorElement())->surface((*seed)[0]->identify());
1472 if(firstphi>0) firstphi=-firstphi;
1481 std::array<Trk::DefinedParameter,5> defPar = {dp0,dp1,dp2,dp3,dp4};
1488 double universal=1.15;
1493 universal*universal, 0., 0., 0., 0.,
1494 0. , 160000./12., 0., 0., 0.,
1495 0. , 0., 0.1, 0., 0.,
1496 0. , 0., 0., 1., 0.,
1497 0. , 0., 0., 0., 1.;
1556 if (not trtcontainer.
isValid()) {
1565 if (!prd_to_track_map.
isValid()) {
1568 prd_to_track_map_cptr = prd_to_track_map.
cptr();
1572 if (!mjo_trtcontainer)
return;
1574 InDet::TRT_DriftCircleContainer::const_iterator
1575 w = trtcontainer->begin(),we = trtcontainer->end();
1581 InDet::TRT_DriftCircleCollection::const_iterator
1582 c = (*w)->begin(), ce = (*w)->end();
1586 if(prd_to_track_map_cptr && prd_to_track_map_cptr->
isUsed(*(*
c)))
continue;
1588 if( (*c)->isNoise() || (*c)->timeOverThreshold()<
m_cutToTLoose) {
1616 const double dPhiLimit{0.03};
1617 const double dzLimit{31.};
1695 std::list<const InDet::TRT_DriftCircle*> & container,
1696 double dPhiLimit,
double dzLimit)
const{
1698 const auto id = (*it)->identify();
1705 for(
auto it2=container.begin();it2!=container.end();++it2){
1706 if(
it==it2)
continue;
1711 if(endcap!=endcap2 || wheel2!=wheel)
continue;
1713 const Amg::Vector3D& sc2 = (*it2)->detectorElement()->strawCenter(ns2);
1714 double dz=
sc.z()-sc2.z();
1715 double dphi=
phidiff(
sc.phi(),sc2.phi());
1718 dphi=std::abs(dphi/300.0);
1719 if(dphi<dPhiLimit and dz<dzLimit){
1738 std::lock_guard<std::mutex> lock(
s_fitMutex);
1739 TVirtualFitter::SetMaxIterations(100000);
1753 <<event_data.
m_fitf_ztanphi->GetParameter(2)<<
" fitresult="<<fitresult1);
1778 if(fitresult1!=0 && fitresult2!=0)
return nullptr;
1780 if(fitresult1==0 && fitresult2!=0)
return event_data.
m_fitf_zphi;
1798 ATH_MSG_VERBOSE(
"Number of hits matching this fit: tan="<<match_tan<<
" pol2="<<match_phi);
1800 if(match_tan>match_phi){
1802 }
else if(match_tan<match_phi){
1821 std::vector<const InDet::TRT_DriftCircle*> *seed)
const
1825 double meanz=0,varz=0;
1835 for(;vit!=vitE;++vit){
1843 varz+=
sc.z()*
sc.z();
1846 if(
count < 2)
return true;
1848 varz=(
count*varz-meanz*meanz);
1853 varz=std::sqrt(varz);
1855 double zmin_sel=100000;
1861 if(std::abs(dcz[
i]-dcz[
i-1])<
zmin)
zmin=std::abs(dcz[
i]-dcz[
i-1]);
1863 if(std::abs(dcz[
i]-dcz[
i+1])<
zmin)
zmin=std::abs(dcz[
i]-dcz[
i+1]);
1881 if(std::abs(zmin_sel-
mean)>2*
var){
1883 ATH_MSG_VERBOSE(
"Hit is suspicious, remove it! "<<std::abs(zmin_sel-
mean)<<
" , "<<2*
var<<
" -> "<<zmin_sel<<
" : "<<dcz[
index]<<
" <-> "<<meanz<<
" ("<<varz<<
")");
1888 ATH_MSG_VERBOSE(
"Hit seems to be ok: "<<std::abs(zmin_sel-
mean)<<
" < "<<2*
var<<
" && "<<std::abs(dcz[
index]-meanz)<<
" < "<<1.5*varz);