Skip to content
Snippets Groups Projects
Commit 1c5608ed authored by Dave Casper's avatar Dave Casper
Browse files

Merge branch 'master-clusterFit-fix' into 'master' 

Clean and fixes of cluster fitter

See merge request faser/calypso!115
parents 190083b2 10e5b599
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 122 additions and 31 deletions
Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx
+ 113
− 25
View file @ 1c5608ed
......@@ -40,6 +40,7 @@ ClusterFitAlg::ClusterFitAlg(const std::string& name, ISvcLocator* pSvcLocator)
// Initialize method:
StatusCode ClusterFitAlg::initialize() {
ATH_MSG_INFO("ClusterFitAlg::initialize()");
ATH_MSG_ALWAYS("zCenter: " << m_zCenter[0] << " / " << m_zCenter[1] << " / " << m_zCenter[2] << " / " << m_zCenter[3]);
ATH_CHECK(m_faserTriggerDataKey.initialize());
ATH_CHECK(m_clusterContainerKey.initialize());
......@@ -76,9 +77,11 @@ StatusCode ClusterFitAlg::execute(const EventContext& ctx) const
SG::ReadHandle<xAOD::FaserTriggerData> triggerContainer{m_faserTriggerDataKey, ctx};
ATH_CHECK(triggerContainer.isValid());
uint32_t eventNumber = triggerContainer->eventId();
uint8_t triggerBits = triggerContainer->tap();
// if ((triggerBits & 0x4) == 0)
if (m_triggerMask != 0 && ((triggerBits & m_triggerMask) == 0))
// if (eventNumber != 10841)
{
ATH_CHECK(trackContainer.record(std::move(outputTracks)));
return StatusCode::SUCCESS;
......@@ -111,22 +114,65 @@ StatusCode ClusterFitAlg::execute(const EventContext& ctx) const
const TrackerDD::SiDetectorElement* elem = (*clusters)->detectorElement();
if (elem != nullptr)
{
ATH_MSG_VERBOSE("( " << m_idHelper->wafer_hash(elem->identify()) << " = " << m_idHelper->layer(elem->identify()) << "/" << m_idHelper->eta_module(elem->identify()) << "/" << m_idHelper->phi_module(elem->identify()) << "|" << m_idHelper->side(elem->identify()) <<
ATH_MSG_VERBOSE("( " << m_idHelper->wafer_hash(elem->identify()) << " = " << m_idHelper->station(elem->identify()) << "/" << m_idHelper->layer(elem->identify()) << "/" << m_idHelper->eta_module(elem->identify()) << "/" << m_idHelper->phi_module(elem->identify()) << "|" << m_idHelper->side(elem->identify()) <<
" ): isStereo? " << elem->isStereo() << " sinStereo: " << elem->sinStereo() << " Swap phi? " << elem->swapPhiReadoutDirection() << " Depth, Eta, Phi directions: " <<
elem->hitDepthDirection() << " : " << elem->hitEtaDirection() << " : " << elem->hitPhiDirection() << " / zGlobal = " << (*clusters)->globalPosition().z());
stations.emplace(m_idHelper->station(elem->identify()));
Identifier id = elem->identify();
stations.emplace(m_idHelper->station(id));
clusterInfo* info = new clusterInfo { };
info->cluster = *clusters;
info->sinAlpha = elem->sinStereo();
info->cosAlpha = sqrt(1.0 - info->sinAlpha * info->sinAlpha);
double alpha = std::abs(asin(elem->sinStereo()));
int layer = m_idHelper->layer(id);
int eta = m_idHelper->eta_module(id); // -1 or +1
int etaIndex = (eta + 1)/2; // 0 or 1
int phi = m_idHelper->phi_module(id); // 0 to 3
int moduleIndex = ((etaIndex + phi)%2 == 1 ? phi : phi + 4); // 0 to 7
switch (moduleIndex)
{
case 0:
case 2:
{
// info->sinAlpha = sin(alpha);
// info->cosAlpha = -cos(alpha);
info->sinAlpha = -sin(alpha);
info->cosAlpha = cos(alpha);
break;
}
case 1:
case 3:
{
info->sinAlpha = -sin(alpha);
info->cosAlpha = cos(alpha);
break;
}
case 4:
case 6:
{
// info->sinAlpha = -sin(alpha);
// info->cosAlpha = -cos(alpha);
info->sinAlpha = sin(alpha);
info->cosAlpha = cos(alpha);
break;
}
case 5:
case 7:
{
info->sinAlpha = sin(alpha);
info->cosAlpha = cos(alpha);
break;
}
}
int side = m_idHelper->side(id);
if (side > 0) info->sinAlpha = - info->sinAlpha; // always reverse angle for back side wafer
info->u = (*clusters)->globalPosition().y() * info->cosAlpha + (*clusters)->globalPosition().x() * info->sinAlpha;
info->z = (*clusters)->globalPosition().z();
ATH_MSG_VERBOSE("Event: " << eventNumber << " Adding cluster " << layer << "/" << eta << "/" << phi << "/" << side << ", u= " << info->u << ", z= " << info->z << ", sin= " << info->sinAlpha << ", cos= " << info->cosAlpha );
info->sigmaSq = (*clusters)->localCovariance()(0,0);
if (info->sinAlpha < 0)
if (info->sinAlpha * info->cosAlpha < 0)
{
clusterMap[hash/2].first.push_back(info);
}
else if (info->sinAlpha > 0)
else if (info->sinAlpha * info->cosAlpha > 0)
{
clusterMap[hash/2].second.push_back(info);
}
......@@ -162,18 +208,21 @@ StatusCode ClusterFitAlg::execute(const EventContext& ctx) const
std::vector<clusterInfo*>& stereoMinus = entry.second.first;
std::vector<clusterInfo*>& stereoPlus = entry.second.second;
ATH_MSG_VERBOSE("Module entry in layer " << layer << " has " << stereoMinus.size() << " negative and " << stereoPlus.size() << " positive stereo clusters");
if (stereoMinus.size() > 0 && stereoPlus.size() > 0 && !layerGood[layer])
if (stereoMinus.size() > 0 && stereoPlus.size() > 0)
{
nGoodLayers++;
layerGood[layer] = true;
if (!layerGood[layer])
{
nGoodLayers++;
layerGood[layer] = true;
}
for (clusterInfo* minus : stereoMinus)
{
for (clusterInfo* plus : stereoPlus)
{
layerCombo* combo = new layerCombo { };
combo->initialize();
combo->addCluster(plus, m_zCenter);
combo->addCluster(minus, m_zCenter);
combo->addCluster(plus, m_zCenter[thisStation]);
combo->addCluster(minus, m_zCenter[thisStation]);
layerCombos[layer].push_back( combo );
}
}
......@@ -218,11 +267,40 @@ StatusCode ClusterFitAlg::execute(const EventContext& ctx) const
}
}
}
ATH_MSG_VERBOSE("best fit: (" << bestFit(0) << "+/-" << sqrt(bestCov(0,0)) << "," <<
ATH_MSG_VERBOSE("Event #" << eventNumber << " best fit: (" << bestFit(0) << "+/-" << sqrt(bestCov(0,0)) << "," <<
bestFit(1) << "+/-" << sqrt(bestCov(1,1)) << "," <<
bestFit(2) << "+/-" << sqrt(bestCov(2,2)) << "," <<
bestFit(3) << "+/-" << sqrt(bestCov(3,3)) << "); chi2 = " << bestChi2 );
Residuals(bestClusters);
if (eventNumber == 10841) std::cout << "fit = ParametricPlot3D[{" << bestFit(0) << " + (" << bestFit(2) << ") z, " << bestFit(1) << " + ( " << bestFit(3) << ") z, z}, {z, -50, 50}, Boxed->False];" << std::endl;
int iC = 0;
for (auto info : bestClusters)
{
if (info->sinAlpha * info->cosAlpha < 0) continue;
double dz = info->z - m_zCenter[thisStation];
double xFit = bestFit(0) + bestFit(2) * dz;
double yFit = bestFit(1) + bestFit(3) * dz;
double uFit = yFit * info->cosAlpha + xFit * info->sinAlpha;
ATH_MSG_VERBOSE("Event: " << eventNumber << " z: " << info->z << " z-zC: " << dz<< " uFit: " << uFit << " u: " << info->u);
auto stripEnds = info->cluster->detectorElement()->endsOfStrip(info->cluster->localPosition());
if (eventNumber == 10841) std::cout << "l" << iC << " = Line[{{" << stripEnds.first.x() << ", " << stripEnds.first.y() << ", " << stripEnds.first.z() - m_zCenter[thisStation] << "}, {" <<
stripEnds.second.x() << ", " << stripEnds.second.y() << ", " << stripEnds.second.z() - m_zCenter[thisStation] << "}}];" << std::endl;
iC++;
}
for (auto info : bestClusters)
{
if (info->sinAlpha * info->cosAlpha > 0) continue;
double dz = info->z - m_zCenter[thisStation];
double xFit = bestFit(0) + bestFit(2) * dz;
double yFit = bestFit(1) + bestFit(3) * dz;
double uFit = yFit * info->cosAlpha + xFit * info->sinAlpha;
ATH_MSG_VERBOSE("Event: " << eventNumber << " z: " << info->z << " z-zC: " << dz<< " uFit: " << uFit << " u: " << info->u);
auto stripEnds = info->cluster->detectorElement()->endsOfStrip(info->cluster->localPosition());
if (eventNumber == 10841) std::cout << "l" << iC << " = Line[{{" << stripEnds.first.x() << ", " << stripEnds.first.y() << ", " << stripEnds.first.z() - m_zCenter[thisStation] << "}, {" <<
stripEnds.second.x() << ", " << stripEnds.second.y() << ", " << stripEnds.second.z() - m_zCenter[thisStation] << "}}];" << std::endl;
iC++;
}
if (eventNumber == 10841) std::cout << "Show[Graphics3D[{l1, l2, l3, l4, l5, l6}, Boxed -> False], fit]" << std::endl;
//Residuals(bestClusters);
ATH_CHECK(AddTrack(outputTracks, bestFit, bestCov, bestClusters, bestChi2, bestClusters.size()-4));
m_chi2->Fill(bestChi2);
setFilterPassed(true, ctx);
......@@ -293,13 +371,15 @@ ClusterFitAlg::AddTrack(std::unique_ptr<TrackCollection>& tracks,
Trk::FitQuality* q = new Trk::FitQuality {chi2, ndof};
DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};
int station = m_idHelper->station(fitClusters[0]->cluster->detectorElement()->identify());
// translate parameters to nominal fit point
s->push_back( GetState(fitResult, fitCovariance, nullptr) );
s->push_back( GetState(fitResult, fitCovariance, nullptr, station) );
for (const clusterInfo* cInfo : fitClusters)
{
s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster) );
}
// for (const clusterInfo* cInfo : fitClusters)
// {
// s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster, station) );
// }
// Create and store track
tracks->push_back(new Trk::Track(i, s , q));
......@@ -309,12 +389,13 @@ ClusterFitAlg::AddTrack(std::unique_ptr<TrackCollection>& tracks,
Trk::TrackStateOnSurface*
ClusterFitAlg::GetState( const Eigen::Matrix< double, 4, 1 >& fitResult,
const Eigen::Matrix< double, 4, 4 >& fitCovariance,
const FaserSCT_Cluster* fitCluster ) const
const FaserSCT_Cluster* fitCluster, int station ) const
{
// position of fit point:
double zFit = m_zCenter;
// int station = m_idHelper->station(fitCluster->detectorElement()->identify());
double zFit = m_zCenter[station];
if (fitCluster != nullptr) zFit = fitCluster->globalPosition()[2];
Amg::Vector3D pos { fitResult[0] + fitResult[2] * (zFit - m_zCenter), fitResult[1] + fitResult[3] * (zFit - m_zCenter), zFit };
Amg::Vector3D pos { fitResult[0] + fitResult[2] * (zFit - m_zCenter[station]), fitResult[1] + fitResult[3] * (zFit - m_zCenter[station]), zFit };
double phi = atan2( fitResult[3], fitResult[2] );
double theta = atan( sqrt(fitResult[2]*fitResult[2] + fitResult[3]*fitResult[3]) );
double qoverp = 1.0/100000.0;
......@@ -352,17 +433,24 @@ ClusterFitAlg::ClusterFit(layerCombo* c0, layerCombo* c1, layerCombo* c2) const
{
sums[i] = c0->sums[i] + c1->sums[i] + c2->sums[i];
}
Amg::MatrixX s(5,5);
Eigen::Matrix< double, 5, 5 > s;
s << sums[1] , sums[2] , sums[4] , sums[5] , sums[10] ,
sums[2] , sums[3] , sums[5] , sums[6] , sums[11] ,
sums[4] , sums[5] , sums[7] , sums[8] , sums[12] ,
sums[5] , sums[6] , sums[8] , sums[9] , sums[13] ,
sums[10] , sums[11] , sums[12] , sums[13] , sums[14];
Amg::MatrixX s4(4,4);
Eigen::Matrix< double, 4, 4 > s4;
s4 = s.block<4,4>(0,0);
Eigen::Matrix< double, 4, 1 > v;
v << sums[10] , sums[11] , sums[12] , sums[13];
Eigen::Matrix< double, 4, 1 > x = s4.colPivHouseholderQr().solve(v);
ATH_MSG_VERBOSE(s4);
ATH_MSG_VERBOSE(v);
ATH_MSG_VERBOSE(s4.inverse());
// Eigen::Matrix< double, 4, 1 > x = s4.colPivHouseholderQr().solve(v);
Eigen::Matrix< double, 4, 1 > x = s4.fullPivLu().solve(v);
ATH_MSG_VERBOSE("Check solution");
ATH_MSG_VERBOSE(s4*x);
ATH_MSG_VERBOSE(v);
Eigen::Matrix< double, 5, 1 > x5;
x5 << x(0), x(1), x(2), x(3), -1;
double chi2 = x5.transpose()*s*x5;
......
Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.h
+ 3
− 2
View file @ 1c5608ed
......@@ -173,7 +173,8 @@ class ClusterFitAlg : public AthReentrantAlgorithm, AthHistogramming
Trk::TrackStateOnSurface* GetState( const Eigen::Matrix< double, 4, 1 >& fitResult,
const Eigen::Matrix< double, 4, 4 >& fitCovariance,
const FaserSCT_Cluster* fitCluster ) const;
const FaserSCT_Cluster* fitCluster,
int station ) const;
std::tuple<Eigen::Matrix<double, 4, 1>,
Eigen::Matrix<double, 4, 4>,
......@@ -189,7 +190,7 @@ class ClusterFitAlg : public AthReentrantAlgorithm, AthHistogramming
SG::ReadHandleKey<FaserSCT_ClusterContainer> m_clusterContainerKey { this, "ClustersName", "SCT_ClusterContainer", "FaserSCT cluster container" };
SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "ClusterFit", "Output track collection name" };
DoubleProperty m_zCenter { this, "ZCenter", 2300.0, "Global z position at which to reconstruct track parameters"};
DoubleArrayProperty m_zCenter { this, "ZCenter", { -1452.2925, 47.7075 , 1237.7075, 2427.7075 }, "Global z position at which to reconstruct track parameters"};
UnsignedIntegerProperty m_triggerMask { this, "TriggerMask", 0x0, "Trigger mask to analyze (0 = pass all)" };
/// a handle on the Hist/TTree registration service
ServiceHandle<ITHistSvc> m_histSvc;
......
Tracker/TrackerRecAlgs/TrackerClusterFit/test/TrackerClusterFitDbg.py
+ 2
− 2
View file @ 1c5608ed
......@@ -26,7 +26,7 @@ Configurable.configurableRun3Behavior = True
# Configure
ConfigFlags.Input.Files = [
'/eos/project-f/faser-commissioning/TI12Data/Run-001261/Faser-Physics-001261-00000.raw',
'/eos/project-f/faser-commissioning/TI12Data/Run-001332/Faser-Physics-001332-00000.raw',
#'/eos/project-f/faser-commissioning/winter2020CosmicsStand/Run-000608/Faser-Physics-000608-00000.raw',
# '/eos/project-f/faser-commissioning/winter2020CosmicsStand/Run-000608/Faser-Physics-000608-00001.raw',
# '/eos/project-f/faser-commissioning/winter2020CosmicsStand/Run-000608/Faser-Physics-000608-00002.raw',
......@@ -92,7 +92,7 @@ ConfigFlags.Input.Files = [
# '/eos/project-f/faser-commissioning/winter2020CosmicsStand/Run-000608/Faser-Physics-000608-00062.raw'
]
#ConfigFlags.Output.ESDFileName = "run608.ESD.pool.root"
ConfigFlags.Output.ESDFileName = "run001261.ESD.pool.root"
ConfigFlags.Output.ESDFileName = "run001332.ESD.pool.root"
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-XXXX-XXX-XX" # Needed to bypass autoconfig, only the "OFLCOND" matters at the moment
ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
......
Tracker/TrackerRecTools/FaserSiClusterizationTool/src/FaserSCT_ClusteringTool.cxx
+ 2
− 1
View file @ 1c5608ed
......@@ -593,7 +593,8 @@ namespace Tracker
// when this is the case here, and set hitsInThirdTimeBin to zero later on
//
double iphipitch = 1./element->phiPitch();
double shift = m_lorentzAngleTool->getLorentzShift(element->identifyHash());
// double shift = m_lorentzAngleTool->getLorentzShift(element->identifyHash());
double shift = 0;
double stripPitch = design->stripPitch();
bool badStripInClusterOnThisModuleSide = (idGroups.size() != tbinGroups.size());
bool rotate = element->design().shape() == TrackerDD::Trapezoid || element->design().shape() == TrackerDD::Annulus;
......
Tracker/TrackerRecTools/FaserSiClusterizationTool/src/TrackerClusterMakerTool.cxx
+ 2
− 1
View file @ 1c5608ed
......@@ -81,7 +81,8 @@ Tracker::FaserSCT_Cluster* TrackerClusterMakerTool::sctCluster(
const TrackerDD::SiDetectorElement* element,
int errorStrategy) const{
double shift = m_sctLorentzAngleTool->getLorentzShift(element->identifyHash());
// double shift = m_sctLorentzAngleTool->getLorentzShift(element->identifyHash());
double shift = 0;
// const InDetDD::SiLocalPosition& localPosition =
// InDetDD::SiLocalPosition(localPos[Trk::locY),
// localPos[Trk::locX)+shift,0);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment