From f8a8a79b1ccf568f489bc4558604e89a15c533e4 Mon Sep 17 00:00:00 2001
From: John Derek Chapman <chapman@hep.phy.cam.ac.uk>
Date: Mon, 11 Jun 2018 15:24:23 +0000
Subject: [PATCH] Merge branch 'MM_mTPC' into '21.3'
Implementation of the mTPC reconstruction for MM
See merge request atlas/athena!11776
(cherry picked from commit 937f65804e4ff852f9f9475a0b2006ca81251131)
5e0d3047 Implementation of the mTPC reconstruction for MM
3b665187 Merge remote-tracking branch 'upstream/21.3' into MM_mTPC
---
.../src/MM_FastDigitizer.cxx | 99 +++++++++++----
.../src/MM_FastDigitizer.h | 1 +
.../src/MuonClusterOnTrackCreator.cxx | 54 ++++++--
.../src/MuonClusterOnTrackCreator.h | 2 +
.../MuonTrackFindingEvent/MuPatHit.h | 2 +-
.../src/MuonClusterSegmentFinderTool.cxx | 118 +++++++++++++-----
.../src/MuonClusterSegmentFinderTool.h | 5 +-
.../src/MuPatHitTool.cxx | 14 ++-
.../src/CombinedMuonTrackBuilder.cxx | 6 +-
.../TrkGlobalChi2Fitter/GlobalChi2Fitter.h | 1 +
.../TrkGlobalChi2Fitter/src/GXFTrajectory.cxx | 3 +-
.../src/GlobalChi2Fitter.cxx | 46 ++++++-
.../src/RIO_OnTrackCreator.cxx | 4 +-
13 files changed, 283 insertions(+), 72 deletions(-)
diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.cxx b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.cxx
index 8a7718e50bd..aae0346ccde 100644
--- a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.cxx
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.cxx
@@ -66,6 +66,7 @@ using namespace Muon;
declareProperty("CheckIds", m_checkIds = false, "Turn on validity checking of Identifiers" );
declareProperty("MaskMultiplet", m_maskMultiplet = 0, "0: all, 1: first, 2: second, 3: both" );
declareProperty("SDOname", m_sdoName = "MMfast_SDO" );
+ declareProperty("MicroTPC", m_microTPC = true, "Turn on microTPC mode" );
}
/*******************************************************************************/
MM_FastDigitizer::~MM_FastDigitizer() {
@@ -339,7 +340,8 @@ StatusCode MM_FastDigitizer::execute() {
double scale;//, scaletop;
// double gasgap = 5.;
- scale = -slpos.z()/ldir.z();
+ if (m_microTPC) scale = 0;
+ else scale = -slpos.z()/ldir.z();
// scaletop = (gasgap+slpos.z())/ldir.z();
Amg::Vector3D hitOnSurface = slpos + scale*ldir;
@@ -352,6 +354,8 @@ StatusCode MM_FastDigitizer::execute() {
Amg::Vector3D hitAfterTimeShift(hitOnSurface.x(),hitOnSurface.y(),shiftTimeOffset);
Amg::Vector3D hitAfterTimeShiftOnSurface = hitAfterTimeShift - (shiftTimeOffset/ldirTime.z())*ldirTime;
+ double tdrift = 0;
+
// resolution = -.01/3 * angle + .64/3.
double resolution;
if (fabs(inAngle_XZ)<3)
@@ -372,8 +376,8 @@ StatusCode MM_FastDigitizer::execute() {
if( fabs(inAngle_XZ) > 70 ){
posOnSurf[0]=(slpos.x()+sp);
// digiMode = 1;
- // if using timing information use hit position after m_shift
- }else if( m_useTimeShift ){
+ // if using timing information use hit position after shift
+ }else if( m_useTimeShift && !m_microTPC ){
posOnSurf[0]=(hitAfterTimeShiftOnSurface.x()+sp);
// digiMode = 2;
}
@@ -497,23 +501,83 @@ StatusCode MM_FastDigitizer::execute() {
std::vector<Identifier> rdoList;
rdoList.push_back(id);
- Amg::MatrixX* cov = new Amg::MatrixX(1,1);
- cov->setIdentity();
- (*cov)(0,0) = resolution*resolution;
- MMPrepData* prd = new MMPrepData( id,hash,posOnSurf,rdoList,cov,detEl);
- prd->setHashAndIndex(col->identifyHash(), col->size()); // <<< add this line to the MM code as well
- col->push_back(prd);
+
ATH_MSG_VERBOSE("Global hit: r " << hit.globalPosition().perp() << " phi " << hit.globalPosition().phi() << " z " << hit.globalPosition().z());
- ATH_MSG_VERBOSE(" Prd: r " << prd->globalPosition().perp() << " phi " << prd->globalPosition().phi() << " z " << prd->globalPosition().z());
ATH_MSG_VERBOSE(" Calculated truth hitOnSurfaceGlobal: r " << hitOnSurfaceGlobal.perp() << " phi " << hitOnSurfaceGlobal.phi() << " z " << hitOnSurfaceGlobal.z());
ATH_MSG_VERBOSE(" detEl: r " << repos.perp() << " phi " << repos.phi() << " z " << repos.z());
ATH_MSG_VERBOSE(" Surface center: r " << surf.center().perp() << " phi " << surf.center().phi() << " z " << surf.center().z());
-
ATH_MSG_VERBOSE("Local hit in Det Element Frame: x " << hit.localPosition().x() << " y " << hit.localPosition().y() << " z " << hit.localPosition().z());
- ATH_MSG_VERBOSE(" Prd: local posOnSurf.x() " << posOnSurf.x() << " posOnSurf.y() " << posOnSurf.y() );
- // << " detEl: x " << lpos.x() << " y " << lpos.y() << " z " << lpos.z() << " strip " << stripNumber);
ATH_MSG_DEBUG(" hit: " << m_idHelperTool->toString(id) << " hitx " << posOnSurf.x() << " hitOnSurface.x() " << hitOnSurface.x() << " residual " << posOnSurf.x() - hitOnSurface.x()
<< " pull " << (posOnSurf.x() - hitOnSurface.x())/resolution );
+ Amg::Vector3D CurrentHitInDriftGap = slpos;
+ // emulating micro track in the drift volume for microTPC
+ if (!m_microTPC) {
+ Amg::MatrixX* cov = new Amg::MatrixX(1,1);
+ cov->setIdentity();
+ (*cov)(0,0) = resolution*resolution;
+ MMPrepData* prd = new MMPrepData( id,hash,Amg::Vector2D(posOnSurf.x(),0.),rdoList,cov,detEl,(int)tdrift,0);
+ prd->setHashAndIndex(col->identifyHash(), col->size()); // <<< add this line to the MM code as well
+ col->push_back(prd);
+
+ std::vector<MuonSimData::Deposit> deposits;
+ MuonSimData::Deposit deposit(hit.particleLink(), MuonMCData(hitOnSurface.x(),hitOnSurface.y()));
+ deposits.push_back(deposit);
+
+ MuonSimData simdata(deposits,0);
+ simdata.setPosition(hitOnSurfaceGlobal);
+ simdata.setTime(m_globalHitTime);
+ h_sdoContainer->insert ( std::make_pair ( id, simdata ) );
+
+ ATH_MSG_VERBOSE(" Prd: local x " << posOnSurf.x() << " y " << 0 );
+ ATH_MSG_VERBOSE(" Prd: r " << prd->globalPosition().perp() << " phi " << prd->globalPosition().phi() << " z " << prd->globalPosition().z());
+ } else {
+ for (int loop_direction = -1; loop_direction <=1; loop_direction+=2) {
+ Amg::Vector3D stepInDriftGap = loop_direction * ldir * (roParam.stripPitch/std::cos(roParam.stereoAngel.at(m_idHelper->gasGap(layid)-1) ))/abs(ldir.x());
+ if (loop_direction == 1) CurrentHitInDriftGap = slpos + stepInDriftGap;
+ while (std::abs(CurrentHitInDriftGap.z()) <= roParam.gasThickness) {
+ Amg::MatrixX* cov = new Amg::MatrixX(1,1);
+ cov->setIdentity();
+ (*cov)(0,0) = resolution*resolution;
+
+ tdrift = CurrentHitInDriftGap.z() / vdrift + CLHEP::RandGauss::shoot(m_rndmEngine, 0., 5.);
+ Amg::Vector2D CurrenPosOnSurf(CurrentHitInDriftGap.x(),CurrentHitInDriftGap.y());
+
+ stripNumber = detEl->stripNumber(CurrenPosOnSurf,layid);
+ if( (stripNumber >= detEl->numberOfStrips(layid)) || (stripNumber == -1) ) {
+ CurrentHitInDriftGap += stepInDriftGap;
+ continue;
+ }
+ id = m_idHelper->channelID(parentID, m_idHelper->multilayer(layid), m_idHelper->gasGap(layid),stripNumber,m_checkIds);
+ if( stripNumber != m_idHelper->channel(id) ) {
+ CurrentHitInDriftGap += stepInDriftGap;
+ continue;
+ }
+
+ MMPrepData* prd = new MMPrepData( id,hash,Amg::Vector2D(CurrenPosOnSurf.x(),0.),rdoList,cov,detEl,(int)tdrift,0);
+ prd->setHashAndIndex(col->identifyHash(), col->size()); // <<< add this line to the MM code as well
+ col->push_back(prd);
+
+ std::vector<MuonSimData::Deposit> deposits;
+ MuonSimData::Deposit deposit(hit.particleLink(), MuonMCData(CurrenPosOnSurf.x(),CurrenPosOnSurf.y()));
+ deposits.push_back(deposit);
+
+ MuonSimData simdata(deposits,0);
+ Amg::Vector3D SDO_GP = surf.transform()*CurrentHitInDriftGap;
+ simdata.setPosition(SDO_GP);
+ simdata.setTime(m_globalHitTime);
+ h_sdoContainer->insert ( std::make_pair ( id, simdata ) );
+
+ ATH_MSG_VERBOSE(" Local CurrentHitInDriftGap.x() " << CurrentHitInDriftGap.x() << " CurrentHitInDriftGap.y() " << CurrentHitInDriftGap.y() << " CurrentHitInDriftGap.z() " << CurrentHitInDriftGap.z() << " drift time " << (int)tdrift );
+ ATH_MSG_VERBOSE(" Prd: local x " << CurrentHitInDriftGap.x() << " y " << 0 << " drift time " << (int)tdrift << " identifier " << id );
+ ATH_MSG_VERBOSE(" Prd: r " << prd->globalPosition().perp() << " phi " << prd->globalPosition().phi() << " z " << prd->globalPosition().z());
+ ATH_MSG_VERBOSE(" SDO: True Global position: x " << SDO_GP.x() << " y " << SDO_GP.y() << " z " << SDO_GP.z());
+
+ CurrentHitInDriftGap += stepInDriftGap;
+ }
+ }
+ }
+
+
m_err = resolution;
m_ich = m_idHelper->channel(id);
m_istr = stripNumber;
@@ -526,15 +590,6 @@ StatusCode MM_FastDigitizer::execute() {
// }
m_ntuple->Fill();
- // create SDO
- MuonSimData::Deposit deposit(hit.particleLink(), MuonMCData(hitOnSurface.x(),hitOnSurface.y()));
- //Record the SDO collection in StoreGate
- std::vector<MuonSimData::Deposit> deposits;
- deposits.push_back(deposit);
- MuonSimData simData(deposits, 0);
- simData.setPosition(hitOnSurfaceGlobal);
- simData.setTime(m_globalHitTime);
- h_sdoContainer->insert ( std::make_pair ( id, simData ) );
// OLD CODE ENDS HERE
previousHit = &hit;
diff --git a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.h b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.h
index 01837dbc834..d75c16fad84 100644
--- a/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.h
+++ b/MuonSpectrometer/MuonDigitization/MuonFastDigitization/src/MM_FastDigitizer.h
@@ -133,6 +133,7 @@ class MM_FastDigitizer : public AthAlgorithm {
double m_energyThreshold;
bool m_checkIds;
int m_maskMultiplet;
+ bool m_microTPC;
};
#endif // MUONDIGITIZATION_MM_DIGITIZER_H
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRIO_OnTrackCreators/MuonClusterOnTrackCreator/src/MuonClusterOnTrackCreator.cxx b/MuonSpectrometer/MuonReconstruction/MuonRIO_OnTrackCreators/MuonClusterOnTrackCreator/src/MuonClusterOnTrackCreator.cxx
index 09992546baa..51cd718907d 100755
--- a/MuonSpectrometer/MuonReconstruction/MuonRIO_OnTrackCreators/MuonClusterOnTrackCreator/src/MuonClusterOnTrackCreator.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonRIO_OnTrackCreators/MuonClusterOnTrackCreator/src/MuonClusterOnTrackCreator.cxx
@@ -69,6 +69,7 @@ namespace Muon {
declareProperty("FixedErrorTgcPhi", m_fixedErrorTgcPhi = 5. );
declareProperty("FixedErrorRpcPhi", m_fixedErrorRpcPhi = 5. );
declareProperty("FixedErrorCscPhi", m_fixedErrorCscPhi = 5. );
+ declareProperty("eDriftVelocityMM", m_eDriftVelocityMM = 0.047 );
}
@@ -113,7 +114,7 @@ namespace Muon {
const MuonClusterOnTrack* MuonClusterOnTrackCreator::createRIO_OnTrack(const Trk::PrepRawData& RIO,
- const Amg::Vector3D& GP) const
+ const Amg::Vector3D& GP, const Amg::Vector3D& GD) const
{
MuonClusterOnTrack* MClT = 0;
@@ -320,13 +321,52 @@ namespace Muon {
MClT = new CscClusterOnTrack(MClus,locpar,loce,positionAlongStrip,MClus->status(),MClus->timeStatus());
}else if( m_idHelper->isMM(RIO.identify()) ){
- // cast to CscPrepData
+ // cast to MMPrepData
const MMPrepData* MClus = dynamic_cast<const MMPrepData*> (&RIO);
if (!MClus) {
ATH_MSG_WARNING ( "RIO not of type MMPrepData, cannot create ROT" );
return 0;
}
-
+
+ Amg::Vector2D lposMM = MClus->localPosition();
+
+ double tOffset = 0.0;
+ double zshift = m_eDriftVelocityMM * (MClus->time() - tOffset);
+
+ ATH_MSG_DEBUG( " Drift time " << MClus->time() << " zshift " << zshift);
+ if (zshift<-5.05) {
+ ATH_MSG_WARNING( "Invalid drift time: zshift " << zshift << " while volume size 10: [-5;5]. Set zshift = -5.05");
+ zshift = -5.05;
+ } else if (zshift>5.) {
+ ATH_MSG_WARNING( "Invalid drift time: zshift " << zshift << " while volume size 10: [-5;5]. Set zshift = 5.05");
+ zshift = 5.05;
+ }
+ ATH_MSG_DEBUG( " MM Prepdata pos locX " << lposMM[0] << " locY " << lposMM[1] << " and zshift " << zshift);
+
+ Amg::Vector3D lposTrack = RIO.detectorElement()->surface(RIO.identify()).transform().inverse()*GP;
+ Amg::Vector3D GPinGasGap = GP + Amg::Vector3D(0.,0.,zshift);
+ Amg::Vector3D GPshifted = GPinGasGap;
+ if (std::abs(GD.z()) > 0.000001)
+ GPshifted-= zshift*GD/GD.z();
+ else
+ ATH_MSG_WARNING ( "GD.z() is 0, mTPC correction cannot be done" );;
+
+ ATH_MSG_DEBUG( " GD x " << GD.x() << " y " << GD.y() << " z " << GD.z());
+ ATH_MSG_DEBUG( " GP x " << GP.x() << " y " << GP.y() << " z " << GP.z());
+ ATH_MSG_DEBUG( " GPshifted x " << GPshifted.x() << " y " << GPshifted.y() << " z " << GPshifted.z());
+
+ ATH_MSG_DEBUG( " lposTrack x " << lposTrack.x() << " y " << lposTrack.y() << " z " << lposTrack.z());
+ Amg::Vector3D lposTrackShifted = RIO.detectorElement()->surface(RIO.identify()).transform().inverse()*GPshifted;
+ ATH_MSG_DEBUG( " lposTrackShifted x " << lposTrackShifted.x() << " y " << lposTrackShifted.y() << " z " << lposTrackShifted.z());
+ double locXshift = lposTrackShifted.x() - lposTrack.x();
+ double locYshift = lposTrackShifted.y() - lposTrack.y();
+ ATH_MSG_DEBUG( " MM ClusterOnTrack locXshift " << locXshift);
+
+ ATH_MSG_DEBUG( " MM ClusterOnTrack old pos locX " << locpar[Trk::locX] << " Old position along strip " << positionAlongStrip);
+ positionAlongStrip += locYshift;
+ locpar[Trk::locX] = lposMM[0] + locXshift;
+ ATH_MSG_DEBUG( " MM ClusterOnTrack new pos locX " << locpar[Trk::locX] << " New position along strip " << positionAlongStrip);
+
MClT = new MMClusterOnTrack(MClus,locpar,loce,positionAlongStrip);
}else if( m_idHelper->issTgc(RIO.identify()) ){
@@ -345,13 +385,13 @@ namespace Muon {
}
const MuonClusterOnTrack* MuonClusterOnTrackCreator::
- createRIO_OnTrack(const Trk::PrepRawData& RIO, const Amg::Vector3D& GP,
- const Amg::Vector3D&) const {
- return createRIO_OnTrack(RIO,GP);
+ createRIO_OnTrack(const Trk::PrepRawData& RIO, const Amg::Vector3D& GP) const {
+ const Amg::Vector3D GD = GP;
+ return createRIO_OnTrack(RIO,GP,GD);
}
const MuonClusterOnTrack* MuonClusterOnTrackCreator::correct(const Trk::PrepRawData& RIO,const Trk::TrackParameters& TP) const
{
- return createRIO_OnTrack(RIO,TP.position());
+ return createRIO_OnTrack(RIO,TP.position(),TP.momentum());
}
}
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRIO_OnTrackCreators/MuonClusterOnTrackCreator/src/MuonClusterOnTrackCreator.h b/MuonSpectrometer/MuonReconstruction/MuonRIO_OnTrackCreators/MuonClusterOnTrackCreator/src/MuonClusterOnTrackCreator.h
index 9188c2db3c1..f047899d614 100755
--- a/MuonSpectrometer/MuonReconstruction/MuonRIO_OnTrackCreators/MuonClusterOnTrackCreator/src/MuonClusterOnTrackCreator.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonRIO_OnTrackCreators/MuonClusterOnTrackCreator/src/MuonClusterOnTrackCreator.h
@@ -108,6 +108,8 @@ namespace Muon {
double m_fixedErrorTgcPhi;
double m_fixedErrorRpcPhi;
double m_fixedErrorCscPhi;
+
+ double m_eDriftVelocityMM;
};
}
#endif // MuonClusterOnTrackCreator_H
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonTrackFindingEvent/MuonTrackFindingEvent/MuPatHit.h b/MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonTrackFindingEvent/MuonTrackFindingEvent/MuPatHit.h
index 1cad2f89aac..7ad98733021 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonTrackFindingEvent/MuonTrackFindingEvent/MuPatHit.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonTrackFindingEvent/MuonTrackFindingEvent/MuPatHit.h
@@ -19,7 +19,7 @@ namespace Muon {
class MuPatHit {
public:
- enum Type { UnknownType = -1, MDT = 0, RPC = 1, TGC = 2, CSC = 3, PREC = 4, Pseudo = 5, Scatterer = 6 };
+ enum Type { UnknownType = -1, MDT = 0, RPC = 1, TGC = 2, CSC = 3, MM = 4, sTGC = 5, PREC = 6, Pseudo = 7, Scatterer = 8 };
enum Status { UnknownStatus = -1, OnTrack = 0, Outlier, NotOnTrack };
enum HitSelection { UnknownSelection = -1, Precise = 0, Broad = 1 };
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/MuonClusterSegmentFinderTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/MuonClusterSegmentFinderTool.cxx
index f7eb4d7dc60..19d277c9b50 100755
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/MuonClusterSegmentFinderTool.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/MuonClusterSegmentFinderTool.cxx
@@ -8,6 +8,7 @@
#include "MuonReadoutGeometry/sTgcReadoutElement.h"
#include "MuonReadoutGeometry/MuonPadDesign.h"
#include "MuonPrepRawData/sTgcPrepData.h"
+#include "MuonPrepRawData/MMPrepData.h"
#include "TrkTrack/Track.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkEventPrimitives/FitQuality.h"
@@ -60,7 +61,7 @@ namespace Muon {
std::vector<const Muon::MuonSegment*>* MuonClusterSegmentFinderTool::find(std::vector< const Muon::MuonClusterOnTrack* >& muonClusters) {
ATH_MSG_DEBUG("Entering MuonClusterSegmentFinderTool with " << muonClusters.size() << " clusters to be fit" );
- if(muonClusters.size() < 4) return 0;
+ if(belowThreshold(muonClusters,4)) return 0;
std::vector<const Muon::MuonSegment*>* segs = findPrecisionSegments(muonClusters);
if(segs == 0) return 0;
if(segs->size() == 0)
@@ -85,7 +86,7 @@ namespace Muon {
bool selectPhiHits(false);
std::vector< const Muon::MuonClusterOnTrack* > clusters = cleanClusters(muonClusters,selectPhiHits);
ATH_MSG_VERBOSE("After hit cleaning, there are " << clusters.size() << " 2D clusters to be fit" );
- if(clusters.size() < 4){
+ if(belowThreshold(clusters,4)){
return 0;
}
@@ -101,9 +102,9 @@ namespace Muon {
std::vector< const Muon::MuonClusterOnTrack* > rioVecPrevious;
//find all clusters near the seed and try to fit
for(unsigned int i=0; i<seeds.size(); ++i) {
- std::vector< const Muon::MuonClusterOnTrack* > rioVec = getClustersOnSegment(orderedClusters,seeds[i]);
+ std::vector< const Muon::MuonClusterOnTrack* > rioVec = getClustersOnSegment(orderedClusters,seeds[i],false);
// make consistent cut
- if(rioVec.size() <4) continue;
+ if(belowThreshold(rioVec,4)) continue;
// logic to reduce combinatorics
if(rioVec.size() == rioVecPrevious.size()) {
bool sameContent = true;
@@ -264,8 +265,9 @@ namespace Muon {
std::vector<const Muon::MuonSegment*>* etaSegs) const {
bool selectPhiHits(true);
std::vector< const Muon::MuonClusterOnTrack* > clusters = cleanClusters(muonClusters,selectPhiHits);
+ std::vector< const Muon::MuonClusterOnTrack* > etaClusters = cleanClusters(muonClusters,false);
ATH_MSG_DEBUG("After hit cleaning, there are " << clusters.size() << " 3D clusters to be fit" );
- if(clusters.size() < 4) {
+ if(belowThreshold(clusters,4)) {
ATH_MSG_DEBUG("Not enough phi hits present, cannot perform the fit!");
return etaSegs;
}
@@ -273,7 +275,8 @@ namespace Muon {
TrackCollection* segTrkColl = new TrackCollection;
//order the clusters by layer
bool useWires(true);
- std::vector< std::vector<const Muon::MuonClusterOnTrack*> > orderedClusters = orderByLayer(clusters,useWires);
+ std::vector< std::vector<const Muon::MuonClusterOnTrack*> > orderedClusters = orderByLayer(clusters,useWires);
+ std::vector< std::vector<const Muon::MuonClusterOnTrack*> > orderedEtaClusters = orderByLayer(etaClusters,false);
if(orderedClusters.size() == 0) {
ATH_MSG_DEBUG( "No phi wire hits found ... moving to pads" );
useWires = false;
@@ -326,7 +329,8 @@ namespace Muon {
seed3D = seeds[i];
}
- std::vector< const Muon::MuonClusterOnTrack* > phiHits = getClustersOnSegment(orderedClusters,seed3D);
+ std::vector< const Muon::MuonClusterOnTrack* > phiHits = getClustersOnSegment(orderedClusters,seed3D,true);
+ std::vector< const Muon::MuonClusterOnTrack* > etaHitsRedone = getClustersOnSegment(orderedEtaClusters,seed3D,true);
if(phiHits.size() < 2) {
delete startpar;
continue;
@@ -352,8 +356,15 @@ namespace Muon {
//interleave the phi hits
std::vector<const Trk::MeasurementBase*> vec2;
const std::vector<const Trk::RIO_OnTrack*> etaHits = (*sit)->containedROTs();
- if(m_ipConstraint) vec2.reserve(phiHits.size()+etaHits.size()+1);
- else vec2.reserve(phiHits.size()+etaHits.size());
+ bool useEtaHitsRedone = false;
+ if(etaHitsRedone.size()>etaHits.size()) {
+ ATH_MSG_VERBOSE(" Found additional eta hits " << etaHitsRedone.size() - etaHits.size());
+ useEtaHitsRedone = true;
+ }
+ unsigned int netas = etaHits.size();
+ if(useEtaHitsRedone) netas = etaHitsRedone.size();
+ if(m_ipConstraint) vec2.reserve(phiHits.size()+netas+1);
+ else vec2.reserve(phiHits.size()+netas);
// pseudo measurement for vtx
Trk::PseudoMeasurementOnTrack* pseudoVtx = nullptr;
@@ -378,10 +389,14 @@ namespace Muon {
iPhi++;
}
else {
- vec2.push_back(etaHits[iEta]);
+ if(!useEtaHitsRedone) {
+ vec2.push_back(etaHits[iEta]);
+ } else {
+ vec2.push_back(etaHitsRedone[iEta]);
+ }
iEta++;
}
- if( iEta >= etaHits.size() && iPhi >= phiHits.size() ) break;
+ if( iEta >= netas && iPhi >= phiHits.size() ) break;
}
@@ -584,30 +599,62 @@ namespace Muon {
}
std::vector< const Muon::MuonClusterOnTrack* > MuonClusterSegmentFinderTool::getClustersOnSegment(std::vector< std::vector<const Muon::MuonClusterOnTrack*> >& clusters,
- std::pair<Amg::Vector3D,Amg::Vector3D>& seed) const {
+ std::pair<Amg::Vector3D,Amg::Vector3D>& seed, bool tight) const {
ATH_MSG_VERBOSE(" getClustersOnSegment: layers " << clusters.size() );
std::vector< const Muon::MuonClusterOnTrack* > rios;
int layer = 0;
for(std::vector<std::vector<const Muon::MuonClusterOnTrack*> >::const_iterator cvecIt=clusters.begin(); cvecIt!=clusters.end(); ++cvecIt) {
const Muon::MuonClusterOnTrack* rio = 0;
double bestDist(9999.);
- for(std::vector< const Muon::MuonClusterOnTrack* >::const_iterator cit=(*cvecIt).begin(); cit!=(*cvecIt).end(); ++cit) {
- double dist = clusterDistanceToSeed( *cit, seed);
- double error = Amg::error((*cit)->localCovariance(),Trk::locX);
- if(m_idHelperTool->isMM((*cit)->identify())) error += 25.;
- ATH_MSG_VERBOSE(" lay " << layer << " dist " << dist << " pull " << dist/error
+ double bestTimeDist(9999.);
+ if (m_idHelperTool->isMM((*((*cvecIt).begin()))->identify())) { // MM specific algorithm
+ for(std::vector< const Muon::MuonClusterOnTrack* >::const_iterator cit=(*cvecIt).begin(); cit!=(*cvecIt).end(); ++cit) {
+ double tdrift = (dynamic_cast<const MMPrepData *>((*cit)->prepRawData()))->time();
+ double dist = clusterDistanceToSeed( *cit, seed);
+ double timedist = std::abs(clusterDistanceToSeed( *cit, seed)) + std::abs(tdrift*0.015); // std::abs(tdrift*0.015) is an ad hoc penalty factor, to be optimised when time resolution is known
+ double error = Amg::error((*cit)->localCovariance(),Trk::locX);
+ if (!tight) error += 1.;
+ ATH_MSG_VERBOSE(" lay " << layer << " tdrift " << tdrift << " dist " << dist << " timedist " << timedist << " pull " << dist/error
<< " cut " << m_maxClustDist << " " << m_idHelperTool->toString((*cit)->identify()) );
- if( fabs(dist/error) < m_maxClustDist) {
- if(fabs(dist) < bestDist) {
- bestDist = fabs(dist);
- rio = (*cit);
- }
- }
- }
- if(rio) {
- ATH_MSG_VERBOSE(" adding " << bestDist << " " << m_idHelperTool->toString(rio->identify()) );
- rios.push_back( rio );
- }else{
+ if( std::abs(dist/error) < m_maxClustDist) {
+ if(std::abs(timedist) < bestTimeDist) {
+ bestTimeDist = std::abs(timedist);
+ bestDist = dist;
+ }
+ }
+ }
+ if(bestDist<9999.) { // check if at least one cluster present close to seed
+ ATH_MSG_VERBOSE(" Best distance " << bestDist);
+ ATH_MSG_VERBOSE(" Looking for RIOs from mTPC around the best distance");
+ for(std::vector< const Muon::MuonClusterOnTrack* >::const_iterator cit=(*cvecIt).begin(); cit!=(*cvecIt).end(); ++cit) {
+ double dist = clusterDistanceToSeed( *cit, seed);
+ double window = std::abs(2.*5.*0.047 * ((*cit)->globalPosition().perp() / (*cit)->globalPosition().z())); // all hits in the range [bestDist-window;bestDist-window] will be accepted; 2-safety factor; 5-time resolution; 0.047-drift velocity; (hardcoded values to be removed once time resolution model is known)
+ double error = Amg::error((*cit)->localCovariance(),Trk::locX);
+ if (!tight) error += 1.;
+ ATH_MSG_VERBOSE(" Current RIO : distance " << dist << " window " << window << " to be attached " << ( (std::abs(std::abs(dist)-bestDist) < window) && (std::abs(dist/error) < m_maxClustDist) ) );
+ if( (std::abs(dist-bestDist) < window) && (std::abs(dist/error) < m_maxClustDist) ) {
+ rios.push_back( (*cit) );
+ ATH_MSG_VERBOSE(" adding " << dist << " " << m_idHelperTool->toString((*cit)->identify()) );
+ }
+ }
+ }
+ } else { // algorithm for all the technoligies but MM
+ for(std::vector< const Muon::MuonClusterOnTrack* >::const_iterator cit=(*cvecIt).begin(); cit!=(*cvecIt).end(); ++cit) {
+ double dist = clusterDistanceToSeed( *cit, seed);
+ double error = Amg::error((*cit)->localCovariance(),Trk::locX);
+ ATH_MSG_VERBOSE(" lay " << layer << " dist " << dist << " pull " << dist/error
+ << " cut " << m_maxClustDist << " " << m_idHelperTool->toString((*cit)->identify()) );
+ if( std::abs(dist/error) < m_maxClustDist) {
+ if(std::abs(dist) < bestDist) {
+ bestDist = std::abs(dist);
+ rio = (*cit);
+ }
+ }
+ }
+ if(rio) {
+ ATH_MSG_VERBOSE(" adding " << bestDist << " " << m_idHelperTool->toString(rio->identify()) );
+ rios.push_back( rio );
+ }
}
++layer;
}
@@ -866,5 +913,20 @@ namespace Muon {
return phiOverlap;
}
+ bool MuonClusterSegmentFinderTool::belowThreshold(std::vector< const Muon::MuonClusterOnTrack* >& muonClusters, int threshold) const {
+ int n_surf_with_hits = 0;
+ if (muonClusters.size()>1) {
+ for(std::vector< const Muon::MuonClusterOnTrack* >::const_iterator cit=muonClusters.begin()+1; cit!=muonClusters.end(); cit++){
+ if ((*cit) && (*(cit-1))) {
+ if ((*cit)->detectorElement()->center((*cit)->identify()) != (*(cit-1))->detectorElement()->center((*(cit-1))->identify())) {
+ n_surf_with_hits++;
+ }
+ }
+ }
+ } else
+ n_surf_with_hits = muonClusters.size();
+ if(n_surf_with_hits < threshold) return true;
+ else return false;
+ }
}//end namespace
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/MuonClusterSegmentFinderTool.h b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/MuonClusterSegmentFinderTool.h
index e53b30c7f55..289d5f2fa7b 100755
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/MuonClusterSegmentFinderTool.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/MuonClusterSegmentFinderTool.h
@@ -77,11 +77,14 @@ namespace Muon {
std::vector<std::pair<double,double> > getPadPhiOverlap(std::vector< std::vector<const Muon::MuonClusterOnTrack*> >& pads) const;
//associate clusters to the segment seeds
std::vector< const Muon::MuonClusterOnTrack* > getClustersOnSegment(std::vector< std::vector<const Muon::MuonClusterOnTrack*> >& clusters,
- std::pair<Amg::Vector3D,Amg::Vector3D>& seed) const;
+ std::pair<Amg::Vector3D,Amg::Vector3D>& seed, bool tight) const;
//distance of cluster to segment seed
double clusterDistanceToSeed(const Muon::MuonClusterOnTrack* clust, std::pair<Amg::Vector3D,Amg::Vector3D>& seed) const;
Amg::Vector3D intersectPlane( const Trk::PlaneSurface& surf, const Amg::Vector3D& pos, const Amg::Vector3D& dir ) const;
Trk::Track* fit( const std::vector<const Trk::MeasurementBase*>& vec2, const Trk::TrackParameters& startpar ) const;
+
+ //check if enough surfaces are hitted
+ bool belowThreshold(std::vector< const Muon::MuonClusterOnTrack* >& muonClusters, int threshold) const;
};
diff --git a/MuonSpectrometer/MuonReconstruction/MuonTrackMakers/MuonTrackMakerTools/MuonTrackSteeringTools/src/MuPatHitTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonTrackMakers/MuonTrackMakerTools/MuonTrackSteeringTools/src/MuPatHitTool.cxx
index e08230638c1..960ce343624 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonTrackMakers/MuonTrackMakerTools/MuonTrackSteeringTools/src/MuPatHitTool.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonTrackMakers/MuonTrackMakerTools/MuonTrackSteeringTools/src/MuPatHitTool.cxx
@@ -400,6 +400,8 @@ namespace Muon {
else if( m_idHelperTool->isTgc(id) ) return MuPatHit::TGC;
else if( m_idHelperTool->isCsc(id) ) return MuPatHit::CSC;
else if( m_idHelperTool->isRpc(id) ) return MuPatHit::RPC;
+ else if( m_idHelperTool->isMM(id) ) return MuPatHit::MM;
+ else if( m_idHelperTool->issTgc(id) ) return MuPatHit::sTGC;
else if( m_idHelperTool->isMuon(id) ) return MuPatHit::PREC;
return MuPatHit::UnknownType;
}
@@ -414,7 +416,7 @@ namespace Muon {
if( hitInfo.id.is_valid() && m_idHelperTool->isMuon(hitInfo.id) ) {
hitInfo.type = getHitType(hitInfo.id);
hitInfo.measuresPhi = m_idHelperTool->measuresPhi(hitInfo.id);
- if( hitInfo.type != MuPatHit::MDT ) hitInfo.id = m_idHelperTool->layerId(hitInfo.id);
+ if( hitInfo.type != MuPatHit::MDT && hitInfo.type != MuPatHit::MM ) hitInfo.id = m_idHelperTool->layerId(hitInfo.id);
}
}
@@ -479,7 +481,7 @@ namespace Muon {
// if we reached the end of the list, insert the hit at the end
if( pos == list.end() ) {
// check whether hit duplicate of last hit in list
- if( isLargerCal(list.back(),hit) != isLargerCal(hit,list.back()) ){
+ if( isLargerCal(list.back(),hit) != isLargerCal(hit,list.back()) || (hit->info().type == MuPatHit::MM && hit->info().id != list.back()->info().id) ){
ATH_MSG_VERBOSE(" inserting hit at back " << m_idHelperTool->toString(hit->info().id)
<< " " << m_printer->print(hit->parameters()) );
list.push_back(hit);
@@ -502,7 +504,7 @@ namespace Muon {
// if we reached the first list item, check whether current hit is smaller. If so insert before first.
if( pos == list.begin() && !isLarger ){
// check whether hit duplicate of last hit in list
- if( isLargerCal(list.front(),hit) != isLargerCal(hit,list.front()) ){
+ if( isLargerCal(list.front(),hit) != isLargerCal(hit,list.front()) || (hit->info().type == MuPatHit::MM && hit->info().id != list.front()->info().id) ){
ATH_MSG_VERBOSE(" inserting hit at front " << m_idHelperTool->toString(hit->info().id)
<< " " << m_printer->print(hit->parameters()) );
@@ -521,7 +523,7 @@ namespace Muon {
++pos;
if( pos == list.end() ) {
// check whether hit duplicate of last hit in list
- if( isLargerCal(list.back(),hit) != isLargerCal(hit,list.back()) ){
+ if( isLargerCal(list.back(),hit) != isLargerCal(hit,list.back()) || (hit->info().type == MuPatHit::MM && hit->info().id != list.back()->info().id) ){
ATH_MSG_VERBOSE(" inserting hit at back " << m_idHelperTool->toString(hit->info().id)
<< " " << m_printer->print(hit->parameters()) );
list.push_back(hit);
@@ -539,7 +541,7 @@ namespace Muon {
// remove duplicates
// check whether hit and entry at pos are a duplicate
- if( isLarger == isLargerCal(*pos,hit) ){
+ if( isLarger == isLargerCal(*pos,hit) && (hit->info().type != MuPatHit::MM || hit->info().id == (*pos)->info().id) ){
// hit is a duplicate
ATH_MSG_VERBOSE(" NOT inserting duplicate hit " << m_idHelperTool->toString(hit->info().id)
<< " " << m_printer->print(hit->parameters()) );
@@ -551,7 +553,7 @@ namespace Muon {
--pos; // move to previous hit
// check whether hit and entry at pos are a duplicate
- if( isLargerCal(hit,*pos) == isLargerCal(*pos,hit) ){
+ if( isLargerCal(hit,*pos) == isLargerCal(*pos,hit) && (hit->info().type != MuPatHit::MM || hit->info().id == (*pos)->info().id) ){
++pos; // move forward to insert position for pos
// hit is a duplicate
ATH_MSG_VERBOSE(" NOT inserting duplicate hit " << m_idHelperTool->toString(hit->info().id)
diff --git a/Reconstruction/MuonIdentification/MuidTrackBuilder/src/CombinedMuonTrackBuilder.cxx b/Reconstruction/MuonIdentification/MuidTrackBuilder/src/CombinedMuonTrackBuilder.cxx
index 5a8dfd45d4f..f3d786d251f 100755
--- a/Reconstruction/MuonIdentification/MuidTrackBuilder/src/CombinedMuonTrackBuilder.cxx
+++ b/Reconstruction/MuonIdentification/MuidTrackBuilder/src/CombinedMuonTrackBuilder.cxx
@@ -3255,7 +3255,8 @@ CombinedMuonTrackBuilder::appendSelectedTSOS(
if (previousSurface
&& std::find(measurementSurfaces.begin(),
measurementSurfaces.end(),
- surface) != measurementSurfaces.end())
+ surface) != measurementSurfaces.end()
+ && !m_idHelperTool->isMM(surface->associatedDetectorElementIdentifier()))
{
std::string type = "";
if (dynamic_cast<const Trk::CompetingRIOsOnTrack*>((**s).measurementOnTrack()))
@@ -4088,7 +4089,8 @@ CombinedMuonTrackBuilder::createSpectrometerTSOS(const Trk::Track& spectrometerT
if (previousSurface
&& std::find(measurementSurfaces.begin(),
measurementSurfaces.end(),
- surface) != measurementSurfaces.end())
+ surface) != measurementSurfaces.end()
+ && !m_idHelperTool->isMM(surface->associatedDetectorElementIdentifier()))
{
std::string type = "";
if (dynamic_cast<const Trk::CompetingRIOsOnTrack*>((**s).measurementOnTrack()))
diff --git a/Tracking/TrkFitter/TrkGlobalChi2Fitter/TrkGlobalChi2Fitter/GlobalChi2Fitter.h b/Tracking/TrkFitter/TrkGlobalChi2Fitter/TrkGlobalChi2Fitter/GlobalChi2Fitter.h
index ee2674d2bd9..5248b04db46 100755
--- a/Tracking/TrkFitter/TrkGlobalChi2Fitter/TrkGlobalChi2Fitter/GlobalChi2Fitter.h
+++ b/Tracking/TrkFitter/TrkGlobalChi2Fitter/TrkGlobalChi2Fitter/GlobalChi2Fitter.h
@@ -293,6 +293,7 @@ private:
mutable std::vector<const Trk::Layer*> m_barrelcylinders;
mutable std::vector<const Trk::Layer*> m_othercylinders;
mutable bool m_fastmat;
+ mutable int m_MMCorrectionStatus;
};
//std::vector<CLHEP::HepMatrix> Trk::GlobalChi2Fitter::m_derivpool;
diff --git a/Tracking/TrkFitter/TrkGlobalChi2Fitter/src/GXFTrajectory.cxx b/Tracking/TrkFitter/TrkGlobalChi2Fitter/src/GXFTrajectory.cxx
index dd2e8609d3f..242693cd8f7 100644
--- a/Tracking/TrkFitter/TrkGlobalChi2Fitter/src/GXFTrajectory.cxx
+++ b/Tracking/TrkFitter/TrkGlobalChi2Fitter/src/GXFTrajectory.cxx
@@ -173,7 +173,8 @@ namespace Trk {
const MeasurementBase *meas = m_states.back()->measurement();
const MeasurementBase *meas2 = state->measurement();
if (&meas->associatedSurface() == &meas2->associatedSurface() &&
- meas->localParameters().parameterKey() == meas2->localParameters().parameterKey()) {
+ meas->localParameters().parameterKey() == meas2->localParameters().parameterKey() &&
+ state->measurementType() != TrackState::MM) {
delete state;
return false;
}
diff --git a/Tracking/TrkFitter/TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx b/Tracking/TrkFitter/TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx
index 9b1cefe2495..4d19e16425c 100644
--- a/Tracking/TrkFitter/TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx
+++ b/Tracking/TrkFitter/TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx
@@ -143,7 +143,8 @@ namespace Trk {
m_updatescat{},
m_useCaloTG(false),
m_caloMaterialProvider("Trk::TrkMaterialProviderTool/TrkMaterialProviderTool"),
- m_rejectLargeNScat(false){
+ m_rejectLargeNScat(false),
+ m_MMCorrectionStatus(0){
// tools and services
declareProperty("ExtrapolationTool", m_extrapolator);
declareProperty("MeasurementUpdateTool", m_updator);
@@ -2238,7 +2239,48 @@ namespace Trk {
MeasurementSet::const_iterator itSet = rots.begin();
MeasurementSet::const_iterator itSetEnd = rots.end();
GXFTrajectory trajectory;
- for (; itSet != itSetEnd; ++itSet) {
+
+ if (m_MMCorrectionStatus==0) {
+ for (itSet = rots.begin(); itSet != itSetEnd; ++itSet) {
+ if ((*itSet)) {
+ if (m_DetID->is_mm((*itSet)->associatedSurface().associatedDetectorElementIdentifier())) {
+ m_MMCorrectionStatus = 1;
+ break;
+ }
+ }
+ }
+ }
+
+ if (m_MMCorrectionStatus==1) {
+ itSet = rots.begin();
+ MeasurementSet rots_new;
+ MeasurementSet rots_tbd;
+ for (itSet = rots.begin(); itSet != itSetEnd; ++itSet) {
+ if (!(*itSet)) {
+ msg(MSG::WARNING) << "There is an empty MeasurementBase object in the track! Skip this object.." << endmsg;
+ } else {
+ const RIO_OnTrack *rot = dynamic_cast<const RIO_OnTrack *>(*itSet);
+ if (rot && m_DetID->is_mm(rot->identify())) {
+ const PlaneSurface* surf = dynamic_cast<const PlaneSurface *>(&rot->associatedSurface());
+ AtaPlane atapl(surf->center(), param.parameters()[Trk::phi], param.parameters()[Trk::theta], param.parameters()[Trk::qOverP], *surf);
+ rot = m_ROTcreator->correct(*(rot->prepRawData()), atapl);
+ rots_tbd.push_back(rot);
+ rots_new.push_back(rot);
+ } else {
+ rots_new.push_back(*itSet);
+ }
+ }
+ }
+ m_MMCorrectionStatus = -1;
+ Track *track = fit(rots_new, param, runOutlier, matEffects);
+ m_MMCorrectionStatus = 0;
+ for (MeasurementSet::const_iterator it = rots_tbd.begin(); it != rots_tbd.end(); it++) {
+ delete *it;
+ }
+ return track;
+ }
+
+ for (itSet = rots.begin(); itSet != itSetEnd; ++itSet) {
if (!(*itSet)) {
msg(MSG::WARNING) << "There is an empty MeasurementBase object in the track! Skip this object.." << endmsg;
}else {
diff --git a/Tracking/TrkTools/TrkRIO_OnTrackCreator/src/RIO_OnTrackCreator.cxx b/Tracking/TrkTools/TrkRIO_OnTrackCreator/src/RIO_OnTrackCreator.cxx
index f0f9814a584..b67d506645a 100755
--- a/Tracking/TrkTools/TrkRIO_OnTrackCreator/src/RIO_OnTrackCreator.cxx
+++ b/Tracking/TrkTools/TrkRIO_OnTrackCreator/src/RIO_OnTrackCreator.cxx
@@ -213,8 +213,8 @@ Trk::RIO_OnTrackCreator::correct(const Trk::PrepRawData& rio,
return m_MuonDriftCircleCor->correct(rio, trk);
}
}
- if ( (m_idHelper->is_csc(id)) || (m_idHelper->is_rpc(id))
- || (m_idHelper->is_tgc(id)) ) {
+ if ( (m_idHelper->is_csc(id)) || (m_idHelper->is_rpc(id))
+ || (m_idHelper->is_tgc(id)) || (m_idHelper->is_mm(id)) || (m_idHelper->is_stgc(id)) ) {
if (m_mode == "indet") {
msg(MSG::WARNING)<<"I have no tool to correct a CSC/RPC/TGC hit! - Giving back nil."<<endmsg;
return 0;
--
GitLab