SurveyInputCSCfromPins.cc

Go to the documentation of this file.

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Alignment/CommonAlignment/interface/SurveyDet.h"
#include "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "Alignment/MuonAlignment/interface/AlignableMuon.h"
#include "FWCore/Framework/interface/ESHandle.h" 
#include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "DataFormats/DetId/interface/DetId.h"

#include "Alignment/MuonAlignment/interface/AlignableCSCChamber.h"
#include "Alignment/MuonAlignment/interface/AlignableCSCStation.h"

#include <iostream>
#include <fstream>
#include "TFile.h"
#include "TTree.h"
#include "TRandom3.h"

#include "SurveyInputCSCfromPins.h"

#define SQR(x) ((x)*(x))

SurveyInputCSCfromPins::SurveyInputCSCfromPins(const edm::ParameterSet& cfg)
  : m_pinPositions(cfg.getParameter<std::string>("pinPositions"))
  , m_rootFile(cfg.getParameter<std::string>("rootFile"))
  , m_verbose(cfg.getParameter<bool>("verbose"))
  , m_errorX(cfg.getParameter<double>("errorX"))
  , m_errorY(cfg.getParameter<double>("errorY"))
  , m_errorZ(cfg.getParameter<double>("errorZ"))
  , m_missingErrorTranslation(cfg.getParameter<double>("missingErrorTranslation"))
  , m_missingErrorAngle(cfg.getParameter<double>("missingErrorAngle"))
  , m_stationErrorX(cfg.getParameter<double>("stationErrorX"))
  , m_stationErrorY(cfg.getParameter<double>("stationErrorY"))
  , m_stationErrorZ(cfg.getParameter<double>("stationErrorZ"))
  , m_stationErrorPhiX(cfg.getParameter<double>("stationErrorPhiX"))
  , m_stationErrorPhiY(cfg.getParameter<double>("stationErrorPhiY"))
  , m_stationErrorPhiZ(cfg.getParameter<double>("stationErrorPhiZ"))
{}

void SurveyInputCSCfromPins::orient(LocalVector LC1, LocalVector LC2, double a, double b, double &T, double &dx, double &dy, double &dz, double &PhX, double &PhZ) {
   double cosPhX, sinPhX, cosPhZ, sinPhZ;

   LocalPoint LP1(LC1.x(), LC1.y() + a, LC1.z() + b);
   LocalPoint LP2(LC2.x(), LC2.y() - a, LC2.z() + b);
  
   LocalPoint P((LP1.x() - LP2.x())/(2.*a), (LP1.y() - LP2.y())/(2.*a), (LP1.z() - LP2.z())/(2.*a));
   LocalPoint Pp((LP1.x() + LP2.x())/(2.), (LP1.y() + LP2.y())/(2.), (LP1.z() + LP2.z())/(2.));
    
   T=P.mag();
    
   sinPhX=P.z()/T;
   cosPhX=sqrt(1-SQR(sinPhX));
   cosPhZ=P.y()/(T*cosPhX);
   sinPhZ=-P.x()/(T*cosPhZ);
        
   PhX=atan2(sinPhX,cosPhX);
    
   PhZ=atan2(sinPhZ,cosPhZ);
    
   dx=Pp.x()-sinPhZ*sinPhX*b;
   dy=Pp.y()+cosPhZ*sinPhX*b;
   dz=Pp.z()-cosPhX*b;
}

void SurveyInputCSCfromPins::errors(double a, double b, bool missing1, bool missing2, double &dx_dx, double &dy_dy, double &dz_dz, double &phix_phix, double &phiz_phiz, double &dy_phix) {
   dx_dx = 0.;
   dy_dy = 0.;
   dz_dz = 0.;
   phix_phix = 0.;
   phiz_phiz = 0.;
   dy_phix = 0.;

   const double trials = 10000.;  // two significant digits

   for (int i = 0;  i < trials;  i++) {
      LocalVector LC1, LC2;
      if (missing1) {
     LC1 = LocalVector(gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation));
      }
      else {
         LC1 = LocalVector(gRandom->Gaus(0., m_errorX), gRandom->Gaus(0., m_errorY), gRandom->Gaus(0., m_errorZ));
      }

      if (missing2) {
     LC2 = LocalVector(gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation));
      }
      else {
     LC2 = LocalVector(gRandom->Gaus(0., m_errorX), gRandom->Gaus(0., m_errorY), gRandom->Gaus(0., m_errorZ));
      }

      double dx, dy, dz, PhX, PhZ, T;
      orient(LC1, LC2, a, b, T, dx, dy, dz, PhX, PhZ);

      dx_dx += dx * dx;
      dy_dy += dy * dy;
      dz_dz += dz * dz;
      phix_phix += PhX * PhX;
      phiz_phiz += PhZ * PhZ;
      dy_phix += dy * PhX;      // the only non-zero off-diagonal element
   }

   dx_dx /= trials;
   dy_dy /= trials;
   dz_dz /= trials;
   phix_phix /= trials;
   phiz_phiz /= trials;
   dy_phix /= trials;
}

