Back to home page

Project CMSSW displayed by LXR
 
 

    

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 // The following looks generic enough to use
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 // Database
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   // Input is XML
0067   _inputXMLCurrent = cfg.getUntrackedParameter<std::string>("inputXMLCurrent");
0068   _inputXMLReference = cfg.getUntrackedParameter<std::string>("inputXMLReference");
0069 
0070   //input is ROOT
0071   _inputFilename1 = cfg.getUntrackedParameter<std::string>("inputROOTFile1");
0072   _inputFilename2 = cfg.getUntrackedParameter<std::string>("inputROOTFile2");
0073   _inputTreename = cfg.getUntrackedParameter<std::string>("treeName");
0074 
0075   //output file
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   // if want to use, make id cut list
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   // turn weightByIdFile into weightByIdVector
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   //root configuration
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   // Unpack the list and create ntuples here.
0166 
0167   int size = _mgacollection.size();
0168   if (size <= 0)
0169     return;  // nothing to do here.
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   // for position plots:
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;                     // can't do anything?
0193   xp[size] = xp[size - 1] + 1;  // wraparound point
0194 
0195   if (1 < minI)
0196     minI = 1;
0197   if (size > maxI)
0198     maxI = size;
0199   maxI++;  // allow for wraparound to show neighbors
0200   int sizeI = maxI + 1 - minI;
0201   float smi = minI - 1;
0202   float sma = maxI + 1;
0203 
0204   // Dx plot
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];  // wraparound point
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   // Dy plot
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];  // wraparound point
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   // Dz plot
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];  // wraparound point
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   // Dphi plot
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];  // wraparound point
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   // Dr plot
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];  // wraparound point
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   // Drphi plot
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];  // wraparound point
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   // Dalpha plot
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];  // wraparound point
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   // Dbeta plot
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];  // wraparound point
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   // Dgamma plot
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];  // wraparound point
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   // Drotx plot
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];  // wraparound point
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   // Droty plot
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];  // wraparound point
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   // Drotz plot
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];  // wraparound point
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   // Vector plots
0528   // First find the maximum length of sqrt(dx*dx+dy*dy):  we'll have to
0529   // scale these for visibility
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   // Don't try to scale rediculously small values
0543   float smallestVcm = .001;  // 10 microns
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   }  // units of pad length!
0553   char scalename[50];
0554   int ret = snprintf(scalename, 50, "#delta #bar{x}   length =%f cm", maxV);
0555   // If ret<=0 we don't want to print the scale!
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     //     ttemp=ttemp*scale;
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);  // ttemp*.3*.05, "->");
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     //     arrow->Draw();
0580     //     arrow->Write();
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   // out of bounds, bail
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     //setting the levels being used in the geometry comparator
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     //compare the goemetries
0678     compare(inputGeometry1, inputGeometry2, inputGeometry2Copy2);
0679 
0680     //write out ntuple
0681     //might be better to do within output module
0682     _theFile->cd();
0683     _alignTree->Write();
0684     endHist();
0685     //   _theFile->Close();
0686 
0687     firstEvent_ = false;
0688   }
0689 }
0690 
0691 /////////////////////////////////////////////////
0692 void MuonGeometryArrange::compare(Alignable* refAli, Alignable* curAli, Alignable* curAliCopy2) {
0693   // First sanity
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   // First sanity
0716   if (refAli == nullptr) {
0717     return;
0718   }
0719   if (curAli == nullptr) {
0720     return;
0721   }
0722   // Is this the Ring we want to align?  If so it will contain the
0723   // chambers specified in the configuration file
0724   if (!isMother(refAli))
0725     return;  // Not the desired alignable object
0726              // But... There are granddaughters involved--and I don't want to monkey with
0727              // the layers of the chambers.  So, if the mother of this is also an approved
0728              // mother, bail.
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   // GlobalVectors is a vector of GlobalVector which is a 3D vector
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   // Use the total displacements here:
0745   CLHEP::Hep3Vector TotalX, TotalL;
0746   TotalX.set(0., 0., 0.);
0747   TotalL.set(0., 0., 0.);
0748   //  CLHEP::Hep3Vector* Rsubtotal, Wsubtotal, DRsubtotal, DWsubtotal;
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   // Create the "center" for the reference alignment chambers, and
0761   // load a vector of their centers
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   // Create the "center" for the current alignment chambers, and
0784   // load a vector of their centers
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   // OK, now load the <very approximate> vectors from the ring "centers"
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   // All right.  Now let the hacking begin.
0816   // First out of the gate let's try using the raw values and see what
0817   // diffRot does for us.
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         //   edm::LogInfo("CompareGeoms") << "Tolerance Exceeded!(alObjId: "
0851         //       << refAli->alignableObjectId()
0852         //   << ", rawId: " << refComp[ich]->geomDetId().rawId()
0853         //   << ", subdetId: "<< detid.subdetId() << "): " << diff;
0854       } else {
0855         TotalX += Rtotal;
0856         break;
0857       }  // end of else
0858     }    // end of for on int i
0859   }      // end of for on ich
0860 
0861   // At this point we should have a total displacement and total L
0862   TotalX = TotalX / nUsed;
0863 
0864   // Now start again!
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);  // move as a chunk
0874 
0875   // Now get the components again.  They should be in new locations
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     }  // end of for on int i
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   }  // end of for on ich
0917 }
0918 
0919 //////////////////////////////////////////////////
0920 
0921 void MuonGeometryArrange::fillTree(Alignable* refAli, const AlgebraicVector& diff) {
0922   _id = refAli->id();
0923   _level = refAli->alignableObjectId();
0924   //need if ali has no mother
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   //getting dR and dPhi
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   //    align::RotationType prot = rot*drot;
1003   float protx;  //, proty, protz;
1004   protx = atan2(prot.yz(), prot.zz());
1005   _drotxVal = protx;  //_rotxVal-protx; //atan2(drot.yz(), drot.zz());
1006   ttt = -prot.xz();
1007   if (ttt > 1.)
1008     ttt = 1.;
1009   if (ttt < -1.)
1010     ttt = -1.;
1011   _drotyVal = asin(ttt);                    // -_rotyVal;
1012   _drotzVal = atan2(prot.xy(), prot.xx());  // - _rotzVal;
1013                                             // Above does not account for 2Pi wraparounds!
1014                                             // Prior knowledge:  these are supposed to be small rotations.  Therefore:
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;  // set fake
1033 
1034   //    if(refAli->alignableObjectId() == align::AlignableDetUnit){
1035   align::GlobalVector dV(_dxVal, _dyVal, _dzVal);
1036   align::LocalVector pointL = refAli->surface().toLocal(dV);
1037   //align::LocalVector pointL = (refAli->mother())->surface().toLocal(dV);
1038   _ldxVal = pointL.x();
1039   _ldyVal = pointL.y();
1040   _ldzVal = pointL.z();
1041   _ldphiVal = pointL.phi();
1042   _ldrVal = pointL.perp();
1043   //    }
1044   //detIdFlag
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   // det module dimension
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   //Fill
1122   _alignTree->Fill();
1123 }
1124 
1125 //////////////////////////////////////////////////
1126 bool MuonGeometryArrange::isMother(Alignable* ali) {
1127   // Is this the mother ring?
1128   if (ali == nullptr)
1129     return false;  // elementary sanity
1130   const auto& aliComp = ali->components();
1131 
1132   int size = aliComp.size();
1133   if (size <= 0)
1134     return false;  // no subcomponents
1135 
1136   for (int i = 0; i < size; i++) {
1137     if (checkChosen(aliComp[i]))
1138       return true;  // A ring has CSC chambers
1139   }                 // as subcomponents
1140   return false;     // 1'st layer of subcomponents weren't CSC chambers
1141 }
1142 //////////////////////////////////////////////////
1143 
1144 bool MuonGeometryArrange::checkChosen(Alignable* ali) {
1145   // Check whether the item passed satisfies the criteria given.
1146   if (ali == nullptr)
1147     return false;  // elementary sanity
1148                    // Is this in the CSC section?  If not, bail.  Later may extend.
1149   if (ali->geomDetId().det() != DetId::Muon || ali->geomDetId().subdetId() != MuonSubdetId::CSC)
1150     return false;
1151   // If it is a CSC alignable, then check that the station, etc are
1152   // those requested.
1153   // One might think of aligning more than a single ring at a time,
1154   // by using a vector of ring numbers.  I don't see the sense in
1155   // trying to align more than one station at a time for comparison.
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   // Check to see if this contains CSC components of the appropriate ring
1171   // Ring will contain N Alignables which represent chambers, each of which
1172   // in turn contains M planes.  For our purposes we don't care about the
1173   // planes.
1174   // Hmm.  Interesting question:  Do I want to try to fit the chamber as
1175   // such, or use the geometry?
1176   // I want to fit the chamber, so I'll try to use its presence as the marker.
1177   // What specifically identifies a chamber as a chamber, and not as a layer?
1178   // The fact that it has layers as sub components, or the fact that it is
1179   // the first item with a non-zero ID breakdown?  Pick the latter.
1180   //
1181   if (ali == nullptr)
1182     return false;
1183   if (checkChosen(ali))
1184     return true;  // If this is one of the desired
1185                   // CSC chambers, accept it
1186   const auto& aliComp = ali->components();
1187 
1188   int size = aliComp.size();
1189   if (size <= 0)
1190     return false;  // no subcomponents
1191 
1192   for (int i = 0; i < size; i++) {
1193     if (checkChosen(aliComp[i]))
1194       return true;  // A ring has CSC chambers
1195   }                 // as subcomponents
1196   return false;     // 1'st layer of subcomponents weren't CSC chambers
1197 }
1198 //////////////////////////////////////////////////
1199 bool MuonGeometryArrange::passIdCut(uint32_t id) {
1200   bool pass = false;
1201   DetId detid(id);
1202   //    if(detid.det()==DetId::Muon && detid.subdetId()== MuonSubdetId::CSC){
1203   //       CSCDetId cscId(refAli->geomDetId());
1204   //       if(cscId.layer()!=1) return false;       // ONLY FIRST LAYER!
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!