CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DTCombinatorialPatternReco4D.cc
Go to the documentation of this file.
1 
7 
10 
12 // For the 2D reco I use thei reconstructor!
15 
22 
24 
25 using namespace std;
26 using namespace edm;
27 
28 // TODO
29 // Throw an exception if a theta segment container is requested and in the event
30 // there isn't it. (Or launch a "lazy" reco on demand)
31 
33  DTRecSegment4DBaseAlgo(pset), theAlgoName("DTCombinatorialPatternReco4D"){
34 
35  // debug parameter
36  debug = pset.getUntrackedParameter<bool>("debug");
37 
38  //do you want the T0 correction?
39  applyT0corr = pset.getParameter<bool>("performT0SegCorrection");
40  computeT0corr = pset.getUntrackedParameter<bool>("computeT0Seg",true);
41 
42  // the updator
43  theUpdator = new DTSegmentUpdator(pset);
44 
45  // the input type.
46  // If true the instructions in setDTRecSegment2DContainer will be schipped and the
47  // theta segment will be recomputed from the 1D rechits
48  // If false the theta segment will be taken from the Event. Caveat: in this case the
49  // event must contain the 2D segments!
50  allDTRecHits = pset.getParameter<bool>("AllDTRecHits");
51 
52  // Get the concrete 2D-segments reconstruction algo from the factory
53  // For the 2D reco I use this reconstructor!
54  the2DAlgo = new DTCombinatorialPatternReco(pset.getParameter<ParameterSet>("Reco2DAlgoConfig"));
55  }
56 
58  delete the2DAlgo;
59  delete theUpdator;
60 }
61 
63  setup.get<MuonGeometryRecord>().get(theDTGeometry);
64  the2DAlgo->setES(setup);
65  theUpdator->setES(setup);
66 }
67 
69  // Set the chamber
70  theChamber = theDTGeometry->chamber(chId);
71 }
72 
74  theHitsFromPhi1.clear();
75  theHitsFromPhi2.clear();
76  theHitsFromTheta.clear();
77 
78  DTRecHitCollection::range rangeHitsFromPhi1 =
80  DTRecHitCollection::range rangeHitsFromPhi2 =
82 
83  vector<DTRecHit1DPair> hitsFromPhi1(rangeHitsFromPhi1.first,rangeHitsFromPhi1.second);
84  vector<DTRecHit1DPair> hitsFromPhi2(rangeHitsFromPhi2.first,rangeHitsFromPhi2.second);
85  if(debug)
86  cout<< "Number of DTRecHit1DPair in the SL 1 (Phi 1): " << hitsFromPhi1.size()<<endl
87  << "Number of DTRecHit1DPair in the SL 3 (Phi 2): " << hitsFromPhi2.size()<<endl;
88 
89  theHitsFromPhi1 = hitsFromPhi1;
90  theHitsFromPhi2 = hitsFromPhi2;
91 
92  if(allDTRecHits){
93  DTRecHitCollection::range rangeHitsFromTheta =
95 
96  vector<DTRecHit1DPair> hitsFromTheta(rangeHitsFromTheta.first,rangeHitsFromTheta.second);
97  if(debug)
98  cout<< "Number of DTRecHit1DPair in the SL 2 (Theta): " << hitsFromTheta.size()<<endl;
99  theHitsFromTheta = hitsFromTheta;
100  }
101 
102 }
103 
105  theSegments2DTheta.clear();
106 
107  if(!allDTRecHits){
108 
109  //Extract the DTRecSegment2DCollection range for the theta SL
110  DTRecSegment2DCollection::range rangeTheta =
111  all2DSegments->get(DTSuperLayerId(theChamber->id(),2));
112 
113  // Fill the DTRecSegment2D container for the theta SL
114  vector<DTSLRecSegment2D> segments2DTheta(rangeTheta.first,rangeTheta.second);
115 
116  if(debug)
117  cout << "Number of 2D-segments in the second SL (Theta): " << segments2DTheta.size() << endl;
118  theSegments2DTheta = segments2DTheta;
119  }
120 
121 }
122 
123 
126 
128 
129  if (debug){
130  cout << "Segments in " << theChamber->id() << endl;
131  cout << "Reconstructing of the Phi segments" << endl;
132  }
133 
134  vector<DTHitPairForFit*> pairPhiOwned;
135  vector<DTSegmentCand*> resultPhi = buildPhiSuperSegmentsCandidates(pairPhiOwned);
136 
137  if (debug) cout << "There are " << resultPhi.size() << " Phi cand" << endl;
138 
139  if (allDTRecHits){
140  // take the theta SL of this chamber
141  const DTSuperLayer* sl = theChamber->superLayer(2);
142  // sl points to 0 if the theta SL was not found
143  if(sl){
144  // reconstruct the theta segments
145  if(debug) cout<<"Reconstructing of the Theta segments"<<endl;
147  vector<DTSLRecSegment2D> segments2DTheta(thetaSegs.begin(),thetaSegs.end());
148  theSegments2DTheta = segments2DTheta;
149  }
150  }
151 
152  bool hasZed = false;
153 
154  // has this chamber the Z-superlayer?
155  if (theSegments2DTheta.size()){
156  hasZed = theSegments2DTheta.size() > 0;
157  if (debug) cout << "There are " << theSegments2DTheta.size() << " Theta cand" << endl;
158  } else {
159  if (debug) cout << "No Theta SL" << endl;
160  }
161 
162  // Now I want to build the concrete DTRecSegment4D.
163  if(debug) cout<<"Building of the concrete DTRecSegment4D"<<endl;
164  if (resultPhi.size()) {
165  for (vector<DTSegmentCand*>::const_iterator phi=resultPhi.begin();
166  phi!=resultPhi.end(); ++phi) {
167 
168  std::auto_ptr<DTChamberRecSegment2D> superPhi(**phi);
169 
170  theUpdator->update(superPhi.get());
171  if(debug) cout << "superPhi: " << *superPhi << endl;
172 
173  if (hasZed) {
174 
175  // Create all the 4D-segment combining the Z view with the Phi one
176  // loop over the Z segments
177  for(vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin();
178  zed != theSegments2DTheta.end(); ++zed){
179 
180  if(debug) cout << "Theta: " << *zed << endl;
181  // Important!!
182  DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
183 
184  // Put the theta segment poistion in its 3D place.
185  // note: (superPhi is in the CHAMBER local frame)
186  const DTSuperLayer* zSL = theChamber->superLayer(ZedSegSLId);
187 
188  // FIXME: should rather extrapolate for Y!
189  LocalPoint zPos(zed->localPosition().x(),
190  (zSL->toLocal(theChamber->toGlobal(superPhi->localPosition()))).y(),
191  0.);
192 
193  const LocalPoint posZInCh = theChamber->toLocal( zSL->toGlobal(zPos));
194  // FIXME: zed->localDirection() is in 2D. Should add the phi direction in the orthogonal plane as well!!
195  const LocalVector dirZInCh = theChamber->toLocal( zSL->toGlobal(zed->localDirection()));
196 
197  DTRecSegment4D* newSeg = new DTRecSegment4D(*superPhi,*zed,posZInCh,dirZInCh);
198 
199  if (debug) cout << "Created a 4D seg " << *newSeg << endl;
200 
202  theUpdator->update(newSeg);
203  if (debug) cout << " seg updated " << *newSeg << endl;
204 
205 
207  if(applyT0corr) theUpdator->update(newSeg,true);
208 
209  result.push_back(newSeg);
210  }
211  } else {
212  // Only phi
213 
214  DTRecSegment4D* newSeg = new DTRecSegment4D(*superPhi);
215 
216  if (debug) cout << "Created a 4D segment using only the 2D Phi segment " << *newSeg << endl;
217 
218  //update the segment with the t0 and possibly vdrift correction
220  if(applyT0corr) theUpdator->update(newSeg,true);
221 
222  result.push_back(newSeg);
223  }
224  }
225  } else {
226  // DTRecSegment4D from zed projection only (unlikely, not so useful, but...)
227  if (hasZed) {
228  for(vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin();
229  zed != theSegments2DTheta.end(); ++zed){
230  if(debug) cout << "Theta: " << *zed << endl;
231 
232  // Important!!
233  DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
234 
235  const LocalPoint posZInCh = theChamber->toLocal( theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localPosition() )) ;
236  const LocalVector dirZInCh = theChamber->toLocal( theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localDirection() )) ;
237 
238  DTRecSegment4D* newSeg = new DTRecSegment4D( *zed,posZInCh,dirZInCh);
239 
240  if (debug) cout << "Created a 4D segment using only the 2D Theta segment " <<
241  *newSeg << endl;
242 
244  if(applyT0corr) theUpdator->update(newSeg,true);
245 
246  result.push_back(newSeg);
247  }
248  }
249  }
250  // finally delete the candidates!
251  for (vector<DTSegmentCand*>::iterator phi=resultPhi.begin();
252  phi!=resultPhi.end(); ++phi) delete *phi;
253  for (vector<DTHitPairForFit*>::iterator phiPair = pairPhiOwned.begin();
254  phiPair!=pairPhiOwned.end(); ++phiPair) delete *phiPair;
255 
256  return result;
257 }
258 
259 
260 
261 vector<DTSegmentCand*> DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandidates(vector<DTHitPairForFit*> &pairPhiOwned){
262 
263  DTSuperLayerId slId;
264 
265  if(theHitsFromPhi1.size())
266  slId = theHitsFromPhi1.front().wireId().superlayerId();
267  else
268  if(theHitsFromPhi2.size())
269  slId = theHitsFromPhi2.front().wireId().superlayerId();
270  else{
271  if(debug) cout<<"DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandidates: "
272  <<"No Hits in the two Phi SL"<<endl;
273  return vector<DTSegmentCand*>();
274  }
275 
276  const DTSuperLayer *sl = theDTGeometry->superLayer(slId);
277 
278  vector<DTHitPairForFit*> pairPhi1 = the2DAlgo->initHits(sl,theHitsFromPhi1);
279  // same sl!! Since the fit will be in the sl phi 1!
280  vector<DTHitPairForFit*> pairPhi2 = the2DAlgo->initHits(sl,theHitsFromPhi2);
281  // copy the pairPhi2 in the pairPhi1 vector
282  copy(pairPhi2.begin(),pairPhi2.end(),back_inserter(pairPhi1));
283 
284  pairPhiOwned.swap(pairPhi1);
285  // Build the segment candidate
286  return the2DAlgo->buildSegments(sl,pairPhiOwned);
287 }
288 
291  // Get the zed projection
292  //if (!seg->hasZed()) return seg;
293  const DTSLRecSegment2D* zedSeg=seg->zSegment();
294  std::vector<DTRecHit1D> hits = zedSeg->specificRecHits();
295 
296  // pick up a hit "in the middle", where the single hit will be put.
297  int nHits=hits.size();
298  DTRecHit1D middle=hits[static_cast<int>(nHits/2.)];
299 
300  // Need to extrapolate pos to the middle layer z
301  LocalPoint posInSL = zedSeg->localPosition();
302  LocalVector dirInSL = zedSeg->localDirection();
303  LocalPoint posInMiddleLayer = posInSL+dirInSL*(-posInSL.z())/cos(dirInSL.theta());
304 
305  // create a hit with position and error as the Zed projection one's
306  auto_ptr<DTRecHit1D> hit(new DTRecHit1D( middle.wireId(),
307  middle.lrSide(),
308  middle.digiTime(),
309  posInMiddleLayer,
310  zedSeg->localPositionError()));
311 
312  std::vector<DTRecHit1D> newHits(1,*hit);
313 
314  // create a new zed segment with that single hits, but same dir and pos
315  LocalPoint pos(zedSeg->localPosition());
316  LocalVector dir(zedSeg->localDirection());
317  AlgebraicSymMatrix cov(zedSeg->covMatrix());
318  double chi2(zedSeg->chi2());
319  //cout << "zed " << *zedSeg << endl;
320  auto_ptr<DTSLRecSegment2D> newZed(new DTSLRecSegment2D(zedSeg->superLayerId(),
321  pos,
322  dir,
323  cov,
324  chi2,
325  newHits));
326  //cout << "newZed " << *newZed << endl;
327 
328  // create a 4d segment with the special zed
330  *newZed,
331  seg->localPosition(),
332  seg->localDirection());
333  // delete the input segment
334  delete seg;
335 
336  // return it
337  return result;
338 }
DTCombinatorialPatternReco * the2DAlgo
T getParameter(std::string const &) const
virtual LocalError localPositionError() const
local position error in SL frame
T getUntrackedParameter(std::string const &, T const &) const
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:50
std::vector< DTHitPairForFit * > initHits(const DTSuperLayer *sl, const std::vector< DTRecHit1DPair > &hits)
virtual void setES(const edm::EventSetup &setup)
void update(DTRecSegment4D *seg, const bool calcT0=false) const
recompute hits position and refit the segment4D
void calculateT0corr(DTRecSegment2D *seg) const
double zPos
const DTChamberRecSegment2D * phiSegment() const
The superPhi segment: 0 if no phi projection available.
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:47
virtual double chi2() const
the chi2 of the fit
static std::pair< DTLayerId, DTSuperLayerIdComparator > layersBySuperLayer(DTSuperLayerId slId)
Access by SL objects written into a RangeMap by layer.
DTSuperLayerId
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:62
std::vector< DTRecHit1DPair > theHitsFromPhi2
DTRecSegment4D * segmentSpecialZed(const DTRecSegment4D *seg)
Build a 4d segment with a zed component made by just one hits.
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
iterator begin()
Definition: OwnVector.h:227
virtual LocalVector localDirection() const
Local direction in Chamber frame.
void push_back(D *&d)
Definition: OwnVector.h:273
void setES(const edm::EventSetup &setup)
set the setup
DTChamberId id() const
Return the DTChamberId of this chamber.
Definition: DTChamber.cc:33
T z() const
Definition: PV3DBase.h:64
tuple result
Definition: query.py:137
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
float digiTime() const
Return the time (ns) of the digi used to build the rechit.
Definition: DTRecHit1D.h:113
DTCombinatorialPatternReco4D(const edm::ParameterSet &pset)
Constructor.
virtual LocalPoint localPosition() const
Local position in Chamber frame.
virtual void setES(const edm::EventSetup &setup)
std::vector< DTRecHit1D > specificRecHits() const
Access to specific components.
DTSuperLayerId superLayerId() const
The id of the superlayer on which reside the segment.
AlgebraicSymMatrix covMatrix() const
the Covariance Matrix
virtual LocalPoint localPosition() const
local position in SL frame
virtual edm::OwnVector< DTSLRecSegment2D > reconstruct(const DTSuperLayer *sl, const std::vector< DTRecHit1DPair > &hits)
this function is called in the producer
std::vector< DTSegmentCand * > buildPhiSuperSegmentsCandidates(std::vector< DTHitPairForFit * > &pairPhiOwned)
iterator end()
Definition: OwnVector.h:232
const DTSLRecSegment2D * zSegment() const
The Z segment: 0 if not zed projection available.
virtual ~DTCombinatorialPatternReco4D()
Destructor.
std::vector< DTRecHit1DPair > theHitsFromPhi1
std::vector< DTSLRecSegment2D > theSegments2DTheta
const T & get() const
Definition: EventSetup.h:55
virtual void setChamber(const DTChamberId &chId)
virtual void setDTRecHit1DContainer(edm::Handle< DTRecHitCollection > all1DHits)
std::vector< DTRecHit1DPair > theHitsFromTheta
edm::ESHandle< DTGeometry > theDTGeometry
std::vector< DTSegmentCand * > buildSegments(const DTSuperLayer *sl, const std::vector< DTHitPairForFit * > &hits)
virtual edm::OwnVector< DTRecSegment4D > reconstruct()
Operations.
virtual LocalVector localDirection() const
the local direction in SL frame
CLHEP::HepSymMatrix AlgebraicSymMatrix
virtual void setDTRecSegment2DContainer(edm::Handle< DTRecSegment2DCollection > all2DSegments)
tuple cout
Definition: gather_cfg.py:121
dbl *** dir
Definition: mlp_gen.cc:35
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
DTEnums::DTCellSide lrSide() const
The side of the cell.
Definition: DTRecHit1D.h:82
const DTSuperLayer * superLayer(DTSuperLayerId id) const
Return the superlayer corresponding to the given id.
Definition: DTChamber.cc:65
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:107
Definition: DDAxes.h:10