Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-08-13 05:00:12 

0001 #include "CondFormats/Alignment/interface/Alignments.h"
0002 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0003 #include "CondFormats/Alignment/interface/AlignmentSurfaceDeformations.h"
0004 #include "CondFormats/Alignment/interface/Definitions.h"
0005 #include "CLHEP/Vector/RotationInterfaces.h"
0006 #include "CondFormats/AlignmentRecord/interface/TrackerSurveyRcd.h"
0007 #include "CondFormats/AlignmentRecord/interface/TrackerSurveyErrorExtendedRcd.h"
0008 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0009 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ESTransientHandle.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 
0018 #include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
0019 #include "CondFormats/DataRecord/interface/SiPixelQualityRcd.h"
0020 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0021 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0022 #include "CondCore/SiPixelPlugins/interface/SiPixelPayloadInspectorHelper.h"
0023 
0024 #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
0025 #include "Geometry/CommonTopologies/interface/SurfaceDeformationFactory.h"
0026 #include "Geometry/CommonTopologies/interface/SurfaceDeformation.h"
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
0029 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0030 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0031 #include "Geometry/GeometryAligner/interface/GeometryAligner.h"
0032 #include "Alignment/CommonAlignment/interface/Utilities.h"
0033 #include "Alignment/CommonAlignment/interface/SurveyDet.h"
0034 #include "Alignment/CommonAlignment/interface/Alignable.h"
0035 #include "DataFormats/DetId/interface/DetId.h"
0036 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
0037 #include "CondFormats/Alignment/interface/DetectorGlobalPosition.h"
0038 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0039 
0040 #include "TrackerGeometryCompare.h"
0041 #include "TFile.h"
0042 #include "CLHEP/Vector/ThreeVector.h"
0043 
0044 // Database
0045 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0046 #include "FWCore/ServiceRegistry/interface/Service.h"
0047 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0048 
0049 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0050 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0051 
0052 #include <iostream>
0053 #include <fstream>
0054 #include <sstream>
0055 
0056 TrackerGeometryCompare::TrackerGeometryCompare(const edm::ParameterSet& cfg)
0057     : cpvTokenDDD_(esConsumes()),
0058       cpvTokenDD4hep_(esConsumes()),
0059       topoToken_(esConsumes()),
0060       geomDetToken_(esConsumes()),
0061       ptpToken_(esConsumes()),
0062       ptitpToken_(esConsumes()),
0063       pixQualityToken_(esConsumes()),
0064       stripQualityToken_(esConsumes()),
0065       referenceTracker(nullptr),
0066       dummyTracker(nullptr),
0067       currentTracker(nullptr),
0068       theSurveyIndex(0),
0069       theSurveyValues(nullptr),
0070       theSurveyErrors(nullptr),
0071       levelStrings_(cfg.getUntrackedParameter<std::vector<std::string> >("levels")),
0072       fromDD4hep_(cfg.getUntrackedParameter<bool>("fromDD4hep")),
0073       writeToDB_(cfg.getUntrackedParameter<bool>("writeToDB")),
0074       commonTrackerLevel_(align::invalid),
0075       moduleListFile_(),
0076       moduleList_(0),
0077       inputRootFile1_(nullptr),
0078       inputRootFile2_(nullptr),
0079       inputTree01_(nullptr),
0080       inputTree02_(nullptr),
0081       inputTree11_(nullptr),
0082       inputTree12_(nullptr),
0083       m_nBins_(10000),
0084       m_rangeLow_(-.1),
0085       m_rangeHigh_(.1),
0086       firstEvent_(true),
0087       m_vtkmap_(13) {
0088   moduleListName_ = cfg.getUntrackedParameter<std::string>("moduleList");
0089 
0090   //input is ROOT
0091   inputFilename1_ = cfg.getUntrackedParameter<std::string>("inputROOTFile1");
0092   inputFilename2_ = cfg.getUntrackedParameter<std::string>("inputROOTFile2");
0093   inputTreenameAlign_ = cfg.getUntrackedParameter<std::string>("treeNameAlign");
0094   inputTreenameDeform_ = cfg.getUntrackedParameter<std::string>("treeNameDeform");
0095 
0096   //output file
0097   filename_ = cfg.getUntrackedParameter<std::string>("outputFile");
0098   //output dir for surface deformations
0099   surfdir_ = cfg.getUntrackedParameter<std::string>("surfDir");
0100 
0101   weightBy_ = cfg.getUntrackedParameter<std::string>("weightBy");
0102   setCommonTrackerSystem_ = cfg.getUntrackedParameter<std::string>("setCommonTrackerSystem");
0103   detIdFlag_ = cfg.getUntrackedParameter<bool>("detIdFlag");
0104   detIdFlagFile_ = cfg.getUntrackedParameter<std::string>("detIdFlagFile");
0105   weightById_ = cfg.getUntrackedParameter<bool>("weightById");
0106   weightByIdFile_ = cfg.getUntrackedParameter<std::string>("weightByIdFile");
0107 
0108   // if want to use, make id cut list
0109   if (detIdFlag_) {
0110     std::ifstream fin;
0111     fin.open(detIdFlagFile_.c_str());
0112 
0113     while (!fin.eof() && fin.good()) {
0114       uint32_t id;
0115       fin >> id;
0116       detIdFlagVector_.push_back(id);
0117     }
0118     fin.close();
0119   }
0120 
0121   // turn weightByIdFile into weightByIdVector
0122   if (weightById_) {
0123     std::ifstream inFile;
0124     inFile.open(weightByIdFile_.c_str());
0125     while (!inFile.eof()) {
0126       unsigned int listId;
0127       inFile >> listId;
0128       inFile.ignore(256, '\n');
0129 
0130       weightByIdVector_.push_back(listId);
0131     }
0132     inFile.close();
0133   }
0134 
0135   //root configuration
0136   theFile_ = new TFile(filename_.c_str(), "RECREATE");
0137   alignTree_ = new TTree("alignTree",
0138                          "alignTree");  //,"id:level:mid:mlevel:sublevel:x:y:z:r:phi:a:b:c:dx:dy:dz:dr:dphi:da:db:dc");
0139   alignTree_->Branch("id", &id_, "id/I");
0140   alignTree_->Branch("badModuleQuality", &badModuleQuality_, "badModuleQuality/I");
0141   alignTree_->Branch("inModuleList", &inModuleList_, "inModuleList/I");
0142   alignTree_->Branch("level", &level_, "level/I");
0143   alignTree_->Branch("mid", &mid_, "mid/I");
0144   alignTree_->Branch("mlevel", &mlevel_, "mlevel/I");
0145   alignTree_->Branch("sublevel", &sublevel_, "sublevel/I");
0146   alignTree_->Branch("x", &xVal_, "x/F");
0147   alignTree_->Branch("y", &yVal_, "y/F");
0148   alignTree_->Branch("z", &zVal_, "z/F");
0149   alignTree_->Branch("r", &rVal_, "r/F");
0150   alignTree_->Branch("phi", &phiVal_, "phi/F");
0151   alignTree_->Branch("eta", &etaVal_, "eta/F");
0152   alignTree_->Branch("alpha", &alphaVal_, "alpha/F");
0153   alignTree_->Branch("beta", &betaVal_, "beta/F");
0154   alignTree_->Branch("gamma", &gammaVal_, "gamma/F");
0155   alignTree_->Branch("dx", &dxVal_, "dx/F");
0156   alignTree_->Branch("dy", &dyVal_, "dy/F");
0157   alignTree_->Branch("dz", &dzVal_, "dz/F");
0158   alignTree_->Branch("dr", &drVal_, "dr/F");
0159   alignTree_->Branch("dphi", &dphiVal_, "dphi/F");
0160   alignTree_->Branch("dalpha", &dalphaVal_, "dalpha/F");
0161   alignTree_->Branch("dbeta", &dbetaVal_, "dbeta/F");
0162   alignTree_->Branch("dgamma", &dgammaVal_, "dgamma/F");
0163   alignTree_->Branch("du", &duVal_, "du/F");
0164   alignTree_->Branch("dv", &dvVal_, "dv/F");
0165   alignTree_->Branch("dw", &dwVal_, "dw/F");
0166   alignTree_->Branch("da", &daVal_, "da/F");
0167   alignTree_->Branch("db", &dbVal_, "db/F");
0168   alignTree_->Branch("dg", &dgVal_, "dg/F");
0169   alignTree_->Branch("useDetId", &useDetId_, "useDetId/I");
0170   alignTree_->Branch("detDim", &detDim_, "detDim/I");
0171   alignTree_->Branch("surW", &surWidth_, "surW/F");
0172   alignTree_->Branch("surL", &surLength_, "surL/F");
0173   alignTree_->Branch("surRot", &surRot_, "surRot[9]/D");
0174   alignTree_->Branch("identifiers", &identifiers_, "identifiers[6]/I");
0175   alignTree_->Branch("type", &type_, "type/I");
0176   alignTree_->Branch("surfDeform", &surfDeform_, "surfDeform[13]/D");
0177 
0178   for (std::vector<TrackerMap>::iterator it = m_vtkmap_.begin(); it != m_vtkmap_.end(); ++it) {
0179     it->setPalette(1);
0180     it->addPixel(true);
0181   }
0182 
0183   edm::Service<TFileService> fs;
0184   TFileDirectory subDir_All = fs->mkdir("AllSubdetectors");
0185   TFileDirectory subDir_PXB = fs->mkdir("PixelBarrel");
0186   TFileDirectory subDir_PXF = fs->mkdir("PixelEndcap");
0187   for (int ii = 0; ii < 13; ++ii) {
0188     std::stringstream histname0;
0189     histname0 << "SurfDeform_Par_" << ii;
0190     m_h1_[histname0.str()] = subDir_All.make<TH1D>(
0191         (histname0.str()).c_str(), (histname0.str()).c_str(), m_nBins_, m_rangeLow_, m_rangeHigh_);
0192 
0193     std::stringstream histname1;
0194     histname1 << "SurfDeform_PixelBarrel_Par_" << ii;
0195     m_h1_[histname1.str()] = subDir_PXB.make<TH1D>(
0196         (histname1.str()).c_str(), (histname1.str()).c_str(), m_nBins_, m_rangeLow_, m_rangeHigh_);
0197 
0198     std::stringstream histname2;
0199     histname2 << "SurfDeform_PixelEndcap_Par_" << ii;
0200     m_h1_[histname2.str()] = subDir_PXF.make<TH1D>(
0201         (histname2.str()).c_str(), (histname2.str()).c_str(), m_nBins_, m_rangeLow_, m_rangeHigh_);
0202   }
0203 }
0204 
0205 void TrackerGeometryCompare::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0206   edm::ParameterSetDescription desc;
0207   desc.setComment("Validates alignment payloads by comparing positions of tracker modules positiona and orientations");
0208   desc.addUntracked<std::vector<std::string> >("levels", {});
0209   desc.addUntracked<bool>("fromDD4hep", false);
0210   desc.addUntracked<bool>("writeToDB", false);
0211   desc.addUntracked<std::string>("moduleList", "moduleList.txt");
0212   desc.addUntracked<std::string>("inputROOTFile1", "IDEAL");
0213   desc.addUntracked<std::string>("inputROOTFile2", "idealtracker2.root");
0214   desc.addUntracked<std::string>("treeNameAlign", "alignTree");
0215   desc.addUntracked<std::string>("treeNameDeform", "alignTreeDeformations");
0216   desc.addUntracked<std::string>("outputFile", "output.root");
0217   desc.addUntracked<std::string>("surfDir", ".");
0218   desc.addUntracked<std::string>("weightBy", "DetUnit");
0219   desc.addUntracked<std::string>("setCommonTrackerSystem", "NONE");
0220   desc.addUntracked<bool>("detIdFlag", false);
0221   desc.addUntracked<std::string>("detIdFlagFile", "blah.txt");
0222   desc.addUntracked<bool>("weightById", false);
0223   desc.addUntracked<std::string>("weightByIdFile", "blah2.txt");
0224   descriptions.addWithDefaultLabel(desc);
0225 }
0226 
0227 void TrackerGeometryCompare::beginJob() { firstEvent_ = true; }
0228 
0229 void TrackerGeometryCompare::endJob() {
0230   int iname(0);
0231   for (std::vector<TrackerMap>::iterator it = m_vtkmap_.begin(); it != m_vtkmap_.end(); ++it) {
0232     std::stringstream mapname;
0233     mapname << surfdir_ << "/TkMap_SurfDeform_0" << iname << ".png";
0234     it->save(true, 0, 0, mapname.str());
0235     mapname.str(std::string());
0236     mapname.clear();
0237     mapname << surfdir_ << "/TkMap_SurfDeform_1" << iname << ".pdf";
0238     it->save(true, 0, 0, mapname.str());
0239     ++iname;
0240   }
0241 
0242   theFile_->cd();
0243   alignTree_->Write();
0244   theFile_->Close();
0245 }
0246 
0247 void TrackerGeometryCompare::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0248   if (firstEvent_) {
0249     //Retrieve tracker topology from geometry
0250     const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0251 
0252     //upload the ROOT geometries
0253     createROOTGeometry(iSetup);
0254 
0255     //setting the levels being used in the geometry comparator
0256     edm::LogInfo("TrackerGeometryCompare") << "levels: " << levelStrings_.size();
0257     for (const auto& level : levelStrings_) {
0258       m_theLevels.push_back(currentTracker->objectIdProvider().stringToId(level));
0259       edm::LogInfo("TrackerGeometryCompare") << "level: " << level;
0260       edm::LogInfo("TrackerGeometryCompare")
0261           << "structure type: " << currentTracker->objectIdProvider().stringToId(level);
0262     }
0263 
0264     //set common tracker system first
0265     // if setting the tracker common system
0266     if (setCommonTrackerSystem_ != "NONE") {
0267       setCommonTrackerSystem();
0268     }
0269 
0270     //compare the goemetries
0271     compareGeometries(referenceTracker, currentTracker, tTopo, iSetup);
0272     compareSurfaceDeformations(inputTree11_, inputTree12_);
0273 
0274     //write out ntuple
0275     //might be better to do within output module
0276 
0277     if (writeToDB_) {
0278       Alignments* myAlignments = currentTracker->alignments();
0279       AlignmentErrorsExtended* myAlignmentErrorsExtended = currentTracker->alignmentErrors();
0280 
0281       // 2. Store alignment[Error]s to DB
0282       edm::Service<cond::service::PoolDBOutputService> poolDbService;
0283       // Call service
0284       if (!poolDbService.isAvailable())  // Die if not available
0285         throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
0286 
0287       poolDbService->writeOneIOV<Alignments>(*myAlignments, poolDbService->beginOfTime(), "TrackerAlignmentRcd");
0288       poolDbService->writeOneIOV<AlignmentErrorsExtended>(
0289           *myAlignmentErrorsExtended, poolDbService->beginOfTime(), "TrackerAlignmentErrorExtendedRcd");
0290     }
0291 
0292     firstEvent_ = false;
0293   }
0294 }
0295 
0296 void TrackerGeometryCompare::createROOTGeometry(const edm::EventSetup& iSetup) {
0297   int inputRawId1, inputRawId2;
0298   double inputX1, inputY1, inputZ1, inputX2, inputY2, inputZ2;
0299   double inputAlpha1, inputBeta1, inputGamma1, inputAlpha2, inputBeta2, inputGamma2;
0300 
0301   //Retrieve tracker topology from geometry
0302   const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0303 
0304   // Fill module IDs from file into a list
0305   moduleListFile_.open(moduleListName_);
0306   if (moduleListFile_.is_open()) {
0307     std::string line;
0308     while (!moduleListFile_.eof()) {
0309       std::getline(moduleListFile_, line);
0310       moduleList_.push_back(std::atoi(line.c_str()));
0311     }
0312   } else {
0313     edm::LogInfo("TrackerGeometryCompare") << "Error: Module list not found! Please verify that given list exists!";
0314   }
0315 
0316   //declare alignments
0317   Alignments* alignments1 = new Alignments();
0318   AlignmentErrorsExtended* alignmentErrors1 = new AlignmentErrorsExtended();
0319   if (inputFilename1_ != "IDEAL") {
0320     inputRootFile1_ = new TFile(inputFilename1_.c_str());
0321     TTree* inputTree01_ = (TTree*)inputRootFile1_->Get(inputTreenameAlign_.c_str());
0322     inputTree01_->SetBranchAddress("rawid", &inputRawId1);
0323     inputTree01_->SetBranchAddress("x", &inputX1);
0324     inputTree01_->SetBranchAddress("y", &inputY1);
0325     inputTree01_->SetBranchAddress("z", &inputZ1);
0326     inputTree01_->SetBranchAddress("alpha", &inputAlpha1);
0327     inputTree01_->SetBranchAddress("beta", &inputBeta1);
0328     inputTree01_->SetBranchAddress("gamma", &inputGamma1);
0329 
0330     int nEntries1 = inputTree01_->GetEntries();
0331     //fill alignments
0332     for (int i = 0; i < nEntries1; ++i) {
0333       inputTree01_->GetEntry(i);
0334       CLHEP::Hep3Vector translation1(inputX1, inputY1, inputZ1);
0335       CLHEP::HepEulerAngles eulerangles1(inputAlpha1, inputBeta1, inputGamma1);
0336       uint32_t detid1 = inputRawId1;
0337       AlignTransform transform1(translation1, eulerangles1, detid1);
0338       alignments1->m_align.push_back(transform1);
0339 
0340       // dummy errors
0341       // APE matrix was 3x3, now it's 6x6 (because of muons), see PR #6483
0342       CLHEP::HepSymMatrix clhepSymMatrix(6, 0);
0343       AlignTransformErrorExtended transformError(clhepSymMatrix, detid1);
0344       alignmentErrors1->m_alignError.push_back(transformError);
0345     }
0346 
0347     // to get the right order, sort by rawId
0348     std::sort(alignments1->m_align.begin(), alignments1->m_align.end());
0349     std::sort(alignmentErrors1->m_alignError.begin(), alignmentErrors1->m_alignError.end());
0350   }
0351   //------------------
0352   Alignments* alignments2 = new Alignments();
0353   AlignmentErrorsExtended* alignmentErrors2 = new AlignmentErrorsExtended();
0354   if (inputFilename2_ != "IDEAL") {
0355     inputRootFile2_ = new TFile(inputFilename2_.c_str());
0356     TTree* inputTree02_ = (TTree*)inputRootFile2_->Get(inputTreenameAlign_.c_str());
0357     inputTree02_->SetBranchAddress("rawid", &inputRawId2);
0358     inputTree02_->SetBranchAddress("x", &inputX2);
0359     inputTree02_->SetBranchAddress("y", &inputY2);
0360     inputTree02_->SetBranchAddress("z", &inputZ2);
0361     inputTree02_->SetBranchAddress("alpha", &inputAlpha2);
0362     inputTree02_->SetBranchAddress("beta", &inputBeta2);
0363     inputTree02_->SetBranchAddress("gamma", &inputGamma2);
0364 
0365     int nEntries2 = inputTree02_->GetEntries();
0366     //fill alignments
0367     for (int i = 0; i < nEntries2; ++i) {
0368       inputTree02_->GetEntry(i);
0369       CLHEP::Hep3Vector translation2(inputX2, inputY2, inputZ2);
0370       CLHEP::HepEulerAngles eulerangles2(inputAlpha2, inputBeta2, inputGamma2);
0371       uint32_t detid2 = inputRawId2;
0372       AlignTransform transform2(translation2, eulerangles2, detid2);
0373       alignments2->m_align.push_back(transform2);
0374 
0375       // dummy errors
0376       // APE matrix was 3x3, now it's 6x6 (because of muons), see PR #6483
0377       CLHEP::HepSymMatrix clhepSymMatrix(6, 0);
0378       AlignTransformErrorExtended transformError(clhepSymMatrix, detid2);
0379       alignmentErrors2->m_alignError.push_back(transformError);
0380     }
0381 
0382     // to get the right order, sort by rawId
0383     std::sort(alignments2->m_align.begin(), alignments2->m_align.end());
0384     std::sort(alignmentErrors2->m_alignError.begin(), alignmentErrors2->m_alignError.end());
0385   }
0386 
0387   //accessing the initial geometry
0388   if (!fromDD4hep_) {
0389     edm::ESTransientHandle<DDCompactView> cpv = iSetup.getTransientHandle(cpvTokenDDD_);
0390   } else {
0391     edm::ESTransientHandle<cms::DDCompactView> cpv = iSetup.getTransientHandle(cpvTokenDD4hep_);
0392   }
0393 
0394   const GeometricDet* theGeometricDet = &iSetup.getData(geomDetToken_);
0395   const PTrackerParameters* ptp = &iSetup.getData(ptpToken_);
0396   const PTrackerAdditionalParametersPerDet* ptitp = &iSetup.getData(ptitpToken_);
0397   TrackerGeomBuilderFromGeometricDet trackerBuilder;
0398 
0399   //reference tracker
0400   TrackerGeometry* theRefTracker = trackerBuilder.build(theGeometricDet, ptitp, *ptp, tTopo);
0401   if (inputFilename1_ != "IDEAL") {
0402     GeometryAligner aligner1;
0403     aligner1.applyAlignments<TrackerGeometry>(
0404         &(*theRefTracker), &(*alignments1), &(*alignmentErrors1), AlignTransform());
0405   }
0406   referenceTracker = new AlignableTracker(&(*theRefTracker), tTopo);
0407   //referenceTracker->setSurfaceDeformation(surfDef1, true) ;
0408 
0409   //set tracker pixel Phase, assume current tracker is the same Phase
0410   if (theRefTracker->isThere(GeomDetEnumerators::P2PXB) || theRefTracker->isThere(GeomDetEnumerators::P2PXEC)) {
0411     phase_ = SiPixelPI::phase::two;
0412   } else if (theRefTracker->isThere(GeomDetEnumerators::P1PXB) || theRefTracker->isThere(GeomDetEnumerators::P1PXEC)) {
0413     phase_ = SiPixelPI::phase::one;
0414   } else {
0415     phase_ = SiPixelPI::phase::zero;
0416   }
0417 
0418   int inputRawid1;
0419   int inputRawid2;
0420   int inputDtype1, inputDtype2;
0421   std::vector<double> inputDpar1;
0422   std::vector<double> inputDpar2;
0423   std::vector<double>* p_inputDpar1 = &inputDpar1;
0424   std::vector<double>* p_inputDpar2 = &inputDpar2;
0425 
0426   const auto& comp1 = referenceTracker->deepComponents();
0427 
0428   SurfaceDeformation* surfDef1;
0429   if (inputFilename1_ != "IDEAL") {
0430     TTree* inputTree11_ = (TTree*)inputRootFile1_->Get(inputTreenameDeform_.c_str());
0431     inputTree11_->SetBranchAddress("irawid", &inputRawid1);
0432     inputTree11_->SetBranchAddress("dtype", &inputDtype1);
0433     inputTree11_->SetBranchAddress("dpar", &p_inputDpar1);
0434 
0435     unsigned int nEntries11 = inputTree11_->GetEntries();
0436     edm::LogInfo("TrackerGeometryCompare") << " nentries11 = " << nEntries11;
0437     for (unsigned int iEntry = 0; iEntry < nEntries11; ++iEntry) {
0438       inputTree11_->GetEntry(iEntry);
0439 
0440       surfDef1 = SurfaceDeformationFactory::create(inputDtype1, inputDpar1);
0441 
0442       if (int(comp1[iEntry]->id()) == inputRawid1) {
0443         comp1[iEntry]->setSurfaceDeformation(surfDef1, true);
0444       }
0445     }
0446   }
0447 
0448   //currernt tracker
0449   TrackerGeometry* theCurTracker = trackerBuilder.build(&*theGeometricDet, ptitp, *ptp, tTopo);
0450   if (inputFilename2_ != "IDEAL") {
0451     GeometryAligner aligner2;
0452     aligner2.applyAlignments<TrackerGeometry>(
0453         &(*theCurTracker), &(*alignments2), &(*alignmentErrors2), AlignTransform());
0454   }
0455   currentTracker = new AlignableTracker(&(*theCurTracker), tTopo);
0456 
0457   const auto& comp2 = currentTracker->deepComponents();
0458 
0459   SurfaceDeformation* surfDef2;
0460   if (inputFilename2_ != "IDEAL") {
0461     TTree* inputTree12_ = (TTree*)inputRootFile2_->Get(inputTreenameDeform_.c_str());
0462     inputTree12_->SetBranchAddress("irawid", &inputRawid2);
0463     inputTree12_->SetBranchAddress("dtype", &inputDtype2);
0464     inputTree12_->SetBranchAddress("dpar", &p_inputDpar2);
0465 
0466     unsigned int nEntries12 = inputTree12_->GetEntries();
0467     edm::LogInfo("TrackerGeometryCompare") << " nentries12 = " << nEntries12;
0468     for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
0469       inputTree12_->GetEntry(iEntry);
0470 
0471       surfDef2 = SurfaceDeformationFactory::create(inputDtype2, inputDpar2);
0472 
0473       if (int(comp2[iEntry]->id()) == inputRawid2) {
0474         comp2[iEntry]->setSurfaceDeformation(surfDef2, true);
0475       }
0476     }
0477   }
0478 
0479   delete alignments1;
0480   delete alignmentErrors1;
0481   delete alignments2;
0482   delete alignmentErrors2;
0483 }
0484 
0485 void TrackerGeometryCompare::compareSurfaceDeformations(TTree* refTree, TTree* curTree) {
0486   if (inputFilename1_ != "IDEAL" && inputFilename2_ != "IDEAL") {
0487     int inputRawid1;
0488     int inputRawid2;
0489     int inputSubdetid1, inputSubdetid2;
0490     int inputDtype1, inputDtype2;
0491     std::vector<double> inputDpar1;
0492     std::vector<double> inputDpar2;
0493     std::vector<double>* p_inputDpar1 = &inputDpar1;
0494     std::vector<double>* p_inputDpar2 = &inputDpar2;
0495 
0496     TTree* refTree = (TTree*)inputRootFile1_->Get(inputTreenameDeform_.c_str());
0497     refTree->SetBranchAddress("irawid", &inputRawid1);
0498     refTree->SetBranchAddress("subdetid", &inputSubdetid1);
0499     refTree->SetBranchAddress("dtype", &inputDtype1);
0500     refTree->SetBranchAddress("dpar", &p_inputDpar1);
0501 
0502     TTree* curTree = (TTree*)inputRootFile2_->Get(inputTreenameDeform_.c_str());
0503     curTree->SetBranchAddress("irawid", &inputRawid2);
0504     curTree->SetBranchAddress("subdetid", &inputSubdetid2);
0505     curTree->SetBranchAddress("dtype", &inputDtype2);
0506     curTree->SetBranchAddress("dpar", &p_inputDpar2);
0507 
0508     unsigned int nEntries11 = refTree->GetEntries();
0509     unsigned int nEntries12 = curTree->GetEntries();
0510 
0511     if (nEntries11 != nEntries12) {
0512       edm::LogError("TrackerGeometryCompare") << " Surface deformation parameters in two geometries differ!\n";
0513       return;
0514     }
0515 
0516     for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
0517       refTree->GetEntry(iEntry);
0518       curTree->GetEntry(iEntry);
0519       for (int ii = 0; ii < 13; ++ii) {
0520         surfDeform_[ii] = -1.0;
0521       }
0522       for (int npar = 0; npar < int(inputDpar2.size()); ++npar) {
0523         if (inputRawid1 == inputRawid2) {
0524           surfDeform_[npar] = inputDpar2.at(npar) - inputDpar1.at(npar);
0525           std::stringstream histname0;
0526           histname0 << "SurfDeform_Par_" << npar;
0527           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0528             m_h1_[histname0.str()]->Fill(surfDeform_[npar]);
0529           if (inputSubdetid1 == 1 && inputSubdetid2 == 1) {
0530             std::stringstream histname1;
0531             histname1 << "SurfDeform_PixelBarrel_Par_" << npar;
0532             if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0533               m_h1_[histname1.str()]->Fill(surfDeform_[npar]);
0534           }
0535           if (inputSubdetid1 == 2 && inputSubdetid2 == 2) {
0536             std::stringstream histname2;
0537             histname2 << "SurfDeform_PixelEndcap_Par_" << npar;
0538             if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0539               m_h1_[histname2.str()]->Fill(surfDeform_[npar]);
0540           }
0541           if (phase_ < SiPixelPI::phase::two) {
0542             if (phase_ == SiPixelPI::phase::zero) {
0543               (m_vtkmap_.at(npar)).fill_current_val(inputRawid1, surfDeform_[npar]);
0544             } else {
0545               // old-style tracker map does not support Phase1 Pixel
0546               if (!(inputSubdetid1 == 1 && inputSubdetid2 == 1) && !(inputSubdetid1 == 2 && inputSubdetid2 == 2)) {
0547                 (m_vtkmap_.at(npar)).fill_current_val(inputRawid1, surfDeform_[npar]);
0548               }
0549             }
0550           }
0551         }
0552       }
0553     }
0554 
0555   } else if (inputFilename1_ == "IDEAL" && inputFilename2_ != "IDEAL") {
0556     int inputRawid2;
0557     int inputSubdetid2;
0558     int inputDtype2;
0559     std::vector<double> inputDpar2;
0560     std::vector<double>* p_inputDpar2 = &inputDpar2;
0561 
0562     TTree* curTree = (TTree*)inputRootFile2_->Get(inputTreenameDeform_.c_str());
0563     curTree->SetBranchAddress("irawid", &inputRawid2);
0564     curTree->SetBranchAddress("subdetid", &inputSubdetid2);
0565     curTree->SetBranchAddress("dtype", &inputDtype2);
0566     curTree->SetBranchAddress("dpar", &p_inputDpar2);
0567 
0568     unsigned int nEntries12 = curTree->GetEntries();
0569     for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
0570       curTree->GetEntry(iEntry);
0571       for (int ii = 0; ii < 12; ++ii) {
0572         surfDeform_[ii] = -1.0;
0573       }
0574       for (int npar = 0; npar < int(inputDpar2.size()); ++npar) {
0575         surfDeform_[npar] = inputDpar2.at(npar);
0576         std::stringstream histname0;
0577         histname0 << "SurfDeform_Par_" << npar;
0578         if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0579           m_h1_[histname0.str()]->Fill(surfDeform_[npar]);
0580         if (inputSubdetid2 == 1) {
0581           std::stringstream histname1;
0582           histname1 << "SurfDeform_PixelBarrel_Par_" << npar;
0583           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0584             m_h1_[histname1.str()]->Fill(surfDeform_[npar]);
0585         }
0586         if (inputSubdetid2 == 2) {
0587           std::stringstream histname2;
0588           histname2 << "SurfDeform_PixelEndcap_Par_" << npar;
0589           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0590             m_h1_[histname2.str()]->Fill(surfDeform_[npar]);
0591         }
0592         if (phase_ < SiPixelPI::phase::two) {
0593           if (phase_ == SiPixelPI::phase::zero) {
0594             (m_vtkmap_.at(npar)).fill_current_val(inputRawid2, surfDeform_[npar]);
0595           } else {
0596             // old-style tracker map does not support Phase1 Pixel
0597             if (inputSubdetid2 != 1 && inputSubdetid2 != 2) {
0598               (m_vtkmap_.at(npar)).fill_current_val(inputRawid2, surfDeform_[npar]);
0599             }
0600           }
0601         }
0602       }
0603     }
0604 
0605   } else if (inputFilename1_ != "IDEAL" && inputFilename2_ == "IDEAL") {
0606     int inputRawid1;
0607     int inputSubdetid1;
0608     int inputDtype1;
0609     std::vector<double> inputDpar1;
0610     std::vector<double>* p_inputDpar1 = &inputDpar1;
0611 
0612     TTree* refTree = (TTree*)inputRootFile1_->Get(inputTreenameDeform_.c_str());
0613     refTree->SetBranchAddress("irawid", &inputRawid1);
0614     refTree->SetBranchAddress("subdetid", &inputSubdetid1);
0615     refTree->SetBranchAddress("dtype", &inputDtype1);
0616     refTree->SetBranchAddress("dpar", &p_inputDpar1);
0617     unsigned int nEntries11 = refTree->GetEntries();
0618 
0619     for (unsigned int iEntry = 0; iEntry < nEntries11; ++iEntry) {
0620       refTree->GetEntry(iEntry);
0621       for (int ii = 0; ii < 12; ++ii) {
0622         surfDeform_[ii] = -1.0;
0623       }
0624       for (int npar = 0; npar < int(inputDpar1.size()); ++npar) {
0625         surfDeform_[npar] = -inputDpar1.at(npar);
0626         std::stringstream histname0;
0627         histname0 << "SurfDeform_Par_" << npar;
0628         if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0629           m_h1_[histname0.str()]->Fill(surfDeform_[npar]);
0630         if (inputSubdetid1 == 1) {
0631           std::stringstream histname1;
0632           histname1 << "SurfDeform_PixelBarrel_Par_" << npar;
0633           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0634             m_h1_[histname1.str()]->Fill(surfDeform_[npar]);
0635         }
0636         if (inputSubdetid1 == 2) {
0637           std::stringstream histname2;
0638           histname2 << "SurfDeform_PixelEndcap_Par_" << npar;
0639           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0640             m_h1_[histname2.str()]->Fill(surfDeform_[npar]);
0641         }
0642         if (phase_ < SiPixelPI::phase::two) {
0643           if (phase_ == SiPixelPI::phase::zero) {
0644             (m_vtkmap_.at(npar)).fill_current_val(inputRawid1, surfDeform_[npar]);
0645           } else {
0646             // old-style tracker map does not support Phase1 Pixel
0647             if (inputSubdetid1 != 1 && inputSubdetid1 != 2) {
0648               (m_vtkmap_.at(npar)).fill_current_val(inputRawid1, surfDeform_[npar]);
0649             }
0650           }
0651         }
0652       }
0653     }
0654 
0655   } else if (inputFilename1_ == "IDEAL" && inputFilename2_ == "IDEAL") {
0656     edm::LogInfo("TrackerGeometryCompare") << ">>>> Comparing IDEAL with IDEAL: nothing to do! <<<<\n";
0657   }
0658   return;
0659 }
0660 
0661 void TrackerGeometryCompare::compareGeometries(Alignable* refAli,
0662                                                Alignable* curAli,
0663                                                const TrackerTopology* tTopo,
0664                                                const edm::EventSetup& iSetup) {
0665   using namespace align;
0666 
0667   const auto& refComp = refAli->components();
0668   const auto& curComp = curAli->components();
0669 
0670   unsigned int nComp = refComp.size();
0671   //only perform for designate levels
0672   bool useLevel = false;
0673   for (unsigned int i = 0; i < m_theLevels.size(); ++i) {
0674     if (refAli->alignableObjectId() == m_theLevels[i])
0675       useLevel = true;
0676   }
0677 
0678   //another added level for difference between det and detunit
0679   //if ((refAli->alignableObjectId()==2)&&(nComp == 1)) useLevel = false;
0680 
0681   //coordinate matching, etc etc
0682   if (useLevel) {
0683     DetId detid(refAli->id());
0684 
0685     CLHEP::Hep3Vector Rtotal, Wtotal, lRtotal, lWtotal;
0686     Rtotal.set(0., 0., 0.);
0687     Wtotal.set(0., 0., 0.);
0688     lRtotal.set(0., 0., 0.);
0689     lWtotal.set(0., 0., 0.);
0690 
0691     bool converged = false;
0692 
0693     AlgebraicVector diff, check;
0694 
0695     for (int i = 0; i < 100; i++) {
0696       // Get differences between alignments for rotations and translations
0697       // both local and global
0698       diff = align::diffAlignables(refAli, curAli, weightBy_, weightById_, weightByIdVector_);
0699 
0700       // 'diffAlignables' returns 'refAli - curAli' for translations and 'curAli - refAli' for rotations.
0701       // The plan is to unify this at some point, but a simple change of the sign for one of them was postponed
0702       // to do some further checks to understand the rotations better
0703       //Updated July 2018: as requested the sign in the translations has been changed to match the one in rotations. A test was done to change the diffAlignables function and solve the issue there, but proved quite time consuming. To unify the sign convention in the least amount of time the choice was made to change the sign here.
0704       CLHEP::Hep3Vector dR(-diff[0], -diff[1], -diff[2]);
0705       CLHEP::Hep3Vector dW(diff[3], diff[4], diff[5]);
0706       CLHEP::Hep3Vector dRLocal(-diff[6], -diff[7], -diff[8]);
0707       CLHEP::Hep3Vector dWLocal(diff[9], diff[10], diff[11]);
0708 
0709       // Translations
0710       Rtotal += dR;
0711       lRtotal += dRLocal;
0712 
0713       //Rotations
0714       CLHEP::HepRotation rot(Wtotal.unit(), Wtotal.mag());
0715       CLHEP::HepRotation drot(dW.unit(), dW.mag());
0716       rot *= drot;
0717       Wtotal.set(rot.axis().x() * rot.delta(), rot.axis().y() * rot.delta(), rot.axis().z() * rot.delta());
0718 
0719       CLHEP::HepRotation rotLocal(lWtotal.unit(), lWtotal.mag());
0720       CLHEP::HepRotation drotLocal(dWLocal.unit(), dWLocal.mag());
0721       rotLocal *= drotLocal;
0722       lWtotal.set(rotLocal.axis().x() * rotLocal.delta(),
0723                   rotLocal.axis().y() * rotLocal.delta(),
0724                   rotLocal.axis().z() * rotLocal.delta());
0725 
0726       // Move current alignable by shift and check if difference
0727       // is smaller than tolerance value
0728       // if true, break the loop
0729       align::moveAlignable(curAli, diff);
0730       float tolerance = 1e-7;
0731       check = align::diffAlignables(refAli, curAli, weightBy_, weightById_, weightByIdVector_);
0732       align::GlobalVector checkR(check[0], check[1], check[2]);
0733       align::GlobalVector checkW(check[3], check[4], check[5]);
0734       if ((checkR.mag() < tolerance) && (checkW.mag() < tolerance)) {
0735         converged = true;
0736         break;
0737       }
0738     }
0739 
0740     // give an exception if difference has not fallen below tolerance level
0741     // i.e. method has not converged
0742     if (!converged) {
0743       edm::LogInfo("TrackerGeometryCompare")
0744           << "Tolerance Exceeded!(alObjId: " << refAli->alignableObjectId()
0745           << ", rawId: " << refAli->geomDetId().rawId() << ", subdetId: " << detid.subdetId() << "): " << diff << check;
0746       throw cms::Exception("Tolerance in TrackerGeometryCompare exceeded");
0747     }
0748     AlgebraicVector TRtot(12);
0749     // global
0750     TRtot(1) = Rtotal.x();
0751     TRtot(2) = Rtotal.y();
0752     TRtot(3) = Rtotal.z();
0753     TRtot(4) = Wtotal.x();
0754     TRtot(5) = Wtotal.y();
0755     TRtot(6) = Wtotal.z();
0756     // local
0757     TRtot(7) = lRtotal.x();
0758     TRtot(8) = lRtotal.y();
0759     TRtot(9) = lRtotal.z();
0760     TRtot(10) = lWtotal.x();
0761     TRtot(11) = lWtotal.y();
0762     TRtot(12) = lWtotal.z();
0763 
0764     fillTree(refAli, TRtot, tTopo, iSetup);
0765   }
0766 
0767   // another added level for difference between det and detunit
0768   for (unsigned int i = 0; i < nComp; ++i)
0769     compareGeometries(refComp[i], curComp[i], tTopo, iSetup);
0770 }
0771 
0772 void TrackerGeometryCompare::setCommonTrackerSystem() {
0773   edm::LogInfo("TrackerGeometryCompare") << "Setting Common Tracker System....";
0774 
0775   // DM_534??AlignableObjectId dummy;
0776   // DM_534??_commonTrackerLevel = dummy.nameToType(_setCommonTrackerSystem);
0777   commonTrackerLevel_ = currentTracker->objectIdProvider().stringToId(setCommonTrackerSystem_);
0778 
0779   diffCommonTrackerSystem(referenceTracker, currentTracker);
0780 
0781   align::EulerAngles dOmega(3);
0782   dOmega[0] = TrackerCommonR_.x();
0783   dOmega[1] = TrackerCommonR_.y();
0784   dOmega[2] = TrackerCommonR_.z();
0785   align::RotationType rot = align::toMatrix(dOmega);
0786   align::GlobalVector theR = TrackerCommonT_;
0787 
0788   edm::LogInfo("TrackerGeometryCompare") << "what we get from overlaying the pixels..." << theR << ", " << rot;
0789 
0790   //transform to the Tracker System
0791   align::PositionType trackerCM = currentTracker->globalPosition();
0792   align::GlobalVector cmDiff(
0793       trackerCM.x() - TrackerCommonCM_.x(), trackerCM.y() - TrackerCommonCM_.y(), trackerCM.z() - TrackerCommonCM_.z());
0794 
0795   edm::LogInfo("TrackerGeometryCompare") << "Pixel CM: " << TrackerCommonCM_ << ", tracker CM: " << trackerCM;
0796 
0797   //adjust translational difference factoring in different rotational CM
0798   //needed because rotateInGlobalFrame is about CM of alignable, not Tracker
0799   const align::GlobalVector::BasicVectorType& lpvgf = cmDiff.basicVector();
0800   align::GlobalVector moveV(rot.multiplyInverse(lpvgf) - lpvgf);
0801   align::GlobalVector theRprime(theR + moveV);
0802 
0803   AlgebraicVector TrackerCommonTR(6);
0804   TrackerCommonTR(1) = theRprime.x();
0805   TrackerCommonTR(2) = theRprime.y();
0806   TrackerCommonTR(3) = theRprime.z();
0807   TrackerCommonTR(4) = TrackerCommonR_.x();
0808   TrackerCommonTR(5) = TrackerCommonR_.y();
0809   TrackerCommonTR(6) = TrackerCommonR_.z();
0810 
0811   edm::LogInfo("TrackerGeometryCompare") << "and after the transformation: " << TrackerCommonTR;
0812   align::moveAlignable(currentTracker, TrackerCommonTR);
0813 }
0814 
0815 void TrackerGeometryCompare::diffCommonTrackerSystem(Alignable* refAli, Alignable* curAli) {
0816   const auto& refComp = refAli->components();
0817   const auto& curComp = curAli->components();
0818   unsigned int nComp = refComp.size();
0819   //only perform for designate levels
0820   bool useLevel = false;
0821   if (refAli->alignableObjectId() == commonTrackerLevel_)
0822     useLevel = true;
0823 
0824   //useLevel = false;
0825   if (useLevel) {
0826     CLHEP::Hep3Vector Rtotal, Wtotal;
0827     Rtotal.set(0., 0., 0.);
0828     Wtotal.set(0., 0., 0.);
0829 
0830     AlgebraicVector diff = align::diffAlignables(refAli, curAli, weightBy_, weightById_, weightByIdVector_);
0831     CLHEP::Hep3Vector dR(diff[0], diff[1], diff[2]);
0832     Rtotal += dR;
0833     CLHEP::Hep3Vector dW(diff[3], diff[4], diff[5]);
0834     CLHEP::HepRotation rot(Wtotal.unit(), Wtotal.mag());
0835     CLHEP::HepRotation drot(dW.unit(), dW.mag());
0836     rot *= drot;
0837     Wtotal.set(rot.axis().x() * rot.delta(), rot.axis().y() * rot.delta(), rot.axis().z() * rot.delta());
0838 
0839     TrackerCommonT_ = align::GlobalVector(Rtotal.x(), Rtotal.y(), Rtotal.z());
0840     TrackerCommonR_ = align::GlobalVector(Wtotal.x(), Wtotal.y(), Wtotal.z());
0841     TrackerCommonCM_ = curAli->globalPosition();
0842 
0843   } else {
0844     for (unsigned int i = 0; i < nComp; ++i)
0845       diffCommonTrackerSystem(refComp[i], curComp[i]);
0846   }
0847 }
0848 
0849 void TrackerGeometryCompare::fillTree(Alignable* refAli,
0850                                       const AlgebraicVector& diff,
0851                                       const TrackerTopology* tTopo,
0852                                       const edm::EventSetup& iSetup) {
0853   //Get bad modules
0854   const SiPixelQuality* SiPixelModules = &iSetup.getData(pixQualityToken_);
0855   const SiStripQuality* SiStripModules = &iSetup.getData(stripQualityToken_);
0856 
0857   id_ = refAli->id();
0858 
0859   badModuleQuality_ = 0;
0860   //check if module has a bad quality tag
0861   if (SiPixelModules->IsModuleBad(id_)) {
0862     badModuleQuality_ = 1;
0863   }
0864   if (SiStripModules->IsModuleBad(id_)) {
0865     badModuleQuality_ = 1;
0866   }
0867 
0868   //check if module is in a given list of bad/untouched etc. modules
0869   inModuleList_ = 0;
0870   for (unsigned int i = 0; i < moduleList_.size(); i++) {
0871     if (moduleList_[i] == id_) {
0872       inModuleList_ = 1;
0873       break;
0874     }
0875   }
0876 
0877   level_ = refAli->alignableObjectId();
0878   //need if ali has no mother
0879   if (refAli->mother()) {
0880     mid_ = refAli->mother()->geomDetId().rawId();
0881     mlevel_ = refAli->mother()->alignableObjectId();
0882   } else {
0883     mid_ = -1;
0884     mlevel_ = -1;
0885   }
0886   DetId detid(id_);
0887   sublevel_ = detid.subdetId();
0888   fillIdentifiers(sublevel_, id_, tTopo);
0889   xVal_ = refAli->globalPosition().x();
0890   yVal_ = refAli->globalPosition().y();
0891   zVal_ = refAli->globalPosition().z();
0892   align::GlobalVector vec(xVal_, yVal_, zVal_);
0893   rVal_ = vec.perp();
0894   phiVal_ = vec.phi();
0895   etaVal_ = vec.eta();
0896   align::RotationType rot = refAli->globalRotation();
0897   align::EulerAngles eulerAngles = align::toAngles(rot);
0898   alphaVal_ = eulerAngles[0];
0899   betaVal_ = eulerAngles[1];
0900   gammaVal_ = eulerAngles[2];
0901   // global
0902   dxVal_ = diff[0];
0903   dyVal_ = diff[1];
0904   dzVal_ = diff[2];
0905   // local
0906   duVal_ = diff[6];
0907   dvVal_ = diff[7];
0908   dwVal_ = diff[8];
0909   //...TODO...
0910   align::GlobalVector g(dxVal_, dyVal_, dzVal_);
0911   //getting dR and dPhi
0912   align::GlobalVector vRef(xVal_, yVal_, zVal_);
0913   align::GlobalVector vCur(xVal_ + dxVal_, yVal_ + dyVal_, zVal_ + dzVal_);
0914   drVal_ = vCur.perp() - vRef.perp();
0915   dphiVal_ = vCur.phi() - vRef.phi();
0916   // global
0917   dalphaVal_ = diff[3];
0918   dbetaVal_ = diff[4];
0919   dgammaVal_ = diff[5];
0920   // local
0921   daVal_ = diff[9];
0922   dbVal_ = diff[10];
0923   dgVal_ = diff[11];
0924   //detIdFlag
0925   if (refAli->alignableObjectId() == align::AlignableDetUnit) {
0926     if (detIdFlag_) {
0927       if ((passIdCut(refAli->id())) || (passIdCut(refAli->mother()->id()))) {
0928         useDetId_ = 1;
0929       } else {
0930         useDetId_ = 0;
0931       }
0932     }
0933   }
0934   // det module dimension
0935   if (refAli->alignableObjectId() == align::AlignableDetUnit) {
0936     if (refAli->mother()->alignableObjectId() != align::AlignableDet)
0937       detDim_ = 1;
0938     else if (refAli->mother()->alignableObjectId() == align::AlignableDet)
0939       detDim_ = 2;
0940   } else
0941     detDim_ = 0;
0942 
0943   surWidth_ = refAli->surface().width();
0944   surLength_ = refAli->surface().length();
0945   align::RotationType rt = refAli->globalRotation();
0946   surRot_[0] = rt.xx();
0947   surRot_[1] = rt.xy();
0948   surRot_[2] = rt.xz();
0949   surRot_[3] = rt.yx();
0950   surRot_[4] = rt.yy();
0951   surRot_[5] = rt.yz();
0952   surRot_[6] = rt.zx();
0953   surRot_[7] = rt.zy();
0954   surRot_[8] = rt.zz();
0955 
0956   //Fill
0957   alignTree_->Fill();
0958 }
0959 
0960 void TrackerGeometryCompare::surveyToTracker(AlignableTracker* ali,
0961                                              Alignments* alignVals,
0962                                              AlignmentErrorsExtended* alignErrors) {
0963   //getting the right alignables for the alignment record
0964   auto detPB = ali->pixelHalfBarrelGeomDets();
0965   auto detPEC = ali->pixelEndcapGeomDets();
0966   auto detTIB = ali->innerBarrelGeomDets();
0967   auto detTID = ali->TIDGeomDets();
0968   auto detTOB = ali->outerBarrelGeomDets();
0969   auto detTEC = ali->endcapGeomDets();
0970 
0971   align::Alignables allGeomDets;
0972   std::copy(detPB.begin(), detPB.end(), std::back_inserter(allGeomDets));
0973   std::copy(detPEC.begin(), detPEC.end(), std::back_inserter(allGeomDets));
0974   std::copy(detTIB.begin(), detTIB.end(), std::back_inserter(allGeomDets));
0975   std::copy(detTID.begin(), detTID.end(), std::back_inserter(allGeomDets));
0976   std::copy(detTOB.begin(), detTOB.end(), std::back_inserter(allGeomDets));
0977   std::copy(detTEC.begin(), detTEC.end(), std::back_inserter(allGeomDets));
0978 
0979   align::Alignables rcdAlis;
0980   for (const auto& i : allGeomDets) {
0981     if (i->components().size() == 1) {
0982       rcdAlis.push_back(i);
0983     } else if (i->components().size() > 1) {
0984       rcdAlis.push_back(i);
0985       const auto& comp = i->components();
0986       for (const auto& j : comp)
0987         rcdAlis.push_back(j);
0988     }
0989   }
0990 
0991   //turning them into alignments
0992   for (const auto& k : rcdAlis) {
0993     const SurveyDet* surveyInfo = k->survey();
0994     const align::PositionType& pos(surveyInfo->position());
0995     align::RotationType rot(surveyInfo->rotation());
0996     CLHEP::Hep3Vector clhepVector(pos.x(), pos.y(), pos.z());
0997     CLHEP::HepRotation clhepRotation(
0998         CLHEP::HepRep3x3(rot.xx(), rot.xy(), rot.xz(), rot.yx(), rot.yy(), rot.yz(), rot.zx(), rot.zy(), rot.zz()));
0999     AlignTransform transform(clhepVector, clhepRotation, k->id());
1000     AlignTransformErrorExtended transformError(CLHEP::HepSymMatrix(3, 1), k->id());
1001     alignVals->m_align.push_back(transform);
1002     alignErrors->m_alignError.push_back(transformError);
1003   }
1004 
1005   // to get the right order, sort by rawId
1006   std::sort(alignVals->m_align.begin(), alignVals->m_align.end());
1007   std::sort(alignErrors->m_alignError.begin(), alignErrors->m_alignError.end());
1008 }
1009 
1010 void TrackerGeometryCompare::addSurveyInfo(Alignable* ali) {
1011   const auto& comp = ali->components();
1012 
1013   unsigned int nComp = comp.size();
1014 
1015   for (unsigned int i = 0; i < nComp; ++i)
1016     addSurveyInfo(comp[i]);
1017 
1018   const SurveyError& error = theSurveyErrors->m_surveyErrors[theSurveyIndex];
1019 
1020   if (ali->geomDetId().rawId() != error.rawId() || ali->alignableObjectId() != error.structureType()) {
1021     throw cms::Exception("DatabaseError") << "Error reading survey info from DB. Mismatched id!";
1022   }
1023 
1024   const CLHEP::Hep3Vector& pos = theSurveyValues->m_align[theSurveyIndex].translation();
1025   const CLHEP::HepRotation& rot = theSurveyValues->m_align[theSurveyIndex].rotation();
1026 
1027   AlignableSurface surf(
1028       align::PositionType(pos.x(), pos.y(), pos.z()),
1029       align::RotationType(rot.xx(), rot.xy(), rot.xz(), rot.yx(), rot.yy(), rot.yz(), rot.zx(), rot.zy(), rot.zz()));
1030 
1031   surf.setWidth(ali->surface().width());
1032   surf.setLength(ali->surface().length());
1033 
1034   ali->setSurvey(new SurveyDet(surf, error.matrix()));
1035 
1036   ++theSurveyIndex;
1037 }
1038 
1039 bool TrackerGeometryCompare::passIdCut(uint32_t id) {
1040   bool pass = false;
1041   int nEntries = detIdFlagVector_.size();
1042 
1043   for (int i = 0; i < nEntries; i++) {
1044     if (detIdFlagVector_[i] == id)
1045       pass = true;
1046   }
1047 
1048   return pass;
1049 }
1050 
1051 void TrackerGeometryCompare::fillIdentifiers(int subdetlevel, int rawid, const TrackerTopology* tTopo) {
1052   switch (subdetlevel) {
1053     case 1: {
1054       identifiers_[0] = tTopo->pxbModule(rawid);
1055       identifiers_[1] = tTopo->pxbLadder(rawid);
1056       identifiers_[2] = tTopo->pxbLayer(rawid);
1057       identifiers_[3] = 999;
1058       identifiers_[4] = 999;
1059       identifiers_[5] = 999;
1060       break;
1061     }
1062     case 2: {
1063       identifiers_[0] = tTopo->pxfModule(rawid);
1064       identifiers_[1] = tTopo->pxfPanel(rawid);
1065       identifiers_[2] = tTopo->pxfBlade(rawid);
1066       identifiers_[3] = tTopo->pxfDisk(rawid);
1067       identifiers_[4] = tTopo->pxfSide(rawid);
1068       identifiers_[5] = 999;
1069       break;
1070     }
1071     case 3: {
1072       identifiers_[0] = tTopo->tibModule(rawid);
1073       identifiers_[1] = tTopo->tibStringInfo(rawid)[0];
1074       identifiers_[2] = tTopo->tibStringInfo(rawid)[1];
1075       identifiers_[3] = tTopo->tibStringInfo(rawid)[2];
1076       identifiers_[4] = tTopo->tibLayer(rawid);
1077       identifiers_[5] = 999;
1078       break;
1079     }
1080     case 4: {
1081       identifiers_[0] = tTopo->tidModuleInfo(rawid)[0];
1082       identifiers_[1] = tTopo->tidModuleInfo(rawid)[1];
1083       identifiers_[2] = tTopo->tidRing(rawid);
1084       identifiers_[3] = tTopo->tidWheel(rawid);
1085       identifiers_[4] = tTopo->tidSide(rawid);
1086       identifiers_[5] = 999;
1087       break;
1088     }
1089     case 5: {
1090       identifiers_[0] = tTopo->tobModule(rawid);
1091       identifiers_[1] = tTopo->tobRodInfo(rawid)[0];
1092       identifiers_[2] = tTopo->tobRodInfo(rawid)[1];
1093       identifiers_[3] = tTopo->tobLayer(rawid);
1094       identifiers_[4] = 999;
1095       identifiers_[5] = 999;
1096       break;
1097     }
1098     case 6: {
1099       identifiers_[0] = tTopo->tecModule(rawid);
1100       identifiers_[1] = tTopo->tecRing(rawid);
1101       identifiers_[2] = tTopo->tecPetalInfo(rawid)[0];
1102       identifiers_[3] = tTopo->tecPetalInfo(rawid)[1];
1103       identifiers_[4] = tTopo->tecWheel(rawid);
1104       identifiers_[5] = tTopo->tecSide(rawid);
1105       break;
1106     }
1107     default: {
1108       edm::LogInfo("TrackerGeometryCompare") << "Error: bad subdetid!!";
1109       break;
1110     }
1111   }
1112 }
1113 
1114 DEFINE_FWK_MODULE(TrackerGeometryCompare);
[ 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!