Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:34 

0001 /****************************************************************************
0002 * Authors: 
0003 *  Jan Kašpar (jan.kaspar@gmail.com) 
0004 ****************************************************************************/
0005 
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 
0011 #include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsMethods.h"
0012 #include "CondFormats/AlignmentRecord/interface/RPRealAlignmentRecord.h"
0013 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0014 #include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
0015 
0016 #include "CalibPPS/AlignmentRelative/interface/StraightTrackAlignment.h"
0017 #include "CalibPPS/AlignmentRelative/interface/Utilities.h"
0018 
0019 #include "CalibPPS/AlignmentRelative/interface/IdealResult.h"
0020 #include "CalibPPS/AlignmentRelative/interface/JanAlignmentAlgorithm.h"
0021 
0022 #include <set>
0023 #include <unordered_set>
0024 #include <vector>
0025 #include <string>
0026 
0027 #include "TDecompLU.h"
0028 #include "TH1D.h"
0029 #include "TFile.h"
0030 #include "TGraph.h"
0031 #include "TCanvas.h"
0032 
0033 //#define DEBUG
0034 
0035 using namespace edm;
0036 using namespace std;
0037 
0038 //----------------------------------------------------------------------------------------------------
0039 
0040 TH1D *StraightTrackAlignment::newResiduaHist(const char *name) {
0041   return new TH1D(name, ";residual   (mm)", 2000, -3., +3.);  // in mm
0042 }
0043 
0044 //----------------------------------------------------------------------------------------------------
0045 
0046 TGraph *newGraph(const string &name, const string &title) {
0047   TGraph *g = new TGraph();
0048   g->SetName(name.c_str());
0049   g->SetTitle(title.c_str());
0050   return g;
0051 }
0052 
0053 //----------------------------------------------------------------------------------------------------
0054 
0055 StraightTrackAlignment::RPSetPlots::RPSetPlots(const string &_name) : name(_name) {
0056   chisqn_lin_fitted = new TH1D("chi^2 norm, lin, fitted", ";#chi^{2}/ndf;", 5000, 0., 500.);
0057   chisqn_lin_selected = new TH1D("chi^2 norm, lin, selected", ";#chi^{2}/ndf;", 5000, 0., 500.);
0058   chisqn_log_fitted = new TH1D("chi^2 norm, log, fitted", ";log_{10}(#chi^{2}/ndf);", 700, -1., 6.);
0059   chisqn_log_selected = new TH1D("chi^2 norm, log, selected", ";log_{10}(#chi^{2}/ndf);", 700, -1., 6.);
0060 
0061   fitAxVsAyGraph_fitted = newGraph("ax vs. ay, fitted", ";a_{x}   (rad);a_{y}   (rad)");
0062   fitAxVsAyGraph_selected = newGraph("ax vs. ay, selected", ";a_{x}   (rad);a_{y}   (rad)");
0063   fitBxVsByGraph_fitted = newGraph("bx vs. by, fitted", ";b_{x}   (mm);b_{y}   (mm)");
0064   fitBxVsByGraph_selected = newGraph("bx vs. by, selected", ";b_{x}   (mm);b_{y}   (mm)");
0065 }
0066 //----------------------------------------------------------------------------------------------------
0067 
0068 void StraightTrackAlignment::RPSetPlots::free() {
0069   delete chisqn_lin_fitted;
0070   delete chisqn_lin_selected;
0071   delete chisqn_log_fitted;
0072   delete chisqn_log_selected;
0073 
0074   delete fitAxVsAyGraph_fitted;
0075   delete fitAxVsAyGraph_selected;
0076   delete fitBxVsByGraph_fitted;
0077   delete fitBxVsByGraph_selected;
0078 }
0079 
0080 //----------------------------------------------------------------------------------------------------
0081 
0082 void StraightTrackAlignment::RPSetPlots::write() const {
0083   chisqn_lin_fitted->Write();
0084   chisqn_lin_selected->Write();
0085   chisqn_log_fitted->Write();
0086   chisqn_log_selected->Write();
0087 
0088   fitAxVsAyGraph_fitted->Write();
0089   fitAxVsAyGraph_selected->Write();
0090   fitBxVsByGraph_fitted->Write();
0091   fitBxVsByGraph_selected->Write();
0092 }
0093 
0094 //----------------------------------------------------------------------------------------------------
0095 //----------------------------------------------------------------------------------------------------
0096 
0097 StraightTrackAlignment::StraightTrackAlignment(const ParameterSet &ps)
0098     : verbosity(ps.getUntrackedParameter<unsigned int>("verbosity", 0)),
0099 
0100       rpIds(ps.getParameter<vector<unsigned int> >("rpIds")),
0101       excludePlanes(ps.getParameter<vector<unsigned int> >("excludePlanes")),
0102       z0(ps.getParameter<double>("z0")),
0103 
0104       maxEvents(ps.getParameter<signed int>("maxEvents")),
0105 
0106       removeImpossible(ps.getParameter<bool>("removeImpossible")),
0107       requireNumberOfUnits(ps.getParameter<unsigned int>("requireNumberOfUnits")),
0108       requireAtLeast3PotsInOverlap(ps.getParameter<bool>("requireAtLeast3PotsInOverlap")),
0109       requireOverlap(ps.getParameter<bool>("requireOverlap")),
0110       cutOnChiSqPerNdf(ps.getParameter<bool>("cutOnChiSqPerNdf")),
0111       chiSqPerNdfCut(ps.getParameter<double>("chiSqPerNdfCut")),
0112       maxTrackAx(ps.getParameter<double>("maxTrackAx")),
0113       maxTrackAy(ps.getParameter<double>("maxTrackAy")),
0114 
0115       fileNamePrefix(ps.getParameter<string>("fileNamePrefix")),
0116       expandedFileNamePrefix(ps.getParameter<string>("expandedFileNamePrefix")),
0117       factoredFileNamePrefix(ps.getParameter<string>("factoredFileNamePrefix")),
0118       preciseXMLFormat(ps.getParameter<bool>("preciseXMLFormat")),
0119       saveXMLUncertainties(ps.getParameter<bool>("saveXMLUncertainties")),
0120 
0121       saveIntermediateResults(ps.getParameter<bool>("saveIntermediateResults")),
0122       taskDataFileName(ps.getParameter<string>("taskDataFileName")),
0123       taskDataFile(nullptr),
0124 
0125       task(ps),
0126       fitter(ps),
0127 
0128       buildDiagnosticPlots(ps.getParameter<bool>("buildDiagnosticPlots")),
0129       diagnosticsFile(ps.getParameter<string>("diagnosticsFile")),
0130       fitNdfHist_fitted(new TH1D("ndf_fitted", ";ndf;", 41, -4.5, 36.5)),
0131       fitNdfHist_selected(new TH1D("ndf_selected", ";ndf;", 41, -4.5, 36.5)),
0132       fitPHist_fitted(new TH1D("p_fitted", ";p value;", 100, 0., 1.)),
0133       fitPHist_selected(new TH1D("p_selected", ";p value;", 100, 0., 1.)),
0134       fitAxHist_fitted(new TH1D("ax_fitted", ";a_{x}   (rad);", 10000, -0.1, 0.1)),
0135       fitAxHist_selected(new TH1D("ax_selected", ";a_{x}   (rad);", 10000, -0.1, 0.1)),
0136       fitAyHist_fitted(new TH1D("ay_fitted", ";a_{y}   (rad);", 10000, -0.1, 0.1)),
0137       fitAyHist_selected(new TH1D("ay_selected", ";a_{y}   (rad);", 10000, -0.1, 0.1)),
0138       fitBxHist_fitted(new TH1D("bx_fitted", ";b_{x}   (mm);", 500, -30., 30.)),
0139       fitBxHist_selected(new TH1D("bx_selected", ";b_{x}   (mm);", 500, -30., 30.)),
0140       fitByHist_fitted(new TH1D("by_fitted", ";b_{y}   (mm);", 500, -30., 30.)),
0141       fitByHist_selected(new TH1D("by_selected", ";b_{y}   (mm);", 500, -30., 30.)),
0142 
0143       globalPlots("global") {
0144   // open task data file
0145   if (!taskDataFileName.empty())
0146     taskDataFile = new TFile(taskDataFileName.c_str(), "recreate");
0147 
0148   // instantiate algorithm objects
0149   // (and save them)
0150   vector<string> alNames(ps.getParameter<vector<string> >("algorithms"));
0151   for (unsigned int i = 0; i < alNames.size(); i++) {
0152     AlignmentAlgorithm *a = nullptr;
0153 
0154     if (alNames[i] == "Ideal") {
0155       IdealResult *ir = new IdealResult(ps, &task);
0156       a = ir;
0157     } else if (alNames[i] == "Jan") {
0158       JanAlignmentAlgorithm *jaa = new JanAlignmentAlgorithm(ps, &task);
0159       a = jaa;
0160     }
0161 
0162     if (a)
0163       algorithms.push_back(a);
0164     else
0165       throw cms::Exception("PPS") << "Unknown alignment algorithm `" << alNames[i] << "'.";
0166   }
0167 
0168   // get constraints type
0169   string ct = ps.getParameter<string>("constraintsType");
0170   if (ct.compare("fixedDetectors") == 0)
0171     constraintsType = ctFixedDetectors;
0172   else if (ct.compare("standard") == 0)
0173     constraintsType = ctStandard;
0174   else
0175     throw cms::Exception("PPS") << "Unknown constraints type `" << ct << "'.";
0176 
0177   // parse additional accepted RP sets
0178   string aars_str = ps.getParameter<string>("additionalAcceptedRPSets");
0179 
0180   size_t idx_b = 0, idx_e = string::npos;
0181   while (idx_b != string::npos) {
0182     // get one block - portion between successive ";"
0183     idx_e = aars_str.find(';', idx_b);
0184     size_t len = (idx_e == string::npos) ? string::npos : idx_e - idx_b;
0185     string block = aars_str.substr(idx_b, len);
0186 
0187     // process the block
0188     if (!block.empty()) {
0189       set<unsigned int> rpSet;
0190 
0191       // isolate bits (= RP ids)
0192       size_t bi_b = 0, bi_e = string::npos;
0193       while (bi_b != string::npos) {
0194         bi_e = block.find(',', bi_b);
0195         size_t bit_len = (bi_e == string::npos) ? string::npos : bi_e - bi_b;
0196         const string &bit = block.substr(bi_b, bit_len);
0197 
0198         unsigned int rp = atoi(bit.c_str());
0199         rpSet.insert(rp);
0200 
0201         bi_b = (bi_e == string::npos) ? string::npos : bi_e + 1;
0202       }
0203 
0204       additionalAcceptedRPSets.push_back(rpSet);
0205     }
0206 
0207     // move to next block
0208     idx_b = (idx_e == string::npos) ? string::npos : idx_e + 1;
0209   }
0210 }
0211 
0212 //----------------------------------------------------------------------------------------------------
0213 
0214 StraightTrackAlignment::~StraightTrackAlignment() {
0215   if (taskDataFile)
0216     delete taskDataFile;
0217 
0218   for (vector<AlignmentAlgorithm *>::iterator it = algorithms.begin(); it != algorithms.end(); ++it)
0219     delete (*it);
0220 
0221   delete fitNdfHist_fitted;
0222   delete fitNdfHist_selected;
0223   delete fitPHist_fitted;
0224   delete fitPHist_selected;
0225   delete fitAxHist_fitted;
0226   delete fitAyHist_fitted;
0227   delete fitAxHist_selected;
0228   delete fitAyHist_selected;
0229   delete fitBxHist_fitted;
0230   delete fitByHist_fitted;
0231   delete fitBxHist_selected;
0232   delete fitByHist_selected;
0233 
0234   globalPlots.free();
0235 
0236   for (auto &p : rpSetPlots)
0237     p.second.free();
0238 
0239   for (auto it = residuaHistograms.begin(); it != residuaHistograms.end(); ++it) {
0240     delete it->second.total_fitted;
0241     delete it->second.total_selected;
0242     delete it->second.selected_vs_chiSq;
0243     for (auto sit = it->second.perRPSet_fitted.begin(); sit != it->second.perRPSet_fitted.end(); ++sit)
0244       delete sit->second;
0245     for (auto sit = it->second.perRPSet_selected.begin(); sit != it->second.perRPSet_selected.end(); ++sit)
0246       delete sit->second;
0247   }
0248 }
0249 
0250 //----------------------------------------------------------------------------------------------------
0251 
0252 void StraightTrackAlignment::begin(edm::ESHandle<CTPPSRPAlignmentCorrectionsData> hRealAlignment,
0253                                    edm::ESHandle<CTPPSGeometry> hRealGeometry,
0254                                    edm::ESHandle<CTPPSGeometry> hMisalignedGeometry) {
0255   // reset counters
0256   eventsTotal = 0;
0257   eventsFitted = 0;
0258   eventsSelected = 0;
0259   fittedTracksPerRPSet.clear();
0260   selectedTracksPerRPSet.clear();
0261   selectedHitsPerPlane.clear();
0262 
0263   // prepare geometry
0264   task.buildGeometry(rpIds, excludePlanes, hRealGeometry.product(), z0, task.geometry);
0265 
0266   // build matrix index mappings
0267   task.buildIndexMaps();
0268 
0269   // print geometry info
0270   if (verbosity > 1)
0271     task.geometry.print();
0272 
0273   // save task (including geometry) and fitter
0274   if (taskDataFile) {
0275     taskDataFile->WriteObject(&task, "task");
0276     taskDataFile->WriteObject(&fitter, "fitter");
0277   }
0278 
0279   // initiate the algorithms
0280   for (const auto &a : algorithms)
0281     a->begin(hRealGeometry.product(), hMisalignedGeometry.product());
0282 
0283   // get initial alignments
0284   initialAlignments = *hRealAlignment;
0285 }
0286 
0287 //----------------------------------------------------------------------------------------------------
0288 
0289 void StraightTrackAlignment::processEvent(const edm::EventID &eventId,
0290                                           const DetSetVector<TotemRPUVPattern> &uvPatternsStrip,
0291                                           const DetSetVector<CTPPSDiamondRecHit> &hitsDiamond,
0292                                           const edm::DetSetVector<CTPPSPixelRecHit> &hitsPixel,
0293                                           const DetSetVector<CTPPSPixelLocalTrack> &tracksPixel) {
0294   eventsTotal++;
0295 
0296   if (verbosity > 9)
0297     printf("\n---------- StraightTrackAlignment::ProcessEvent (%u:%llu)\n", eventId.run(), eventId.event());
0298 
0299   // -------------------- STEP 1: get hits from selected RPs
0300 
0301   HitCollection hitSelection;
0302 
0303   // strips
0304   for (const auto &pv : uvPatternsStrip) {
0305     const CTPPSDetId detId(pv.detId());
0306     const unsigned int rpDecId = detId.arm() * 100 + detId.station() * 10 + detId.rp();
0307 
0308     // skip if RP not selected
0309     if (find(rpIds.begin(), rpIds.end(), rpDecId) == rpIds.end())
0310       continue;
0311 
0312     // require exactly 1 U and 1 V pattern, i.e. unique U-V association
0313     unsigned int n_U = 0, n_V = 0;
0314     unsigned int idx_U = 0, idx_V = 0;
0315     for (unsigned int pi = 0; pi < pv.size(); pi++) {
0316       const TotemRPUVPattern &pattern = pv[pi];
0317 
0318       switch (pattern.projection()) {
0319         case TotemRPUVPattern::projU:
0320           n_U++;
0321           idx_U = pi;
0322           break;
0323 
0324         case TotemRPUVPattern::projV:
0325           n_V++;
0326           idx_V = pi;
0327           break;
0328 
0329         default:
0330           break;
0331       }
0332     }
0333 
0334     if (n_U != 1 || n_V != 1)
0335       continue;
0336 
0337     // skip if patterns not reasonable
0338     if (!pv[idx_U].fittable() || !pv[idx_V].fittable())
0339       continue;
0340 
0341     // add hits from the U and V pattern
0342     for (const auto &pattern : {pv[idx_U], pv[idx_V]}) {
0343       for (const auto &hitsDetSet : pattern.hits()) {
0344         // skip if sensor not in geometry
0345         if (!task.geometry.isValidSensorId(hitsDetSet.detId()))
0346           continue;
0347 
0348         const double &z = task.geometry.get(hitsDetSet.detId()).z;
0349         for (auto &hit : hitsDetSet)
0350           hitSelection.emplace_back(Hit(hitsDetSet.detId(), 2, hit.position(), hit.sigma(), z));
0351       }
0352     }
0353   }
0354 
0355   // diamonds
0356   for (const auto &pv : hitsDiamond) {
0357     // skip if RP not selected
0358     const CTPPSDetId detId(pv.detId());
0359     const unsigned int rpDecId = detId.arm() * 100 + detId.station() * 10 + detId.rp();
0360 
0361     // skip if RP not selected
0362     if (find(rpIds.begin(), rpIds.end(), rpDecId) == rpIds.end())
0363       continue;
0364 
0365     // skip if sensor not in geometry
0366     if (!task.geometry.isValidSensorId(detId))
0367       continue;
0368 
0369     const double &z = task.geometry.get(pv.detId()).z;
0370 
0371     for (const auto &h : pv)
0372       hitSelection.emplace_back(Hit(pv.detId(), 1, h.x(), h.xWidth() / sqrt(12.), z));
0373   }
0374 
0375   // pixels: data from rec hits
0376   map<unsigned int, unsigned int> pixelPlaneMultiplicity;
0377   for (const auto &pv : hitsPixel)
0378     pixelPlaneMultiplicity[pv.detId()] += pv.size();
0379 
0380   for (const auto &pv : hitsPixel) {
0381     // skip if RP not selected
0382     const CTPPSDetId detId(pv.detId());
0383     const unsigned int rpDecId = detId.arm() * 100 + detId.station() * 10 + detId.rp();
0384 
0385     // skip if RP not selected
0386     if (find(rpIds.begin(), rpIds.end(), rpDecId) == rpIds.end())
0387       continue;
0388 
0389     // skip if sensor not in geometry
0390     if (!task.geometry.isValidSensorId(detId))
0391       continue;
0392 
0393     // skip if plane multiplicity greater than 1
0394     if (pixelPlaneMultiplicity[pv.detId()] > 1)
0395       continue;
0396 
0397     for (const auto &h : pv) {
0398       const auto &dg = task.geometry.get(pv.detId());
0399       const double dz1 = dg.getDirectionData(1).dz;
0400       const double dz2 = dg.getDirectionData(2).dz;
0401       const double z = dg.z + h.point().x() * dz1 + h.point().y() * dz2;
0402 
0403       double x_unc = sqrt(h.error().xx());
0404       double y_unc = sqrt(h.error().yy());
0405 
0406       // TODO: Check this
0407 
0408       if (x_unc <= 0.)
0409         x_unc = 10E-3;  // mm
0410       if (y_unc <= 0.)
0411         y_unc = 10E-3;  // mm
0412 
0413       hitSelection.emplace_back(Hit(pv.detId(), 1, h.point().x(), x_unc, z));
0414       hitSelection.emplace_back(Hit(pv.detId(), 2, h.point().y(), y_unc, z));
0415     }
0416   }
0417 
0418   // pixels: data from tracks
0419   for (const auto &pv : tracksPixel) {
0420     const CTPPSDetId rpId(pv.detId());
0421     const unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0422 
0423     // skip if RP not selected
0424     if (find(rpIds.begin(), rpIds.end(), rpDecId) == rpIds.end())
0425       continue;
0426 
0427     // skip if more than 1 (valid) track in the RP
0428     unsigned int n_valid_tracks = 0;
0429     for (const auto &tr : pv) {
0430       if (tr.isValid())
0431         n_valid_tracks++;
0432     }
0433 
0434     if (n_valid_tracks > 1)
0435       continue;
0436 
0437     // go through all valid tracks
0438     for (const auto &tr : pv) {
0439       if (!tr.isValid())
0440         continue;
0441 
0442       // go through all associated rec hits
0443       for (const auto &hv : tr.hits()) {
0444         const CTPPSPixelDetId senId(hv.detId());
0445 
0446         // skip if sensor not in geometry
0447         if (!task.geometry.isValidSensorId(senId))
0448           continue;
0449 
0450         for (const auto &h : hv) {
0451           // skip hit if not used for fit
0452           if (!h.isUsedForFit())
0453             continue;
0454 
0455           const auto &dg = task.geometry.get(senId);
0456           const double dz1 = dg.getDirectionData(1).dz;
0457           const double dz2 = dg.getDirectionData(2).dz;
0458           const double z = dg.z + h.point().x() * dz1 + h.point().y() * dz2;
0459 
0460           double x_unc = sqrt(h.error().xx());
0461           double y_unc = sqrt(h.error().yy());
0462 
0463           if (x_unc <= 0.)
0464             x_unc = 10E-3;  // mm
0465           if (y_unc <= 0.)
0466             y_unc = 10E-3;  // mm
0467 
0468           hitSelection.emplace_back(Hit(senId, 1, h.point().x(), x_unc, z));
0469           hitSelection.emplace_back(Hit(senId, 2, h.point().y(), y_unc, z));
0470         }
0471       }
0472     }
0473   }
0474 
0475   if (hitSelection.empty())
0476     return;
0477 
0478   // -------------------- STEP 2: fit + outlier rejection
0479 
0480   LocalTrackFit trackFit;
0481   if (!fitter.fit(hitSelection, task.geometry, trackFit))
0482     return;
0483 
0484   set<unsigned int> selectedRPs;
0485   for (const auto &hit : hitSelection) {
0486     CTPPSDetId detId(hit.id);
0487     const unsigned int decRPId = detId.arm() * 100 + detId.station() * 10 + detId.rp();
0488     selectedRPs.insert(decRPId);
0489   }
0490 
0491   eventsFitted++;
0492   fittedTracksPerRPSet[selectedRPs]++;
0493 
0494   // -------------------- STEP 3: quality checks
0495 
0496   bool top = false, bottom = false, horizontal = false;
0497   unordered_set<unsigned int> units;
0498   for (const auto &rp : selectedRPs) {
0499     unsigned int rpIdx = rp % 10;
0500     unsigned int stId = rp / 10;
0501     unsigned int unitId = stId * 10;
0502     if (rpIdx > 2)
0503       unitId++;
0504 
0505     if (rpIdx == 0 || rpIdx == 4)
0506       top = true;
0507     if (rpIdx == 1 || rpIdx == 5)
0508       bottom = true;
0509     if (rpIdx == 2 || rpIdx == 3)
0510       horizontal = true;
0511 
0512     units.insert(unitId);
0513   }
0514 
0515   bool overlap = (top && horizontal) || (bottom && horizontal);
0516 
0517   bool rp_set_accepted = true;
0518 
0519   // impossible signature
0520   if (removeImpossible && top && bottom)
0521     rp_set_accepted = false;
0522 
0523   // cleanliness cuts
0524   if (units.size() < requireNumberOfUnits)
0525     rp_set_accepted = false;
0526 
0527   if (requireOverlap && !overlap)
0528     rp_set_accepted = false;
0529 
0530   if (requireAtLeast3PotsInOverlap && overlap && selectedRPs.size() < 3)
0531     rp_set_accepted = false;
0532 
0533   // is it an additional accepted RP set?
0534   if (find(additionalAcceptedRPSets.begin(), additionalAcceptedRPSets.end(), selectedRPs) !=
0535       additionalAcceptedRPSets.end())
0536     rp_set_accepted = true;
0537 
0538   if (verbosity > 5)
0539     printf("* rp set accepted: %u\n", rp_set_accepted);
0540 
0541   bool selected = rp_set_accepted;
0542 
0543   // too bad chisq
0544   if (cutOnChiSqPerNdf && trackFit.chiSqPerNdf() > chiSqPerNdfCut)
0545     selected = false;
0546 
0547   // parallelity cut
0548   if (fabs(trackFit.ax) > maxTrackAx || fabs(trackFit.ay) > maxTrackAy)
0549     selected = false;
0550 
0551   updateDiagnosticHistograms(hitSelection, selectedRPs, trackFit, selected);
0552 
0553   if (verbosity > 5)
0554     printf("* event selected: %u\n", selected);
0555 
0556   if (!selected)
0557     return;
0558 
0559   // update counters
0560   eventsSelected++;
0561   selectedTracksPerRPSet[selectedRPs]++;
0562 
0563   for (const auto &h : hitSelection)
0564     selectedHitsPerPlane[h.id]++;
0565 
0566   // -------------------- STEP 4: FEED ALGORITHMS
0567 
0568   for (auto &a : algorithms)
0569     a->feed(hitSelection, trackFit);
0570 
0571   // -------------------- STEP 5: ENOUGH TRACKS?
0572 
0573   if (eventsSelected == maxEvents)
0574     throw "StraightTrackAlignment: Number of tracks processed reached maximum";
0575 }
0576 
0577 //----------------------------------------------------------------------------------------------------
0578 
0579 void StraightTrackAlignment::updateDiagnosticHistograms(const HitCollection &selection,
0580                                                         const set<unsigned int> &selectedRPs,
0581                                                         const LocalTrackFit &trackFit,
0582                                                         bool trackSelected) {
0583   if (!buildDiagnosticPlots)
0584     return;
0585 
0586   fitNdfHist_fitted->Fill(trackFit.ndf);
0587   fitPHist_fitted->Fill(trackFit.pValue());
0588   fitAxHist_fitted->Fill(trackFit.ax);
0589   fitAyHist_fitted->Fill(trackFit.ay);
0590   fitBxHist_fitted->Fill(trackFit.bx);
0591   fitByHist_fitted->Fill(trackFit.by);
0592 
0593   globalPlots.chisqn_lin_fitted->Fill(trackFit.chiSqPerNdf());
0594   globalPlots.chisqn_log_fitted->Fill(log10(trackFit.chiSqPerNdf()));
0595   globalPlots.fitAxVsAyGraph_fitted->SetPoint(globalPlots.fitAxVsAyGraph_fitted->GetN(), trackFit.ax, trackFit.ay);
0596   globalPlots.fitBxVsByGraph_fitted->SetPoint(globalPlots.fitBxVsByGraph_fitted->GetN(), trackFit.bx, trackFit.by);
0597 
0598   if (trackSelected) {
0599     fitNdfHist_selected->Fill(trackFit.ndf);
0600     fitPHist_selected->Fill(trackFit.pValue());
0601     fitAxHist_selected->Fill(trackFit.ax);
0602     fitAyHist_selected->Fill(trackFit.ay);
0603     fitBxHist_selected->Fill(trackFit.bx);
0604     fitByHist_selected->Fill(trackFit.by);
0605 
0606     globalPlots.chisqn_lin_selected->Fill(trackFit.chiSqPerNdf());
0607     globalPlots.chisqn_log_selected->Fill(log10(trackFit.chiSqPerNdf()));
0608     globalPlots.fitAxVsAyGraph_selected->SetPoint(
0609         globalPlots.fitAxVsAyGraph_selected->GetN(), trackFit.ax, trackFit.ay);
0610     globalPlots.fitBxVsByGraph_selected->SetPoint(
0611         globalPlots.fitBxVsByGraph_selected->GetN(), trackFit.bx, trackFit.by);
0612   }
0613 
0614   auto it = rpSetPlots.find(selectedRPs);
0615   if (it == rpSetPlots.end())
0616     it = rpSetPlots.insert({selectedRPs, RPSetPlots(setToString(selectedRPs))}).first;
0617 
0618   it->second.chisqn_lin_fitted->Fill(trackFit.chiSqPerNdf());
0619   it->second.chisqn_log_fitted->Fill(log10(trackFit.chiSqPerNdf()));
0620   it->second.fitAxVsAyGraph_fitted->SetPoint(it->second.fitAxVsAyGraph_fitted->GetN(), trackFit.ax, trackFit.ay);
0621   it->second.fitBxVsByGraph_fitted->SetPoint(it->second.fitBxVsByGraph_fitted->GetN(), trackFit.bx, trackFit.by);
0622 
0623   if (trackSelected) {
0624     it->second.chisqn_lin_selected->Fill(trackFit.chiSqPerNdf());
0625     it->second.chisqn_log_selected->Fill(log10(trackFit.chiSqPerNdf()));
0626     it->second.fitAxVsAyGraph_selected->SetPoint(it->second.fitAxVsAyGraph_selected->GetN(), trackFit.ax, trackFit.ay);
0627     it->second.fitBxVsByGraph_selected->SetPoint(it->second.fitBxVsByGraph_selected->GetN(), trackFit.bx, trackFit.by);
0628   }
0629 
0630   for (const auto &hit : selection) {
0631     unsigned int id = hit.id;
0632 
0633     const DetGeometry &geom = task.geometry.get(id);
0634     const auto dirData = geom.getDirectionData(hit.dirIdx);
0635 
0636     double m = hit.position + dirData.s - (hit.z - geom.z) * dirData.dz;
0637     double x = trackFit.ax * hit.z + trackFit.bx;
0638     double y = trackFit.ay * hit.z + trackFit.by;
0639     double f = x * dirData.dx + y * dirData.dy;
0640     double R = m - f;
0641 
0642     auto it = residuaHistograms.find(id);
0643     if (it == residuaHistograms.end()) {
0644       it = residuaHistograms.insert(pair<unsigned int, ResiduaHistogramSet>(id, ResiduaHistogramSet())).first;
0645       char buf[30];
0646       sprintf(buf, "%u: total_fitted", id);
0647       it->second.total_fitted = newResiduaHist(buf);
0648       sprintf(buf, "%u: total_selected", id);
0649       it->second.total_selected = newResiduaHist(buf);
0650       it->second.selected_vs_chiSq = new TGraph();
0651       sprintf(buf, "%u: selected_vs_chiSq", id);
0652       it->second.selected_vs_chiSq->SetName(buf);
0653     }
0654 
0655     it->second.total_fitted->Fill(R);
0656 
0657     if (trackSelected) {
0658       it->second.total_selected->Fill(R);
0659       it->second.selected_vs_chiSq->SetPoint(it->second.selected_vs_chiSq->GetN(), trackFit.chiSqPerNdf(), R);
0660     }
0661 
0662     auto sit = it->second.perRPSet_fitted.find(selectedRPs);
0663     if (sit == it->second.perRPSet_fitted.end()) {
0664       char buf[10];
0665       sprintf(buf, "%u: ", id);
0666       string label = buf;
0667       label += setToString(selectedRPs);
0668       sit =
0669           it->second.perRPSet_fitted.insert(pair<set<unsigned int>, TH1D *>(selectedRPs, newResiduaHist(label.c_str())))
0670               .first;
0671     }
0672 
0673     sit->second->Fill(R);
0674 
0675     if (trackSelected) {
0676       sit = it->second.perRPSet_selected.find(selectedRPs);
0677       if (sit == it->second.perRPSet_selected.end()) {
0678         char buf[10];
0679         sprintf(buf, "%u: ", id);
0680         string label = buf;
0681         label += setToString(selectedRPs);
0682         sit = it->second.perRPSet_selected
0683                   .insert(pair<set<unsigned int>, TH1D *>(selectedRPs, newResiduaHist(label.c_str())))
0684                   .first;
0685       }
0686 
0687       sit->second->Fill(R);
0688     }
0689   }
0690 }
0691 
0692 //----------------------------------------------------------------------------------------------------
0693 
0694 void StraightTrackAlignment::buildConstraints(vector<AlignmentConstraint> &constraints) {
0695   constraints.clear();
0696 
0697   switch (constraintsType) {
0698     case ctFixedDetectors:
0699       task.buildFixedDetectorsConstraints(constraints);
0700       return;
0701 
0702     case ctStandard:
0703       task.buildStandardConstraints(constraints);
0704       return;
0705   }
0706 }
0707 
0708 //----------------------------------------------------------------------------------------------------
0709 
0710 void StraightTrackAlignment::finish() {
0711   // print statistics
0712   if (verbosity) {
0713     printf("----------------------------------------------------------------------------------------------------\n");
0714     printf("\n>> StraightTrackAlignment::Finish\n");
0715     printf("\tevents total = %i\n", eventsTotal);
0716     printf("\tevents fitted = %i\n", eventsFitted);
0717     printf("\tevents selected = %i\n", eventsSelected);
0718 
0719     printf("\n* events per RP set:\n");
0720     printf("%30s  %10s%10s\n", "set of RPs", "fitted", "selected");
0721 
0722     for (auto it = fittedTracksPerRPSet.begin(); it != fittedTracksPerRPSet.end(); ++it) {
0723       const string &label = setToString(it->first);
0724 
0725       auto sit = selectedTracksPerRPSet.find(it->first);
0726       unsigned long sv = (sit == selectedTracksPerRPSet.end()) ? 0 : sit->second;
0727 
0728       printf("%30s :%10lu%10lu\n", label.c_str(), it->second, sv);
0729     }
0730 
0731     if (verbosity >= 2) {
0732       printf("\n* hits per plane:\n");
0733       for (const auto it : selectedHitsPerPlane) {
0734         printf("    ");
0735         printId(it.first);
0736         printf(" : %u\n", it.second);
0737       }
0738     }
0739   }
0740 
0741   // write diagnostics plots
0742   saveDiagnostics();
0743 
0744   // run analysis
0745   for (auto a : algorithms)
0746     a->analyze();
0747 
0748   // build constraints
0749   vector<AlignmentConstraint> constraints;
0750   buildConstraints(constraints);
0751 
0752   // save constraints
0753   if (taskDataFile)
0754     taskDataFile->WriteObject(&constraints, "constraints");
0755 
0756   if (verbosity) {
0757     printf("\n>> StraightTrackAlignment::Finish > %lu constraints built\n", constraints.size());
0758     for (unsigned int i = 0; i < constraints.size(); i++)
0759       printf("\t%s\n", constraints[i].name.c_str());
0760   }
0761 
0762   // solve
0763   vector<map<unsigned int, AlignmentResult> > results;
0764   for (auto algorithm : algorithms) {
0765     TDirectory *dir = nullptr;
0766     if (taskDataFile && saveIntermediateResults)
0767       dir = taskDataFile->mkdir((algorithm->getName() + "_data").c_str());
0768 
0769     results.resize(results.size() + 1);
0770     unsigned int rf = algorithm->solve(constraints, results.back(), dir);
0771 
0772     if (rf)
0773       throw cms::Exception("PPS") << "The Solve method of `" << algorithm->getName()
0774                                   << "' algorithm has failed (return value " << rf << ").";
0775   }
0776 
0777   // print results
0778   if (verbosity) {
0779     printf("\n>> StraightTrackAlignment::Finish > Print\n");
0780 
0781     printLineSeparator(results);
0782     printQuantitiesLine(results);
0783     printAlgorithmsLine(results);
0784 
0785     signed int prevRPId = -1;
0786 
0787     for (const auto &dit : task.geometry.getSensorMap()) {
0788       signed int rpId = CTPPSDetId(dit.first).rpId();
0789       if (rpId != prevRPId)
0790         printLineSeparator(results);
0791       prevRPId = rpId;
0792 
0793       printId(dit.first);
0794       printf(" ║");
0795 
0796       for (unsigned int q = 0; q < task.quantityClasses.size(); q++) {
0797         for (unsigned int a = 0; a < results.size(); a++) {
0798           const auto it = results[a].find(dit.first);
0799           if (it == results[a].end()) {
0800             if (algorithms[a]->hasErrorEstimate())
0801               printf("%19s", "----");
0802             else
0803               printf("%8s", "----");
0804 
0805             if (a + 1 == results.size())
0806               printf("║");
0807             else
0808               printf("│");
0809 
0810             continue;
0811           }
0812 
0813           const auto &ac = it->second;
0814           double v = 0., e = 0.;
0815           switch (task.quantityClasses[q]) {
0816             case AlignmentTask::qcShR1:
0817               v = ac.getShR1();
0818               e = ac.getShR1Unc();
0819               break;
0820             case AlignmentTask::qcShR2:
0821               v = ac.getShR2();
0822               e = ac.getShR2Unc();
0823               break;
0824             case AlignmentTask::qcShZ:
0825               v = ac.getShZ();
0826               e = ac.getShZUnc();
0827               break;
0828             case AlignmentTask::qcRotZ:
0829               v = ac.getRotZ();
0830               e = ac.getRotZUnc();
0831               break;
0832           }
0833 
0834           if (algorithms[a]->hasErrorEstimate())
0835             printf("%+8.1f ± %7.1f", v * 1E3, e * 1E3);
0836           else
0837             printf("%+8.1f", v * 1E3);
0838 
0839           if (a + 1 == results.size())
0840             printf("║");
0841           else
0842             printf("│");
0843         }
0844       }
0845 
0846       printf("\n");
0847     }
0848 
0849     printLineSeparator(results);
0850     printAlgorithmsLine(results);
0851     printQuantitiesLine(results);
0852     printLineSeparator(results);
0853   }
0854 
0855   // save results as alignment XML files
0856   for (unsigned int a = 0; a < results.size(); a++) {
0857     // convert readout-direction corrections to X and Y
0858     CTPPSRPAlignmentCorrectionsData convertedAlignments;
0859     for (const auto &p : results[a]) {
0860       const DetGeometry &d = task.geometry.get(p.first);
0861       const auto dir1 = d.getDirectionData(1);
0862       const auto dir2 = d.getDirectionData(2);
0863 
0864       const double det = dir1.dx * dir2.dy - dir1.dy * dir2.dx;
0865       const double sh_x = (dir2.dy * p.second.getShR1() - dir1.dy * p.second.getShR2()) / det;
0866       const double sh_y = (-dir2.dx * p.second.getShR1() + dir1.dx * p.second.getShR2()) / det;
0867 
0868       const double sh_x_e =
0869           sqrt(pow(dir2.dy / det * p.second.getShR1Unc(), 2) + pow(dir1.dy / det * p.second.getShR2Unc(), 2));
0870       const double sh_y_e =
0871           sqrt(pow(dir2.dx / det * p.second.getShR1Unc(), 2) + pow(dir1.dx / det * p.second.getShR2Unc(), 2));
0872 
0873       const CTPPSRPAlignmentCorrectionData corr(sh_x,
0874                                                 sh_x_e,
0875                                                 sh_y,
0876                                                 sh_y_e,
0877                                                 p.second.getShZ(),
0878                                                 p.second.getShZUnc(),
0879                                                 0.,
0880                                                 0.,
0881                                                 0.,
0882                                                 0.,
0883                                                 p.second.getRotZ(),
0884                                                 p.second.getRotZUnc());
0885       convertedAlignments.setSensorCorrection(p.first, corr);
0886     }
0887 
0888     // write non-cumulative alignments
0889     if (!fileNamePrefix.empty()) {
0890       CTPPSRPAlignmentCorrectionsMethods::writeToXML(convertedAlignments,
0891                                                      fileNamePrefix + algorithms[a]->getName() + ".xml",
0892                                                      preciseXMLFormat,
0893                                                      saveXMLUncertainties,
0894                                                      true,
0895                                                      true,
0896                                                      true,
0897                                                      true);
0898     }
0899 
0900     // merge alignments
0901     CTPPSRPAlignmentCorrectionsData cumulativeAlignments;
0902 
0903     cumulativeAlignments.addCorrections(initialAlignments, false, true, true);
0904     cumulativeAlignments.addCorrections(convertedAlignments, false, true, true);
0905 
0906     // write expanded and factored results
0907     if (!expandedFileNamePrefix.empty() || !factoredFileNamePrefix.empty()) {
0908       CTPPSRPAlignmentCorrectionsData expandedAlignments;
0909       CTPPSRPAlignmentCorrectionsData factoredAlignments;
0910 
0911       if (verbosity)
0912         printf(">> Factorizing results of %s algorithm\n", algorithms[a]->getName().c_str());
0913 
0914       const bool equalWeights = false;
0915       factorRPFromSensorCorrections(
0916           cumulativeAlignments, expandedAlignments, factoredAlignments, task.geometry, equalWeights, verbosity);
0917 
0918       if (!expandedFileNamePrefix.empty()) {
0919         CTPPSRPAlignmentCorrectionsMethods::writeToXML(expandedAlignments,
0920                                                        expandedFileNamePrefix + algorithms[a]->getName() + ".xml",
0921                                                        preciseXMLFormat,
0922                                                        saveXMLUncertainties,
0923                                                        true,
0924                                                        true,
0925                                                        true,
0926                                                        true);
0927       }
0928 
0929       if (!factoredFileNamePrefix.empty()) {
0930         CTPPSRPAlignmentCorrectionsMethods::writeToXML(factoredAlignments,
0931                                                        factoredFileNamePrefix + algorithms[a]->getName() + ".xml",
0932                                                        preciseXMLFormat,
0933                                                        saveXMLUncertainties,
0934                                                        true,
0935                                                        true,
0936                                                        true,
0937                                                        true);
0938       }
0939     }
0940   }
0941 
0942   // prepare algorithms for destructions
0943   for (const auto &algorithm : algorithms)
0944     algorithm->end();
0945 }
0946 
0947 //----------------------------------------------------------------------------------------------------
0948 
0949 string StraightTrackAlignment::setToString(const set<unsigned int> &s) {
0950   unsigned int N = s.size();
0951   if (N == 0)
0952     return "empty";
0953 
0954   string str;
0955   char buf[10];
0956   unsigned int i = 0;
0957   for (set<unsigned int>::iterator it = s.begin(); it != s.end(); ++it, ++i) {
0958     sprintf(buf, "%u", *it);
0959     str += buf;
0960     if (i < N - 1)
0961       str += ", ";
0962   }
0963 
0964   return str;
0965 }
0966 
0967 //----------------------------------------------------------------------------------------------------
0968 
0969 void StraightTrackAlignment::printN(const char *str, unsigned int N) {
0970   for (unsigned int i = 0; i < N; i++)
0971     printf("%s", str);
0972 }
0973 
0974 //----------------------------------------------------------------------------------------------------
0975 
0976 void StraightTrackAlignment::printLineSeparator(const std::vector<std::map<unsigned int, AlignmentResult> > &results) {
0977   printf("═════════════════════════╬");
0978   for (unsigned int q = 0; q < task.quantityClasses.size(); q++) {
0979     for (unsigned int a = 0; a < results.size(); a++) {
0980       printN("═", algorithms[a]->hasErrorEstimate() ? 18 : 8);
0981       if (a + 1 != results.size())
0982         printf("═");
0983     }
0984     printf("╬");
0985   }
0986   printf("\n");
0987 }
0988 
0989 //----------------------------------------------------------------------------------------------------
0990 
0991 void StraightTrackAlignment::printQuantitiesLine(const std::vector<std::map<unsigned int, AlignmentResult> > &results) {
0992   printf("                         ║");
0993 
0994   for (unsigned int q = 0; q < task.quantityClasses.size(); q++) {
0995     unsigned int size = 0;
0996     for (unsigned int a = 0; a < results.size(); a++)
0997       size += (algorithms[a]->hasErrorEstimate()) ? 18 : 8;
0998     size += algorithms.size() - 1;
0999 
1000     const string &tag = task.quantityClassTag(task.quantityClasses[q]);
1001     unsigned int space = (size - tag.size()) / 2;
1002     printN(" ", space);
1003     printf("%s", tag.c_str());
1004     printN(" ", size - space - tag.size());
1005     printf("║");
1006   }
1007   printf("\n");
1008 }
1009 
1010 //----------------------------------------------------------------------------------------------------
1011 
1012 void StraightTrackAlignment::printAlgorithmsLine(const std::vector<std::map<unsigned int, AlignmentResult> > &results) {
1013   printf("                         ║");
1014 
1015   for (unsigned int q = 0; q < task.quantityClasses.size(); q++) {
1016     for (unsigned int a = 0; a < results.size(); a++) {
1017       printf((algorithms[a]->hasErrorEstimate()) ? "%18s" : "%8s", algorithms[a]->getName().substr(0, 8).c_str());
1018 
1019       if (a + 1 == results.size())
1020         printf("║");
1021       else
1022         printf("│");
1023     }
1024   }
1025   printf("\n");
1026 }
1027 
1028 //----------------------------------------------------------------------------------------------------
1029 
1030 void StraightTrackAlignment::saveDiagnostics() const {
1031   if (diagnosticsFile.empty())
1032     return;
1033 
1034   TFile *df = new TFile(diagnosticsFile.c_str(), "recreate");
1035   if (df->IsZombie())
1036     throw cms::Exception("PPS") << "Cannot open file `" << diagnosticsFile << "' for writing.";
1037 
1038   if (buildDiagnosticPlots) {
1039     TDirectory *commonDir = df->mkdir("common");
1040     gDirectory = commonDir;
1041 
1042     fitNdfHist_fitted->Write();
1043     fitNdfHist_selected->Write();
1044     fitAxHist_fitted->Write();
1045     fitAyHist_fitted->Write();
1046     fitAxHist_selected->Write();
1047     fitAyHist_selected->Write();
1048     fitBxHist_fitted->Write();
1049     fitByHist_fitted->Write();
1050     fitBxHist_selected->Write();
1051     fitByHist_selected->Write();
1052     fitPHist_fitted->Write();
1053     fitPHist_selected->Write();
1054 
1055     gDirectory = commonDir->mkdir("plots global");
1056     globalPlots.write();
1057 
1058     TDirectory *ppsDir = commonDir->mkdir("plots per RP set");
1059     for (map<set<unsigned int>, RPSetPlots>::const_iterator it = rpSetPlots.begin(); it != rpSetPlots.end(); ++it) {
1060       gDirectory = ppsDir->mkdir(setToString(it->first).c_str());
1061 
1062       it->second.write();
1063     }
1064 
1065     TDirectory *resDir = commonDir->mkdir("residuals");
1066     for (map<unsigned int, ResiduaHistogramSet>::const_iterator it = residuaHistograms.begin();
1067          it != residuaHistograms.end();
1068          ++it) {
1069       char buf[10];
1070       sprintf(buf, "%u", it->first);
1071       gDirectory = resDir->mkdir(buf);
1072       it->second.total_fitted->Write();
1073       it->second.total_selected->Write();
1074       it->second.selected_vs_chiSq->Write();
1075 
1076       /*
1077       gDirectory = gDirectory->mkdir("fitted per RP set");
1078       for (map< set<unsigned int>, TH1D* >::iterator sit = it->second.perRPSet_fitted.begin();
1079           sit != it->second.perRPSet_fitted.end(); ++sit)
1080         sit->second->Write();
1081       gDirectory->cd("..");
1082 */
1083 
1084       gDirectory = gDirectory->mkdir("selected per RP set");
1085       TCanvas *c = new TCanvas;
1086       c->SetName("alltogether");
1087       unsigned int idx = 0;
1088       for (map<set<unsigned int>, TH1D *>::const_iterator sit = it->second.perRPSet_selected.begin();
1089            sit != it->second.perRPSet_selected.end();
1090            ++sit, ++idx) {
1091         sit->second->SetLineColor(idx + 1);
1092         sit->second->Draw((idx == 0) ? "" : "same");
1093         sit->second->Write();
1094       }
1095       c->Write();
1096     }
1097   }
1098 
1099   // save diagnostics of algorithms
1100   for (vector<AlignmentAlgorithm *>::const_iterator it = algorithms.begin(); it != algorithms.end(); ++it) {
1101     TDirectory *algDir = df->mkdir((*it)->getName().c_str());
1102     (*it)->saveDiagnostics(algDir);
1103   }
1104 
1105   delete df;
1106 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!