Merge pull request #43935 from mmusich/mm_dev_ZMuMuMassConstraintPara…

…meterFinder

Miscellaneous improvements to `Alignment/MillePedeAlignmentAlgorithm` unit tests

Alignment/MillePedeAlignmentAlgorithm/plugins/ZMuMuMassConstraintParameterFinder.cc

@@ -47,7 +47,7 @@
 class ZMuMuMassConstraintParameterFinder : public edm::one::EDAnalyzer<edm::one::SharedResources> {
 public:
   explicit ZMuMuMassConstraintParameterFinder(const edm::ParameterSet&);
-  ~ZMuMuMassConstraintParameterFinder() override;
+  ~ZMuMuMassConstraintParameterFinder() override = default;
 
   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
 
@@ -117,8 +117,6 @@ ZMuMuMassConstraintParameterFinder ::ZMuMuMassConstraintParameterFinder(const ed
   muonInfoFromZ_.setupTree("di_muon_from_Z", fs);
 }
 
-ZMuMuMassConstraintParameterFinder ::~ZMuMuMassConstraintParameterFinder() {}
-
 //
 // member functions
 //
@@ -158,8 +156,8 @@ void ZMuMuMassConstraintParameterFinder ::fillDescriptions(edm::ConfigurationDes
   desc.add<double>("ptMin", 15.0);
   desc.add<double>("etaMin", -3.0);
   desc.add<double>("etaMax", 3.0);
-  desc.add<double>("phiMin", -3.1416);
-  desc.add<double>("phiMax", 3.1416);
+  desc.add<double>("phiMin", -M_PI);
+  desc.add<double>("phiMax", M_PI);
   desc.add<double>("minMassPair", 85.8);
   desc.add<double>("maxMassPair", 95.8);
   descriptions.add("zMuMuMassConstraintParameterFinder", desc);

Alignment/MillePedeAlignmentAlgorithm/test/AlignmentRcdChecker.cpp

@@ -8,7 +8,10 @@
 // system includes
 #include <string>
 #include <map>
+#include <ranges>
 #include <vector>
+#include <tuple>
+#include <iterator>
 
 // user includes
 #include "CondFormats/Alignment/interface/AlignTransform.h"
@@ -27,6 +30,20 @@
 // ROOT includes
 #include <TMath.h>
 
+namespace checker {
+  template <typename T1, typename T2>
+  auto zip(const std::vector<T1>& v1, const std::vector<T2>& v2) {
+    std::vector<std::tuple<T1, T2>> result;
+
+    auto minSize = std::min(v1.size(), v2.size());
+    for (size_t i = 0; i < minSize; ++i) {
+      result.emplace_back(v1[i], v2[i]);
+    }
+
+    return result;
+  }
+}  // namespace checker
+
 class AlignmentRcdChecker : public edm::one::EDAnalyzer<> {
 public:
   explicit AlignmentRcdChecker(const edm::ParameterSet& iConfig);
@@ -40,12 +57,14 @@ class AlignmentRcdChecker : public edm::one::EDAnalyzer<> {
                      const Alignments* refAlignments,
                      const Alignments* alignments);
 
-  bool verbose_;
-  bool compareStrict_;
-  std::string label_;
+  const bool verbose_;
+  const bool compareStrict_;
+  const std::string label_;
 
   const edm::ESGetToken<Alignments, TrackerAlignmentRcd> tkAliTokenRef_;
   const edm::ESGetToken<Alignments, TrackerAlignmentRcd> tkAliTokenNew_;
+
+  static constexpr double strictTolerance_ = 1e-4;  // if in cm (i.e. 1 micron)
 };
 
 AlignmentRcdChecker::AlignmentRcdChecker(const edm::ParameterSet& iConfig)
@@ -69,41 +88,35 @@ void AlignmentRcdChecker::inspectRecord(const std::string& rcdname,
   edm::LogPrint("inspectRecord") << " with " << alignments->m_align.size() << " entries";
 
   if (refAlignments && alignments) {
-    double meanX = 0;
-    double rmsX = 0;
-    double meanY = 0;
-    double rmsY = 0;
-    double meanZ = 0;
-    double rmsZ = 0;
-    double meanR = 0;
-    double rmsR = 0;
-    double dPhi;
-    double meanPhi = 0;
-    double rmsPhi = 0;
-
-    std::vector<AlignTransform>::const_iterator iref = refAlignments->m_align.begin();
-    for (std::vector<AlignTransform>::const_iterator i = alignments->m_align.begin(); i != alignments->m_align.end();
-         ++i, ++iref) {
-      meanX += i->translation().x() - iref->translation().x();
-      rmsX += pow(i->translation().x() - iref->translation().x(), 2);
-
-      meanY += i->translation().y() - iref->translation().y();
-      rmsY += pow(i->translation().y() - iref->translation().y(), 2);
-
-      meanZ += i->translation().z() - iref->translation().z();
-      rmsZ += pow(i->translation().z() - iref->translation().z(), 2);
-
-      meanR += i->translation().perp() - iref->translation().perp();
-      rmsR += pow(i->translation().perp() - iref->translation().perp(), 2);
-
-      dPhi = i->translation().phi() - iref->translation().phi();
+    double meanX = 0, rmsX = 0;
+    double meanY = 0, rmsY = 0;
+    double meanZ = 0, rmsZ = 0;
+    double meanR = 0, rmsR = 0;
+    double dPhi, meanPhi = 0, rmsPhi = 0;
+
+    for (const auto& [alignment, refAlignment] : checker::zip(alignments->m_align, refAlignments->m_align)) {
+      auto delta = alignment.translation() - refAlignment.translation();
+
+      meanX += delta.x();
+      rmsX += pow(delta.x(), 2);
+
+      meanY += delta.y();
+      rmsY += pow(delta.y(), 2);
+
+      meanZ += delta.z();
+      rmsZ += pow(delta.z(), 2);
+
+      meanR += delta.perp();
+      rmsR += pow(delta.perp(), 2);
+
+      dPhi = alignment.translation().phi() - refAlignment.translation().phi();
       if (dPhi > M_PI)
         dPhi -= 2.0 * M_PI;
       if (dPhi < -M_PI)
         dPhi += 2.0 * M_PI;
 
       meanPhi += dPhi;
-      rmsPhi += dPhi * dPhi;
+      rmsPhi += std::pow(dPhi, 2);
     }
 
     meanX /= alignments->m_align.size();
@@ -120,24 +133,25 @@ void AlignmentRcdChecker::inspectRecord(const std::string& rcdname,
     if (verbose_) {
       edm::LogPrint("inspectRecord") << "  Compared to previous record:";
       edm::LogPrint("inspectRecord") << "    mean X shift:   " << std::setw(12) << std::scientific
-                                     << std::setprecision(3) << meanX << " (RMS = " << sqrt(rmsX) << ")";
+                                     << std::setprecision(3) << meanX << "  [cm] (RMS = " << sqrt(rmsX) << " [cm])";
       edm::LogPrint("inspectRecord") << "    mean Y shift:   " << std::setw(12) << std::scientific
-                                     << std::setprecision(3) << meanY << " (RMS = " << sqrt(rmsY) << ")";
+                                     << std::setprecision(3) << meanY << "  [cm] (RMS = " << sqrt(rmsY) << " [cm])";
       edm::LogPrint("inspectRecord") << "    mean Z shift:   " << std::setw(12) << std::scientific
-                                     << std::setprecision(3) << meanZ << " (RMS = " << sqrt(rmsZ) << ")";
+                                     << std::setprecision(3) << meanZ << "  [cm] (RMS = " << sqrt(rmsZ) << " [cm])";
       edm::LogPrint("inspectRecord") << "    mean R shift:   " << std::setw(12) << std::scientific
-                                     << std::setprecision(3) << meanR << " (RMS = " << sqrt(rmsR) << ")";
+                                     << std::setprecision(3) << meanR << "  [cm] (RMS = " << sqrt(rmsR) << " [cm])";
       edm::LogPrint("inspectRecord") << "    mean Phi shift: " << std::setw(12) << std::scientific
-                                     << std::setprecision(3) << meanPhi << " (RMS = " << sqrt(rmsPhi) << ")";
+                                     << std::setprecision(3) << meanPhi << " [rad] (RMS = " << sqrt(rmsPhi)
+                                     << " [rad])";
     }  // verbose
 
     if (compareStrict_) {
       // do not let any of the coordinates to fluctuate less then 1um
-      assert(meanX < 1e-4);
-      assert(meanY < 1e-4);
-      assert(meanZ < 1e-4);
-      assert(meanR < 1e-4);
-      assert(meanPhi < 1e-4);
+      assert(meanX < strictTolerance_);    // 1 micron
+      assert(meanY < strictTolerance_);    // 1 micron
+      assert(meanZ < strictTolerance_);    // 1 micron
+      assert(meanR < strictTolerance_);    // 1 micron
+      assert(meanPhi < strictTolerance_);  // 10 micro-rad
     }
 
   } else {

Alignment/MillePedeAlignmentAlgorithm/test/BuildFile.xml

@@ -10,7 +10,12 @@
 <test name="test-pede" command="pede -t">
   <use name="millepede"/>
 </test>
-<test name="testPedeCampaign" command="test_pede.sh"/>
-<test name="testPayloadSanity" command="test_payload_sanity.sh">
-  <flags PRE_TEST="testPedeCampaign"/>
+<test name="test_PedeCampaign" command="test_pede.sh"/>
+<test name="test_PayloadSanity" command="test_payload_sanity.sh">
+  <flags PRE_TEST="test_PedeCampaign"/>
 </test>
+<test name="test_ZMuMuMassConstraintParameterFinder" command="test_ZMuMuMassConstraintParameterFinder.sh"/>
+<bin file="testFindMassContraintParameters.cpp" name="test_ExtractMassConstraint">
+  <flags PRE_TEST="test_ZMuMuMassConstraintParameterFinder"/>
+  <use name="rootcore"/>
+</bin>

Alignment/MillePedeAlignmentAlgorithm/test/ZMuMuMassConstraintParameterFinder/README.md

@@ -0,0 +1,8 @@
+# ZMuMuMassConstraintParameterFinder
+
+To obtain the 'Z → µµ' mass constraint parameters do the following:
+
+- modify the value of `dataset` in `submit_jobs.sh`, if needed
+- modify cuts in `zmumudistribution_cfg.py`, if needed (see possible cuts in `fillDescriptions` in [`Alignment/MillePedeAlignmentAlgorithm/plugins/ZMuMuMassConstraintParameterFinder.cc`](https://github.com/cms-sw/cmssw/blob/master/Alignment/MillePedeAlignmentAlgorithm/plugins/ZMuMuMassConstraintParameterFinder.cc))
+- execute `./submit_jobs.sh`
+- when all jobs are finished you will find the parameters in `submit_${dataset}/zMuMuMassConstraintParameters.txt`

Alignment/MillePedeAlignmentAlgorithm/test/ZMuMuMassConstraintParameterFinder/printParameters.C

@@ -1,4 +1,23 @@
+#include <TFile.h>
+#include <TTree.h>
+#include <TPad.h>
+#include <TH1F.h>
+#include <iostream>
+
+bool isValidFile(const TString& fileName) {
+  TFile* file = TFile::Open(fileName, "read");
+  if (!file || file->IsZombie()) {
+    std::cout << "Error: Invalid file or file is a zombie.\n";
+    return false;
+  }
+  return true;
+}
+
 void printParameters(const TString& fileName) {
+  if (!isValidFile(fileName)) {
+    exit(EXIT_FAILURE);
+  }
+
   TFile* file = TFile::Open(fileName, "read");
   TTree* tree = static_cast<TTree*>(file->Get("zMuMuMassConstraintParameterFinder/di_muon_from_Z"));
   tree->Draw("di_muon_mass>>htemp", "in_mass_window");
@@ -10,4 +29,6 @@ void printParameters(const TString& fileName) {
   std::cout << "  PrimaryMass  = " << htemp->GetMean() << "\n";
   std::cout << "  PrimaryWidth = " << htemp->GetRMS() << "\n";
   std::cout << "========================================\n";
+
+  exit(EXIT_SUCCESS);
 }

Alignment/MillePedeAlignmentAlgorithm/test/testFindMassContraintParameters.cpp

@@ -0,0 +1,5 @@
+#include <iostream>
+#include <sstream>
+#include "Alignment/MillePedeAlignmentAlgorithm/test/ZMuMuMassConstraintParameterFinder/printParameters.C"
+
+int main(int argc, char** argv) { printParameters("output.root"); }

Alignment/MillePedeAlignmentAlgorithm/test/test_ZMuMuMassConstraintParameterFinder.sh

@@ -0,0 +1,4 @@
+#!/bin/bash
+function die { echo $1: status $2; exit $2; }
+
+(cmsRun ${SCRAM_TEST_PATH}/ZMuMuMassConstraintParameterFinder/zmumudistribution_cfg.py inputFiles=/store/relval/CMSSW_14_0_0_pre2/RelValZMM_14/GEN-SIM/133X_mcRun3_2024_realistic_v5_STD_2024_PU-v1/2590000/c38cee3f-99d7-48aa-b236-86f6bbc869b3.root,/store/relval/CMSSW_14_0_0_pre2/RelValZMM_14/GEN-SIM/133X_mcRun3_2024_realistic_v5_STD_2024_PU-v1/2590000/5bf98cca-d491-4e95-98b0-d3acb6ea0807.root,/store/relval/CMSSW_14_0_0_pre2/RelValZMM_14/GEN-SIM/133X_mcRun3_2024_realistic_v5_STD_2024_PU-v1/2590000/1e362cc1-235b-4c32-bb24-178ccac4659f.root) || die 'failed running ZMuMuMassConstraintParameterFinder' $?

Alignment/MillePedeAlignmentAlgorithm/test/test_payload_sanity.sh

@@ -3,5 +3,12 @@ function die { echo $1: status $2; exit $2; }
 
 echo -e "Content of the current directory is: "`ls .`
 INPUTFILE=alignments_MP.db
+
+# Check if the file exists
+if [ ! -f "$INPUTFILE" ]; then
+    echo "Error: $INPUTFILE does not exist."
+    exit 1
+fi
+
 (cmsRun ${SCRAM_TEST_PATH}/AlignmentRcdChecker_cfg.py inputSqliteFile=${INPUTFILE}) || die 'failed running AlignmentRcdChecker' $?
 rm $INPUTFILE