CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DTParametrizedDriftAlgo.cc
Go to the documentation of this file.
1 /*
2  * See header file for a description of this class.
3  *
4  * $Date: 2009/04/30 09:30:07 $
5  * $Revision: 1.3 $
6  * \author G. Cerminara - INFN Torino
7  */
8 
12 
15 
18 
22 
23 #include <iostream>
24 #include <iomanip>
25 
26 using namespace std;
27 using namespace edm;
28 
29 
31  DTRecHitBaseAlgo(config) {
32  interpolate = config.getParameter<bool>("interpolate");
33 
34  minTime = config.getParameter<double>("minTime"); // FIXME: Default was -3 ns
35 
36  maxTime = config.getParameter<double>("maxTime"); // FIXME: Default was 415 ns
37 
38  // Set verbose output
39  debug = config.getUntrackedParameter<bool>("debug","false");
40 
41  }
42 
43 
44 
46 
47 
48 
50  theSync->setES(setup);
51  // Access the magnetic field
52  ESHandle<MagneticField> magneticField;
53  setup.get<IdealMagneticFieldRecord>().get(magneticField);
54  magField = &*magneticField;
55 }
56 
57 
58 
59 // First Step
61  const DTDigi& digi,
62  LocalPoint& leftPoint,
63  LocalPoint& rightPoint,
64  LocalError& error) const {
65  // Get the layerId
66  DTLayerId layerId = layer->id();//FIXME: pass it instead of get it from layer
67  const DTWireId wireId(layerId, digi.wire());
68 
69  // Get Wire position
70  if(!layer->specificTopology().isWireValid(wireId.wire())) return false;
71  LocalPoint locWirePos(layer->specificTopology().wirePosition(wireId.wire()), 0, 0);
72  const GlobalPoint globWirePos = layer->toGlobal(locWirePos);
73 
74  // impact angle on DT chamber
75  // compute the impact angle using the centre of the wire,
76  // only for RZ superlayer (eta). In the other cases theta is very close to 0.
77  float angle=0.0;
78 
79  if (layerId.superlayer() == 2 ) {
80  //GlobalPoint lPos=layer->position();
81  GlobalVector lDir=(GlobalPoint()-globWirePos).unit();
82  LocalVector lDirLoc=layer->toLocal(lDir);
83 
84  angle = atan(lDirLoc.x()/-lDirLoc.z());
85  }
86 
87  return compute(layer, wireId, digi.time(), angle, globWirePos,
88  leftPoint, rightPoint, error, 1);
89 }
90 
91 
92 
93 // Second step
95  const DTRecHit1D& recHit1D,
96  const float& angle,
97  DTRecHit1D& newHit1D) const {
98  const DTWireId wireId = recHit1D.wireId();
99 
100  // Get Wire position
101  if(!layer->specificTopology().isWireValid(wireId.wire())) return false;
102  LocalPoint locWirePos(layer->specificTopology().wirePosition(wireId.wire()), 0, 0);
103  const GlobalPoint globWirePos = layer->toGlobal(locWirePos);
104 
105  return compute(layer, wireId, recHit1D.digiTime(), angle, globWirePos,
106  newHit1D, 2);
107 }
108 
109 
110 
111 // Third step
112 bool DTParametrizedDriftAlgo::compute(const DTLayer* layer,
113  const DTRecHit1D& recHit1D,
114  const float& angle,
115  const GlobalPoint& globPos,
116  DTRecHit1D& newHit1D) const {
117  return compute(layer, recHit1D.wireId(), recHit1D.digiTime(), angle,
118  globPos, newHit1D, 3);
119 }
120 
121 
122 
123 // Do the actual work.
124 bool DTParametrizedDriftAlgo::compute(const DTLayer* layer,
125  const DTWireId& wireId,
126  const float digiTime,
127  const float& angle,
128  const GlobalPoint& globPos,
129  LocalPoint& leftPoint,
130  LocalPoint& rightPoint,
131  LocalError& error,
132  int step) const {
133  // Subtract Offset according to DTTTrigBaseSync concrete instance
134  // chosen with the 'tZeroMode' parameter
135  float driftTime = digiTime - theSync->offset(layer, wireId, globPos);
136 
137  // check for out-of-time only at step 1
138  if (step==1 && (driftTime < minTime || driftTime > maxTime)) {
139  if (debug)
140  cout << "*** Drift time out of window for in-time hits "
141  << driftTime << endl;
142 
143  if(step == 1) { //FIXME: protection against failure at 2nd and 3rd steps, must be checked!!!
144  // Hits are interpreted as coming from out-of-time pile-up and recHit
145  // is ignored.
146  return false;
147  }
148  }
149 
150  // Small negative times interpreted as hits close to the wire.
151  if (driftTime<0.) driftTime=0;
152 
153  //----------------------------------------------------------------------
154  // Magnetic Field in layer frame
155  const LocalVector BLoc =
156  layer->toLocal(magField->inTesla(globPos));
157 
158  float By = BLoc.y();
159  float Bz = BLoc.z();
160 
161  //--------------------------------------------------------------------
162  // Calculate the drift distance and the resolution from the parametrization
163 
165  static DTTime2DriftParametrization par;
166 
167  bool parStatus =
168  par.computeDriftDistance_mean(driftTime, angle, By, Bz, interpolate, &DX);
169 
170  if (!parStatus) {
171  if (debug)
172  cout << "[DTParametrizedDriftAlgo]*** WARNING: call to parametrization failed" << endl;
173  return false;
174  }
175 
176  float sigma_l = DX.x_width_m;
177  float sigma_r = DX.x_width_p;
178  float drift = DX.x_drift;
179 
180 
181  float reso = 0;
182 
183  bool variableSigma = false;
184  // By defualt the errors are obtained from a fit of the residuals in the various
185  // stations/SL.
186  // The error returned by DTTime2DriftParametrization can not be used easily
187  // to determine the hit error due to the way the parametrization of the error
188  // is done (please contcat the authors for details).
189  // Anyhow changing variableSigma==true, an attempt is done to set a variable error
190  // according to the sigma calculated by DTTime2DriftParametrization.
191  // Additionally, contributions from uncertaionties in TOF and signal propagation
192  // corrections are added.
193  // Uncertainty in the determination of incident angle and hit position
194  // (ie B value) are NOT accounted.
195  // This is not the default since it does not give good results...
196 
197  int wheel = abs(wireId.wheel());
198  if (variableSigma) {
199  float sTDelays=0;
200  if (step==1) { // 1. step
201  reso = (sigma_l+sigma_r)/2.; // FIXME: theta/B are not yet known...
202  if (wireId.superlayer()==2) { // RZ
203  sTDelays = 2.92;
204  } else { // RPhi
205  if (wheel==0) {
206  sTDelays = 2.74;
207  } else if (wheel==1) {
208  sTDelays = 1.83;
209  } else if (wheel==2){
210  sTDelays = 1.25;
211  }
212  }
213  } else if (step==2) { // 2. step
214  reso = (sigma_l+sigma_r)/2.; // FIXME: B is not yet known...
215  if (wireId.superlayer()==2) { // RZ
216  sTDelays = 0.096;
217  } else { // RPhi
218  if (wheel==0) {
219  sTDelays = 0.022;
220  } else if (wheel==1) {
221  sTDelays = 0.079;
222  } else if (wheel==2){
223  sTDelays = 0.124;
224  }
225  }
226  } else if (step==3) { // 3. step
227  reso = (sigma_l+sigma_r)/2.;
228  if (wireId.superlayer()==2) { // RZ
229  sTDelays = 0.096;
230  } else { // RPhi
231  if (wheel==0) {
232  sTDelays = 0.022;
233  } else if (wheel==1) {
234  sTDelays = 0.079;
235  } else if (wheel==2){
236  sTDelays = 0.124;
237  }
238  }
239  }
240  float sXDelays = sTDelays*DX.v_drift/10.;
241  reso = sqrt(reso*reso + sXDelays*sXDelays);
242  } else { // Use a fixed sigma, from fit of residuals.
243  if (step==1) { // 1. step
244  if (wireId.superlayer()==2) {
245  if (wheel==0) {
246  reso = 0.0250;
247  } else if (wheel==1) {
248  reso = 0.0271;
249  } else if (wheel==2){
250  reso = 0.0308;
251  }
252  } else {
253  reso = 0.0237;
254  }
255  } else if (step==2) { // 2. step //FIXME
256  if (wireId.superlayer()==2) {
257  if (wheel==0) {
258  reso = 0.0250;
259  } else if (wheel==1) {
260  reso = 0.0271;
261  } else if (wheel==2){
262  reso = 0.0305;
263  }
264  } else {
265  reso = 0.0231;
266  }
267  } else if (step==3) { // 3. step
268  if (wireId.superlayer()==2) {
269  if (wheel==0) {
270  reso = 0.0196;
271  } else if (wheel==1) {
272  reso = 0.0210;
273  } else if (wheel==2){
274  reso = 0.0228;
275  }
276  } else {
277  reso = 0.0207;
278  }
279  }
280  }
281  //--------------------------------------------------------------------
282 
283  error = LocalError(reso*reso,0.,0.);
284 
285  // Get Wire position
286  if(!layer->specificTopology().isWireValid(wireId.wire())) return false;
287  LocalPoint locWirePos(layer->specificTopology().wirePosition(wireId.wire()), 0, 0);
288 
289  //Build the two possible points and the error on the position
290  leftPoint = LocalPoint(locWirePos.x()-drift,
291  locWirePos.y(),
292  locWirePos.z());
293  rightPoint = LocalPoint(locWirePos.x()+drift,
294  locWirePos.y(),
295  locWirePos.z());
296 
297  if (debug) {
298  int prevW = cout.width();
299  cout << setiosflags(ios::showpoint | ios::fixed) << setw(3)
300  << "[DTParametrizedDriftAlgo]: step " << step << endl
301  << " Global Position " << globPos << endl
302  << " Local Position " << layer->toLocal(globPos) << endl
303  // << " y along Wire " << wireCoord << endl
304  << " Digi time " << digiTime << endl
305  // << " dpropDelay " << propDelay << endl
306  // << " deltaTOF " << deltaTOF << endl
307  << " >Drif time " << driftTime << endl
308  << " >Angle " << angle * 180./M_PI << endl
309  << " <Drift distance " << drift << endl
310  << " <sigma_l " << sigma_l << endl
311  << " <sigma_r " << sigma_r << endl
312  << " reso " << reso << endl
313  << " leftPoint " << leftPoint << endl
314  << " rightPoint " << rightPoint << endl
315  << " error " << error
316  << resetiosflags(ios::showpoint | ios::fixed) << setw(prevW) << endl
317  << endl;
318  }
319 
320  return true;
321 }
322 
323 
324 // Interface to the method which does the actual work suited for 2nd and 3rd steps
325 bool DTParametrizedDriftAlgo::compute(const DTLayer* layer,
326  const DTWireId& wireId,
327  const float digiTime,
328  const float& angle,
329  const GlobalPoint& globPos,
330  DTRecHit1D& newHit1D,
331  int step) const {
332  LocalPoint leftPoint;
333  LocalPoint rightPoint;
335 
336  if(compute(layer, wireId, digiTime, angle, globPos, leftPoint, rightPoint, error, step)) {
337  // Set the position and the error of the rechit which is being updated
338  switch(newHit1D.lrSide()) {
339 
340  case DTEnums::Left:
341  {
342  // Keep the original y position of newHit1D: for step==3, it's the
343  // position along the wire. Needed for rotation alignment
344  LocalPoint leftPoint3D(leftPoint.x(), newHit1D.localPosition().y(), leftPoint.z());
345  newHit1D.setPositionAndError(leftPoint3D, error);
346  break;
347  }
348 
349  case DTEnums::Right:
350  {
351  // as above: 3d position
352  LocalPoint rightPoint3D(rightPoint.x(), newHit1D.localPosition().y(), rightPoint.z());
353  newHit1D.setPositionAndError(rightPoint3D, error);
354  break;
355  }
356 
357  default:
358  throw cms::Exception("InvalidDTCellSide") << "[DTParametrizedDriftAlgo] Compute at Step "
359  << step << ", Hit side "
360  << newHit1D.lrSide()
361  << " is invalid!" << endl;
362  return false;
363  }
364 
365  return true;
366  }else {
367  return false;
368  }
369 }
370 
371 
373 
374 
376 
377 
379 
380 
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
float wirePosition(int wireNumber) const
Returns the x position in the layer of a given wire number.
Definition: DTTopology.cc:88
virtual void setES(const edm::EventSetup &setup)
Pass the Event Setup to the algo at each event.
bool computeDriftDistance_mean(double time, double alpha, double by, double bz, short interpolate, drift_distance *dx) const
DTLayerId id() const
Return the DetId of this SL.
Definition: DTLayer.cc:48
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:47
LocalVector drift(const StripGeomDetUnit *, const MagneticField &, const SiStripLorentzAngle &)
Definition: ShallowTools.cc:39
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:63
#define abs(x)
Definition: mlp_lapack.h:159
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:62
virtual bool compute(const DTLayer *layer, const DTDigi &digi, LocalPoint &leftPoint, LocalPoint &rightPoint, LocalError &error) const
virtual ~DTParametrizedDriftAlgo()
Destructor.
int wire() const
Return wire number.
Definition: DTDigi.cc:69
const DTTopology & specificTopology() const
Definition: DTLayer.cc:44
double time() const
Get time in ns.
Definition: DTDigi.cc:65
string unit
Definition: csvLumiCalc.py:46
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
float digiTime() const
Return the time (ns) of the digi used to build the rechit.
Definition: DTRecHit1D.h:115
Definition: DTDigi.h:19
DTParametrizedDriftAlgo(const edm::ParameterSet &config)
Constructor.
virtual LocalPoint localPosition() const
Return the 3-dimensional local position.
Definition: DTRecHit1D.h:62
int wire() const
Return the wire number.
Definition: DTWireId.h:58
int superlayer() const
Return the superlayer number (deprecated method name)
#define M_PI
Definition: BFit3D.cc:3
const T & get() const
Definition: EventSetup.h:55
bool isWireValid(const int wireNumber) const
Definition: DTTopology.h:67
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
tuple cout
Definition: gather_cfg.py:121
void setPositionAndError(LocalPoint pos, LocalError err)
Set the local position and its error.
Definition: DTRecHit1D.h:102
#define debug
Definition: MEtoEDMFormat.h:34
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:47
T x() const
Definition: PV3DBase.h:62
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
Structure used to return output values.
DTEnums::DTCellSide lrSide() const
The side of the cell.
Definition: DTRecHit1D.h:84
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:109