/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/Alignment/SurveyAnalysis/src/DTSurvey.cc

Go to the documentation of this file.

#include <fstream>

#include "Alignment/SurveyAnalysis/interface/DTSurveyChamber.h"
#include "DataFormats/MuonDetId/interface/DTChamberId.h" 
#include "FWCore/Framework/interface/ESHandle.h"
#include "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/DTGeometry/interface/DTChamber.h" 

#include "Alignment/SurveyAnalysis/interface/DTSurvey.h"
#include <iostream>

DTSurvey::DTSurvey(const std::string& Wheel, const std::string& Chambers, int n) {
  
  nameOfWheelInfoFile = Wheel;
  nameOfChamberInfoFile = Chambers;
  id = n;   
  
  FillWheelInfo();
}


DTSurvey::~DTSurvey() {
  delete [] chambers;
}


void DTSurvey::CalculateChambers() {
  for(int stationCounter = 0; stationCounter < 4; stationCounter++) {
    for(int sectorCounter = 0; sectorCounter < 14; sectorCounter++) {
      if(chambers[stationCounter][sectorCounter]->getNumberPoints() > 2) {
        chambers[stationCounter][sectorCounter]->compute();
      }
    }
  }
}


const DTSurveyChamber * DTSurvey::getChamber(int station, int sector) const {return chambers[station][sector];}

void DTSurvey::ReadChambers(edm::ESHandle<DTGeometry> pDD) {
  
  //Create the chambers
  chambers = new DTSurveyChamber ** [4];
  for (int cont_stat = 0; cont_stat < 4; cont_stat++) {
    chambers[cont_stat] = new DTSurveyChamber * [14];
    for(int cont_sect = 0; cont_sect < 14; cont_sect++) {
      DTChamberId mId(id, cont_stat+1, cont_sect+1);
      chambers[cont_stat][cont_sect] = new DTSurveyChamber(id, cont_stat+1, cont_sect+1, mId.rawId());
    }
  }

  std::cout << nameOfChamberInfoFile << std::endl;
  std::ifstream file(nameOfChamberInfoFile.c_str());
  while(!file.eof()) {
    int code, station, sector;
    double x, y, z, rms, dx, dy, dz;
    file >> code >> x >> y >> z >> rms >> dx >> dy >> dz;
    if(file.eof()) break;
    x = x/10.0; y=y/10.0; z=z/10.0; dx=dx/10.0; dy=dy/10.0; dz=dz/10.0;rms=rms/10.0; 
    station = code/10000 - 1;
    sector = (code-(station+1)*10000)/100 - 1;
    //De momento vamos a actuar como si no hubiera otra forma de resolver esto
    TMatrixD r(3,1);
    r(0,0) = x; r(1,0) = y;r(2,0) = z+OffsetZ;
    TMatrixD disp(3,1);
    disp(0,0) = dx; disp(1,0) = dy; disp(2,0) = dz;
    TMatrixD rp = Rot*r-delta;
    disp = disp-r+rp;
    
    GlobalPoint rg(r(0,0), r(1,0), r(2,0));
    GlobalPoint rt(r(0,0)-disp(0,0), r(1,0)-disp(1,0), r(2,0)-disp(2,0));
    DTChamberId mId(id, station+1, sector+1);
    const DTChamber *mChamber = static_cast<const DTChamber *>(pDD->idToDet(mId));
    LocalPoint rl = mChamber->toLocal(rg);
    LocalPoint rtl = mChamber->toLocal(rt);
    TMatrixD rLocal(3,1);
    rLocal(0,0) = rl.x(); rLocal(1,0) = rl.y(); rLocal(2,0) = rl.z();
    TMatrixD rTeo(3,1);
    rTeo(0,0) = rtl.x(); rTeo(1,0) = rtl.y(); rTeo(2,0) = rtl.z();
    TMatrixD diff = rLocal-rTeo;
    TMatrixD errors(3,1);
    errors(0,0) = rms; errors(1,0) = rms; errors(2,0) = rms;
    chambers[station][sector]->addPoint(code, rLocal, diff, errors);
  }
  file.close();
}

/*
void DTSurvey::ToDB(MuonAlignment *myMuonAlignment) {

  for(int station = 0; station < 4; station++) {
    for(int sector = 0; sector < 14; sector++) {
      if(chambers[station][sector]->getNumberPoints() > 2) {
        std::vector<float> displacements;
        std::vector<float> rotations;
        displacements.push_back(chambers[station][sector]->getDeltaX());
        displacements.push_back(chambers[station][sector]->getDeltaY());
        displacements.push_back(chambers[station][sector]->getDeltaZ());
        rotations.push_back(chambers[station][sector]->getAlpha());
        rotations.push_back(chambers[station][sector]->getBeta());
        rotations.push_back(chambers[station][sector]->getGamma());
        DTChamberId mId(id, station+1, sector+1);
        myMuonAlignment->moveAlignableLocalCoord(mId, displacements, rotations);
      }
    }
  }
}
*/

void DTSurvey::FillWheelInfo() {

  std::ifstream wheeltowheel(nameOfWheelInfoFile.c_str());
  float zOffset, deltax, deltay, deltaz, alpha, beta, gamma;
  wheeltowheel >> zOffset >> deltax >> deltay >> deltaz >> alpha >> beta >> gamma;
  wheeltowheel.close();
  
  OffsetZ = zOffset;

  //Build displacement vector
  delta.ResizeTo(3,1);
  delta(0,0) = deltax/10.0;
  delta(1,0) = deltay/10.0;
  delta(2,0) = deltaz/10.0;
  
  //Build rotation matrix
  Rot.ResizeTo(3,3);
  TMatrixD alpha_m(3,3);
  TMatrixD beta_m(3,3);
  TMatrixD gamma_m(3,3);
  alpha_m.Zero();
  beta_m.Zero();
  gamma_m.Zero();
  for(int k = 0; k < 3; k++) {
    alpha_m(k,k) = 1.0;
    beta_m(k,k) = 1.0;
    gamma_m(k,k) = 1.0;
  }
  alpha /= 1000.0; //New scale: angles in radians
  beta /= 1000.0;
  gamma /= 1000.0;
  alpha_m(1,1) = cos(alpha);
  alpha_m(1,2) = sin(alpha);
  alpha_m(2,1) = -sin(alpha);
  alpha_m(2,2) = cos(alpha);
  beta_m(0,0) = cos(beta);
  beta_m(0,2) = -sin(beta);
  beta_m(2,0) = sin(beta);
  beta_m(2,2) = cos(beta);
  gamma_m(0,0) = cos(gamma);
  gamma_m(0,1) = sin(gamma);
  gamma_m(1,0) = -sin(gamma);
  gamma_m(1,1) = cos(gamma);
  Rot = alpha_m*beta_m*gamma_m;
}


std::ostream &operator<<(std::ostream & flux, const DTSurvey& obj) {

  for(int stationCounter = 0; stationCounter < 4; stationCounter++) {
    for(int sectorCounter = 0; sectorCounter < 14; sectorCounter++) {
      if(obj.getChamber(stationCounter,sectorCounter)->getNumberPoints() > 2) {
        const DTSurveyChamber *m_chamber = obj.getChamber(stationCounter, sectorCounter);
        flux << *m_chamber;
      }
    }
  }
  return flux;
}