File indexing completed on 2021-10-22 02:30:52
0001 #include "CondFormats/Alignment/interface/Alignments.h"
0002 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0003 #include "CLHEP/Vector/RotationInterfaces.h"
0004
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
0011
0012 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
0013 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0014 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0015
0016 #include "Alignment/CommonAlignment/interface/Utilities.h"
0017 #include "Alignment/CommonAlignment/interface/SurveyDet.h"
0018 #include "Alignment/CommonAlignment/interface/Alignable.h"
0019 #include "DataFormats/DetId/interface/DetId.h"
0020 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0021 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
0022 #include "CondFormats/Alignment/interface/DetectorGlobalPosition.h"
0023
0024 #include "Alignment/MuonAlignment/interface/MuonAlignmentInputXML.h"
0025 #include "Alignment/MuonAlignment/interface/MuonAlignment.h"
0026
0027 #include "MuonGeometryArrange.h"
0028 #include "TFile.h"
0029 #include "TLatex.h"
0030 #include "TArrow.h"
0031 #include "TGraph.h"
0032 #include "TH1F.h"
0033 #include "TH2F.h"
0034 #include "CLHEP/Vector/ThreeVector.h"
0035
0036
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038
0039 #include <iostream>
0040 #include <fstream>
0041
0042 MuonGeometryArrange::MuonGeometryArrange(const edm::ParameterSet& cfg)
0043 : theSurveyIndex(0),
0044 _levelStrings(cfg.getUntrackedParameter<std::vector<std::string> >("levels")),
0045 _writeToDB(false),
0046 _commonMuonLevel(align::invalid),
0047 firstEvent_(true),
0048 idealInputLabel1("MuonGeometryArrangeLabel1"),
0049 idealInputLabel2("MuonGeometryArrangeLabel2"),
0050 idealInputLabel2a("MuonGeometryArrangeLabel2a"),
0051 geomIdeal("MuonGeometryArrangeGeomIdeal"),
0052 dtGeomToken1_(esConsumes(edm::ESInputTag("", idealInputLabel1))),
0053 cscGeomToken1_(esConsumes(edm::ESInputTag("", idealInputLabel1))),
0054 gemGeomToken1_(esConsumes(edm::ESInputTag("", idealInputLabel1))),
0055 dtGeomToken2_(esConsumes(edm::ESInputTag("", idealInputLabel2))),
0056 cscGeomToken2_(esConsumes(edm::ESInputTag("", idealInputLabel2))),
0057 gemGeomToken2_(esConsumes(edm::ESInputTag("", idealInputLabel2))),
0058 dtGeomToken3_(esConsumes(edm::ESInputTag("", idealInputLabel2a))),
0059 cscGeomToken3_(esConsumes(edm::ESInputTag("", idealInputLabel2a))),
0060 gemGeomToken3_(esConsumes(edm::ESInputTag("", idealInputLabel2a))),
0061 dtGeomIdealToken_(esConsumes(edm::ESInputTag("", geomIdeal))),
0062 cscGeomIdealToken_(esConsumes(edm::ESInputTag("", geomIdeal))),
0063 gemGeomIdealToken_(esConsumes(edm::ESInputTag("", geomIdeal))) {
0064 referenceMuon = nullptr;
0065 currentMuon = nullptr;
0066
0067 _inputXMLCurrent = cfg.getUntrackedParameter<std::string>("inputXMLCurrent");
0068 _inputXMLReference = cfg.getUntrackedParameter<std::string>("inputXMLReference");
0069
0070
0071 _inputFilename1 = cfg.getUntrackedParameter<std::string>("inputROOTFile1");
0072 _inputFilename2 = cfg.getUntrackedParameter<std::string>("inputROOTFile2");
0073 _inputTreename = cfg.getUntrackedParameter<std::string>("treeName");
0074
0075
0076 _filename = cfg.getUntrackedParameter<std::string>("outputFile");
0077
0078 _weightBy = cfg.getUntrackedParameter<std::string>("weightBy");
0079 _detIdFlag = cfg.getUntrackedParameter<bool>("detIdFlag");
0080 _detIdFlagFile = cfg.getUntrackedParameter<std::string>("detIdFlagFile");
0081 _weightById = cfg.getUntrackedParameter<bool>("weightById");
0082 _weightByIdFile = cfg.getUntrackedParameter<std::string>("weightByIdFile");
0083 _endcap = cfg.getUntrackedParameter<int>("endcapNumber");
0084 _station = cfg.getUntrackedParameter<int>("stationNumber");
0085 _ring = cfg.getUntrackedParameter<int>("ringNumber");
0086
0087
0088 if (_detIdFlag) {
0089 std::ifstream fin;
0090 fin.open(_detIdFlagFile.c_str());
0091
0092 while (!fin.eof() && fin.good()) {
0093 uint32_t id;
0094 fin >> id;
0095 _detIdFlagVector.push_back(id);
0096 }
0097 fin.close();
0098 }
0099
0100
0101 unsigned int lastID = 999999999;
0102 if (_weightById) {
0103 std::ifstream inFile;
0104 inFile.open(_weightByIdFile.c_str());
0105 int ctr = 0;
0106 while (!inFile.eof()) {
0107 ctr++;
0108 unsigned int listId;
0109 inFile >> listId;
0110 inFile.ignore(256, '\n');
0111 if (listId != lastID) {
0112 _weightByIdVector.push_back(listId);
0113 }
0114 lastID = listId;
0115 }
0116 inFile.close();
0117 }
0118
0119
0120 _theFile = new TFile(_filename.c_str(), "RECREATE");
0121 _alignTree = new TTree("alignTree", "alignTree");
0122 _alignTree->Branch("id", &_id, "id/I");
0123 _alignTree->Branch("level", &_level, "level/I");
0124 _alignTree->Branch("mid", &_mid, "mid/I");
0125 _alignTree->Branch("mlevel", &_mlevel, "mlevel/I");
0126 _alignTree->Branch("sublevel", &_sublevel, "sublevel/I");
0127 _alignTree->Branch("x", &_xVal, "x/F");
0128 _alignTree->Branch("y", &_yVal, "y/F");
0129 _alignTree->Branch("z", &_zVal, "z/F");
0130 _alignTree->Branch("r", &_rVal, "r/F");
0131 _alignTree->Branch("phi", &_phiVal, "phi/F");
0132 _alignTree->Branch("eta", &_etaVal, "eta/F");
0133 _alignTree->Branch("alpha", &_alphaVal, "alpha/F");
0134 _alignTree->Branch("beta", &_betaVal, "beta/F");
0135 _alignTree->Branch("gamma", &_gammaVal, "gamma/F");
0136 _alignTree->Branch("dx", &_dxVal, "dx/F");
0137 _alignTree->Branch("dy", &_dyVal, "dy/F");
0138 _alignTree->Branch("dz", &_dzVal, "dz/F");
0139 _alignTree->Branch("dr", &_drVal, "dr/F");
0140 _alignTree->Branch("dphi", &_dphiVal, "dphi/F");
0141 _alignTree->Branch("dalpha", &_dalphaVal, "dalpha/F");
0142 _alignTree->Branch("dbeta", &_dbetaVal, "dbeta/F");
0143 _alignTree->Branch("dgamma", &_dgammaVal, "dgamma/F");
0144 _alignTree->Branch("ldx", &_ldxVal, "ldx/F");
0145 _alignTree->Branch("ldy", &_ldyVal, "ldy/F");
0146 _alignTree->Branch("ldz", &_ldzVal, "ldz/F");
0147 _alignTree->Branch("ldr", &_ldrVal, "ldr/F");
0148 _alignTree->Branch("ldphi", &_ldphiVal, "ldphi/F");
0149 _alignTree->Branch("useDetId", &_useDetId, "useDetId/I");
0150 _alignTree->Branch("detDim", &_detDim, "detDim/I");
0151 _alignTree->Branch("rotx", &_rotxVal, "rotx/F");
0152 _alignTree->Branch("roty", &_rotyVal, "roty/F");
0153 _alignTree->Branch("rotz", &_rotzVal, "rotz/F");
0154 _alignTree->Branch("drotx", &_drotxVal, "drotx/F");
0155 _alignTree->Branch("droty", &_drotyVal, "droty/F");
0156 _alignTree->Branch("drotz", &_drotzVal, "drotz/F");
0157 _alignTree->Branch("surW", &_surWidth, "surW/F");
0158 _alignTree->Branch("surL", &_surLength, "surL/F");
0159 _alignTree->Branch("surRot", &_surRot, "surRot[9]/D");
0160
0161 _mgacollection.clear();
0162 }
0163
0164 void MuonGeometryArrange::endHist() {
0165
0166
0167 int size = _mgacollection.size();
0168 if (size <= 0)
0169 return;
0170 std::vector<float> xp(size + 1);
0171 std::vector<float> yp(size + 1);
0172 int i;
0173 float minV, maxV;
0174 int minI, maxI;
0175
0176 minV = 99999999.;
0177 maxV = -minV;
0178 minI = 9999999;
0179 maxI = -minI;
0180 TGraph* grx = nullptr;
0181 TH2F* dxh = nullptr;
0182
0183
0184 for (i = 0; i < size; i++) {
0185 if (_mgacollection[i].phipos < minI)
0186 minI = _mgacollection[i].phipos;
0187 if (_mgacollection[i].phipos > maxI)
0188 maxI = _mgacollection[i].phipos;
0189 xp[i] = _mgacollection[i].phipos;
0190 }
0191 if (minI >= maxI)
0192 return;
0193 xp[size] = xp[size - 1] + 1;
0194
0195 if (1 < minI)
0196 minI = 1;
0197 if (size > maxI)
0198 maxI = size;
0199 maxI++;
0200 int sizeI = maxI + 1 - minI;
0201 float smi = minI - 1;
0202 float sma = maxI + 1;
0203
0204
0205
0206 for (i = 0; i < size; i++) {
0207 if (_mgacollection[i].ldx < minV)
0208 minV = _mgacollection[i].ldx;
0209 if (_mgacollection[i].ldx > maxV)
0210 maxV = _mgacollection[i].ldx;
0211 yp[i] = _mgacollection[i].ldx;
0212 }
0213 yp[size] = yp[0];
0214
0215 makeGraph(sizeI,
0216 smi,
0217 sma,
0218 minV,
0219 maxV,
0220 dxh,
0221 grx,
0222 "delX_vs_position",
0223 "Local #delta X vs position",
0224 "GdelX_vs_position",
0225 "#delta x in cm",
0226 xp.data(),
0227 yp.data(),
0228 size);
0229
0230 minV = 99999999.;
0231 maxV = -minV;
0232 for (i = 0; i < size; i++) {
0233 if (_mgacollection[i].ldy < minV)
0234 minV = _mgacollection[i].ldy;
0235 if (_mgacollection[i].ldy > maxV)
0236 maxV = _mgacollection[i].ldy;
0237 yp[i] = _mgacollection[i].ldy;
0238 }
0239 yp[size] = yp[0];
0240
0241 makeGraph(sizeI,
0242 smi,
0243 sma,
0244 minV,
0245 maxV,
0246 dxh,
0247 grx,
0248 "delY_vs_position",
0249 "Local #delta Y vs position",
0250 "GdelY_vs_position",
0251 "#delta y in cm",
0252 xp.data(),
0253 yp.data(),
0254 size);
0255
0256
0257 minV = 99999999.;
0258 maxV = -minV;
0259 for (i = 0; i < size; i++) {
0260 if (_mgacollection[i].dz < minV)
0261 minV = _mgacollection[i].dz;
0262 if (_mgacollection[i].dz > maxV)
0263 maxV = _mgacollection[i].dz;
0264 yp[i] = _mgacollection[i].dz;
0265 }
0266 yp[size] = yp[0];
0267
0268 makeGraph(sizeI,
0269 smi,
0270 sma,
0271 minV,
0272 maxV,
0273 dxh,
0274 grx,
0275 "delZ_vs_position",
0276 "Local #delta Z vs position",
0277 "GdelZ_vs_position",
0278 "#delta z in cm",
0279 xp.data(),
0280 yp.data(),
0281 size);
0282
0283
0284 minV = 99999999.;
0285 maxV = -minV;
0286 for (i = 0; i < size; i++) {
0287 if (_mgacollection[i].dphi < minV)
0288 minV = _mgacollection[i].dphi;
0289 if (_mgacollection[i].dphi > maxV)
0290 maxV = _mgacollection[i].dphi;
0291 yp[i] = _mgacollection[i].dphi;
0292 }
0293 yp[size] = yp[0];
0294
0295 makeGraph(sizeI,
0296 smi,
0297 sma,
0298 minV,
0299 maxV,
0300 dxh,
0301 grx,
0302 "delphi_vs_position",
0303 "#delta #phi vs position",
0304 "Gdelphi_vs_position",
0305 "#delta #phi in radians",
0306 xp.data(),
0307 yp.data(),
0308 size);
0309
0310
0311 minV = 99999999.;
0312 maxV = -minV;
0313 for (i = 0; i < size; i++) {
0314 if (_mgacollection[i].dr < minV)
0315 minV = _mgacollection[i].dr;
0316 if (_mgacollection[i].dr > maxV)
0317 maxV = _mgacollection[i].dr;
0318 yp[i] = _mgacollection[i].dr;
0319 }
0320 yp[size] = yp[0];
0321
0322 makeGraph(sizeI,
0323 smi,
0324 sma,
0325 minV,
0326 maxV,
0327 dxh,
0328 grx,
0329 "delR_vs_position",
0330 "#delta R vs position",
0331 "GdelR_vs_position",
0332 "#delta R in cm",
0333 xp.data(),
0334 yp.data(),
0335 size);
0336
0337
0338 minV = 99999999.;
0339 maxV = -minV;
0340 for (i = 0; i < size; i++) {
0341 float ttemp = _mgacollection[i].r * _mgacollection[i].dphi;
0342 if (ttemp < minV)
0343 minV = ttemp;
0344 if (ttemp > maxV)
0345 maxV = ttemp;
0346 yp[i] = ttemp;
0347 }
0348 yp[size] = yp[0];
0349
0350 makeGraph(sizeI,
0351 smi,
0352 sma,
0353 minV,
0354 maxV,
0355 dxh,
0356 grx,
0357 "delRphi_vs_position",
0358 "R #delta #phi vs position",
0359 "GdelRphi_vs_position",
0360 "R #delta #phi in cm",
0361 xp.data(),
0362 yp.data(),
0363 size);
0364
0365
0366 minV = 99999999.;
0367 maxV = -minV;
0368 for (i = 0; i < size; i++) {
0369 if (_mgacollection[i].dalpha < minV)
0370 minV = _mgacollection[i].dalpha;
0371 if (_mgacollection[i].dalpha > maxV)
0372 maxV = _mgacollection[i].dalpha;
0373 yp[i] = _mgacollection[i].dalpha;
0374 }
0375 yp[size] = yp[0];
0376
0377 makeGraph(sizeI,
0378 smi,
0379 sma,
0380 minV,
0381 maxV,
0382 dxh,
0383 grx,
0384 "delalpha_vs_position",
0385 "#delta #alpha vs position",
0386 "Gdelalpha_vs_position",
0387 "#delta #alpha in rad",
0388 xp.data(),
0389 yp.data(),
0390 size);
0391
0392
0393 minV = 99999999.;
0394 maxV = -minV;
0395 for (i = 0; i < size; i++) {
0396 if (_mgacollection[i].dbeta < minV)
0397 minV = _mgacollection[i].dbeta;
0398 if (_mgacollection[i].dbeta > maxV)
0399 maxV = _mgacollection[i].dbeta;
0400 yp[i] = _mgacollection[i].dbeta;
0401 }
0402 yp[size] = yp[0];
0403
0404 makeGraph(sizeI,
0405 smi,
0406 sma,
0407 minV,
0408 maxV,
0409 dxh,
0410 grx,
0411 "delbeta_vs_position",
0412 "#delta #beta vs position",
0413 "Gdelbeta_vs_position",
0414 "#delta #beta in rad",
0415 xp.data(),
0416 yp.data(),
0417 size);
0418
0419
0420 minV = 99999999.;
0421 maxV = -minV;
0422 for (i = 0; i < size; i++) {
0423 if (_mgacollection[i].dgamma < minV)
0424 minV = _mgacollection[i].dgamma;
0425 if (_mgacollection[i].dgamma > maxV)
0426 maxV = _mgacollection[i].dgamma;
0427 yp[i] = _mgacollection[i].dgamma;
0428 }
0429 yp[size] = yp[0];
0430
0431 makeGraph(sizeI,
0432 smi,
0433 sma,
0434 minV,
0435 maxV,
0436 dxh,
0437 grx,
0438 "delgamma_vs_position",
0439 "#delta #gamma vs position",
0440 "Gdelgamma_vs_position",
0441 "#delta #gamma in rad",
0442 xp.data(),
0443 yp.data(),
0444 size);
0445
0446
0447 minV = 99999999.;
0448 maxV = -minV;
0449 for (i = 0; i < size; i++) {
0450 if (_mgacollection[i].drotx < minV)
0451 minV = _mgacollection[i].drotx;
0452 if (_mgacollection[i].drotx > maxV)
0453 maxV = _mgacollection[i].drotx;
0454 yp[i] = _mgacollection[i].drotx;
0455 }
0456 yp[size] = yp[0];
0457
0458 makeGraph(sizeI,
0459 smi,
0460 sma,
0461 minV,
0462 maxV,
0463 dxh,
0464 grx,
0465 "delrotX_vs_position",
0466 "#delta rotX vs position",
0467 "GdelrotX_vs_position",
0468 "#delta rotX in rad",
0469 xp.data(),
0470 yp.data(),
0471 size);
0472
0473
0474 minV = 99999999.;
0475 maxV = -minV;
0476 for (i = 0; i < size; i++) {
0477 if (_mgacollection[i].droty < minV)
0478 minV = _mgacollection[i].droty;
0479 if (_mgacollection[i].droty > maxV)
0480 maxV = _mgacollection[i].droty;
0481 yp[i] = _mgacollection[i].droty;
0482 }
0483 yp[size] = yp[0];
0484
0485 makeGraph(sizeI,
0486 smi,
0487 sma,
0488 minV,
0489 maxV,
0490 dxh,
0491 grx,
0492 "delrotY_vs_position",
0493 "#delta rotY vs position",
0494 "GdelrotY_vs_position",
0495 "#delta rotY in rad",
0496 xp.data(),
0497 yp.data(),
0498 size);
0499
0500
0501 minV = 99999999.;
0502 maxV = -minV;
0503 for (i = 0; i < size; i++) {
0504 if (_mgacollection[i].drotz < minV)
0505 minV = _mgacollection[i].drotz;
0506 if (_mgacollection[i].drotz > maxV)
0507 maxV = _mgacollection[i].drotz;
0508 yp[i] = _mgacollection[i].drotz;
0509 }
0510 yp[size] = yp[0];
0511
0512 makeGraph(sizeI,
0513 smi,
0514 sma,
0515 minV,
0516 maxV,
0517 dxh,
0518 grx,
0519 "delrotZ_vs_position",
0520 "#delta rotZ vs position",
0521 "GdelrotZ_vs_position",
0522 "#delta rotZ in rad",
0523 xp.data(),
0524 yp.data(),
0525 size);
0526
0527
0528
0529
0530 maxV = -99999999.;
0531 float ttemp, rtemp;
0532 float maxR = -9999999.;
0533 for (i = 0; i < size; i++) {
0534 ttemp = sqrt(_mgacollection[i].dx * _mgacollection[i].dx + _mgacollection[i].dy * _mgacollection[i].dy);
0535 rtemp = sqrt(_mgacollection[i].x * _mgacollection[i].x + _mgacollection[i].y * _mgacollection[i].y);
0536 if (ttemp > maxV)
0537 maxV = ttemp;
0538 if (rtemp > maxR)
0539 maxR = rtemp;
0540 }
0541
0542
0543 float smallestVcm = .001;
0544 if (maxV < smallestVcm)
0545 maxV = smallestVcm;
0546 float scale = 0.;
0547 float lside = 1.1 * maxR;
0548 if (lside <= 0)
0549 lside = 100.;
0550 if (maxV > 0) {
0551 scale = .09 * lside / maxV;
0552 }
0553 char scalename[50];
0554 int ret = snprintf(scalename, 50, "#delta #bar{x} length =%f cm", maxV);
0555
0556
0557 if (ret > 0) {
0558 dxh = new TH2F("vecdrplot", scalename, 80, -lside, lside, 80, -lside, lside);
0559 } else {
0560 dxh = new TH2F("vecdrplot", "delta #bar{x} Bad scale", 80, -lside, lside, 80, -lside, lside);
0561 }
0562 dxh->GetXaxis()->SetTitle("x in cm");
0563 dxh->GetYaxis()->SetTitle("y in cm");
0564 dxh->SetStats(kFALSE);
0565 dxh->Draw();
0566 TArrow* arrow;
0567 for (i = 0; i < size; i++) {
0568 ttemp = sqrt(_mgacollection[i].dx * _mgacollection[i].dx + _mgacollection[i].dy * _mgacollection[i].dy);
0569
0570 float nx = _mgacollection[i].x + scale * _mgacollection[i].dx;
0571 float ny = _mgacollection[i].y + scale * _mgacollection[i].dy;
0572 arrow = new TArrow(_mgacollection[i].x, _mgacollection[i].y, nx, ny);
0573 arrow->SetLineWidth(2);
0574 arrow->SetArrowSize(ttemp * .2 * .05 / maxV);
0575 arrow->SetLineColor(1);
0576 arrow->SetLineStyle(1);
0577 arrow->Paint();
0578 dxh->GetListOfFunctions()->Add(static_cast<TObject*>(arrow));
0579
0580
0581 }
0582 dxh->Write();
0583
0584 _theFile->Write();
0585 _theFile->Close();
0586 }
0587
0588 void MuonGeometryArrange::makeGraph(int sizeI,
0589 float smi,
0590 float sma,
0591 float minV,
0592 float maxV,
0593 TH2F* dxh,
0594 TGraph* grx,
0595 const char* name,
0596 const char* title,
0597 const char* titleg,
0598 const char* axis,
0599 const float* xp,
0600 const float* yp,
0601 int size) {
0602 if (minV >= maxV || smi >= sma || sizeI <= 1 || xp == nullptr || yp == nullptr)
0603 return;
0604
0605 float diff = maxV - minV;
0606 float over = .05 * diff;
0607 double ylo = minV - over;
0608 double yhi = maxV + over;
0609 double dsmi, dsma;
0610 dsmi = smi;
0611 dsma = sma;
0612 dxh = new TH2F(name, title, sizeI + 2, dsmi, dsma, 50, ylo, yhi);
0613 dxh->GetXaxis()->SetTitle("Position around ring");
0614 dxh->GetYaxis()->SetTitle(axis);
0615 dxh->SetStats(kFALSE);
0616 dxh->Draw();
0617 grx = new TGraph(size, xp, yp);
0618 grx->SetName(titleg);
0619 grx->SetTitle(title);
0620 grx->SetMarkerColor(2);
0621 grx->SetMarkerStyle(3);
0622 grx->GetXaxis()->SetLimits(dsmi, dsma);
0623 grx->GetXaxis()->SetTitle("position number");
0624 grx->GetYaxis()->SetLimits(ylo, yhi);
0625 grx->GetYaxis()->SetTitle(axis);
0626 grx->Draw("A*");
0627 grx->Write();
0628 return;
0629 }
0630
0631 void MuonGeometryArrange::beginJob() { firstEvent_ = true; }
0632
0633
0634 void MuonGeometryArrange::createROOTGeometry(const edm::EventSetup& iSetup) {}
0635
0636 void MuonGeometryArrange::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0637 if (firstEvent_) {
0638 MuonAlignmentInputXML inputMethod1(_inputXMLCurrent,
0639 &iSetup.getData(dtGeomToken1_),
0640 &iSetup.getData(cscGeomToken1_),
0641 &iSetup.getData(gemGeomToken1_),
0642 &iSetup.getData(dtGeomToken1_),
0643 &iSetup.getData(cscGeomToken1_),
0644 &iSetup.getData(gemGeomToken1_));
0645 inputAlign1 = new MuonAlignment(iSetup, inputMethod1);
0646 inputAlign1->fillGapsInSurvey(0, 0);
0647 MuonAlignmentInputXML inputMethod2(_inputXMLReference,
0648 &iSetup.getData(dtGeomToken2_),
0649 &iSetup.getData(cscGeomToken2_),
0650 &iSetup.getData(gemGeomToken2_),
0651 &iSetup.getData(dtGeomToken1_),
0652 &iSetup.getData(cscGeomToken1_),
0653 &iSetup.getData(gemGeomToken1_));
0654 inputAlign2 = new MuonAlignment(iSetup, inputMethod2);
0655 inputAlign2->fillGapsInSurvey(0, 0);
0656 MuonAlignmentInputXML inputMethod2a(_inputXMLReference,
0657 &iSetup.getData(dtGeomToken3_),
0658 &iSetup.getData(cscGeomToken3_),
0659 &iSetup.getData(gemGeomToken3_),
0660 &iSetup.getData(dtGeomToken1_),
0661 &iSetup.getData(cscGeomToken1_),
0662 &iSetup.getData(gemGeomToken1_));
0663 inputAlign2a = new MuonAlignment(iSetup, inputMethod2a);
0664 inputAlign2a->fillGapsInSurvey(0, 0);
0665
0666 inputGeometry1 = static_cast<Alignable*>(inputAlign1->getAlignableMuon());
0667 inputGeometry2 = static_cast<Alignable*>(inputAlign2->getAlignableMuon());
0668 auto inputGeometry2Copy2 = inputAlign2a->getAlignableMuon();
0669
0670
0671 edm::LogInfo("MuonGeometryArrange") << "levels: " << _levelStrings.size();
0672 for (const auto& level : _levelStrings) {
0673 theLevels.push_back(inputGeometry2Copy2->objectIdProvider().stringToId(level));
0674 edm::LogInfo("MuonGeometryArrange") << "level: " << level;
0675 }
0676
0677
0678 compare(inputGeometry1, inputGeometry2, inputGeometry2Copy2);
0679
0680
0681
0682 _theFile->cd();
0683 _alignTree->Write();
0684 endHist();
0685
0686
0687 firstEvent_ = false;
0688 }
0689 }
0690
0691
0692 void MuonGeometryArrange::compare(Alignable* refAli, Alignable* curAli, Alignable* curAliCopy2) {
0693
0694 if (refAli == nullptr) {
0695 return;
0696 }
0697 if (curAli == nullptr) {
0698 return;
0699 }
0700
0701 const auto& refComp = refAli->components();
0702 const auto& curComp = curAli->components();
0703 const auto& curComp2 = curAliCopy2->components();
0704 compareGeometries(refAli, curAli, curAliCopy2);
0705
0706 int nComp = refComp.size();
0707 for (int i = 0; i < nComp; i++) {
0708 compare(refComp[i], curComp[i], curComp2[i]);
0709 }
0710 return;
0711 }
0712
0713
0714 void MuonGeometryArrange::compareGeometries(Alignable* refAli, Alignable* curAli, Alignable* curCopy) {
0715
0716 if (refAli == nullptr) {
0717 return;
0718 }
0719 if (curAli == nullptr) {
0720 return;
0721 }
0722
0723
0724 if (!isMother(refAli))
0725 return;
0726
0727
0728
0729 if (isMother(refAli->mother()))
0730 return;
0731 const auto& refComp = refAli->components();
0732 const auto& curComp = curCopy->components();
0733 if (refComp.size() != curComp.size()) {
0734 return;
0735 }
0736
0737 align::GlobalVectors originalVectors;
0738 align::GlobalVectors currentVectors;
0739 align::GlobalVectors originalRelativeVectors;
0740 align::GlobalVectors currentRelativeVectors;
0741
0742 int nComp = refComp.size();
0743 int nUsed = 0;
0744
0745 CLHEP::Hep3Vector TotalX, TotalL;
0746 TotalX.set(0., 0., 0.);
0747 TotalL.set(0., 0., 0.);
0748
0749 std::vector<CLHEP::Hep3Vector> Positions;
0750 std::vector<CLHEP::Hep3Vector> DelPositions;
0751
0752 double xrcenter = 0.;
0753 double yrcenter = 0.;
0754 double zrcenter = 0.;
0755 double xccenter = 0.;
0756 double yccenter = 0.;
0757 double zccenter = 0.;
0758
0759 bool useIt;
0760
0761
0762 for (int ich = 0; ich < nComp; ich++) {
0763 useIt = true;
0764 if (_weightById) {
0765 if (!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
0766 useIt = false;
0767 }
0768 if (!useIt)
0769 continue;
0770 align::GlobalVectors curVs;
0771 align::createPoints(&curVs, refComp[ich], _weightBy, _weightById, _weightByIdVector);
0772 align::GlobalVector pointsCM = align::centerOfMass(curVs);
0773 originalVectors.push_back(pointsCM);
0774 nUsed++;
0775 xrcenter += pointsCM.x();
0776 yrcenter += pointsCM.y();
0777 zrcenter += pointsCM.z();
0778 }
0779 xrcenter = xrcenter / nUsed;
0780 yrcenter = yrcenter / nUsed;
0781 zrcenter = zrcenter / nUsed;
0782
0783
0784
0785 for (int ich = 0; ich < nComp; ich++) {
0786 useIt = true;
0787 if (_weightById) {
0788 if (!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
0789 useIt = false;
0790 }
0791 if (!useIt)
0792 continue;
0793 align::GlobalVectors curVs;
0794 align::createPoints(&curVs, curComp[ich], _weightBy, _weightById, _weightByIdVector);
0795 align::GlobalVector pointsCM = align::centerOfMass(curVs);
0796 currentVectors.push_back(pointsCM);
0797
0798 xccenter += pointsCM.x();
0799 yccenter += pointsCM.y();
0800 zccenter += pointsCM.z();
0801 }
0802 xccenter = xccenter / nUsed;
0803 yccenter = yccenter / nUsed;
0804 zccenter = zccenter / nUsed;
0805
0806
0807 align::GlobalVector CCur(xccenter, yccenter, zccenter);
0808 align::GlobalVector CRef(xrcenter, yrcenter, zrcenter);
0809 int nCompR = currentVectors.size();
0810 for (int ich = 0; ich < nCompR; ich++) {
0811 originalRelativeVectors.push_back(originalVectors[ich] - CRef);
0812 currentRelativeVectors.push_back(currentVectors[ich] - CCur);
0813 }
0814
0815
0816
0817
0818
0819 align::RotationType rtype3 = align::diffRot(currentRelativeVectors, originalRelativeVectors);
0820
0821 align::EulerAngles angles(3);
0822 angles = align::toAngles(rtype3);
0823
0824 for (int ich = 0; ich < nComp; ich++) {
0825 if (_weightById) {
0826 if (!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
0827 continue;
0828 }
0829 CLHEP::Hep3Vector Rtotal, Wtotal;
0830 Rtotal.set(0., 0., 0.);
0831 Wtotal.set(0., 0., 0.);
0832 for (int i = 0; i < 100; i++) {
0833 AlgebraicVector diff =
0834 align::diffAlignables(refComp[ich], curComp[ich], _weightBy, _weightById, _weightByIdVector);
0835 CLHEP::Hep3Vector dR(diff[0], diff[1], diff[2]);
0836 Rtotal += dR;
0837 CLHEP::Hep3Vector dW(diff[3], diff[4], diff[5]);
0838 CLHEP::HepRotation rot(Wtotal.unit(), Wtotal.mag());
0839 CLHEP::HepRotation drot(dW.unit(), dW.mag());
0840 rot *= drot;
0841 Wtotal.set(rot.axis().x() * rot.delta(), rot.axis().y() * rot.delta(), rot.axis().z() * rot.delta());
0842 align::moveAlignable(curComp[ich], diff);
0843 float tolerance = 1e-7;
0844 AlgebraicVector check =
0845 align::diffAlignables(refComp[ich], curComp[ich], _weightBy, _weightById, _weightByIdVector);
0846 align::GlobalVector checkR(check[0], check[1], check[2]);
0847 align::GlobalVector checkW(check[3], check[4], check[5]);
0848 DetId detid(refComp[ich]->id());
0849 if ((checkR.mag() > tolerance) || (checkW.mag() > tolerance)) {
0850
0851
0852
0853
0854 } else {
0855 TotalX += Rtotal;
0856 break;
0857 }
0858 }
0859 }
0860
0861
0862 TotalX = TotalX / nUsed;
0863
0864
0865 AlgebraicVector change(6);
0866 change(1) = TotalX.x();
0867 change(2) = TotalX.y();
0868 change(3) = TotalX.z();
0869
0870 change(4) = angles[0];
0871 change(5) = angles[1];
0872 change(6) = angles[2];
0873 align::moveAlignable(curAli, change);
0874
0875
0876 const auto& curComp2 = curAli->components();
0877
0878 for (int ich = 0; ich < nComp; ich++) {
0879 CLHEP::Hep3Vector Rtotal, Wtotal;
0880 Rtotal.set(0., 0., 0.);
0881 Wtotal.set(0., 0., 0.);
0882 if (_weightById) {
0883 if (!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
0884 continue;
0885 }
0886
0887 for (int i = 0; i < 100; i++) {
0888 AlgebraicVector diff =
0889 align::diffAlignables(refComp[ich], curComp2[ich], _weightBy, _weightById, _weightByIdVector);
0890 CLHEP::Hep3Vector dR(diff[0], diff[1], diff[2]);
0891 Rtotal += dR;
0892 CLHEP::Hep3Vector dW(diff[3], diff[4], diff[5]);
0893 CLHEP::HepRotation rot(Wtotal.unit(), Wtotal.mag());
0894 CLHEP::HepRotation drot(dW.unit(), dW.mag());
0895 rot *= drot;
0896 Wtotal.set(rot.axis().x() * rot.delta(), rot.axis().y() * rot.delta(), rot.axis().z() * rot.delta());
0897 align::moveAlignable(curComp2[ich], diff);
0898 float tolerance = 1e-7;
0899 AlgebraicVector check =
0900 align::diffAlignables(refComp[ich], curComp2[ich], _weightBy, _weightById, _weightByIdVector);
0901 align::GlobalVector checkR(check[0], check[1], check[2]);
0902 align::GlobalVector checkW(check[3], check[4], check[5]);
0903 if ((checkR.mag() > tolerance) || (checkW.mag() > tolerance)) {
0904 } else {
0905 break;
0906 }
0907 }
0908 AlgebraicVector TRtot(6);
0909 TRtot(1) = Rtotal.x();
0910 TRtot(2) = Rtotal.y();
0911 TRtot(3) = Rtotal.z();
0912 TRtot(4) = Wtotal.x();
0913 TRtot(5) = Wtotal.y();
0914 TRtot(6) = Wtotal.z();
0915 fillTree(refComp[ich], TRtot);
0916 }
0917 }
0918
0919
0920
0921 void MuonGeometryArrange::fillTree(Alignable* refAli, const AlgebraicVector& diff) {
0922 _id = refAli->id();
0923 _level = refAli->alignableObjectId();
0924
0925 if (refAli->mother()) {
0926 _mid = refAli->mother()->geomDetId().rawId();
0927 _mlevel = refAli->mother()->alignableObjectId();
0928 } else {
0929 _mid = -1;
0930 _mlevel = -1;
0931 }
0932 DetId detid(_id);
0933 _sublevel = detid.subdetId();
0934 int ringPhiPos = -99;
0935 if (detid.det() == DetId::Muon && detid.subdetId() == MuonSubdetId::CSC) {
0936 CSCDetId cscId(refAli->geomDetId());
0937 ringPhiPos = cscId.chamber();
0938 }
0939 _xVal = refAli->globalPosition().x();
0940 _yVal = refAli->globalPosition().y();
0941 _zVal = refAli->globalPosition().z();
0942 align::GlobalVector vec(_xVal, _yVal, _zVal);
0943 _rVal = vec.perp();
0944 _phiVal = vec.phi();
0945 _etaVal = vec.eta();
0946 align::RotationType rot = refAli->globalRotation();
0947 align::EulerAngles eulerAngles = align::toAngles(rot);
0948 _rotxVal = atan2(rot.yz(), rot.zz());
0949 float ttt = -rot.xz();
0950 if (ttt > 1.)
0951 ttt = 1.;
0952 if (ttt < -1.)
0953 ttt = -1.;
0954 _rotyVal = asin(ttt);
0955 _rotzVal = atan2(rot.xy(), rot.xx());
0956 _alphaVal = eulerAngles[0];
0957 _betaVal = eulerAngles[1];
0958 _gammaVal = eulerAngles[2];
0959 _dxVal = diff[0];
0960 _dyVal = diff[1];
0961 _dzVal = diff[2];
0962
0963 align::GlobalVector vRef(_xVal, _yVal, _zVal);
0964 align::GlobalVector vCur(_xVal - _dxVal, _yVal - _dyVal, _zVal - _dzVal);
0965 _drVal = vCur.perp() - vRef.perp();
0966 _dphiVal = vCur.phi() - vRef.phi();
0967
0968 _dalphaVal = diff[3];
0969 _dbetaVal = diff[4];
0970 _dgammaVal = diff[5];
0971 _drotxVal = -999.;
0972 _drotyVal = -999.;
0973 _drotzVal = -999.;
0974
0975 align::EulerAngles deuler(3);
0976 deuler(1) = _dalphaVal;
0977 deuler(2) = _dbetaVal;
0978 deuler(3) = _dgammaVal;
0979 align::RotationType drot = align::toMatrix(deuler);
0980 double xx = rot.xx();
0981 double xy = rot.xy();
0982 double xz = rot.xz();
0983 double yx = rot.yx();
0984 double yy = rot.yy();
0985 double yz = rot.yz();
0986 double zx = rot.zx();
0987 double zy = rot.zy();
0988 double zz = rot.zz();
0989 double detrot = (zz * yy - zy * yz) * xx + (-zz * yx + zx * yz) * xy + (zy * yx - zx * yy) * xz;
0990 detrot = 1 / detrot;
0991 double ixx = (zz * yy - zy * yz) * detrot;
0992 double ixy = (-zz * xy + zy * xz) * detrot;
0993 double ixz = (yz * xy - yy * xz) * detrot;
0994 double iyx = (-zz * yx + zx * yz) * detrot;
0995 double iyy = (zz * xx - zx * xz) * detrot;
0996 double iyz = (-yz * xx + yx * xz) * detrot;
0997 double izx = (zy * yx - zx * yy) * detrot;
0998 double izy = (-zy * xx + zx * xy) * detrot;
0999 double izz = (yy * xx - yx * xy) * detrot;
1000 align::RotationType invrot(ixx, ixy, ixz, iyx, iyy, iyz, izx, izy, izz);
1001 align::RotationType prot = rot * drot * invrot;
1002
1003 float protx;
1004 protx = atan2(prot.yz(), prot.zz());
1005 _drotxVal = protx;
1006 ttt = -prot.xz();
1007 if (ttt > 1.)
1008 ttt = 1.;
1009 if (ttt < -1.)
1010 ttt = -1.;
1011 _drotyVal = asin(ttt);
1012 _drotzVal = atan2(prot.xy(), prot.xx());
1013
1014
1015 if (_drotxVal > 3.141592656)
1016 _drotxVal = -6.2831853072 + _drotxVal;
1017 if (_drotxVal < -3.141592656)
1018 _drotxVal = 6.2831853072 + _drotxVal;
1019 if (_drotyVal > 3.141592656)
1020 _drotyVal = -6.2831853072 + _drotyVal;
1021 if (_drotyVal < -3.141592656)
1022 _drotyVal = 6.2831853072 + _drotyVal;
1023 if (_drotzVal > 3.141592656)
1024 _drotzVal = -6.2831853072 + _drotzVal;
1025 if (_drotzVal < -3.141592656)
1026 _drotzVal = 6.2831853072 + _drotzVal;
1027
1028 _ldxVal = -999.;
1029 _ldyVal = -999.;
1030 _ldxVal = -999.;
1031 _ldrVal = -999.;
1032 _ldphiVal = -999;
1033
1034
1035 align::GlobalVector dV(_dxVal, _dyVal, _dzVal);
1036 align::LocalVector pointL = refAli->surface().toLocal(dV);
1037
1038 _ldxVal = pointL.x();
1039 _ldyVal = pointL.y();
1040 _ldzVal = pointL.z();
1041 _ldphiVal = pointL.phi();
1042 _ldrVal = pointL.perp();
1043
1044
1045 if (refAli->alignableObjectId() == align::AlignableDetUnit) {
1046 if (_detIdFlag) {
1047 if ((passIdCut(refAli->id())) || (passIdCut(refAli->mother()->id()))) {
1048 _useDetId = 1;
1049 } else {
1050 _useDetId = 0;
1051 }
1052 }
1053 }
1054
1055 if (refAli->alignableObjectId() == align::AlignableDetUnit) {
1056 if (refAli->mother()->alignableObjectId() != align::AlignableDet) {
1057 _detDim = 1;
1058 } else if (refAli->mother()->alignableObjectId() == align::AlignableDet) {
1059 _detDim = 2;
1060 }
1061 } else
1062 _detDim = 0;
1063
1064 _surWidth = refAli->surface().width();
1065 _surLength = refAli->surface().length();
1066 align::RotationType rt = refAli->globalRotation();
1067 _surRot[0] = rt.xx();
1068 _surRot[1] = rt.xy();
1069 _surRot[2] = rt.xz();
1070 _surRot[3] = rt.yx();
1071 _surRot[4] = rt.yy();
1072 _surRot[5] = rt.yz();
1073 _surRot[6] = rt.zx();
1074 _surRot[7] = rt.zy();
1075 _surRot[8] = rt.zz();
1076
1077 MGACollection holdit;
1078 holdit.id = _id;
1079 holdit.level = _level;
1080 holdit.mid = _mid;
1081 holdit.mlevel = _mlevel;
1082 holdit.sublevel = _sublevel;
1083 holdit.x = _xVal;
1084 holdit.y = _yVal;
1085 holdit.z = _zVal;
1086 holdit.r = _rVal;
1087 holdit.phi = _phiVal;
1088 holdit.eta = _etaVal;
1089 holdit.alpha = _alphaVal;
1090 holdit.beta = _betaVal;
1091 holdit.gamma = _gammaVal;
1092 holdit.dx = _dxVal;
1093 holdit.dy = _dyVal;
1094 holdit.dz = _dzVal;
1095 holdit.dr = _drVal;
1096 holdit.dphi = _dphiVal;
1097 holdit.dalpha = _dalphaVal;
1098 holdit.dbeta = _dbetaVal;
1099 holdit.dgamma = _dgammaVal;
1100 holdit.useDetId = _useDetId;
1101 holdit.detDim = _detDim;
1102 holdit.surW = _surWidth;
1103 holdit.surL = _surLength;
1104 holdit.ldx = _ldxVal;
1105 holdit.ldy = _ldyVal;
1106 holdit.ldz = _ldzVal;
1107 holdit.ldr = _ldrVal;
1108 holdit.ldphi = _ldphiVal;
1109 holdit.rotx = _rotxVal;
1110 holdit.roty = _rotyVal;
1111 holdit.rotz = _rotzVal;
1112 holdit.drotx = _drotxVal;
1113 holdit.droty = _drotyVal;
1114 holdit.drotz = _drotzVal;
1115 for (int i = 0; i < 9; i++) {
1116 holdit.surRot[i] = _surRot[i];
1117 }
1118 holdit.phipos = ringPhiPos;
1119 _mgacollection.push_back(holdit);
1120
1121
1122 _alignTree->Fill();
1123 }
1124
1125
1126 bool MuonGeometryArrange::isMother(Alignable* ali) {
1127
1128 if (ali == nullptr)
1129 return false;
1130 const auto& aliComp = ali->components();
1131
1132 int size = aliComp.size();
1133 if (size <= 0)
1134 return false;
1135
1136 for (int i = 0; i < size; i++) {
1137 if (checkChosen(aliComp[i]))
1138 return true;
1139 }
1140 return false;
1141 }
1142
1143
1144 bool MuonGeometryArrange::checkChosen(Alignable* ali) {
1145
1146 if (ali == nullptr)
1147 return false;
1148
1149 if (ali->geomDetId().det() != DetId::Muon || ali->geomDetId().subdetId() != MuonSubdetId::CSC)
1150 return false;
1151
1152
1153
1154
1155
1156 CSCDetId cscId(ali->geomDetId());
1157 #ifdef jnbdebug
1158 std::cout << "JNB " << ali->id() << " " << cscId.endcap() << " " << cscId.station() << " " << cscId.ring() << " "
1159 << cscId.chamber() << " " << _endcap << " " << _station << " " << _ring << "\n"
1160 << std::flush;
1161 #endif
1162 if (cscId.endcap() == _endcap && cscId.station() == _station && cscId.ring() == _ring) {
1163 return true;
1164 }
1165 return false;
1166 }
1167
1168
1169 bool MuonGeometryArrange::passChosen(Alignable* ali) {
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181 if (ali == nullptr)
1182 return false;
1183 if (checkChosen(ali))
1184 return true;
1185
1186 const auto& aliComp = ali->components();
1187
1188 int size = aliComp.size();
1189 if (size <= 0)
1190 return false;
1191
1192 for (int i = 0; i < size; i++) {
1193 if (checkChosen(aliComp[i]))
1194 return true;
1195 }
1196 return false;
1197 }
1198
1199 bool MuonGeometryArrange::passIdCut(uint32_t id) {
1200 bool pass = false;
1201 DetId detid(id);
1202
1203
1204
1205
1206 int nEntries = _detIdFlagVector.size();
1207
1208 for (int i = 0; i < nEntries; i++) {
1209 if (_detIdFlagVector[i] == id)
1210 pass = true;
1211 }
1212
1213 return pass;
1214 }
1215
1216
1217 DEFINE_FWK_MODULE(MuonGeometryArrange);