void SurveyInputCSCfromPins::analyze(const edm::Event&, const edm::EventSetup& iSetup)
{
  if (theFirstEvent) {

    edm::LogInfo("SurveyInputCSCfromPins") << "***************ENTERING INITIALIZATION******************" << "  \n";
    
    std::ifstream in;
    in.open(m_pinPositions.c_str());
  
    Double_t x1, y1, z1, x2, y2, z2, a, b, tot=0.0, maxErr=0.0, h, s1, dx, dy, dz, PhX, PhZ, T;
   
    int ID1, ID2, ID3, ID4, ID5, i=1, ii=0;
    
    TFile *file1 = new TFile(m_rootFile.c_str(),"recreate");
    TTree *tree1 = new TTree("tree1","alignment pins");
      
    if (m_verbose) {
   
        tree1->Branch("displacement_x_pin1_cm", &x1, "x1/D");
        tree1->Branch("displacement_y_pin1_cm", &y1, "y1/D");
        tree1->Branch("displacement_z_pin1_cm", &z1, "z1/D");
        tree1->Branch("displacement_x_pin2_cm", &x2, "x2/D");
        tree1->Branch("displacement_y_pin2_cm", &y2, "y2/D");
        tree1->Branch("displacement_z_pin2_cm", &z2, "z2/D");     
        tree1->Branch("error_vector_length_cm",&h, "h/D"); 
        tree1->Branch("stretch_diff_cm",&s1, "s1/D");
        tree1->Branch("stretch_factor",&T, "T/D");
        tree1->Branch("chamber_displacement_x_cm",&dx, "dx/D");
        tree1->Branch("chamber_displacement_y_cm",&dy, "dy/D");
        tree1->Branch("chamber_displacement_z_cm",&dz, "dz/D");
        tree1->Branch("chamber_rotation_x_rad",&PhX, "PhX/D");
        tree1->Branch("chamber_rotation_z_rad",&PhZ, "PhZ/D");
    }
  
    edm::ESHandle<DTGeometry> dtGeometry;
    edm::ESHandle<CSCGeometry> cscGeometry;
    iSetup.get<MuonGeometryRecord>().get( dtGeometry );     
    iSetup.get<MuonGeometryRecord>().get( cscGeometry );
 
    AlignableMuon* theAlignableMuon = new AlignableMuon( &(*dtGeometry) , &(*cscGeometry) );
    AlignableNavigator* theAlignableNavigator = new AlignableNavigator( theAlignableMuon );

        const auto& theEndcaps = theAlignableMuon->CSCEndcaps();

        for (const auto& aliiter: theEndcaps) {
                addComponent(aliiter);
        }
    
    
    while (in.good())
    {
       bool missing1 = false;
       bool missing2 = false;
    
            in >> ID1 >> ID2 >> ID3 >> ID4 >> ID5 >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> a >> b;

        if (fabs(x1 - 1000.) < 1e5  &&  fabs(y1 - 1000.) < 1e5  &&  fabs(z1 - 1000.) < 1e5) {
           missing1 = true;
           x1 = x2;
           y1 = y2;
           z1 = z2;
        }

        if (fabs(x2 - 1000.) < 1e5  &&  fabs(y2 - 1000.) < 1e5  &&  fabs(z2 - 1000.) < 1e5) {
           missing2 = true;
           x2 = x1;
           y2 = y1;
           z2 = z1;
        }

        x1=x1/10.0;
        y1=y1/10.0;
        z1=z1/10.0;
        x2=x2/10.0;
        y2=y2/10.0;
        z2=z2/10.0;
   
        CSCDetId layerID(ID1, ID2, ID3, ID4, 1);
        
        // We cannot use chamber ID (when ID5=0), because AlignableNavigator gives the error (aliDet and aliDetUnit are undefined for chambers)
        
        
        Alignable* theAlignable1 = theAlignableNavigator->alignableFromDetId( layerID ); 
        Alignable* chamberAli=theAlignable1->mother();

        LocalVector LC1 = chamberAli->surface().toLocal(GlobalVector(x1,y1,z1));
        LocalVector LC2 = chamberAli->surface().toLocal(GlobalVector(x2,y2,z2));
  
        orient(LC1, LC2, a, b, T, dx, dy, dz, PhX, PhZ);    
   
        GlobalPoint PG1 = chamberAli->surface().toGlobal(LocalPoint(LC1.x(), LC1.y() + a, LC1.z() + b));
        chamberAli->surface().toGlobal(LocalPoint(LC2.x(), LC2.y() - a, LC2.z() + b));
 
    
        LocalVector lvector( dx, dy, dz);
        GlobalVector gvector = ( chamberAli->surface()).toGlobal( lvector );
        chamberAli->move( gvector );
  
        chamberAli->rotateAroundLocalX( PhX );
        chamberAli->rotateAroundLocalZ( PhZ );
        
        double dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix;
        errors(a, b, missing1, missing2, dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix);
        align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
        error(0,0) = dx_dx;
        error(1,1) = dy_dy;
        error(2,2) = dz_dz;
        error(3,3) = phix_phix;
        error(4,4) = m_missingErrorAngle;
        error(5,5) = phiz_phiz;
        error(1,3) = dy_phix;
        error(3,1) = dy_phix;  // just in case

        chamberAli->setSurvey( new SurveyDet(chamberAli->surface(), error) );
    
            if (m_verbose) {
                
            edm::LogInfo("SurveyInputCSCfromPins") << " survey information = " << chamberAli->survey() << "  \n";
            
            LocalPoint LP1n = chamberAli->surface().toLocal(PG1);
    
            LocalPoint hiP(LP1n.x(), LP1n.y() - a*T, LP1n.z() - b);

            h=hiP.mag();
            s1=LP1n.y() - a;

            if (h>maxErr) { maxErr=h; 
            
                ii=i;
            }
    
            edm::LogInfo("SurveyInputCSCfromPins") << "  \n" 
             << "i " << i++ << " " << ID1 << " " << ID2 << " " << ID3 << " " << ID4 << " " << ID5 << " error  " << h  << " \n"  
             <<" x1 " << x1 << " y1 " << y1 << " z1 " << z1 << " x2 " << x2 << " y2 " << y2 << " z2 " << z2  << " \n"
             << " error " << h <<  " S1 " << s1 << " \n"
             <<" dx " << dx << " dy " << dy << " dz " << dz << " PhX " << PhX << " PhZ " << PhZ  << " \n";

            tot+=h;
    
            tree1->Fill();
            }
    
    } 
 
    in.close();

    if (m_verbose) {
    
        file1->Write();
        edm::LogInfo("SurveyInputCSCfromPins") << " Total error  " << tot << "  Max Error " << maxErr << " N " << ii << "  \n";
    }

    file1->Close();
   
        for (const auto& aliiter: theEndcaps) {
                fillAllRecords(aliiter);
        }

    delete theAlignableMuon;
    delete theAlignableNavigator;

    edm::LogInfo("SurveyInputCSCfromPins") << "*************END INITIALIZATION***************" << "  \n";

    theFirstEvent = false;
  }
}


