SurveyAnalysis/test/TestConverter2.cpp

0001 // -*- C++ -*-
0002 //
0003 // Package:    TestConverter2
0004 // Class:      TestConverter2
0005 //
0006 //
0007 // Description: Module to test the SurveyConverter software
0008 //
0009 //
0010 // Original Author:  Roberto Covarelli
0011 //         Created:  March 16, 2006
0012 //
0013
0014 // system include files
0015 #include <string>
0016 #include "TTree.h"
0017 #include "TFile.h"
0018 // #include "TRotMatrix.h"
0019
0020 // user include files
0021 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028
0029 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0030 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0031 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0032
0033 #include "CondFormats/Alignment/interface/Alignments.h"
0034 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0035 #include "CLHEP/Vector/RotationInterfaces.h"
0036 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0037 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
0038
0039 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0040 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0041 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0042 //
0043 //
0044 // class declaration
0045 //
0046 using namespace std;
0047
0048 class TestConverter2 : public edm::one::EDAnalyzer<> {
0049 public:
0050   explicit TestConverter2(const edm::ParameterSet&);
0051   ~TestConverter2();
0052
0053   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0054
0055 private:
0056   // ----------member data ---------------------------
0057   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0058   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0059   const edm::ESGetToken<Alignments, TrackerAlignmentRcd> aliToken_;
0060   const edm::ESGetToken<AlignmentErrorsExtended, TrackerAlignmentErrorExtendedRcd> aliErrToken_;
0061
0062   TTree* theTree;
0063   TFile* theFile;
0064   edm::ParameterSet theParameterSet;
0065   float dx_, dy_, dz_;
0066   float dtx_, dty_, dtz_;
0067   float dkx_, dky_, dkz_;
0068   float dnx_, dny_, dnz_;
0069   float errx_, erry_, errz_;
0070   int subdid_, fwbw_, layerdisk_, frontback_, stringrod_, petal_, module_;
0071   // TRotMatrix* rot_;
0072 };
0073
0074 //
0075 // constructors and destructor
0076 //
0077 TestConverter2::TestConverter2(const edm::ParameterSet& iConfig)
0078     : tTopoToken_(esConsumes()),
0079       tkGeomToken_(esConsumes()),
0080       aliToken_(esConsumes()),
0081       aliErrToken_(esConsumes()),
0082       theParameterSet(iConfig) {
0083   // Open root file and define tree
0084   std::string fileName = theParameterSet.getUntrackedParameter<std::string>("fileName", "test.root");
0085   theFile = new TFile(fileName.c_str(), "RECREATE");
0086   theTree = new TTree("theTree", "Detector units positions");
0087
0088   theTree->Branch("subDetId", &subdid_, "subDetId/I");
0089   theTree->Branch("Fw_Bw", &fwbw_, "Fw_Bw/I");
0090   theTree->Branch("FrBa_InEx", &frontback_, "FrBa_inEx/I");
0091   theTree->Branch("LayerDisk", &layerdisk_, "LayerDisk/I");
0092   theTree->Branch("Petal", &petal_, "Petal/I");
0093   theTree->Branch("StRingRod", &stringrod_, "StRingRod/I");
0094   theTree->Branch("Module", &module_, "Module/I");
0095   theTree->Branch("dx", &dx_, "dx/F");
0096   theTree->Branch("dy", &dy_, "dy/F");
0097   theTree->Branch("dz", &dz_, "dz/F");
0098   theTree->Branch("dtx", &dtx_, "dtx/F");
0099   theTree->Branch("dty", &dty_, "dty/F");
0100   theTree->Branch("dtz", &dtz_, "dtz/F");
0101   theTree->Branch("dkx", &dkx_, "dkx/F");
0102   theTree->Branch("dky", &dky_, "dky/F");
0103   theTree->Branch("dkz", &dkz_, "dkz/F");
0104   theTree->Branch("dnx", &dnx_, "dnx/F");
0105   theTree->Branch("dny", &dny_, "dny/F");
0106   theTree->Branch("dnz", &dnz_, "dnz/F");
0107   theTree->Branch("errx", &errx_, "errx/F");
0108   theTree->Branch("erry", &erry_, "erry/F");
0109   theTree->Branch("errz", &errz_, "errz/F");
0110 }
0111
0112 TestConverter2::~TestConverter2() {
0113   theTree->Write();
0114   theFile->Close();
0115 }
0116
0117 void TestConverter2::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0118   //Retrieve tracker topology from geometry
0119
0120   edm::LogInfo("TrackerAlignment") << "Starting!";
0121   const TrackerTopology* const tTopo = &iSetup.getData(tTopoToken_);
0122   //
0123   // Retrieve tracker geometry from event setup
0124   const TrackerGeometry* trackerGeometry = &iSetup.getData(tkGeomToken_);
0125
0126   // Retrieve alignment[Error]s from DBase
0127   const Alignments* alignments = &iSetup.getData(aliToken_);
0128   const AlignmentErrorsExtended* alignmentErrors = &iSetup.getData(aliErrToken_);
0129   auto alignErrors = alignmentErrors->m_alignError;
0130
0131   // Now loop on detector units, and store difference position and orientation w.r.t. survey
0132
0133   for (auto iGeomDet = alignments->m_align.begin(); iGeomDet != alignments->m_align.end(); iGeomDet++) {
0134     for (auto iDet = trackerGeometry->dets().begin(); iDet != trackerGeometry->dets().end(); iDet++) {
0135       // std::cout << (*iDet)->geographicalId().rawId() << " " << (*iGeomDet).rawId() << std::endl;
0136       if ((*iDet)->geographicalId().rawId() == (*iGeomDet).rawId()) {
0137         for (auto it = alignErrors.begin(); it != alignErrors.end(); it++) {
0138           if ((*it).rawId() == (*iGeomDet).rawId()) {
0139             DetId thisId((*iGeomDet).rawId());
0140
0141             if (thisId.subdetId() == int(StripSubdetector::TIB)) {
0142               subdid_ = 3;
0143               layerdisk_ = tTopo->tibLayer(thisId);
0144               std::vector<unsigned int> theString = tTopo->tibStringInfo(thisId);
0145               fwbw_ = theString[0];
0146               frontback_ = theString[1];
0147               stringrod_ = theString[2];
0148               petal_ = 0;
0149               module_ = tTopo->tibModule(thisId);
0150
0151             } else if (thisId.subdetId() == int(StripSubdetector::TID)) {
0152               subdid_ = 4;
0153               layerdisk_ = tTopo->tidWheel(thisId);
0154               std::vector<unsigned int> theModule = tTopo->tidModuleInfo(thisId);
0155               frontback_ = theModule[0];
0156               module_ = theModule[1];
0157               petal_ = 0;
0158               fwbw_ = tTopo->tidSide(thisId);
0159               stringrod_ = tTopo->tidRing(thisId);
0160
0161             } else if (thisId.subdetId() == int(StripSubdetector::TOB)) {
0162               subdid_ = 5;
0163               layerdisk_ = tTopo->tobLayer(thisId);
0164               std::vector<unsigned int> theRod = tTopo->tobRodInfo(thisId);
0165               fwbw_ = theRod[0];
0166               stringrod_ = theRod[1];
0167               petal_ = 0;
0168               frontback_ = 0;
0169               module_ = tTopo->tobModule(thisId);
0170
0171             } else if (thisId.subdetId() == int(StripSubdetector::TEC)) {
0172               subdid_ = 6;
0173               layerdisk_ = tTopo->tecWheel(thisId);
0174               std::vector<unsigned int> thePetal = tTopo->tecPetalInfo(thisId);
0175               frontback_ = thePetal[0];
0176               petal_ = thePetal[1];
0177               fwbw_ = tTopo->tecSide(thisId);
0178               stringrod_ = tTopo->tecRing(thisId);
0179               module_ = tTopo->tecModule(thisId);
0180
0181             } else {
0182               std::cout << "WARNING!!! this DetId (" << thisId.rawId() << ") does not belong to SiStrip tracker"
0183                         << std::endl;
0184               break;
0185             }
0186
0187             // CLHEP::HepRotation fromAngles( (*iGeomDet).eulerAngles()  );
0188             CLHEP::HepRotation fromAngles((*iGeomDet).rotation());
0189             align::RotationType rotation(fromAngles.xx(),
0190                                          fromAngles.xy(),
0191                                          fromAngles.xz(),
0192                                          fromAngles.yx(),
0193                                          fromAngles.yy(),
0194                                          fromAngles.yz(),
0195                                          fromAngles.zx(),
0196                                          fromAngles.zy(),
0197                                          fromAngles.zz());
0198
0199             dx_ = (*iGeomDet).translation().x() - (*iDet)->position().x();
0200             dy_ = (*iGeomDet).translation().y() - (*iDet)->position().y();
0201             dz_ = (*iGeomDet).translation().z() - (*iDet)->position().z();
0202             dtx_ = rotation.xx() - (*iDet)->rotation().xx();
0203             dty_ = rotation.xy() - (*iDet)->rotation().xy();
0204             dtz_ = rotation.xz() - (*iDet)->rotation().xz();
0205             dkx_ = rotation.yx() - (*iDet)->rotation().yx();
0206             dky_ = rotation.yy() - (*iDet)->rotation().yy();
0207             dkz_ = rotation.yz() - (*iDet)->rotation().yz();
0208             dnx_ = rotation.zx() - (*iDet)->rotation().zx();
0209             dny_ = rotation.zy() - (*iDet)->rotation().zy();
0210             dnz_ = rotation.zz() - (*iDet)->rotation().zz();
0211             CLHEP::HepSymMatrix errMat = (*it).matrix();
0212             errx_ = sqrt(errMat[0][0]);
0213             erry_ = sqrt(errMat[1][1]);
0214             errz_ = sqrt(errMat[2][2]);
0215
0216             theTree->Fill();
0217
0218             cout << "DetId = " << (*iGeomDet).rawId() << " " << endl;
0219             cout << "DetId decodified = " << subdid_ << " " << fwbw_ << " " << frontback_ << " " << layerdisk_ << " "
0220                  << stringrod_ << " " << petal_ << " " << module_ << endl;
0221             cout << "X pos : TRACKER_MOVED = " << std::fixed << std::setprecision(2) << (*iGeomDet).translation().x()
0222                  << " / TRACKER_ORIGINAL = " << (*iDet)->position().x() << endl;
0223             cout << "Y pos : TRACKER_MOVED = " << std::fixed << std::setprecision(2) << (*iGeomDet).translation().y()
0224                  << " / TRACKER_ORIGINAL = " << (*iDet)->position().y() << endl;
0225             cout << "Z pos : TRACKER_MOVED = " << std::fixed << std::setprecision(2) << (*iGeomDet).translation().z()
0226                  << " / TRACKER_ORIGINAL = " << (*iDet)->position().z() << endl;
0227             cout << "SPATIAL DISTANCE = " << std::fixed << std::setprecision(3)
0228                  << sqrt(pow(dx_, 2) + pow(dy_, 2) + pow(dz_, 2)) << endl;
0229             cout << "Trans vect X : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.xx()
0230                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().xx() << endl;
0231             cout << "Trans vect Y : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.xy()
0232                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().xy() << endl;
0233             cout << "Trans vect Z : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.xz()
0234                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().xz() << endl;
0235             cout << "Long vect X : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.yx()
0236                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().yx() << endl;
0237             cout << "Long vect Y : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.yy()
0238                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().yy() << endl;
0239             cout << "Long vect Z : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.yz()
0240                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().yz() << endl;
0241             cout << "Norm vect X : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.zx()
0242                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().zx() << endl;
0243             cout << "Norm vect Y : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.zy()
0244                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().zy() << endl;
0245             cout << "Norm vect Z : TRACKER_MOVED = " << std::fixed << std::setprecision(5) << rotation.zz()
0246                  << " / TRACKER_ORIGINAL = " << (*iDet)->rotation().zz() << endl;
0247           }
0248         }
0249       }
0250     }
0251   }
0252
0253   edm::LogInfo("TrackerAlignment") << "Done!";
0254 }
0255
0256 //define this as a plug-in
0257 DEFINE_FWK_MODULE(TestConverter2);