SurveyAnalysis/test/TestIdealGeometry2.cpp

0001 // -*- C++ -*-
0002 //
0003 // Package:    TestIdealGeometry2
0004 // Class:      TestIdealGeometry2
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 "CondFormats/Alignment/interface/Alignments.h"
0030 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0031
0032 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0033 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0034 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0035
0036 #include "Alignment/SurveyAnalysis/interface/SurveyDataReader.h"
0037 //
0038 //
0039 // class declaration
0040 //
0041
0042 class TestIdealGeometry2 : public edm::one::EDAnalyzer<> {
0043   typedef SurveyDataReader::MapType MapType;
0044   typedef SurveyDataReader::PairType PairType;
0045   typedef SurveyDataReader::MapTypeOr MapTypeOr;
0046   typedef SurveyDataReader::PairTypeOr PairTypeOr;
0047
0048 public:
0049   explicit TestIdealGeometry2(const edm::ParameterSet&);
0050   ~TestIdealGeometry2();
0051
0052   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0053
0054 private:
0055   // ----------member data ---------------------------
0056   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0057   const edm::ESGetToken<Alignments, TrackerAlignmentRcd> aliToken_;
0058
0059   TTree* theTree;
0060   TFile* theFile;
0061   edm::ParameterSet theParameterSet;
0062   float dx_, dy_, dz_;
0063   float dtx_, dty_, dtz_;
0064   float dkx_, dky_, dkz_;
0065   float dnx_, dny_, dnz_;
0066   int Id_;
0067   // TRotMatrix* rot_;
0068
0069   static const int NFILES = 2;
0070 };
0071
0072 //
0073 // constructors and destructor
0074 //
0075 TestIdealGeometry2::TestIdealGeometry2(const edm::ParameterSet& iConfig)
0076     : tTopoToken_(esConsumes()), aliToken_(esConsumes()), theParameterSet(iConfig) {
0077   // Open root file and define tree
0078   std::string fileName = theParameterSet.getUntrackedParameter<std::string>("fileName", "testideal.root");
0079   theFile = new TFile(fileName.c_str(), "RECREATE");
0080   theTree = new TTree("theTree", "Detector units positions");
0081
0082   theTree->Branch("Id", &Id_, "Id/I");
0083   theTree->Branch("dx", &dx_, "dx/F");
0084   theTree->Branch("dy", &dy_, "dy/F");
0085   theTree->Branch("dz", &dz_, "dz/F");
0086   theTree->Branch("dtx", &dtx_, "dtx/F");
0087   theTree->Branch("dty", &dty_, "dty/F");
0088   theTree->Branch("dtz", &dtz_, "dtz/F");
0089   theTree->Branch("dkx", &dkx_, "dkx/F");
0090   theTree->Branch("dky", &dky_, "dky/F");
0091   theTree->Branch("dkz", &dkz_, "dkz/F");
0092   theTree->Branch("dnx", &dnx_, "dnx/F");
0093   theTree->Branch("dny", &dny_, "dny/F");
0094   theTree->Branch("dnz", &dnz_, "dnz/F");
0095 }
0096
0097 TestIdealGeometry2::~TestIdealGeometry2() {
0098   theTree->Write();
0099   theFile->Close();
0100 }
0101
0102 void TestIdealGeometry2::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0103   //Retrieve tracker topology from geometry
0104   const TrackerTopology* const tTopo = &iSetup.getData(tTopoToken_);
0105   edm::LogInfo("TrackerAlignment") << "Starting!";
0106
0107   //
0108   // Read in the survey information from the text files
0109   //
0110   edm::ParameterSet textFiles = theParameterSet.getParameter<edm::ParameterSet>("textFileNames");
0111   std::string textFileNames[NFILES];
0112   std::string fileType[NFILES];
0113   textFileNames[0] = textFiles.getUntrackedParameter<std::string>("forTIB", "NONE");
0114   fileType[0] = "TIB";
0115   textFileNames[1] = textFiles.getUntrackedParameter<std::string>("forTID", "NONE");
0116   fileType[1] = "TID";
0117
0118   SurveyDataReader dataReader;
0119   for (int ii = 0; ii < NFILES; ii++) {
0120     if (textFileNames[ii] == "NONE")
0121       throw cms::Exception("BadConfig") << fileType[ii] << " input file not found in configuration";
0122     dataReader.readFile(textFileNames[ii], fileType[ii], tTopo);
0123   }
0124
0125   edm::LogInfo("TrackerAlignment") << "Files read";
0126
0127   const MapTypeOr& theSurveyMap = dataReader.surveyMap();
0128
0129   edm::LogInfo("TrackerAlignment") << "Map written";
0130
0131   //
0132   // Retrieve tracker geometry from event setup
0133   //
0134   // edm::ESHandle<TrackerGeometry> trackerGeometry;
0135   // iSetup.get<TrackerDigiGeometryRecord>().get( trackerGeometry );
0136
0137   // Retrieve alignment[Error]s from DBase
0138   const Alignments* alignments = &iSetup.getData(aliToken_);
0139
0140   int countDet = 0;
0141
0142   // Now loop on detector units, and store difference position and orientation w.r.t. survey
0143
0144   for (std::vector<AlignTransform>::const_iterator iGeomDet = alignments->m_align.begin();
0145        iGeomDet != alignments->m_align.end();
0146        iGeomDet++) {
0147     if (countDet == 0)
0148       countDet = countDet + 3;
0149     unsigned int comparisonVect[6] = {0, 0, 0, 0, 0, 0};
0150
0151     DetId thisId((*iGeomDet).rawId());
0152     if (thisId.subdetId() == int(StripSubdetector::TIB)) {
0153       comparisonVect[0] = int(StripSubdetector::TIB);
0154
0155       comparisonVect[1] = tTopo->tibLayer(thisId);
0156       if (comparisonVect[1] < 3) {
0157         countDet--;
0158       } else {
0159         countDet = countDet - 3;
0160       }
0161       std::vector<unsigned int> theString = tTopo->tibStringInfo(thisId);
0162       comparisonVect[2] = theString[0];
0163       comparisonVect[3] = theString[1];
0164       comparisonVect[4] = theString[2];
0165       comparisonVect[5] = tTopo->tibModule(thisId);
0166
0167     } else if (thisId.subdetId() == int(StripSubdetector::TID)) {
0168       comparisonVect[0] = int(StripSubdetector::TID);
0169
0170       comparisonVect[1] = tTopo->tidSide(thisId);
0171       comparisonVect[2] = tTopo->tidWheel(thisId);
0172       comparisonVect[3] = tTopo->tidRing(thisId);
0173       if (comparisonVect[3] < 3) {
0174         countDet--;
0175       } else {
0176         countDet = countDet - 3;
0177       }
0178       std::vector<unsigned int> theModule = tTopo->tidModuleInfo(thisId);
0179       comparisonVect[4] = theModule[0];
0180       comparisonVect[5] = theModule[1];
0181     }
0182
0183     if (countDet == 0) {  // Store only r-phi for double-sided modules
0184
0185       for (MapTypeOr::const_iterator it = theSurveyMap.begin(); it != theSurveyMap.end(); it++) {
0186         const std::vector<int>& locPos = (it)->first;
0187         const align::Scalars& align_params = (it)->second;
0188
0189         if (locPos[0] == int(comparisonVect[0]) && locPos[1] == int(comparisonVect[1]) &&
0190             locPos[2] == int(comparisonVect[2]) && locPos[3] == int(comparisonVect[3]) &&
0191             locPos[4] == int(comparisonVect[4]) && locPos[5] == int(comparisonVect[5])) {
0192           const CLHEP::HepRotation& rot = (*iGeomDet).rotation();
0193           align::RotationType rotation(
0194               rot.xx(), rot.xy(), rot.xz(), rot.yx(), rot.yy(), rot.yz(), rot.zx(), rot.zy(), rot.zz());
0195
0196           Id_ = (*iGeomDet).rawId();
0197           // cout << "DetId = " << Id_ << " " << endl;
0198           // cout << "DetId decodified = " << comparisonVect[0] << " " << comparisonVect[1] << " " << comparisonVect[2] << " " << comparisonVect[3] << " " << comparisonVect[4] << " " << comparisonVect[5] << endl;
0199           dx_ = (*iGeomDet).translation().x() - align_params[0];
0200           // cout << "X pos : TRACKER_ALIGN = " << std::fixed << std::setprecision(2) << (*iGeomDet).translation().x() << " / IDEAL RICCARDO = " << align_params[0] << endl;
0201           dy_ = (*iGeomDet).translation().y() - align_params[1];
0202           // cout << "Y pos : TRACKER_ALIGN = " << std::fixed << std::setprecision(2) << (*iGeomDet).translation().y() << " / IDEAL RICCARDO = " << align_params[1] << endl;
0203           dz_ = (*iGeomDet).translation().z() - align_params[2];
0204           // cout << "Z pos : TRACKER_ALIGN = " << std::fixed << std::setprecision(2) << (*iGeomDet).translation().z() << " / IDEAL RICCARDO = " << align_params[2] << endl;
0205           // cout << "SPATIAL DISTANCE = " << std::fixed << std::setprecision(3) << sqrt(pow(dx_,2)+pow(dy_,2)+pow(dz_,2)) << endl;
0206           dtx_ = rotation.xx() - align_params[6];
0207           // cout << "Trans vect X : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.xx() << " / IDEAL RICCARDO = " << align_params[6] << endl;
0208           dty_ = rotation.xy() - align_params[7];
0209           // cout << "Trans vect Y : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.xy() << " / IDEAL RICCARDO = " << align_params[7] << endl;
0210           dtz_ = rotation.xz() - align_params[8];
0211           // cout << "Trans vect Z : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.xz() << " / IDEAL RICCARDO = " << align_params[8] << endl;
0212           dkx_ = rotation.yx() - align_params[9];
0213           // cout << "Long vect X : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.yx() << " / IDEAL RICCARDO = " << align_params[9] << endl;
0214           dky_ = rotation.yy() - align_params[10];
0215           // cout << "Long vect Y : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.yy() << " / IDEAL RICCARDO = " << align_params[10] << endl;
0216           dkz_ = rotation.yz() - align_params[11];
0217           // cout << "Long vect Z : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.yz() << " / IDEAL RICCARDO = " << align_params[11] << endl;
0218           dnx_ = rotation.zx() - align_params[3];
0219           // cout << "Norm vect X : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.zx() << " / IDEAL RICCARDO = " << align_params[3] << endl;
0220           dny_ = rotation.zy() - align_params[4];
0221           // cout << "Norm vect Y : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.zy() << " / IDEAL RICCARDO = " << align_params[4] << endl;
0222           dnz_ = rotation.zz() - align_params[5];
0223           // cout << "Norm vect Z : TRACKER_ALIGN = " << std::fixed << std::setprecision(3) << rotation.zz() << " / IDEAL RICCARDO = " << align_params[5] << endl;
0224           theTree->Fill();
0225         }
0226       }
0227     }
0228   }
0229
0230   edm::LogInfo("TrackerAlignment") << "Done!";
0231 }
0232
0233 //define this as a plug-in
0234 DEFINE_FWK_MODULE(TestIdealGeometry2);