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