Back to home page

Project CMSSW displayed by LXR
 
 

    

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