SurveyAnalysis/test/TestConverter.cpp

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