diff --git a/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx b/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx
index e8daf6cea653aac8a9793fb493d7e0d40ebd01b2..7dc182d8b6a8dd0c4083105d436875ddc2cd903b 100644
--- a/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx
+++ b/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx
@@ -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;
diff --git a/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.h b/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.h
index a6bc1bf981d443b466ac9450de3ca8118bdbc5f2..54868f69ebbfabce756163e0c7180e00e794747f 100644
--- a/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.h
+++ b/Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.h
@@ -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;
diff --git a/Tracker/TrackerRecAlgs/TrackerClusterFit/test/TrackerClusterFitDbg.py b/Tracker/TrackerRecAlgs/TrackerClusterFit/test/TrackerClusterFitDbg.py
index 77182425696ee70043778ba15026b8a167ebfad5..8c0309ca4f3495d42e87e6be720f63ec7db20688 100644
--- a/Tracker/TrackerRecAlgs/TrackerClusterFit/test/TrackerClusterFitDbg.py
+++ b/Tracker/TrackerRecAlgs/TrackerClusterFit/test/TrackerClusterFitDbg.py
@@ -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
diff --git a/Tracker/TrackerRecTools/FaserSiClusterizationTool/src/FaserSCT_ClusteringTool.cxx b/Tracker/TrackerRecTools/FaserSiClusterizationTool/src/FaserSCT_ClusteringTool.cxx
index b8c6cfbae285b0470f32d511df4a8b220613d4f7..dece2904daf8acc0ac4ac9e2d71c681713df8353 100644
--- a/Tracker/TrackerRecTools/FaserSiClusterizationTool/src/FaserSCT_ClusteringTool.cxx
+++ b/Tracker/TrackerRecTools/FaserSiClusterizationTool/src/FaserSCT_ClusteringTool.cxx
@@ -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;
diff --git a/Tracker/TrackerRecTools/FaserSiClusterizationTool/src/TrackerClusterMakerTool.cxx b/Tracker/TrackerRecTools/FaserSiClusterizationTool/src/TrackerClusterMakerTool.cxx
index 78bcb43a710e8e7866a06dfea1136b819c162bb8..2313a1584a00f63827cbfd9183c5fbf9acd29feb 100644
--- a/Tracker/TrackerRecTools/FaserSiClusterizationTool/src/TrackerClusterMakerTool.cxx
+++ b/Tracker/TrackerRecTools/FaserSiClusterizationTool/src/TrackerClusterMakerTool.cxx
@@ -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);