Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-06 09:39:45

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "Alignment/CommonAlignment/interface/SurveyDet.h"
0003 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
0004 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0007 #include "DataFormats/DetId/interface/DetId.h"
0008 
0009 #include "Alignment/MuonAlignment/interface/AlignableCSCChamber.h"
0010 #include "Alignment/MuonAlignment/interface/AlignableCSCStation.h"
0011 
0012 #include <iostream>
0013 #include <fstream>
0014 #include "TFile.h"
0015 #include "TTree.h"
0016 #include "TRandom3.h"
0017 
0018 #include "SurveyInputCSCfromPins.h"
0019 
0020 #define SQR(x) ((x) * (x))
0021 
0022 SurveyInputCSCfromPins::SurveyInputCSCfromPins(const edm::ParameterSet &cfg)
0023     : DTGeoToken_(esConsumes()),
0024       CSCGeoToken_(esConsumes()),
0025       GEMGeoToken_(esConsumes()),
0026       m_pinPositions(cfg.getParameter<std::string>("pinPositions")),
0027       m_rootFile(cfg.getParameter<std::string>("rootFile")),
0028       m_verbose(cfg.getParameter<bool>("verbose")),
0029       m_errorX(cfg.getParameter<double>("errorX")),
0030       m_errorY(cfg.getParameter<double>("errorY")),
0031       m_errorZ(cfg.getParameter<double>("errorZ")),
0032       m_missingErrorTranslation(cfg.getParameter<double>("missingErrorTranslation")),
0033       m_missingErrorAngle(cfg.getParameter<double>("missingErrorAngle")),
0034       m_stationErrorX(cfg.getParameter<double>("stationErrorX")),
0035       m_stationErrorY(cfg.getParameter<double>("stationErrorY")),
0036       m_stationErrorZ(cfg.getParameter<double>("stationErrorZ")),
0037       m_stationErrorPhiX(cfg.getParameter<double>("stationErrorPhiX")),
0038       m_stationErrorPhiY(cfg.getParameter<double>("stationErrorPhiY")),
0039       m_stationErrorPhiZ(cfg.getParameter<double>("stationErrorPhiZ")) {}
0040 
0041 void SurveyInputCSCfromPins::orient(align::LocalVector LC1,
0042                                     align::LocalVector LC2,
0043                                     double a,
0044                                     double b,
0045                                     double &T,
0046                                     double &dx,
0047                                     double &dy,
0048                                     double &dz,
0049                                     double &PhX,
0050                                     double &PhZ) {
0051   double cosPhX, sinPhX, cosPhZ, sinPhZ;
0052 
0053   LocalPoint LP1(LC1.x(), LC1.y() + a, LC1.z() + b);
0054   LocalPoint LP2(LC2.x(), LC2.y() - a, LC2.z() + b);
0055 
0056   LocalPoint P((LP1.x() - LP2.x()) / (2. * a), (LP1.y() - LP2.y()) / (2. * a), (LP1.z() - LP2.z()) / (2. * a));
0057   LocalPoint Pp((LP1.x() + LP2.x()) / (2.), (LP1.y() + LP2.y()) / (2.), (LP1.z() + LP2.z()) / (2.));
0058 
0059   T = P.mag();
0060 
0061   sinPhX = P.z() / T;
0062   cosPhX = sqrt(1 - SQR(sinPhX));
0063   cosPhZ = P.y() / (T * cosPhX);
0064   sinPhZ = -P.x() / (T * cosPhZ);
0065 
0066   PhX = atan2(sinPhX, cosPhX);
0067 
0068   PhZ = atan2(sinPhZ, cosPhZ);
0069 
0070   dx = Pp.x() - sinPhZ * sinPhX * b;
0071   dy = Pp.y() + cosPhZ * sinPhX * b;
0072   dz = Pp.z() - cosPhX * b;
0073 }
0074 
0075 void SurveyInputCSCfromPins::errors(double a,
0076                                     double b,
0077                                     bool missing1,
0078                                     bool missing2,
0079                                     double &dx_dx,
0080                                     double &dy_dy,
0081                                     double &dz_dz,
0082                                     double &phix_phix,
0083                                     double &phiz_phiz,
0084                                     double &dy_phix) {
0085   dx_dx = 0.;
0086   dy_dy = 0.;
0087   dz_dz = 0.;
0088   phix_phix = 0.;
0089   phiz_phiz = 0.;
0090   dy_phix = 0.;
0091 
0092   const double trials = 10000.;  // two significant digits
0093 
0094   for (int i = 0; i < trials; i++) {
0095     LocalVector LC1, LC2;
0096     if (missing1) {
0097       LC1 = LocalVector(gRandom->Gaus(0., m_missingErrorTranslation),
0098                         gRandom->Gaus(0., m_missingErrorTranslation),
0099                         gRandom->Gaus(0., m_missingErrorTranslation));
0100     } else {
0101       LC1 = LocalVector(gRandom->Gaus(0., m_errorX), gRandom->Gaus(0., m_errorY), gRandom->Gaus(0., m_errorZ));
0102     }
0103 
0104     if (missing2) {
0105       LC2 = LocalVector(gRandom->Gaus(0., m_missingErrorTranslation),
0106                         gRandom->Gaus(0., m_missingErrorTranslation),
0107                         gRandom->Gaus(0., m_missingErrorTranslation));
0108     } else {
0109       LC2 = LocalVector(gRandom->Gaus(0., m_errorX), gRandom->Gaus(0., m_errorY), gRandom->Gaus(0., m_errorZ));
0110     }
0111 
0112     double dx, dy, dz, PhX, PhZ, T;
0113     orient(LC1, LC2, a, b, T, dx, dy, dz, PhX, PhZ);
0114 
0115     dx_dx += dx * dx;
0116     dy_dy += dy * dy;
0117     dz_dz += dz * dz;
0118     phix_phix += PhX * PhX;
0119     phiz_phiz += PhZ * PhZ;
0120     dy_phix += dy * PhX;  // the only non-zero off-diagonal element
0121   }
0122 
0123   dx_dx /= trials;
0124   dy_dy /= trials;
0125   dz_dz /= trials;
0126   phix_phix /= trials;
0127   phiz_phiz /= trials;
0128   dy_phix /= trials;
0129 }
0130 
0131 void SurveyInputCSCfromPins::analyze(const edm::Event &, const edm::EventSetup &iSetup) {
0132   if (theFirstEvent) {
0133     edm::LogInfo("SurveyInputCSCfromPins") << "***************ENTERING INITIALIZATION******************"
0134                                            << "  \n";
0135 
0136     std::ifstream in;
0137     in.open(m_pinPositions.c_str());
0138 
0139     Double_t x1, y1, z1, x2, y2, z2, a, b, tot = 0.0, maxErr = 0.0, h, s1, dx, dy, dz, PhX, PhZ, T;
0140 
0141     int ID1, ID2, ID3, ID4, ID5, i = 1, ii = 0;
0142 
0143     TFile *file1 = new TFile(m_rootFile.c_str(), "recreate");
0144     TTree *tree1 = new TTree("tree1", "alignment pins");
0145 
0146     if (m_verbose) {
0147       tree1->Branch("displacement_x_pin1_cm", &x1, "x1/D");
0148       tree1->Branch("displacement_y_pin1_cm", &y1, "y1/D");
0149       tree1->Branch("displacement_z_pin1_cm", &z1, "z1/D");
0150       tree1->Branch("displacement_x_pin2_cm", &x2, "x2/D");
0151       tree1->Branch("displacement_y_pin2_cm", &y2, "y2/D");
0152       tree1->Branch("displacement_z_pin2_cm", &z2, "z2/D");
0153       tree1->Branch("error_vector_length_cm", &h, "h/D");
0154       tree1->Branch("stretch_diff_cm", &s1, "s1/D");
0155       tree1->Branch("stretch_factor", &T, "T/D");
0156       tree1->Branch("chamber_displacement_x_cm", &dx, "dx/D");
0157       tree1->Branch("chamber_displacement_y_cm", &dy, "dy/D");
0158       tree1->Branch("chamber_displacement_z_cm", &dz, "dz/D");
0159       tree1->Branch("chamber_rotation_x_rad", &PhX, "PhX/D");
0160       tree1->Branch("chamber_rotation_z_rad", &PhZ, "PhZ/D");
0161     }
0162 
0163     const DTGeometry *dtGeometry = &iSetup.getData(DTGeoToken_);
0164     const CSCGeometry *cscGeometry = &iSetup.getData(CSCGeoToken_);
0165     const GEMGeometry *gemGeometry = &iSetup.getData(GEMGeoToken_);
0166 
0167     AlignableMuon *theAlignableMuon = new AlignableMuon(dtGeometry, cscGeometry, gemGeometry);
0168     AlignableNavigator *theAlignableNavigator = new AlignableNavigator(theAlignableMuon);
0169 
0170     const auto &theEndcaps = theAlignableMuon->CSCEndcaps();
0171 
0172     for (const auto &aliiter : theEndcaps) {
0173       addComponent(aliiter);
0174     }
0175 
0176     while (in.good()) {
0177       bool missing1 = false;
0178       bool missing2 = false;
0179 
0180       in >> ID1 >> ID2 >> ID3 >> ID4 >> ID5 >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> a >> b;
0181 
0182       if (fabs(x1 - 1000.) < 1e5 && fabs(y1 - 1000.) < 1e5 && fabs(z1 - 1000.) < 1e5) {
0183         missing1 = true;
0184         x1 = x2;
0185         y1 = y2;
0186         z1 = z2;
0187       }
0188 
0189       if (fabs(x2 - 1000.) < 1e5 && fabs(y2 - 1000.) < 1e5 && fabs(z2 - 1000.) < 1e5) {
0190         missing2 = true;
0191         x2 = x1;
0192         y2 = y1;
0193         z2 = z1;
0194       }
0195 
0196       x1 = x1 / 10.0;
0197       y1 = y1 / 10.0;
0198       z1 = z1 / 10.0;
0199       x2 = x2 / 10.0;
0200       y2 = y2 / 10.0;
0201       z2 = z2 / 10.0;
0202 
0203       CSCDetId layerID(ID1, ID2, ID3, ID4, 1);
0204 
0205       // We cannot use chamber ID (when ID5=0), because AlignableNavigator gives the error (aliDet and aliDetUnit are undefined for chambers)
0206 
0207       Alignable *theAlignable1 = theAlignableNavigator->alignableFromDetId(layerID);
0208       Alignable *chamberAli = theAlignable1->mother();
0209 
0210       LocalVector LC1 = chamberAli->surface().toLocal(GlobalVector(x1, y1, z1));
0211       LocalVector LC2 = chamberAli->surface().toLocal(GlobalVector(x2, y2, z2));
0212 
0213       orient(LC1, LC2, a, b, T, dx, dy, dz, PhX, PhZ);
0214 
0215       GlobalPoint PG1 = chamberAli->surface().toGlobal(LocalPoint(LC1.x(), LC1.y() + a, LC1.z() + b));
0216       chamberAli->surface().toGlobal(LocalPoint(LC2.x(), LC2.y() - a, LC2.z() + b));
0217 
0218       LocalVector lvector(dx, dy, dz);
0219       GlobalVector gvector = (chamberAli->surface()).toGlobal(lvector);
0220       chamberAli->move(gvector);
0221 
0222       chamberAli->rotateAroundLocalX(PhX);
0223       chamberAli->rotateAroundLocalZ(PhZ);
0224 
0225       double dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix;
0226       errors(a, b, missing1, missing2, dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix);
0227       align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
0228       error(0, 0) = dx_dx;
0229       error(1, 1) = dy_dy;
0230       error(2, 2) = dz_dz;
0231       error(3, 3) = phix_phix;
0232       error(4, 4) = m_missingErrorAngle;
0233       error(5, 5) = phiz_phiz;
0234       error(1, 3) = dy_phix;
0235       error(3, 1) = dy_phix;  // just in case
0236 
0237       chamberAli->setSurvey(new SurveyDet(chamberAli->surface(), error));
0238 
0239       if (m_verbose) {
0240         edm::LogInfo("SurveyInputCSCfromPins") << " survey information = " << chamberAli->survey() << "  \n";
0241 
0242         LocalPoint LP1n = chamberAli->surface().toLocal(PG1);
0243 
0244         LocalPoint hiP(LP1n.x(), LP1n.y() - a * T, LP1n.z() - b);
0245 
0246         h = hiP.mag();
0247         s1 = LP1n.y() - a;
0248 
0249         if (h > maxErr) {
0250           maxErr = h;
0251 
0252           ii = i;
0253         }
0254 
0255         edm::LogInfo("SurveyInputCSCfromPins")
0256             << "  \n"
0257             << "i " << i++ << " " << ID1 << " " << ID2 << " " << ID3 << " " << ID4 << " " << ID5 << " error  " << h
0258             << " \n"
0259             << " x1 " << x1 << " y1 " << y1 << " z1 " << z1 << " x2 " << x2 << " y2 " << y2 << " z2 " << z2 << " \n"
0260             << " error " << h << " S1 " << s1 << " \n"
0261             << " dx " << dx << " dy " << dy << " dz " << dz << " PhX " << PhX << " PhZ " << PhZ << " \n";
0262 
0263         tot += h;
0264 
0265         tree1->Fill();
0266       }
0267     }
0268 
0269     in.close();
0270 
0271     if (m_verbose) {
0272       file1->Write();
0273       edm::LogInfo("SurveyInputCSCfromPins")
0274           << " Total error  " << tot << "  Max Error " << maxErr << " N " << ii << "  \n";
0275     }
0276 
0277     file1->Close();
0278 
0279     for (const auto &aliiter : theEndcaps) {
0280       fillAllRecords(aliiter);
0281     }
0282 
0283     delete theAlignableMuon;
0284     delete theAlignableNavigator;
0285 
0286     edm::LogInfo("SurveyInputCSCfromPins") << "*************END INITIALIZATION***************"
0287                                            << "  \n";
0288 
0289     theFirstEvent = false;
0290   }
0291 }
0292 
0293 void SurveyInputCSCfromPins::fillAllRecords(Alignable *ali) {
0294   if (ali->survey() == nullptr) {
0295     AlignableCSCChamber *ali_AlignableCSCChamber = dynamic_cast<AlignableCSCChamber *>(ali);
0296     AlignableCSCStation *ali_AlignableCSCStation = dynamic_cast<AlignableCSCStation *>(ali);
0297 
0298     if (ali_AlignableCSCChamber != nullptr) {
0299       CSCDetId detid(ali->geomDetId());
0300       if (abs(detid.station()) == 1 && (detid.ring() == 1 || detid.ring() == 4)) {
0301         align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
0302         error(0, 0) = m_missingErrorTranslation;
0303         error(1, 1) = m_missingErrorTranslation;
0304         error(2, 2) = m_missingErrorTranslation;
0305         error(3, 3) = m_missingErrorAngle;
0306         error(4, 4) = m_missingErrorAngle;
0307         error(5, 5) = m_missingErrorAngle;
0308 
0309         ali->setSurvey(new SurveyDet(ali->surface(), error));
0310       } else {
0311         double a = 100.;
0312         double b = -9.4034;
0313         if (abs(detid.station()) == 1 && detid.ring() == 2)
0314           a = -90.260;
0315         else if (abs(detid.station()) == 1 && detid.ring() == 3)
0316           a = -85.205;
0317         else if (abs(detid.station()) == 2 && detid.ring() == 1)
0318           a = -97.855;
0319         else if (abs(detid.station()) == 2 && detid.ring() == 2)
0320           a = -164.555;
0321         else if (abs(detid.station()) == 3 && detid.ring() == 1)
0322           a = -87.870;
0323         else if (abs(detid.station()) == 3 && detid.ring() == 2)
0324           a = -164.555;
0325         else if (abs(detid.station()) == 4 && detid.ring() == 1)
0326           a = -77.890;
0327 
0328         double dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix;
0329         errors(a, b, true, true, dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix);
0330         align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
0331         error(0, 0) = dx_dx;
0332         error(1, 1) = dy_dy;
0333         error(2, 2) = dz_dz;
0334         error(3, 3) = phix_phix;
0335         error(4, 4) = m_missingErrorAngle;
0336         error(5, 5) = phiz_phiz;
0337         error(1, 3) = dy_phix;
0338         error(3, 1) = dy_phix;  // just in case
0339 
0340         ali->setSurvey(new SurveyDet(ali->surface(), error));
0341       }
0342     }
0343 
0344     else if (ali_AlignableCSCStation != nullptr) {
0345       align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
0346       error(0, 0) = m_stationErrorX;
0347       error(1, 1) = m_stationErrorY;
0348       error(2, 2) = m_stationErrorZ;
0349       error(3, 3) = m_stationErrorPhiX;
0350       error(4, 4) = m_stationErrorPhiY;
0351       error(5, 5) = m_stationErrorPhiZ;
0352 
0353       ali->setSurvey(new SurveyDet(ali->surface(), error));
0354     }
0355 
0356     else {
0357       align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
0358       ali->setSurvey(new SurveyDet(ali->surface(), error * (1e-10)));
0359     }
0360   }
0361 
0362   for (const auto &iter : ali->components()) {
0363     fillAllRecords(iter);
0364   }
0365 }
0366 
0367 // Plug in to framework
0368 #include "FWCore/Framework/interface/MakerMacros.h"
0369 DEFINE_FWK_MODULE(SurveyInputCSCfromPins);