Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-21 00:19:07 

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