24 (
const std::string&
t,
const std::string&
n,
const IInterface*
p)
26 m_TRTManagerName(
"TRT"),
28 m_maxTotalHits(21000),
30 m_minHitsForSegment(20),
31 m_minHitsAboveTOT(-1),
35 m_magneticField(true),
36 m_mergeSegments(false)
40 declareInterface<ITRT_TrackSegmentsMaker>(
this);
60 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::initialize(), March 2012"
71 return StatusCode::FAILURE;
81 if (m_minSeedTOT<0. || m_minSeedTOT>20.)
84 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::initialize(), jobProperties: "
92 return StatusCode::SUCCESS;
101 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::newEvent()");
106 if (not trtcontainer.
isValid()) {
107 msg(MSG::ERROR) <<
"Could not find TRT_DriftCircles collection!" <<
endmsg;
108 std::stringstream
msg;
110 throw std::runtime_error(
msg.str());
113 std::unique_ptr<TRT_TrackSegmentsMaker_BarrelCosmics::EventData>
114 event_data = std::make_unique<TRT_TrackSegmentsMaker_BarrelCosmics::EventData>(trtcontainer.
cptr());
116 for(
const auto *
it : *trtcontainer) {
119 if (!colNext) {
msg(MSG::WARNING) <<
"newEvent(): !colNext " <<
endmsg;
continue; }
121 for (
const auto *circleit : *colNext){
125 event_data->m_listHits.push_back( circleit );
127 double onMyTime = circleit->timeOverThreshold();
128 if ( circleit->firstBinHigh() && (circleit->trailingEdge() != 24) && (circleit->trailingEdge() != 0) ) onMyTime = ( (
double) circleit->trailingEdge() + 1 ) * 3.125;
131 const Amg::Vector3D ¢er = circleit->detectorElement()->surface(circleit->identify()).center();
132 event_data->m_listHitCenter.push_back( center );
138 ATH_MSG_DEBUG(
"skipping high occupancy event of " << event_data->m_listHits.size() <<
" barrel hits, limit at "
143 ATH_MSG_DEBUG(
"newEvent(): Number of TRT barrel hits: " << event_data->m_listHits.size()
144 <<
" Number of hits with TOT > " <<
m_minSeedTOT <<
": " << event_data->m_listHitCenter.size() );
146 return std::unique_ptr<InDet::ITRT_TrackSegmentsMaker::IEventData>(event_data.release());
149 std::unique_ptr<InDet::ITRT_TrackSegmentsMaker::IEventData>
152 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::newRegion()" );
155 if (not trtcontainer.
isValid()) {
156 msg(MSG::ERROR) <<
"m_trtcontainer is empty!!!" <<
endmsg;
157 std::stringstream
msg;
159 throw std::runtime_error(
msg.str());
162 std::unique_ptr<TRT_TrackSegmentsMaker_BarrelCosmics::EventData>
163 event_data = std::make_unique<TRT_TrackSegmentsMaker_BarrelCosmics::EventData>(trtcontainer.
cptr());
167 for ( InDet::TRT_DriftCircleContainer::const_iterator
w = trtcontainer->indexFind(
d);
w!=trtcontainer->end(); ++
w ) {
168 for(
const auto *circleit : **
w) {
172 event_data->m_listHits.push_back(circleit);
174 double onMyTime = circleit->timeOverThreshold();
175 if ( circleit->firstBinHigh() && (circleit->trailingEdge() != 24) && (circleit->trailingEdge() != 0) ) onMyTime = ( (
double) circleit->trailingEdge() + 1 ) * 3.125;
177 const Amg::Vector3D ¢er = circleit->detectorElement()->surface(circleit->identify()).center();
178 event_data->m_listHitCenter.push_back( center );
184 ATH_MSG_DEBUG(
"skipping high occupancy event of " << event_data->m_listHits.size() <<
" barrel hits, limit at "
189 ATH_MSG_DEBUG(
"newRegion(): Number of TRT barrel hits: " << event_data->m_listHits.size()
190 <<
" Number of hits with TOT > " <<
m_minSeedTOT <<
": " << event_data->m_listHitCenter.size() );
192 return std::unique_ptr<InDet::ITRT_TrackSegmentsMaker::IEventData>(event_data.release());
199 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::endEvent()" );
203 msg(MSG::WARNING) <<
"endEvent() you called the function t create the segments but not retrived them later??" <<
endmsg;
204 msg(MSG::WARNING) <<
"endEvent() deleting remaining elements of m_segments" <<
endmsg;
222 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::find()" );
228 msg(MSG::WARNING) <<
"TRT_TrackSegmentsMaker_BarrelCosmics::find() probably called twice per event? or newEvent / newRegion have not been called. check your program" <<
endmsg;
229 event_data.
clear();
return;
232 std::vector<double> x0,
phi, pivotX, pivotY, Xparabola, cotanParabola, InverseR;
233 std::vector<int> nHitsPosY, nHitsNegY;
236 int nIterations(15), countCalls(0);
237 while (countCalls<150) {
239 double xRange(1000.),
phiRange(M_PI_4);
240 double par[] = {0., M_PI_2, 0., 0., 0., 0., 0., 0.};
243 for (
int i=0;
i<nIterations;
i++) {
247 if ( !foundSeed )
break;
251 if ( xRange < 0.01 ||
phiRange < 0.00001)
break;
253 if (!foundSeed)
break;
258 int countAssociatedHits[] = {0, 0};
260 double sinphi = std::sqrt( 1. - cosphi*cosphi );
264 double measx[200], measy[200];
269 double a = (
par[3]-
it->x())*sinphi+(
it->y()-
par[4])*cosphi;
273 if (
m_magneticField && countMeas<200) { measx[countMeas] =
it->x(); measy[countMeas] =
it->y(); countMeas++; }
275 countAssociatedHits[(
it->y()>0)?0:1]++;
279 if ( countAssociatedHits[0]+countAssociatedHits[1] <
m_minHitsAboveTOT )
continue;
283 ATH_MSG_DEBUG(
"countAssociatedHits " << countAssociatedHits[0] <<
" "
285 x0.push_back(
par[0] );
phi.push_back(
par[1] ); nHitsPosY.push_back( countAssociatedHits[0] ); nHitsNegY.push_back( countAssociatedHits[1] );
287 pivotX.push_back(
par[3] ); pivotY.push_back(
par[4] ); Xparabola.push_back(
par[5] ); cotanParabola.push_back(
par[6] ); InverseR.push_back(
par[7] );
293 int nFoundSegments = x0.size();
if (nFoundSegments>10) nFoundSegments = 10;
295 for (
int i=0;
i<nFoundSegments;
i++) {
303 for (
unsigned int i=0;
i<event_data.
m_listHits.size();
i++) {
306 for (
int j=0; j<nFoundSegments; j++) {
308 if (nHitsPosY[j]<5 && center.y()>0.)
continue;
309 if (nHitsNegY[j]<5 && center.y()<0.)
continue;
312 double sinphi = std::sqrt(1./(1+cotanParabola[j]*cotanParabola[j]));
313 double cosphi = std::sqrt(1.-sinphi*sinphi);
if (cotanParabola[j]<0.) cosphi *= -1.;
314 double a = (pivotX[j]+Xparabola[j]-center.x())*sinphi+(center.y()-pivotY[j])*cosphi;
320 double sinphi = std::sqrt(1.-cosphi*cosphi);
321 residual = (x0[j]-center.x())*sinphi+center.y()*cosphi;
342 double mergeX0(100.), mergePhi(0.1);
343 for (
int i=0;
i<nFoundSegments;
i++) {
344 for (
int j=
i+1; j<nFoundSegments; j++) {
345 if ( std::abs(x0[
i]-x0[j])<mergeX0 && std::abs(
phi[
i]-
phi[j])<mergePhi ) {
348 mergeX0 = std::abs(x0[
i]-x0[j]);
349 mergePhi = std::abs(
phi[
i]-
phi[j]);
353 if (mergeI != mergeJ) {
354 ATH_MSG_DEBUG(
"Merge segments " << mergeI <<
" and " << mergeJ <<
" of size "
356 <<
"; difference in the impact par x: " << mergeX0 <<
" phi: " << mergePhi );
367 msg(
MSG::DEBUG) <<
"find() debug (" << nFoundSegments <<
")" ;
376 double trackpar[] = {x0[
i],
phi[
i]};
394 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::next(): return "
412 double *bestParameters,
418 double searchWindow = (xbin>2.) ? xbin : 2.;
422 int binInX[100];
for (
int &
i : binInX)
i = 0;
426 double phi = phimin+(0.5+j)*phibin;
432 double transformX =
k.x() * sinphi -
k.y() * cosphi;
433 indexmin = (
int) ((transformX-transformXMin-searchWindow)/xbin);
434 indexmax = (
int) ((transformX-transformXMin+searchWindow)/xbin);
443 bestParameters[0] = ( xbin *
index + transformXMin) / sinphi;
444 bestParameters[1] =
phi;
455 double sinphi = std::sqrt( 1. - cosphi*cosphi );
457 double parTransformX =
par[0] * sinphi;
459 double meanTransformY(0.);
460 int countMeanTransformY(0);
463 double transformX =
k.x() * sinphi -
k.y() * cosphi;
464 if ( std::abs(transformX-parTransformX) > 2. )
continue;
465 meanTransformY +=
k.x() * cosphi +
k.y() * sinphi;
466 countMeanTransformY++;
468 if (!countMeanTransformY) {
469 ATH_MSG_WARNING(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::findSeedInverseR(), no hits in the seed region???" );
472 meanTransformY /= (
double) countMeanTransformY;
481 double mean[] = { meanTransformY * cosphi + parTransformX * sinphi, meanTransformY * sinphi - parTransformX * cosphi };
482 int scanInverseR[101];
for (
int &
i : scanInverseR)
i = 0;
488 double transformX =
it.x() * sinphi -
it.y() * cosphi;
489 if ( std::abs(transformX-parTransformX) > 200. )
continue;
498 if (std::abs(
b)>0.0000001) {
499 double i1_tmp = ( window -
a ) /
b + 0.5;
500 double i2_tmp = ( -window -
a ) /
b + 0.5;
501 if (std::abs(i1_tmp)<1000.) { i1 = (
int) (std::ceil( i1_tmp )); i1 += 50; }
502 if (std::abs(i2_tmp)<1000.) { i2 = (
int) (std::ceil( i2_tmp )); i2 += 50; }
504 if (i1>100) { i1 = 100; }
else {
if (i1<0) i1 = 0; }
505 if (i2>100) { i2 = 100; }
else {
if (i2<0) i2 = 0; }
506 if (i1+i2==0 || i1+i2==200)
continue;
507 if (i1<=i2) {
for (
int i=i1;
i<=i2;
i++) scanInverseR[
i]++; }
508 else {
for (
int i=i2;
i<=i1;
i++) scanInverseR[
i]++; }
511 int iMax(-1), hitsMax(0);
512 for (
int i=0;
i<101;
i++)
if (scanInverseR[
i]>hitsMax) { iMax =
i; hitsMax = scanInverseR[
i]; }
514 par[2] = 0.00004 * (iMax-50);
549 double inverseSin3phi = 1. +
a[1]*
a[1]; inverseSin3phi *= std::sqrt(inverseSin3phi);
550 a[2] *= 2. / inverseSin3phi; }
564 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::convert() segment " << event_data.
m_segments.size()
565 <<
", N hits = " <<
hits.size());
567 if (
hits.size()<5) {
msg(MSG::ERROR) <<
"convert(): empty list of hits! size: " <<
hits.size() <<
endmsg;
return; }
573 int countOppositeSide(0);
581 if (countOppositeSide>5) {
582 ATH_MSG_DEBUG(
"convert(): removed " << countOppositeSide <<
" hits from the other side, N remaining hits: " <<
hits.size() );
585 ATH_MSG_WARNING(
"convert(): not enough hits after opposite side removal: " <<
hits.size() <<
", removed: "
586 << countOppositeSide );
591 double mean[10];
for (
double &
i :
mean)
i = 0.;
592 for (
auto & hit :
hits) {
594 mean[0] += (hit->detectorElement())->surface(hit->identify()).center().x();
595 mean[1] += (hit->detectorElement())->surface(hit->identify()).center().y();
597 for (
int i=0;
i<2;
i++)
mean[
i] /= (
double)
hits.size();
598 int iPivot(-1);
double yPivot(10000000.);
599 for (
unsigned int i=0;
i<
hits.size();
i++) {
600 double x = (
hits[
i]->detectorElement())->surface(
hits[
i]->identify()).center().x() -
mean[0];
601 double y = (
hits[
i]->detectorElement())->surface(
hits[
i]->identify()).center().y() -
mean[1];
612 if (
x2+
y2<yPivot) { yPivot =
x2+
y2; iPivot =
i; }
616 msg(MSG::ERROR) <<
"SF pivot index not found!!! " << yPivot <<
" " << iPivot <<
endmsg;
620 double cotanPhi =
mean[2] /
mean[4];
623 ATH_MSG_DEBUG(
"compare parameters phi: " << trackpar[1] <<
" vs. " <<
phi );
645 double Xparabola =
a[0];
646 double cotanParabola =
a[1];
647 double inverseSin3phi = 1. +
a[1]*
a[1]; inverseSin3phi *= std::sqrt(inverseSin3phi);
648 double inverseR = 2. *
a[2] / inverseSin3phi;
651 msg(
MSG::DEBUG) <<
"TRT_TrackSegmentsMaker_BarrelCosmics: parabola fit, X: " << Xparabola <<
endmsg;
652 msg(
MSG::DEBUG) <<
" parabola fit, cotan: " << cotanParabola <<
" compare to " << cotanPhi <<
endmsg;
656 qOverp = inverseR / 0.6;
666 double f33 = 4. / ( 3. * (1.+cotanPhi*cotanPhi) *
mean[4] );
669 0. ,45633., 0., 0., 0.,
670 0. , 0., f33, 0., 0.,
671 0. , 0., 0., 0.2, 0.,
675 (
hits[iPivot]->detectorElement())->surface(
hits[iPivot]->identify()).center().
z() );
685 if (phi<-M_PI || phi>0.)
msg(MSG::ERROR) <<
"phi value problem: " <<
phi <<
endmsg;
697 for (
const auto *DC :
hits) {
734 Amg::Vector3D faketransl(firstcenter.x(),firstcenter.y()-5.,firstcenter.z());
739 std::move(pseudoPar),
740 std::move(pseudocov),
742 rio.push_back( pseudo );
743 delete pseudoSurface;
747 chi2 /= ( 1. + cotanPhi*cotanPhi );
768 if (nhits<10)
return;
770 double shift[] = {0., 0.};
771 for (
int i=0;
i<nhits;
i++) { shift[0] += measx[
i]; shift[1] += measy[
i]; }
772 for (
double &
i : shift)
i /= (
double) nhits;
773 for (
int i=0;
i<nhits;
i++) { measx[
i] -= shift[0]; measy[
i] -= shift[1]; }
775 double mean[10];
for (
double &
i :
mean)
i = 0.;
776 for (
int i=0;
i<nhits;
i++) {
777 double x2 = measx[
i]*measx[
i];
778 double y2 = measy[
i]*measy[
i];
781 mean[2] += measx[
i]*measy[
i];
807 double Xparabola =
a[0];
808 double cotanParabola =
a[1];
809 double inverseSin3phi = 1. +
a[1]*
a[1]; inverseSin3phi *= std::sqrt(inverseSin3phi);
810 double inverseR = 2. *
a[2] / inverseSin3phi;
811 double sinphi = std::sqrt(1./(1+
a[1]*
a[1]));
812 double cosphi = std::sqrt(1.-sinphi*sinphi);
if (cotanParabola<0.) cosphi *= -1.;
817 result[3] = cotanParabola;
822 for (
int i=0;
i<nhits;
i++) {
823 double tmp = (Xparabola-measx[
i])*sinphi+measy[
i]*cosphi;
829 for (
int i=0;
i<nhits;
i++) { measx[
i] += shift[0]; measy[
i] += shift[1]; }
843 ATH_MSG_DEBUG(
"InDet::TRT_TrackSegmentsMaker_BarrelCosmics::find()");
848 ATH_MSG_WARNING(
"find probably called twice per event? or newEvent / newRegion have not been called. check your program" );
849 event_data.
clear();
return;
852 std::vector<double> x0,
phi;
853 std::vector<int> nHitsPosY, nHitsNegY;
856 int nIterations(15), countCalls(0);
857 while (countCalls<150) {
859 double xRange(1000.),
phiRange(M_PI_4);
860 double par[] = {0., M_PI_2};
863 for (
int i=0;
i<nIterations;
i++) {
867 if ( !foundSeed )
break;
871 if ( xRange < 0.01 ||
phiRange < 0.00001)
break;
874 if (!foundSeed)
break;
877 int countAssociatedHits[] = {0, 0};
879 double sinphi = std::sqrt(1.-cosphi*cosphi);
883 if ( std::abs( (
par[0]-
it->x())*sinphi +
it->y()*cosphi ) > 2. ) { ++
it;
continue; }
884 countAssociatedHits[(
it->y()>0)?0:1]++;
888 if ( countAssociatedHits[0]+countAssociatedHits[1] <
m_minHitsAboveTOT )
continue;
890 x0.push_back(
par[0] );
phi.push_back(
par[1] ); nHitsPosY.push_back( countAssociatedHits[0] ); nHitsNegY.push_back( countAssociatedHits[1] );
896 int nFoundSegments = x0.size();
897 if (nFoundSegments>10) nFoundSegments = 10;
898 for (
int i=0;
i<nFoundSegments;
i++) {
902 for (
unsigned int i=0;
i<event_data.
m_listHits.size();
i++) {
904 center = event_data.
m_listHits[
i]->detectorElement()->surface(event_data.
m_listHits[
i]->identify()).center();
906 for (
int j=0; j<nFoundSegments; j++) {
907 if (nHitsPosY[j]<5 && center.y()>0.)
continue;
908 if (nHitsNegY[j]<5 && center.y()<0.)
continue;
910 double sinphi = std::sqrt(1.-cosphi*cosphi);
911 if ( std::abs((x0[j]-center.x())*sinphi+center.y()*cosphi) < 2.) {
921 msg(
MSG::DEBUG) <<
"find() debug (" << nFoundSegments <<
")" ;
931 double trackpar[] = {x0[
i],
phi[
i]};