00001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00002 #include "Alignment/CommonAlignment/interface/SurveyDet.h"
00003 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00004 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00005 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
00006 #include <FWCore/Framework/interface/ESHandle.h>
00007 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
00008 #include "FWCore/Framework/interface/EventSetup.h"
00009 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00010 #include "DataFormats/DetId/interface/DetId.h"
00011
00012 #include "Alignment/MuonAlignment/interface/AlignableCSCChamber.h"
00013 #include "Alignment/MuonAlignment/interface/AlignableCSCStation.h"
00014
00015 #include <iostream>
00016 #include <fstream>
00017 #include "TFile.h"
00018 #include "TTree.h"
00019 #include "TRandom3.h"
00020
00021 #include "SurveyInputCSCfromPins.h"
00022
00023 #define SQR(x) ((x)*(x))
00024
00025 SurveyInputCSCfromPins::SurveyInputCSCfromPins(const edm::ParameterSet& cfg)
00026 : m_pinPositions(cfg.getParameter<std::string>("pinPositions"))
00027 , m_rootFile(cfg.getParameter<std::string>("rootFile"))
00028 , m_verbose(cfg.getParameter<bool>("verbose"))
00029 , m_errorX(cfg.getParameter<double>("errorX"))
00030 , m_errorY(cfg.getParameter<double>("errorY"))
00031 , m_errorZ(cfg.getParameter<double>("errorZ"))
00032 , m_missingErrorTranslation(cfg.getParameter<double>("missingErrorTranslation"))
00033 , m_missingErrorAngle(cfg.getParameter<double>("missingErrorAngle"))
00034 , m_stationErrorX(cfg.getParameter<double>("stationErrorX"))
00035 , m_stationErrorY(cfg.getParameter<double>("stationErrorY"))
00036 , m_stationErrorZ(cfg.getParameter<double>("stationErrorZ"))
00037 , m_stationErrorPhiX(cfg.getParameter<double>("stationErrorPhiX"))
00038 , m_stationErrorPhiY(cfg.getParameter<double>("stationErrorPhiY"))
00039 , m_stationErrorPhiZ(cfg.getParameter<double>("stationErrorPhiZ"))
00040 {}
00041
00042 void SurveyInputCSCfromPins::orient(LocalVector LC1, LocalVector LC2, double a, double b, double &T, double &dx, double &dy, double &dz, double &PhX, double &PhZ) {
00043 double cosPhX, sinPhX, cosPhZ, sinPhZ;
00044
00045 LocalPoint LP1(LC1.x(), LC1.y() + a, LC1.z() + b);
00046 LocalPoint LP2(LC2.x(), LC2.y() - a, LC2.z() + b);
00047
00048 LocalPoint P((LP1.x() - LP2.x())/(2.*a), (LP1.y() - LP2.y())/(2.*a), (LP1.z() - LP2.z())/(2.*a));
00049 LocalPoint Pp((LP1.x() + LP2.x())/(2.), (LP1.y() + LP2.y())/(2.), (LP1.z() + LP2.z())/(2.));
00050
00051 T=P.mag();
00052
00053 sinPhX=P.z()/T;
00054 cosPhX=sqrt(1-SQR(sinPhX));
00055 cosPhZ=P.y()/(T*cosPhX);
00056 sinPhZ=-P.x()/(T*cosPhZ);
00057
00058 PhX=atan2(sinPhX,cosPhX);
00059
00060 PhZ=atan2(sinPhZ,cosPhZ);
00061
00062 dx=Pp.x()-sinPhZ*sinPhX*b;
00063 dy=Pp.y()+cosPhZ*sinPhX*b;
00064 dz=Pp.z()-cosPhX*b;
00065 }
00066
00067 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) {
00068 dx_dx = 0.;
00069 dy_dy = 0.;
00070 dz_dz = 0.;
00071 phix_phix = 0.;
00072 phiz_phiz = 0.;
00073 dy_phix = 0.;
00074
00075 const double trials = 10000.;
00076
00077 for (int i = 0; i < trials; i++) {
00078 LocalVector LC1, LC2;
00079 if (missing1) {
00080 LC1 = LocalVector(gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation));
00081 }
00082 else {
00083 LC1 = LocalVector(gRandom->Gaus(0., m_errorX), gRandom->Gaus(0., m_errorY), gRandom->Gaus(0., m_errorZ));
00084 }
00085
00086 if (missing2) {
00087 LC2 = LocalVector(gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation));
00088 }
00089 else {
00090 LC2 = LocalVector(gRandom->Gaus(0., m_errorX), gRandom->Gaus(0., m_errorY), gRandom->Gaus(0., m_errorZ));
00091 }
00092
00093 double dx, dy, dz, PhX, PhZ, T;
00094 orient(LC1, LC2, a, b, T, dx, dy, dz, PhX, PhZ);
00095
00096 dx_dx += dx * dx;
00097 dy_dy += dy * dy;
00098 dz_dz += dz * dz;
00099 phix_phix += PhX * PhX;
00100 phiz_phiz += PhZ * PhZ;
00101 dy_phix += dy * PhX;
00102 }
00103
00104 dx_dx /= trials;
00105 dy_dy /= trials;
00106 dz_dz /= trials;
00107 phix_phix /= trials;
00108 phiz_phiz /= trials;
00109 dy_phix /= trials;
00110 }
00111
00112 void SurveyInputCSCfromPins::beginJob(const edm::EventSetup& iSetup)
00113 {
00114 edm::LogInfo("SurveyInputCSCfromPins") << "***************ENTERING BEGIN JOB******************" << " \n";
00115
00116 std::ifstream in;
00117 in.open(m_pinPositions.c_str());
00118
00119 Double_t x1, y1, z1, x2, y2, z2, a, b, tot=0.0, maxErr=0.0, h, s1, dx, dy, dz, PhX, PhZ, T;
00120
00121 int ID1, ID2, ID3, ID4, ID5, i=1, ii=0;
00122
00123 TFile *file1 = new TFile(m_rootFile.c_str(),"recreate");
00124 TTree *tree1 = new TTree("tree1","alignment pins");
00125
00126 if (m_verbose) {
00127
00128 tree1->Branch("displacement_x_pin1_cm", &x1, "x1/D");
00129 tree1->Branch("displacement_y_pin1_cm", &y1, "y1/D");
00130 tree1->Branch("displacement_z_pin1_cm", &z1, "z1/D");
00131 tree1->Branch("displacement_x_pin2_cm", &x2, "x2/D");
00132 tree1->Branch("displacement_y_pin2_cm", &y2, "y2/D");
00133 tree1->Branch("displacement_z_pin2_cm", &z2, "z2/D");
00134 tree1->Branch("error_vector_length_cm",&h, "h/D");
00135 tree1->Branch("stretch_diff_cm",&s1, "s1/D");
00136 tree1->Branch("stretch_factor",&T, "T/D");
00137 tree1->Branch("chamber_displacement_x_cm",&dx, "dx/D");
00138 tree1->Branch("chamber_displacement_y_cm",&dy, "dy/D");
00139 tree1->Branch("chamber_displacement_z_cm",&dz, "dz/D");
00140 tree1->Branch("chamber_rotation_x_rad",&PhX, "PhX/D");
00141 tree1->Branch("chamber_rotation_z_rad",&PhZ, "PhZ/D");
00142 }
00143
00144 edm::ESHandle<DTGeometry> dtGeometry;
00145 edm::ESHandle<CSCGeometry> cscGeometry;
00146 iSetup.get<MuonGeometryRecord>().get( dtGeometry );
00147 iSetup.get<MuonGeometryRecord>().get( cscGeometry );
00148
00149 AlignableMuon* theAlignableMuon = new AlignableMuon( &(*dtGeometry) , &(*cscGeometry) );
00150 AlignableNavigator* theAlignableNavigator = new AlignableNavigator( theAlignableMuon );
00151
00152
00153
00154 std::vector<Alignable*> theEndcaps = theAlignableMuon->CSCEndcaps();
00155
00156 for (std::vector<Alignable*>::const_iterator aliiter = theEndcaps.begin(); aliiter != theEndcaps.end(); ++aliiter) {
00157
00158 addComponent(*aliiter);
00159 }
00160
00161
00162 while (in.good())
00163 {
00164 bool missing1 = false;
00165 bool missing2 = false;
00166
00167 in >> ID1 >> ID2 >> ID3 >> ID4 >> ID5 >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> a >> b;
00168
00169 if (fabs(x1 - 1000.) < 1e5 && fabs(y1 - 1000.) < 1e5 && fabs(z1 - 1000.) < 1e5) {
00170 missing1 = true;
00171 x1 = x2;
00172 y1 = y2;
00173 z1 = z2;
00174 }
00175
00176 if (fabs(x2 - 1000.) < 1e5 && fabs(y2 - 1000.) < 1e5 && fabs(z2 - 1000.) < 1e5) {
00177 missing2 = true;
00178 x2 = x1;
00179 y2 = y1;
00180 z2 = z1;
00181 }
00182
00183 x1=x1/10.0;
00184 y1=y1/10.0;
00185 z1=z1/10.0;
00186 x2=x2/10.0;
00187 y2=y2/10.0;
00188 z2=z2/10.0;
00189
00190 CSCDetId layerID(ID1, ID2, ID3, ID4, 1);
00191
00192
00193
00194
00195 Alignable* theAlignable1 = theAlignableNavigator->alignableFromDetId( layerID );
00196 Alignable* chamberAli=theAlignable1->mother();
00197
00198 LocalVector LC1 = chamberAli->surface().toLocal(GlobalVector(x1,y1,z1));
00199 LocalVector LC2 = chamberAli->surface().toLocal(GlobalVector(x2,y2,z2));
00200
00201 orient(LC1, LC2, a, b, T, dx, dy, dz, PhX, PhZ);
00202
00203 GlobalPoint PG1 = chamberAli->surface().toGlobal(LocalPoint(LC1.x(), LC1.y() + a, LC1.z() + b));
00204 GlobalPoint PG2 = chamberAli->surface().toGlobal(LocalPoint(LC2.x(), LC2.y() - a, LC2.z() + b));
00205
00206
00207 LocalVector lvector( dx, dy, dz);
00208 GlobalVector gvector = ( chamberAli->surface()).toGlobal( lvector );
00209 chamberAli->move( gvector );
00210
00211 chamberAli->rotateAroundLocalX( PhX );
00212 chamberAli->rotateAroundLocalZ( PhZ );
00213
00214 double dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix;
00215 errors(a, b, missing1, missing2, dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix);
00216 align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
00217 error(0,0) = dx_dx;
00218 error(1,1) = dy_dy;
00219 error(2,2) = dz_dz;
00220 error(3,3) = phix_phix;
00221 error(4,4) = m_missingErrorAngle;
00222 error(5,5) = phiz_phiz;
00223 error(1,3) = dy_phix;
00224 error(3,1) = dy_phix;
00225
00226 chamberAli->setSurvey( new SurveyDet(chamberAli->surface(), error) );
00227
00228 if (m_verbose) {
00229
00230 edm::LogInfo("SurveyInputCSCfromPins") << " survey information = " << chamberAli->survey() << " \n";
00231
00232 LocalPoint LP1n = chamberAli->surface().toLocal(PG1);
00233
00234 LocalPoint hiP(LP1n.x(), LP1n.y() - a*T, LP1n.z() - b);
00235
00236 h=hiP.mag();
00237 s1=LP1n.y() - a;
00238
00239 if (h>maxErr) { maxErr=h;
00240
00241 ii=i;
00242 }
00243
00244 edm::LogInfo("SurveyInputCSCfromPins") << " \n"
00245 << "i " << i++ << " " << ID1 << " " << ID2 << " " << ID3 << " " << ID4 << " " << ID5 << " error " << h << " \n"
00246 <<" x1 " << x1 << " y1 " << y1 << " z1 " << z1 << " x2 " << x2 << " y2 " << y2 << " z2 " << z2 << " \n"
00247 << " error " << h << " S1 " << s1 << " \n"
00248 <<" dx " << dx << " dy " << dy << " dz " << dz << " PhX " << PhX << " PhZ " << PhZ << " \n";
00249
00250 tot+=h;
00251
00252 tree1->Fill();
00253 }
00254
00255 }
00256
00257 in.close();
00258
00259 if (m_verbose) {
00260
00261 file1->Write();
00262 edm::LogInfo("SurveyInputCSCfromPins") << " Total error " << tot << " Max Error " << maxErr << " N " << ii << " \n";
00263 }
00264
00265 file1->Close();
00266
00267 for (std::vector<Alignable*>::const_iterator aliiter = theEndcaps.begin(); aliiter != theEndcaps.end(); ++aliiter) {
00268
00269 fillAllRecords(*aliiter);
00270 }
00271
00272 edm::LogInfo("SurveyInputCSCfromPins") << "*************END BEGIN JOB***************" << " \n";
00273 }
00274
00275
00276 void SurveyInputCSCfromPins::fillAllRecords(Alignable *ali) {
00277 if (ali->survey() == 0) {
00278
00279 AlignableCSCChamber *ali_AlignableCSCChamber = dynamic_cast<AlignableCSCChamber*>(ali);
00280 AlignableCSCStation *ali_AlignableCSCStation = dynamic_cast<AlignableCSCStation*>(ali);
00281
00282 if (ali_AlignableCSCChamber != 0) {
00283 CSCDetId detid(ali->geomDetId());
00284 if (abs(detid.station()) == 1 && (detid.ring() == 1 || detid.ring() == 4)) {
00285 align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
00286 error(0,0) = m_missingErrorTranslation;
00287 error(1,1) = m_missingErrorTranslation;
00288 error(2,2) = m_missingErrorTranslation;
00289 error(3,3) = m_missingErrorAngle;
00290 error(4,4) = m_missingErrorAngle;
00291 error(5,5) = m_missingErrorAngle;
00292
00293 ali->setSurvey(new SurveyDet(ali->surface(), error));
00294 }
00295 else {
00296 double a = 100.;
00297 double b = -9.4034;
00298 if (abs(detid.station()) == 1 && detid.ring() == 2) a = -90.260;
00299 else if (abs(detid.station()) == 1 && detid.ring() == 3) a = -85.205;
00300 else if (abs(detid.station()) == 2 && detid.ring() == 1) a = -97.855;
00301 else if (abs(detid.station()) == 2 && detid.ring() == 2) a = -164.555;
00302 else if (abs(detid.station()) == 3 && detid.ring() == 1) a = -87.870;
00303 else if (abs(detid.station()) == 3 && detid.ring() == 2) a = -164.555;
00304 else if (abs(detid.station()) == 4 && detid.ring() == 1) a = -77.890;
00305
00306 double dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix;
00307 errors(a, b, true, true, dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix);
00308 align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
00309 error(0,0) = dx_dx;
00310 error(1,1) = dy_dy;
00311 error(2,2) = dz_dz;
00312 error(3,3) = phix_phix;
00313 error(4,4) = m_missingErrorAngle;
00314 error(5,5) = phiz_phiz;
00315 error(1,3) = dy_phix;
00316 error(3,1) = dy_phix;
00317
00318 ali->setSurvey(new SurveyDet(ali->surface(), error));
00319 }
00320 }
00321
00322 else if (ali_AlignableCSCStation != 0) {
00323 align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
00324 error(0,0) = m_stationErrorX;
00325 error(1,1) = m_stationErrorY;
00326 error(2,2) = m_stationErrorZ;
00327 error(3,3) = m_stationErrorPhiX;
00328 error(4,4) = m_stationErrorPhiY;
00329 error(5,5) = m_stationErrorPhiZ;
00330
00331 ali->setSurvey(new SurveyDet(ali->surface(), error));
00332 }
00333
00334 else {
00335 align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
00336 ali->setSurvey(new SurveyDet(ali->surface(), error*(1e-10)));
00337 }
00338 }
00339
00340 std::vector<Alignable*> components = ali->components();
00341 for (std::vector<Alignable*>::const_iterator iter = components.begin(); iter != components.end(); ++iter) {
00342
00343 fillAllRecords(*iter);
00344 }
00345 }
00346
00347
00348
00349 #include "FWCore/Framework/interface/MakerMacros.h"
00350 DEFINE_FWK_MODULE(SurveyInputCSCfromPins);
1.5.4