CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SurveyInputCSCfromPins.cc
Go to the documentation of this file.
1 #include "FWCore/MessageLogger/interface/MessageLogger.h"
2 #include "Alignment/CommonAlignment/interface/SurveyDet.h"
3 #include "Geometry/DTGeometry/interface/DTGeometry.h"
4 #include "Geometry/Records/interface/MuonGeometryRecord.h"
5 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
6 #include <FWCore/Framework/interface/ESHandle.h>
7 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
8 #include "FWCore/Framework/interface/EventSetup.h"
9 #include "Geometry/Records/interface/IdealGeometryRecord.h"
10 #include "DataFormats/DetId/interface/DetId.h"
11 
12 #include "Alignment/MuonAlignment/interface/AlignableCSCChamber.h"
13 #include "Alignment/MuonAlignment/interface/AlignableCSCStation.h"
14 
15 #include <iostream>
16 #include <fstream>
17 #include "TFile.h"
18 #include "TTree.h"
19 #include "TRandom3.h"
20 
21 #include "SurveyInputCSCfromPins.h"
22 
23 #define SQR(x) ((x)*(x))
24 
25 SurveyInputCSCfromPins::SurveyInputCSCfromPins(const edm::ParameterSet& cfg)
26  : m_pinPositions(cfg.getParameter<std::string>("pinPositions"))
27  , m_rootFile(cfg.getParameter<std::string>("rootFile"))
28  , m_verbose(cfg.getParameter<bool>("verbose"))
29  , m_errorX(cfg.getParameter<double>("errorX"))
30  , m_errorY(cfg.getParameter<double>("errorY"))
31  , m_errorZ(cfg.getParameter<double>("errorZ"))
32  , m_missingErrorTranslation(cfg.getParameter<double>("missingErrorTranslation"))
33  , m_missingErrorAngle(cfg.getParameter<double>("missingErrorAngle"))
34  , m_stationErrorX(cfg.getParameter<double>("stationErrorX"))
35  , m_stationErrorY(cfg.getParameter<double>("stationErrorY"))
36  , m_stationErrorZ(cfg.getParameter<double>("stationErrorZ"))
37  , m_stationErrorPhiX(cfg.getParameter<double>("stationErrorPhiX"))
38  , m_stationErrorPhiY(cfg.getParameter<double>("stationErrorPhiY"))
39  , m_stationErrorPhiZ(cfg.getParameter<double>("stationErrorPhiZ"))
40 {}
41 
42 void SurveyInputCSCfromPins::orient(LocalVector LC1, LocalVector LC2, double a, double b, double &T, double &dx, double &dy, double &dz, double &PhX, double &PhZ) {
43  double cosPhX, sinPhX, cosPhZ, sinPhZ;
44 
45  LocalPoint LP1(LC1.x(), LC1.y() + a, LC1.z() + b);
46  LocalPoint LP2(LC2.x(), LC2.y() - a, LC2.z() + b);
47 
48  LocalPoint P((LP1.x() - LP2.x())/(2.*a), (LP1.y() - LP2.y())/(2.*a), (LP1.z() - LP2.z())/(2.*a));
49  LocalPoint Pp((LP1.x() + LP2.x())/(2.), (LP1.y() + LP2.y())/(2.), (LP1.z() + LP2.z())/(2.));
50 
51  T=P.mag();
52 
53  sinPhX=P.z()/T;
54  cosPhX=sqrt(1-SQR(sinPhX));
55  cosPhZ=P.y()/(T*cosPhX);
56  sinPhZ=-P.x()/(T*cosPhZ);
57 
58  PhX=atan2(sinPhX,cosPhX);
59 
60  PhZ=atan2(sinPhZ,cosPhZ);
61 
62  dx=Pp.x()-sinPhZ*sinPhX*b;
63  dy=Pp.y()+cosPhZ*sinPhX*b;
64  dz=Pp.z()-cosPhX*b;
65 }
66 
67 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) {
68  dx_dx = 0.;
69  dy_dy = 0.;
70  dz_dz = 0.;
71  phix_phix = 0.;
72  phiz_phiz = 0.;
73  dy_phix = 0.;
74 
75  const double trials = 10000.; // two significant digits
76 
77  for (int i = 0; i < trials; i++) {
78  LocalVector LC1, LC2;
79  if (missing1) {
80  LC1 = LocalVector(gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation));
81  }
82  else {
83  LC1 = LocalVector(gRandom->Gaus(0., m_errorX), gRandom->Gaus(0., m_errorY), gRandom->Gaus(0., m_errorZ));
84  }
85 
86  if (missing2) {
87  LC2 = LocalVector(gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation), gRandom->Gaus(0., m_missingErrorTranslation));
88  }
89  else {
90  LC2 = LocalVector(gRandom->Gaus(0., m_errorX), gRandom->Gaus(0., m_errorY), gRandom->Gaus(0., m_errorZ));
91  }
92 
93  double dx, dy, dz, PhX, PhZ, T;
94  orient(LC1, LC2, a, b, T, dx, dy, dz, PhX, PhZ);
95 
96  dx_dx += dx * dx;
97  dy_dy += dy * dy;
98  dz_dz += dz * dz;
99  phix_phix += PhX * PhX;
100  phiz_phiz += PhZ * PhZ;
101  dy_phix += dy * PhX; // the only non-zero off-diagonal element
102  }
103 
104  dx_dx /= trials;
105  dy_dy /= trials;
106  dz_dz /= trials;
107  phix_phix /= trials;
108  phiz_phiz /= trials;
109  dy_phix /= trials;
110 }
111 
112 void SurveyInputCSCfromPins::analyze(const edm::Event&, const edm::EventSetup& iSetup)
113 {
114  if (theFirstEvent) {
115 
116  edm::LogInfo("SurveyInputCSCfromPins") << "***************ENTERING INITIALIZATION******************" << " \n";
117 
118  std::ifstream in;
119  in.open(m_pinPositions.c_str());
120 
121  Double_t x1, y1, z1, x2, y2, z2, a, b, tot=0.0, maxErr=0.0, h, s1, dx, dy, dz, PhX, PhZ, T;
122 
123  int ID1, ID2, ID3, ID4, ID5, i=1, ii=0;
124 
125  TFile *file1 = new TFile(m_rootFile.c_str(),"recreate");
126  TTree *tree1 = new TTree("tree1","alignment pins");
127 
128  if (m_verbose) {
129 
130  tree1->Branch("displacement_x_pin1_cm", &x1, "x1/D");
131  tree1->Branch("displacement_y_pin1_cm", &y1, "y1/D");
132  tree1->Branch("displacement_z_pin1_cm", &z1, "z1/D");
133  tree1->Branch("displacement_x_pin2_cm", &x2, "x2/D");
134  tree1->Branch("displacement_y_pin2_cm", &y2, "y2/D");
135  tree1->Branch("displacement_z_pin2_cm", &z2, "z2/D");
136  tree1->Branch("error_vector_length_cm",&h, "h/D");
137  tree1->Branch("stretch_diff_cm",&s1, "s1/D");
138  tree1->Branch("stretch_factor",&T, "T/D");
139  tree1->Branch("chamber_displacement_x_cm",&dx, "dx/D");
140  tree1->Branch("chamber_displacement_y_cm",&dy, "dy/D");
141  tree1->Branch("chamber_displacement_z_cm",&dz, "dz/D");
142  tree1->Branch("chamber_rotation_x_rad",&PhX, "PhX/D");
143  tree1->Branch("chamber_rotation_z_rad",&PhZ, "PhZ/D");
144  }
145 
146  edm::ESHandle<DTGeometry> dtGeometry;
147  edm::ESHandle<CSCGeometry> cscGeometry;
148  iSetup.get<MuonGeometryRecord>().get( dtGeometry );
149  iSetup.get<MuonGeometryRecord>().get( cscGeometry );
150 
151  AlignableMuon* theAlignableMuon = new AlignableMuon( &(*dtGeometry) , &(*cscGeometry) );
152  AlignableNavigator* theAlignableNavigator = new AlignableNavigator( theAlignableMuon );
153 
154  std::vector<Alignable*> theEndcaps = theAlignableMuon->CSCEndcaps();
155 
156  for (std::vector<Alignable*>::const_iterator aliiter = theEndcaps.begin(); aliiter != theEndcaps.end(); ++aliiter) {
157 
158  addComponent(*aliiter);
159  }
160 
161 
162  while (in.good())
163  {
164  bool missing1 = false;
165  bool missing2 = false;
166 
167  in >> ID1 >> ID2 >> ID3 >> ID4 >> ID5 >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> a >> b;
168 
169  if (fabs(x1 - 1000.) < 1e5 && fabs(y1 - 1000.) < 1e5 && fabs(z1 - 1000.) < 1e5) {
170  missing1 = true;
171  x1 = x2;
172  y1 = y2;
173  z1 = z2;
174  }
175 
176  if (fabs(x2 - 1000.) < 1e5 && fabs(y2 - 1000.) < 1e5 && fabs(z2 - 1000.) < 1e5) {
177  missing2 = true;
178  x2 = x1;
179  y2 = y1;
180  z2 = z1;
181  }
182 
183  x1=x1/10.0;
184  y1=y1/10.0;
185  z1=z1/10.0;
186  x2=x2/10.0;
187  y2=y2/10.0;
188  z2=z2/10.0;
189 
190  CSCDetId layerID(ID1, ID2, ID3, ID4, 1);
191 
192  // We cannot use chamber ID (when ID5=0), because AlignableNavigator gives the error (aliDet and aliDetUnit are undefined for chambers)
193 
194 
195  Alignable* theAlignable1 = theAlignableNavigator->alignableFromDetId( layerID );
196  Alignable* chamberAli=theAlignable1->mother();
197 
198  LocalVector LC1 = chamberAli->surface().toLocal(GlobalVector(x1,y1,z1));
199  LocalVector LC2 = chamberAli->surface().toLocal(GlobalVector(x2,y2,z2));
200 
201  orient(LC1, LC2, a, b, T, dx, dy, dz, PhX, PhZ);
202 
203  GlobalPoint PG1 = chamberAli->surface().toGlobal(LocalPoint(LC1.x(), LC1.y() + a, LC1.z() + b));
204  chamberAli->surface().toGlobal(LocalPoint(LC2.x(), LC2.y() - a, LC2.z() + b));
205 
206 
207  LocalVector lvector( dx, dy, dz);
208  GlobalVector gvector = ( chamberAli->surface()).toGlobal( lvector );
209  chamberAli->move( gvector );
210 
211  chamberAli->rotateAroundLocalX( PhX );
212  chamberAli->rotateAroundLocalZ( PhZ );
213 
214  double dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix;
215  errors(a, b, missing1, missing2, dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix);
216  align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
217  error(0,0) = dx_dx;
218  error(1,1) = dy_dy;
219  error(2,2) = dz_dz;
220  error(3,3) = phix_phix;
221  error(4,4) = m_missingErrorAngle;
222  error(5,5) = phiz_phiz;
223  error(1,3) = dy_phix;
224  error(3,1) = dy_phix; // just in case
225 
226  chamberAli->setSurvey( new SurveyDet(chamberAli->surface(), error) );
227 
228  if (m_verbose) {
229 
230  edm::LogInfo("SurveyInputCSCfromPins") << " survey information = " << chamberAli->survey() << " \n";
231 
232  LocalPoint LP1n = chamberAli->surface().toLocal(PG1);
233 
234  LocalPoint hiP(LP1n.x(), LP1n.y() - a*T, LP1n.z() - b);
235 
236  h=hiP.mag();
237  s1=LP1n.y() - a;
238 
239  if (h>maxErr) { maxErr=h;
240 
241  ii=i;
242  }
243 
244  edm::LogInfo("SurveyInputCSCfromPins") << " \n"
245  << "i " << i++ << " " << ID1 << " " << ID2 << " " << ID3 << " " << ID4 << " " << ID5 << " error " << h << " \n"
246  <<" x1 " << x1 << " y1 " << y1 << " z1 " << z1 << " x2 " << x2 << " y2 " << y2 << " z2 " << z2 << " \n"
247  << " error " << h << " S1 " << s1 << " \n"
248  <<" dx " << dx << " dy " << dy << " dz " << dz << " PhX " << PhX << " PhZ " << PhZ << " \n";
249 
250  tot+=h;
251 
252  tree1->Fill();
253  }
254 
255  }
256 
257  in.close();
258 
259  if (m_verbose) {
260 
261  file1->Write();
262  edm::LogInfo("SurveyInputCSCfromPins") << " Total error " << tot << " Max Error " << maxErr << " N " << ii << " \n";
263  }
264 
265  file1->Close();
266 
267  for (std::vector<Alignable*>::const_iterator aliiter = theEndcaps.begin(); aliiter != theEndcaps.end(); ++aliiter) {
268 
269  fillAllRecords(*aliiter);
270  }
271 
272  delete theAlignableMuon;
273  delete theAlignableNavigator;
274 
275  edm::LogInfo("SurveyInputCSCfromPins") << "*************END INITIALIZATION***************" << " \n";
276 
277  theFirstEvent = false;
278  }
279 }
280 
281 
282 void SurveyInputCSCfromPins::fillAllRecords(Alignable *ali) {
283  if (ali->survey() == 0) {
284 
285  AlignableCSCChamber *ali_AlignableCSCChamber = dynamic_cast<AlignableCSCChamber*>(ali);
286  AlignableCSCStation *ali_AlignableCSCStation = dynamic_cast<AlignableCSCStation*>(ali);
287 
288  if (ali_AlignableCSCChamber != 0) {
289  CSCDetId detid(ali->geomDetId());
290  if (abs(detid.station()) == 1 && (detid.ring() == 1 || detid.ring() == 4)) {
291  align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
292  error(0,0) = m_missingErrorTranslation;
293  error(1,1) = m_missingErrorTranslation;
294  error(2,2) = m_missingErrorTranslation;
295  error(3,3) = m_missingErrorAngle;
296  error(4,4) = m_missingErrorAngle;
297  error(5,5) = m_missingErrorAngle;
298 
299  ali->setSurvey(new SurveyDet(ali->surface(), error));
300  }
301  else {
302  double a = 100.;
303  double b = -9.4034;
304  if (abs(detid.station()) == 1 && detid.ring() == 2) a = -90.260;
305  else if (abs(detid.station()) == 1 && detid.ring() == 3) a = -85.205;
306  else if (abs(detid.station()) == 2 && detid.ring() == 1) a = -97.855;
307  else if (abs(detid.station()) == 2 && detid.ring() == 2) a = -164.555;
308  else if (abs(detid.station()) == 3 && detid.ring() == 1) a = -87.870;
309  else if (abs(detid.station()) == 3 && detid.ring() == 2) a = -164.555;
310  else if (abs(detid.station()) == 4 && detid.ring() == 1) a = -77.890;
311 
312  double dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix;
313  errors(a, b, true, true, dx_dx, dy_dy, dz_dz, phix_phix, phiz_phiz, dy_phix);
314  align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
315  error(0,0) = dx_dx;
316  error(1,1) = dy_dy;
317  error(2,2) = dz_dz;
318  error(3,3) = phix_phix;
319  error(4,4) = m_missingErrorAngle;
320  error(5,5) = phiz_phiz;
321  error(1,3) = dy_phix;
322  error(3,1) = dy_phix; // just in case
323 
324  ali->setSurvey(new SurveyDet(ali->surface(), error));
325  }
326  }
327 
328  else if (ali_AlignableCSCStation != 0) {
329  align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
330  error(0,0) = m_stationErrorX;
331  error(1,1) = m_stationErrorY;
332  error(2,2) = m_stationErrorZ;
333  error(3,3) = m_stationErrorPhiX;
334  error(4,4) = m_stationErrorPhiY;
335  error(5,5) = m_stationErrorPhiZ;
336 
337  ali->setSurvey(new SurveyDet(ali->surface(), error));
338  }
339 
340  else {
341  align::ErrorMatrix error = ROOT::Math::SMatrixIdentity();
342  ali->setSurvey(new SurveyDet(ali->surface(), error*(1e-10)));
343  }
344  }
345 
346  std::vector<Alignable*> components = ali->components();
347  for (std::vector<Alignable*>::const_iterator iter = components.begin(); iter != components.end(); ++iter) {
348 
349  fillAllRecords(*iter);
350  }
351 }
352 
353 
354 // Plug in to framework
355 #include "FWCore/Framework/interface/MakerMacros.h"
356 DEFINE_FWK_MODULE(SurveyInputCSCfromPins);
i
int i
Definition: DBlmapReader.cc:9
SurveyInputCSCfromPins::orient
void orient(LocalVector LC1, LocalVector LC2, double a, double b, double &T, double &dx, double &dy, double &dz, double &PhX, double &PhZ)
Definition: SurveyInputCSCfromPins.cc:42
SurveyInputCSCfromPins::errors
void 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)
Definition: SurveyInputCSCfromPins.cc:67
LocalVector
Local3DVector LocalVector
Definition: LocalVector.h:12
Alignable::rotateAroundLocalZ
virtual void rotateAroundLocalZ(Scalar radians)
Rotation around local z-axis.
Definition: Alignable.cc:187
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
Alignable::survey
const SurveyDet * survey() const
Return survey info.
Definition: Alignable.h:205
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
abs
#define abs(x)
Definition: mlp_lapack.h:159
P
#define P
Definition: PFJetBenchmark.cc:33
SurveyInputCSCfromPins::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Read ideal tracker geometry from DB.
Definition: SurveyInputCSCfromPins.cc:112
Alignable::move
virtual void move(const GlobalVector &displacement)=0
Movement with respect to the global reference frame.
SurveyInputBase::addComponent
static void addComponent(Alignable *)
Definition: SurveyInputBase.cc:14
Alignable::components
virtual Alignables components() const =0
Return vector of all direct components.
AlignableMuon::CSCEndcaps
align::Alignables CSCEndcaps()
Definition: AlignableMuon.cc:377
AlignableSurface::toLocal
align::RotationType toLocal(const align::RotationType &) const
Return in local frame a rotation given in global frame.
SurveyInputCSCfromPins::fillAllRecords
void fillAllRecords(Alignable *ali)
Definition: SurveyInputCSCfromPins.cc:282
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
PV3DBase::z
T z() const
Definition: PV3DBase.h:63
SurveyInputCSCfromPins::SurveyInputCSCfromPins
SurveyInputCSCfromPins(const edm::ParameterSet &)
Definition: SurveyInputCSCfromPins.cc:25
Alignable::surface
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:126
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
Alignable::setSurvey
void setSurvey(const SurveyDet *)
Set survey info.
Definition: Alignable.cc:268
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
Alignable::rotateAroundLocalX
virtual void rotateAroundLocalX(Scalar radians)
Rotation around local x-axis.
Definition: Alignable.cc:135
AlignableSurface::toGlobal
align::GlobalPoints toGlobal(const align::LocalPoints &) const
Return in global coord given a set of local points.
Definition: AlignableSurface.cc:22
LocalPoint
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
align::ErrorMatrix
math::Error< 6 >::type ErrorMatrix
Definition: Definitions.h:39
T
long double T
Definition: Basic3DVectorLD.h:59
AlignableMuon
Constructor of the full muon geometry.
Definition: AlignableMuon.h:36
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
AlignableNavigator::alignableFromDetId
AlignableDetOrUnitPtr alignableFromDetId(const DetId &detid)
Returns AlignableDetOrUnitPtr corresponding to given DetId.
Definition: AlignableNavigator.cc:96
Alignable::mother
Alignable * mother() const
Return pointer to container alignable (if any)
Definition: Alignable.h:85
Alignable::geomDetId
const DetId & geomDetId() const
Definition: Alignable.h:177
SQR
#define SQR(x)
Definition: SurveyInputCSCfromPins.cc:23
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
AlignableCSCChamber
A muon CSC Chamber( an AlignableDet )
Definition: AlignableCSCChamber.h:25
Generated for CMSSW Reference Manual by doxygen 1.8.5