31 for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
32 if(dynamic_cast<const RPCChamber* >( *it ) != 0 ){
33 auto ch =
dynamic_cast<const RPCChamber*
>( *it );
34 std::vector< const RPCRoll*> roles = (ch->rolls());
35 for(std::vector<const RPCRoll*>::const_iterator
r = roles.begin();
r != roles.end(); ++
r){
39 int wheel=rpcId.
ring();
43 std::set<RPCDetId> myrolls;
45 myrolls.insert(rpcId);
54 if(sector1==13) sector1=4;
55 if(sector1==14) sector1=10;
57 if(sector2==13) sector2=4;
58 if(sector2==14) sector2=10;
60 int distance =
std::abs(sector1 - sector2);
61 if(distance>6) distance = 12-distance;
66 int distance =
std::abs(wheel1 - wheel2);
102 if(all4DSegments->size()>8){
103 if(debug)
std::cout<<
"Too many segments in this event we are not doing the extrapolation"<<std::endl;
121 std::map<DTChamberId,int> DTSegmentCounter;
124 for (segment = all4DSegments->begin();segment!=all4DSegments->end(); ++segment){
125 DTSegmentCounter[segment->chamberId()]++;
140 for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){
144 if(debug)
std::cout<<
"DT \t \t This Segment is in Chamber id: "<<DTId<<std::endl;
145 if(debug)
std::cout<<
"DT \t \t Number of segments in this DT = "<<DTSegmentCounter[DTId]<<std::endl;
146 if(debug)
std::cout<<
"DT \t \t Is the only one in this DT? and is not in the 4th Station?"<<std::endl;
148 if(DTSegmentCounter[DTId]!=1 || DTId.
station()==4){
149 if(debug)
std::cout<<
"DT \t \t More than one segment in this chamber, or we are in Station 4"<<std::endl;
153 int dtWheel = DTId.
wheel();
154 int dtStation = DTId.
station();
155 int dtSector = DTId.
sector();
157 LocalPoint segmentPosition= segment->localPosition();
158 LocalVector segmentDirection=segment->localDirection();
160 const GeomDet* gdet=dtGeo->idToDet(segment->geographicalId());
165 if(debug)
std::cout<<
"DT \t \t Is the segment 4D?"<<std::endl;
167 if(segment->dimension()!=4){
168 if(debug)
std::cout<<
"DT \t \t no"<<std::endl;
172 if(debug)
std::cout<<
"DT \t \t yes"<<std::endl;
173 if(debug)
std::cout<<
"DT \t \t DT Segment Dimension "<<segment->dimension()<<std::endl;
175 float Xo=segmentPosition.
x();
176 float Yo=segmentPosition.
y();
177 float Zo=segmentPosition.
z();
178 float dx=segmentDirection.
x();
179 float dy=segmentDirection.
y();
180 float dz=segmentDirection.
z();
182 if(debug)
std::cout<<
"Calling to Object Map class"<<std::endl;
184 if(debug)
std::cout<<
"Creating the DTIndex"<<std::endl;
186 if(debug)
std::cout<<
"Getting the Rolls for the given index"<<std::endl;
190 if(debug)
std::cout<<
"DT \t \t Number of rolls for this DT = "<<rollsForThisDT.size()<<std::endl;
192 assert(rollsForThisDT.size()>=1);
194 if(debug)
std::cout<<
"DT \t \t Loop over all the rolls asociated to this DT"<<std::endl;
195 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
196 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
203 if(debug)
std::cout<<
"DT \t \t \t RollName: "<<nameRoll<<std::endl;
204 if(debug)
std::cout<<
"DT \t \t \t Doing the extrapolation to this roll"<<std::endl;
205 if(debug)
std::cout<<
"DT \t \t \t DT Segment Direction in DTLocal "<<segmentDirection<<std::endl;
206 if(debug)
std::cout<<
"DT \t \t \t DT Segment Point in DTLocal "<<segmentPosition<<std::endl;
210 LocalPoint CenterRollinDTFrame = DTSurface.toLocal(CenterPointRollGlobal);
212 if(debug)
std::cout<<
"DT \t \t \t Center (0,0,0) Roll In DTLocal"<<CenterRollinDTFrame<<std::endl;
213 if(debug)
std::cout<<
"DT \t \t \t Center (0,0,0) of the Roll in Global"<<CenterPointRollGlobal<<std::endl;
215 float D=CenterRollinDTFrame.
z();
221 const RectangularStripTopology* top_=
dynamic_cast<const RectangularStripTopology*
> (&(rollasociated->
topology()));
223 if(debug)
std::cout<<
"DT \t \t \t xmin of this Roll "<<xmin<<
"cm"<<std::endl;
225 if(debug)
std::cout<<
"DT \t \t \t xmax of this Roll "<<xmax<<
"cm"<<std::endl;
226 float rsize = fabs( xmax.
x()-xmin.
x() );
227 if(debug)
std::cout<<
"DT \t \t \t Roll Size "<<rsize<<
"cm"<<std::endl;
228 float stripl = top_->stripLength();
230 float stripw = top_->pitch();
232 if(debug)
std::cout<<
"DT \t \t \t Strip Lenght "<<stripl<<
"cm"<<std::endl;
233 if(debug)
std::cout<<
"DT \t \t \t Strip Width "<<stripw<<
"cm"<<std::endl;
234 if(debug)
std::cout<<
"DT \t \t \t X Predicted in DTLocal= "<<X<<
"cm"<<std::endl;
235 if(debug)
std::cout<<
"DT \t \t \t Y Predicted in DTLocal= "<<Y<<
"cm"<<std::endl;
236 if(debug)
std::cout<<
"DT \t \t \t Z Predicted in DTLocal= "<<Z<<
"cm"<<std::endl;
238 float extrapolatedDistance =
sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
240 if(debug)
std::cout<<
"DT \t \t \t Is the distance of extrapolation less than MaxD? ="<<extrapolatedDistance<<
"cm"<<
"MaxD="<<
MaxD<<
"cm"<<std::endl;
242 if(extrapolatedDistance<=
MaxD){
243 if(debug)
std::cout<<
"DT \t \t \t yes"<<std::endl;
245 if(debug)
std::cout<<
"DT \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
246 LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
248 if(debug)
std::cout<<
"DT \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
249 if(debug)
std::cout<<
"DT \t \t \t Corner of the Roll = ("<<rsize*eyr<<
","<<stripl*eyr<<
")"<<std::endl;
250 if(debug)
std::cout<<
"DT \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.
x())<<
","
251 <<fabs(PointExtrapolatedRPCFrame.
y())<<
","<<fabs(PointExtrapolatedRPCFrame.
z())<<
")"<<std::endl;
252 if(debug)
std::cout<<
"DT \t \t \t Does the extrapolation go inside this roll?"<<std::endl;
254 if(fabs(PointExtrapolatedRPCFrame.
z()) < 1. &&
255 fabs(PointExtrapolatedRPCFrame.
x()) < rsize*eyr &&
256 fabs(PointExtrapolatedRPCFrame.
y()) < stripl*eyr){
257 if(debug)
std::cout<<
"DT \t \t \t \t yes"<<std::endl;
258 if(debug)
std::cout<<
"DT \t \t \t \t Creating the RecHit"<<std::endl;
259 RPCRecHit RPCPoint(rpcId,0,PointExtrapolatedRPCFrame);
260 if(debug)
std::cout<<
"DT \t \t \t \t Clearing the vector"<<std::endl;
262 if(debug)
std::cout<<
"DT \t \t \t \t Pushing back"<<std::endl;
264 if(debug)
std::cout<<
"DT \t \t \t \t Putting the vector"<<std::endl;
267 if(debug)
std::cout<<
"DT \t \t \t \t Filling container with "<<nameRoll
268 <<
" Point.x="<<PointExtrapolatedRPCFrame.
x()<<
" Point.y="<<PointExtrapolatedRPCFrame.
y()<<
" size="<<
RPCPointVector.
size()<<std::endl;
271 if(debug)
std::cout<<
"DT \t \t \t \t No the prediction is outside of this roll"<<std::endl;
274 if(debug)
std::cout<<
"DT \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
281 if(all4DSegments->size()>0){
282 if(debug)
std::cout<<
"MB4 \t \t Loop Over all4DSegments "<<all4DSegments->size()<<std::endl;
284 for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){
288 if(debug)
std::cout<<
"MB4 \t \t This Segment is in Chamber id: "<<DTId<<std::endl;
289 if(debug)
std::cout<<
"MB4 \t \t Number of segments in this DT = "<<DTSegmentCounter[DTId]<<std::endl;
290 if(debug)
std::cout<<
"MB4 \t \t \t Is the only one in this DT? and is in the Station 4?"<<std::endl;
292 if(DTSegmentCounter[DTId] == 1 && DTId.
station()==4){
294 if(debug)
std::cout<<
"MB4 \t \t \t yes"<<std::endl;
295 int dtWheel = DTId.
wheel();
296 int dtStation = DTId.
station();
297 int dtSector = DTId.
sector();
299 LocalPoint segmentPosition= segment->localPosition();
300 LocalVector segmentDirection=segment->localDirection();
302 if(debug)
std::cout<<
"MB4 \t \t \t \t The Segment in MB4 is 2D?"<<std::endl;
303 if(segment->dimension()==2){
304 if(debug)
std::cout<<
"MB4 \t \t \t \t yes"<<std::endl;
306 LocalPoint segmentPositionMB4=segmentPosition;
308 const BoundPlane& DTSurface4 = dtGeo->idToDet(DTId)->surface();
312 if(debug)
std::cout<<
"MB4 \t \t \t \t Loop on segments in =sector && MB3 && adjacent sectors && y dim=4"<<std::endl;
313 for(segMB3=all4DSegments->begin();segMB3!=all4DSegments->end();++segMB3){
317 if(debug)
std::cout<<
"MB4 \t \t \t \t Segment in Chamber ="<<dtid3<<std::endl;
323 && DTSegmentCounter[dtid3] == 1
324 && segMB3->dimension()==4){
329 const GeomDet* gdet3=dtGeo->idToDet(segMB3->geographicalId());
330 const BoundPlane & DTSurface3 = gdet3->surface();
332 LocalVector segmentDirectionMB3 = segMB3->localDirection();
333 GlobalPoint segmentPositionMB3inGlobal = DTSurface3.toGlobal(segMB3->localPosition());
334 GlobalPoint segmentPositionMB4inGlobal = DTSurface4.toGlobal(segmentPosition);
337 LocalVector segDirMB3inMB4Frame=DTSurface4.toLocal(DTSurface3.toGlobal(segmentDirectionMB3));
339 GlobalVector segDirMB4inGlobalFrame=DTSurface4.toGlobal(segmentDirectionMB4);
340 GlobalVector segDirMB3inGlobalFrame=DTSurface3.toGlobal(segmentDirectionMB3);
342 float dx=segDirMB4inGlobalFrame.
x();
343 float dy=segDirMB4inGlobalFrame.
y();
345 float dx3=segDirMB3inGlobalFrame.
x();
346 float dy3=segDirMB3inGlobalFrame.
y();
348 double cosAng=fabs(dx*dx3+dy*dy3/
sqrt((dx3*dx3+dy3*dy3)*(dx*dx+dy*dy)));
350 if(debug)
std::cout<<
"MB4 \t \t \t \t cosAng"<<cosAng<<
"Beetween "<<dtid3<<
" and "<<DTId<<std::endl;
353 std::cout<<
"MB4 \t \t \t \t dx="<<dx<<
" dy="<<dy<<std::endl;
354 std::cout<<
"MB4 \t \t \t \t dx3="<<dx3<<
" dy3="<<dy<<std::endl;
355 std::cout<<
"MB4 \t \t \t \t cosAng="<<cosAng<<std::endl;
358 float DistanceBetweenSegments = ((segmentPositionMB3inGlobal) - (segmentPositionMB4inGlobal)).
mag();
362 if(debug)
std::cout<<
"MB4 \t \t \t \t Distance between segments="<<DistanceBetweenSegments<<std::endl;
364 if(debug)
std::cout<<
"MB4 \t \t We found compatible Segments (similar direction and close enough) in "<<dtid3<<
" and "<<DTId<<std::endl;
373 if(debug)
std::cout<<
"Calling to Object Map class"<<std::endl;
375 if(debug)
std::cout<<
"Creating the DTIndex"<<std::endl;
377 if(debug)
std::cout<<
"Getting the Rolls for the given index"<<std::endl;
381 if(debug)
std::cout<<
"MB4 \t \t Number of rolls for this DT = "<<rollsForThisDT.size()<<std::endl;
383 assert(rollsForThisDT.size()>=1);
385 if(debug)
std::cout<<
"MB4 \t \t Loop over all the rolls asociated to this DT"<<std::endl;
386 for (std::set<RPCDetId>::iterator iteraRoll=rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
387 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
394 if(debug)
std::cout<<
"MB4 \t \t \t RollName: "<<nameRoll<<std::endl;
395 if(debug)
std::cout<<
"MB4 \t \t \t Doing the extrapolation to this roll"<<std::endl;
398 LocalPoint CenterRollinMB4Frame = DTSurface4.toLocal(CenterPointRollGlobal);
399 LocalPoint segmentPositionMB3inMB4Frame = DTSurface4.toLocal(segmentPositionMB3inGlobal);
401 LocalVector segmentDirectionMB3inMB4Frame = DTSurface4.toLocal(segDirMB3inGlobalFrame);
404 float Dxz=CenterRollinMB4Frame.
z();
405 float Xo4=segmentPositionMB4.
x();
406 float dxl=segmentDirectionMB4.
x();
407 float dzl=segmentDirectionMB4.
z();
409 float X=Xo4+dxl*Dxz/dzl;
413 float Yo34 = segmentPositionMB3inMB4Frame.
y();
414 float dy34 = segmentDirectionMB3inMB4Frame.
y();
415 float dz34 = segmentDirectionMB3inMB4Frame.
z();
416 float Dy=Dxz-(segmentPositionMB3inMB4Frame.
z());
418 if(debug)
std::cout<<
"MB4 \t \t \t The distance to extrapolate in Y from MB3 is "<<Dy<<
"cm"<<std::endl;
420 float Y=Yo34+dy34*Dy/dz34;
422 const RectangularStripTopology* top_
423 =
dynamic_cast<const RectangularStripTopology*
>(&(rollasociated->
topology()));
426 float rsize = fabs( xmax.
x()-xmin.
x() );
427 float stripl = top_->stripLength();
428 float stripw = top_->pitch();
431 if(debug)
std::cout<<
"MB4 \t \t \t Strip Lenght "<<stripl<<
"cm"<<std::endl;
432 if(debug)
std::cout<<
"MB4 \t \t \t Strip Width "<<stripw<<
"cm"<<std::endl;
434 if(debug)
std::cout<<
"MB4 \t \t \t X Predicted in MB4DTLocal= "<<X<<
"cm"<<std::endl;
435 if(debug)
std::cout<<
"MB4 \t \t \t Y Predicted in MB4DTLocal= "<<Y<<
"cm"<<std::endl;
436 if(debug)
std::cout<<
"MB4 \t \t \t Z Predicted in MB4DTLocal= "<<Z<<
"cm"<<std::endl;
438 float extrapolatedDistance =
sqrt((Y-Yo34)*(Y-Yo34)+Dy*Dy);
440 if(debug)
std::cout<<
"MB4 \t \t \t segmentPositionMB3inMB4Frame"<<segmentPositionMB3inMB4Frame<<std::endl;
441 if(debug)
std::cout<<
"MB4 \t \t \t segmentPositionMB4inMB4Frame"<<segmentPosition<<std::endl;
443 if(debug)
std::cout<<
"MB4 \t \t \t segmentDirMB3inMB4Frame"<<segDirMB3inMB4Frame<<std::endl;
444 if(debug)
std::cout<<
"MB4 \t \t \t segmentDirMB4inMB4Frame"<<segmentDirectionMB4<<std::endl;
446 if(debug)
std::cout<<
"MB4 \t \t \t CenterRB4PositioninMB4Frame"<<CenterRollinMB4Frame<<std::endl;
448 if(debug)
std::cout<<
"MB4 \t \t \t Is the extrapolation distance ="<<extrapolatedDistance<<
"less than "<<
MaxDrb4<<std::endl;
451 if(extrapolatedDistance<=
MaxDrb4){
452 if(debug)
std::cout<<
"MB4 \t \t \t yes"<<std::endl;
456 if(debug)
std::cout<<
"MB4 \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
458 LocalPoint PointExtrapolatedRPCFrame = RPCSurfaceRB4.toLocal(GlobalPointExtrapolated);
460 if(debug)
std::cout<<
"MB4 \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
461 if(debug)
std::cout<<
"MB4 \t \t \t Corner of the Roll = ("<<rsize*eyr<<
","<<stripl*eyr<<
")"<<std::endl;
462 if(debug)
std::cout<<
"MB4 \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.
x())<<
","
463 <<fabs(PointExtrapolatedRPCFrame.
y())<<
","<<fabs(PointExtrapolatedRPCFrame.
z())<<
")"<<std::endl;
465 if(debug)
std::cout<<
"MB4 \t \t \t Does the extrapolation go inside this roll?"<<std::endl;
467 if(fabs(PointExtrapolatedRPCFrame.
z()) < 5. &&
468 fabs(PointExtrapolatedRPCFrame.
x()) < rsize*eyr &&
469 fabs(PointExtrapolatedRPCFrame.
y()) < stripl*eyr){
470 if(debug)
std::cout<<
"MB4 \t \t \t \t yes"<<std::endl;
471 if(debug)
std::cout<<
"MB4 \t \t \t \t Creating the RecHit"<<std::endl;
472 RPCRecHit RPCPointMB4(rpcId,0,PointExtrapolatedRPCFrame);
473 if(debug)
std::cout<<
"MB4 \t \t \t \t Clearing the RPCPointVector"<<std::endl;
475 if(debug)
std::cout<<
"MB4 \t \t \t \t Pushing Back"<<std::endl;
477 if(debug)
std::cout<<
"MB4 \t \t \t \t Putting for "<<rpcId<<std::endl;
478 if(debug)
std::cout<<
"MB4 \t \t \t \t Filling container with "<<nameRoll
479 <<
" Point.x="<<PointExtrapolatedRPCFrame.
x()<<
" Point.y="<<PointExtrapolatedRPCFrame.
y()<<
" size="<<
RPCPointVector.
size()<<std::endl;
485 if(debug)
std::cout<<
"MB4 \t \t \t \t roll already extrapolated "<<rpcId<<std::endl;
490 if(debug)
std::cout<<
"MB4 \t \t \t \t No the prediction is outside of this roll"<<std::endl;
494 if(debug)
std::cout<<
"MB4 \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
498 if(debug)
std::cout<<
"MB4 \t \t \t \t I found segments in MB4 and MB3 adjacent wheel and/or sector but not compatibles, Diferent Directions"<<std::endl;
501 if(debug)
std::cout<<
"MB4 \t \t \t No the same station or same wheel or segment dim in mb3 not 4D"<<std::endl;
505 if(debug)
std::cout<<
"MB4 \t \t \t Is NOT a 2D Segment"<<std::endl;
508 if(debug)
std::cout<<
"MB4 \t \t \t \t There is not just one segment or is not in station 4"<<std::endl;
512 if(debug)
std::cout<<
"MB4 \t This event doesn't have 4D Segment"<<std::endl;
const double Z[kNumberCalorimeter]
double MaxDistanceBetweenSegments
const Topology & topology() const
std::set< RPCDetId > GetRolls(DTStationIndex dtstationindex)
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
int distsector(int sector1, int sector2)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
edm::OwnVector< RPCRecHit > RPCPointVector
static ObjectMap * mapInstance
const Plane & surface() const
The nominal surface of the GeomDet.
uint32_t rawId() const
get the raw id
DTSegtoRPC(edm::Handle< DTRecSegment4DCollection > all4DSegments, const edm::EventSetup &iSetup, const edm::Event &iEvent, bool debug, double eyr)
std::map< DTStationIndex, std::set< RPCDetId > > rollstoreDT
virtual std::string name()
Abs< T >::type abs(const T &t)
std::vector< uint32_t > extrapolatedRolls
static ObjectMap * GetInstance(const edm::EventSetup &iSetup)
edm::RangeMap< RPCDetId, edm::OwnVector< RPCRecHit, edm::ClonePolicy< RPCRecHit > >, edm::ClonePolicy< RPCRecHit > > RPCRecHitCollection
int sector() const
Sector id: the group of chambers at same phi (and increasing r)
int distwheel(int wheel1, int wheel2)
ObjectMap(const edm::EventSetup &iSetup)
int station() const
Return the station number.
int wheel() const
Return the wheel number.
RPCRecHitCollection * _ThePoints
DecomposeProduct< arg, typename Div::arg > D
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.