SurveyAnalysis/src/DTSurvey.cc

0001 #include <fstream>
0002
0003 #include "Alignment/SurveyAnalysis/interface/DTSurveyChamber.h"
0004 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0007 #include "Geometry/DTGeometry/interface/DTChamber.h"
0008
0009 #include "Alignment/SurveyAnalysis/interface/DTSurvey.h"
0010 #include <iostream>
0011
0012 DTSurvey::DTSurvey(const std::string &Wheel, const std::string &Chambers, int n) : chambers(nullptr) {
0013   nameOfWheelInfoFile = Wheel;
0014   nameOfChamberInfoFile = Chambers;
0015   id = n;
0016
0017   FillWheelInfo();
0018 }
0019
0020 DTSurvey::~DTSurvey() { delete[] chambers; }
0021
0022 void DTSurvey::CalculateChambers() {
0023   for (int stationCounter = 0; stationCounter < 4; stationCounter++) {
0024     for (int sectorCounter = 0; sectorCounter < 14; sectorCounter++) {
0025       if (chambers[stationCounter][sectorCounter]->getNumberPoints() > 2) {
0026         chambers[stationCounter][sectorCounter]->compute();
0027       }
0028     }
0029   }
0030 }
0031
0032 const DTSurveyChamber *DTSurvey::getChamber(int station, int sector) const { return chambers[station][sector]; }
0033
0034 void DTSurvey::ReadChambers(edm::ESHandle<DTGeometry> pDD) {
0035   //Create the chambers
0036   chambers = new DTSurveyChamber **[4];
0037   for (int cont_stat = 0; cont_stat < 4; cont_stat++) {
0038     chambers[cont_stat] = new DTSurveyChamber *[14];
0039     for (int cont_sect = 0; cont_sect < 14; cont_sect++) {
0040       DTChamberId mId(id, cont_stat + 1, cont_sect + 1);
0041       chambers[cont_stat][cont_sect] = new DTSurveyChamber(id, cont_stat + 1, cont_sect + 1, mId.rawId());
0042     }
0043   }
0044
0045   std::cout << nameOfChamberInfoFile << std::endl;
0046   std::ifstream file(nameOfChamberInfoFile.c_str());
0047   while (!file.eof()) {
0048     int code, station, sector;
0049     double x, y, z, rms, dx, dy, dz;
0050     file >> code >> x >> y >> z >> rms >> dx >> dy >> dz;
0051     if (file.eof())
0052       break;
0053     x = x / 10.0;
0054     y = y / 10.0;
0055     z = z / 10.0;
0056     dx = dx / 10.0;
0057     dy = dy / 10.0;
0058     dz = dz / 10.0;
0059     rms = rms / 10.0;
0060     station = code / 10000 - 1;
0061     sector = (code - (station + 1) * 10000) / 100 - 1;
0062     //De momento vamos a actuar como si no hubiera otra forma de resolver esto
0063     TMatrixD r(3, 1);
0064     r(0, 0) = x;
0065     r(1, 0) = y;
0066     r(2, 0) = z + OffsetZ;
0067     TMatrixD disp(3, 1);
0068     disp(0, 0) = dx;
0069     disp(1, 0) = dy;
0070     disp(2, 0) = dz;
0071     TMatrixD rp = Rot * r - delta;
0072     disp = disp - r + rp;
0073
0074     GlobalPoint rg(r(0, 0), r(1, 0), r(2, 0));
0075     GlobalPoint rt(r(0, 0) - disp(0, 0), r(1, 0) - disp(1, 0), r(2, 0) - disp(2, 0));
0076     DTChamberId mId(id, station + 1, sector + 1);
0077     const DTChamber *mChamber = static_cast<const DTChamber *>(pDD->idToDet(mId));
0078     LocalPoint rl = mChamber->toLocal(rg);
0079     LocalPoint rtl = mChamber->toLocal(rt);
0080     TMatrixD rLocal(3, 1);
0081     rLocal(0, 0) = rl.x();
0082     rLocal(1, 0) = rl.y();
0083     rLocal(2, 0) = rl.z();
0084     TMatrixD rTeo(3, 1);
0085     rTeo(0, 0) = rtl.x();
0086     rTeo(1, 0) = rtl.y();
0087     rTeo(2, 0) = rtl.z();
0088     TMatrixD diff = rLocal - rTeo;
0089     TMatrixD errors(3, 1);
0090     errors(0, 0) = rms;
0091     errors(1, 0) = rms;
0092     errors(2, 0) = rms;
0093     chambers[station][sector]->addPoint(code, rLocal, diff, errors);
0094   }
0095   file.close();
0096 }
0097
0098 /*
0099 void DTSurvey::ToDB(MuonAlignment *myMuonAlignment) {
0100
0101   for(int station = 0; station < 4; station++) {
0102     for(int sector = 0; sector < 14; sector++) {
0103       if(chambers[station][sector]->getNumberPoints() > 2) {
0104         std::vector<float> displacements;
0105         std::vector<float> rotations;
0106         displacements.push_back(chambers[station][sector]->getDeltaX());
0107         displacements.push_back(chambers[station][sector]->getDeltaY());
0108         displacements.push_back(chambers[station][sector]->getDeltaZ());
0109         rotations.push_back(chambers[station][sector]->getAlpha());
0110         rotations.push_back(chambers[station][sector]->getBeta());
0111         rotations.push_back(chambers[station][sector]->getGamma());
0112         DTChamberId mId(id, station+1, sector+1);
0113         myMuonAlignment->moveAlignableLocalCoord(mId, displacements, rotations);
0114       }
0115     }
0116   }
0117 }
0118 */
0119
0120 void DTSurvey::FillWheelInfo() {
0121   std::ifstream wheeltowheel(nameOfWheelInfoFile.c_str());
0122   float zOffset, deltax, deltay, deltaz, alpha, beta, gamma;
0123   wheeltowheel >> zOffset >> deltax >> deltay >> deltaz >> alpha >> beta >> gamma;
0124   wheeltowheel.close();
0125
0126   OffsetZ = zOffset;
0127
0128   //Build displacement vector
0129   delta.ResizeTo(3, 1);
0130   delta(0, 0) = deltax / 10.0;
0131   delta(1, 0) = deltay / 10.0;
0132   delta(2, 0) = deltaz / 10.0;
0133
0134   //Build rotation matrix
0135   Rot.ResizeTo(3, 3);
0136   TMatrixD alpha_m(3, 3);
0137   TMatrixD beta_m(3, 3);
0138   TMatrixD gamma_m(3, 3);
0139   alpha_m.Zero();
0140   beta_m.Zero();
0141   gamma_m.Zero();
0142   for (int k = 0; k < 3; k++) {
0143     alpha_m(k, k) = 1.0;
0144     beta_m(k, k) = 1.0;
0145     gamma_m(k, k) = 1.0;
0146   }
0147   alpha /= 1000.0;  //New scale: angles in radians
0148   beta /= 1000.0;
0149   gamma /= 1000.0;
0150   alpha_m(1, 1) = cos(alpha);
0151   alpha_m(1, 2) = sin(alpha);
0152   alpha_m(2, 1) = -sin(alpha);
0153   alpha_m(2, 2) = cos(alpha);
0154   beta_m(0, 0) = cos(beta);
0155   beta_m(0, 2) = -sin(beta);
0156   beta_m(2, 0) = sin(beta);
0157   beta_m(2, 2) = cos(beta);
0158   gamma_m(0, 0) = cos(gamma);
0159   gamma_m(0, 1) = sin(gamma);
0160   gamma_m(1, 0) = -sin(gamma);
0161   gamma_m(1, 1) = cos(gamma);
0162   Rot = alpha_m * beta_m * gamma_m;
0163 }
0164
0165 std::ostream &operator<<(std::ostream &flux, const DTSurvey &obj) {
0166   for (int stationCounter = 0; stationCounter < 4; stationCounter++) {
0167     for (int sectorCounter = 0; sectorCounter < 14; sectorCounter++) {
0168       if (obj.getChamber(stationCounter, sectorCounter)->getNumberPoints() > 2) {
0169         const DTSurveyChamber *m_chamber = obj.getChamber(stationCounter, sectorCounter);
0170         flux << *m_chamber;
0171       }
0172     }
0173   }
0174   return flux;
0175 }