void SurveyInputCSCfromPins::fillAllRecords(Alignable *ali) {
    if (ali->survey() == nullptr) {

       AlignableCSCChamber *ali_AlignableCSCChamber = dynamic_cast<AlignableCSCChamber*>(ali);
       AlignableCSCStation *ali_AlignableCSCStation = dynamic_cast<AlignableCSCStation*>(ali);

       if (ali_AlignableCSCChamber != nullptr) {
          CSCDetId detid(ali->geomDetId());
          if (abs(detid.station()) == 1  &&  (detid.ring() == 1  ||  detid.ring() == 4)) {
         align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
         error(0,0) = m_missingErrorTranslation;
         error(1,1) = m_missingErrorTranslation;
         error(2,2) = m_missingErrorTranslation;
         error(3,3) = m_missingErrorAngle;
         error(4,4) = m_missingErrorAngle;
         error(5,5) = m_missingErrorAngle;

         ali->setSurvey(new SurveyDet(ali->surface(), error));
          }
          else {
         double a = 100.;
         double b = -9.4034;
         if      (abs(detid.station()) == 1  &&  detid.ring() == 2) a = -90.260;
         else if (abs(detid.station()) == 1  &&  detid.ring() == 3) a = -85.205;
         else if (abs(detid.station()) == 2  &&  detid.ring() == 1) a = -97.855;
         else if (abs(detid.station()) == 2  &&  detid.ring() == 2) a = -164.555;
         else if (abs(detid.station()) == 3  &&  detid.ring() == 1) a = -87.870;
         else if (abs(detid.station()) == 3  &&  detid.ring() == 2) a = -164.555;
         else if (abs(detid.station()) == 4  &&  detid.ring() == 1) a = -77.890;

        double dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix;
        errors(a, b, true, true, dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix);
        align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
        error(0,0) = dx_dx;
        error(1,1) = dy_dy;
        error(2,2) = dz_dz;
        error(3,3) = phix_phix;
        error(4,4) = m_missingErrorAngle;
        error(5,5) = phiz_phiz;
        error(1,3) = dy_phix;
        error(3,1) = dy_phix;  // just in case

        ali->setSurvey(new SurveyDet(ali->surface(), error));
          }
       }

       else if (ali_AlignableCSCStation != nullptr) {
          align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
          error(0,0) = m_stationErrorX;
          error(1,1) = m_stationErrorY;
          error(2,2) = m_stationErrorZ;
          error(3,3) = m_stationErrorPhiX;
          error(4,4) = m_stationErrorPhiY;
          error(5,5) = m_stationErrorPhiZ;

          ali->setSurvey(new SurveyDet(ali->surface(), error));
       }

       else {
          align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
          ali->setSurvey(new SurveyDet(ali->surface(), error*(1e-10)));
       }
    }

        for (const auto& iter: ali->components()) {
                fillAllRecords(iter);
        }
}


// Plug in to framework
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(SurveyInputCSCfromPins);