File indexing completed on 2024-04-06 11:57:23
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.;
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;
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
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;
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;
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
0368 #include "FWCore/Framework/interface/MakerMacros.h"
0369 DEFINE_FWK_MODULE(SurveyInputCSCfromPins);
