CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
RPCEfficiency.cc
Go to the documentation of this file.
1 /***************************************
2 Original Author: Camilo Carrillo
3 ****************************************/
4 
6 #include <sstream>
15 
16 //FW Core
18 
19 int distsector_tmp(int sector1,int sector2){
20 
21  if(sector1==13) sector1=4;
22  if(sector1==14) sector1=10;
23 
24  if(sector2==13) sector2=4;
25  if(sector2==14) sector2=10;
26 
27  int distance = abs(sector1 - sector2);
28  if(distance>6) distance = 12-distance;
29  return distance;
30 }
31 
32 
34  incldt=iConfig.getUntrackedParameter<bool>("incldt",true);
35  incldtMB4=iConfig.getUntrackedParameter<bool>("incldtMB4",true);
36  inclcsc=iConfig.getUntrackedParameter<bool>("inclcsc",true);
37  debug=iConfig.getUntrackedParameter<bool>("debug",false);
38 
39  rangestrips = iConfig.getUntrackedParameter<double>("rangestrips",4.);
40  rangestripsRB4=iConfig.getUntrackedParameter<double>("rangestripsRB4",4.);
41  dupli = iConfig.getUntrackedParameter<int>("DuplicationCorrection",2);
42  MinCosAng=iConfig.getUntrackedParameter<double>("MinCosAng",0.96);
43  MaxD=iConfig.getUntrackedParameter<double>("MaxD",80.);
44  MaxDrb4=iConfig.getUntrackedParameter<double>("MaxDrb4",150.);
45 
46 
47  cscSegments= consumes<CSCSegmentCollection>(iConfig.getParameter<edm::InputTag>("cscSegments"));
48  dt4DSegments = consumes<DTRecSegment4DCollection>(iConfig.getParameter<edm::InputTag>("dt4DSegments"));
49  RPCRecHitLabel_ = consumes<RPCRecHitCollection>(iConfig.getParameter<edm::InputTag>("RecHitLabel"));
50 
51  folderPath=iConfig.getUntrackedParameter<std::string>("folderPath","RPC/RPCEfficiency/");
52 
53  EffSaveRootFile = iConfig.getUntrackedParameter<bool>("EffSaveRootFile", false);
54  EffRootFileName = iConfig.getUntrackedParameter<std::string>("EffRootFileName", "RPCEfficiency.root");
55 
56 }
57 
58 
59 void RPCEfficiency::bookHistograms(DQMStore::IBooker & ibooker, edm::Run const & run, edm::EventSetup const & iSetup ) {
60 
61  std::string folder;
62  ibooker.cd();
64  statistics = ibooker.book1D("Statistics","All Statistics",33,0.5,33.5);
65 
66 
67  statistics->setBinLabel(1,"Events ",1);
68  statistics->setBinLabel(2,"Events with DT seg",1);
69 
70  std::stringstream sstr;
71  for( int i = 1; i<=15; i++ ){ //DT form bin 3 to bin 17
72  sstr.str("");
73  sstr<<i<<" DT seg";
74  statistics->setBinLabel(i+2, sstr.str() ,1);
75  }
76 
77  statistics->setBinLabel(18,"Events with CSC seg",1);
78  for( int i = 1; i<=15; i++ ){ //CSC form bin 19 to bin 33
79  sstr.str("");
80  sstr<<i<<" CSC seg";
81  statistics->setBinLabel(i+18, sstr.str(),1);
82  }
83 
84 
85  LogDebug("rpcefficiency")<<"booking Global histograms with "<<folderPath;
86 
87  folder = folderPath+"MuonSegEff/Residuals/Barrel";
88  ibooker.setCurrentFolder(folder);
89 
90  //Barrel
91  std::stringstream histoName, histoTitle;
92 
93  for (int layer = 1 ; layer<= 6 ;layer++){
94  histoName.str("");
95  histoTitle.str("");
96  histoName<<"GlobalResidualsClu1La"<<layer;
97  histoTitle<<"RPC Residuals Layer "<<layer<<" Cluster Size 1";
98  hGlobalResClu1La[layer-1] = ibooker.book1D(histoName.str(), histoTitle.str(),101,-10.,10.);
99 
100  histoName.str("");
101  histoTitle.str("");
102  histoName<<"GlobalResidualsClu2La"<<layer;
103  histoTitle<<"RPC Residuals Layer "<<layer<<" Cluster Size 2";
104  hGlobalResClu2La[layer-1] = ibooker.book1D(histoName.str(), histoTitle.str(),101,-10.,10.);
105 
106  histoName.str("");
107  histoTitle.str("");
108  histoName<<"GlobalResidualsClu3La"<<layer;
109  histoTitle<<"RPC Residuals Layer "<<layer<<" Cluster Size 3";
110  hGlobalResClu3La[layer-1] = ibooker.book1D(histoName.str(), histoTitle.str(),101,-10.,10.);
111 
112  }
113 
114  LogDebug("rpcefficiency")<<"Booking Residuals for EndCap";
115  folder = folderPath+"MuonSegEff/Residuals/Endcap";
116  ibooker.setCurrentFolder(folder);
117 
118  //Endcap
119  hGlobalResClu1R3C = ibooker.book1D("GlobalResidualsClu1R3C","RPC Residuals Ring 3 Roll C Cluster Size 1",101,-10.,10.);
120  hGlobalResClu1R3B = ibooker.book1D("GlobalResidualsClu1R3B","RPC Residuals Ring 3 Roll B Cluster Size 1",101,-10.,10.);
121  hGlobalResClu1R3A = ibooker.book1D("GlobalResidualsClu1R3A","RPC Residuals Ring 3 Roll A Cluster Size 1",101,-10.,10.);
122  hGlobalResClu1R2C = ibooker.book1D("GlobalResidualsClu1R2C","RPC Residuals Ring 2 Roll C Cluster Size 1",101,-10.,10.);
123  hGlobalResClu1R2B = ibooker.book1D("GlobalResidualsClu1R2B","RPC Residuals Ring 2 Roll B Cluster Size 1",101,-10.,10.);
124  hGlobalResClu1R2A = ibooker.book1D("GlobalResidualsClu1R2A","RPC Residuals Ring 2 Roll A Cluster Size 1",101,-10.,10.);
125 
126  hGlobalResClu2R3C = ibooker.book1D("GlobalResidualsClu2R3C","RPC Residuals Ring 3 Roll C Cluster Size 2",101,-10.,10.);
127  hGlobalResClu2R3B = ibooker.book1D("GlobalResidualsClu2R3B","RPC Residuals Ring 3 Roll B Cluster Size 2",101,-10.,10.);
128  hGlobalResClu2R3A = ibooker.book1D("GlobalResidualsClu2R3A","RPC Residuals Ring 3 Roll A Cluster Size 2",101,-10.,10.);
129  hGlobalResClu2R2C = ibooker.book1D("GlobalResidualsClu2R2C","RPC Residuals Ring 2 Roll C Cluster Size 2",101,-10.,10.);
130  hGlobalResClu2R2B = ibooker.book1D("GlobalResidualsClu2R2B","RPC Residuals Ring 2 Roll B Cluster Size 2",101,-10.,10.);
131  hGlobalResClu2R2A = ibooker.book1D("GlobalResidualsClu2R2A","RPC Residuals Ring 2 Roll A Cluster Size 2",101,-10.,10.);
132 
133  hGlobalResClu3R3C = ibooker.book1D("GlobalResidualsClu3R3C","RPC Residuals Ring 3 Roll C Cluster Size 3",101,-10.,10.);
134  hGlobalResClu3R3B = ibooker.book1D("GlobalResidualsClu3R3B","RPC Residuals Ring 3 Roll B Cluster Size 3",101,-10.,10.);
135  hGlobalResClu3R3A = ibooker.book1D("GlobalResidualsClu3R3A","RPC Residuals Ring 3 Roll A Cluster Size 3",101,-10.,10.);
136  hGlobalResClu3R2C = ibooker.book1D("GlobalResidualsClu3R2C","RPC Residuals Ring 2 Roll C Cluster Size 3",101,-10.,10.);
137  hGlobalResClu3R2B = ibooker.book1D("GlobalResidualsClu3R2B","RPC Residuals Ring 2 Roll B Cluster Size 3",101,-10.,10.);
138  hGlobalResClu3R2A = ibooker.book1D("GlobalResidualsClu3R2A","RPC Residuals Ring 2 Roll A Cluster Size 3",101,-10.,10.);
139 
140 
142  iSetup.get<MuonGeometryRecord>().get(rpcGeo);
143 
144 
145  for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){ //Loop on all detector units
146 
147  if(dynamic_cast< const RPCChamber* >( *it ) != 0 ){ // check if chamber exists
148  const RPCChamber* ch = dynamic_cast< const RPCChamber* >( *it );
149  std::vector< const RPCRoll*> roles = (ch->rolls()); //get all rolls in a chambers
150 
151  for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){ //Loop on all rolls
152 
153  RPCDetId rpcId = (*r)->id();
154  int region=rpcId.region();
155 
156  LogDebug("rpcefficiency")<<"Booking for "<<rpcId.rawId();
157 
158  bookDetUnitSeg(ibooker, rpcId,(*r)->nstrips(), folderPath+"MuonSegEff/", meCollection[rpcId.rawId()] ); // book histograms
159 
160  if(region==0&&(incldt||incldtMB4)){
161  //LogDebug("rpcefficiency")<<"--Filling the dtstore"<<rpcId;
162  int wheel=rpcId.ring();
163  int sector=rpcId.sector();
164  int station=rpcId.station();
165  DTStationIndex ind(region,wheel,sector,station);
166  std::set<RPCDetId> myrolls;
167  if (rollstoreDT.find(ind)!=rollstoreDT.end()) {myrolls=rollstoreDT[ind];}
168  myrolls.insert(rpcId);
169  rollstoreDT[ind]=myrolls;
170 
171  }else if(region!=0 && inclcsc){
172  int station=rpcId.station();
173  int ring=rpcId.ring();
174  int cscring=ring;
175  int cscstation=station;
176  RPCGeomServ rpcsrv(rpcId);
177  int rpcsegment = rpcsrv.segment();
178  int cscchamber = rpcsegment;
179 
180  if((station==2||station==3||station==4)&&ring==3){
181  cscring = 2;
182  }
183 
184  CSCStationIndex ind(region,cscstation,cscring,cscchamber);
185  std::set<RPCDetId> myrolls;
186  if (rollstoreCSC.find(ind)!=rollstoreCSC.end()){ myrolls=rollstoreCSC[ind];}
187  myrolls.insert(rpcId);
188  rollstoreCSC[ind]=myrolls;
189 
190  if(rpcId.ring()==2 || rpcId.ring()==3){
191 
192  cscchamber = rpcsegment+1;
193  if(cscchamber==37){cscchamber=1;}
194  CSCStationIndex ind2(region,cscstation,cscring,cscchamber);
195  std::set<RPCDetId> myrolls2;
196  if (rollstoreCSC.find(ind2)!=rollstoreCSC.end()){myrolls2=rollstoreCSC[ind2];}
197  myrolls2.insert(rpcId);
198  rollstoreCSC[ind2]=myrolls2;
199 
200  cscchamber = rpcsegment-1;
201  if(cscchamber==0){cscchamber=36;}
202  CSCStationIndex ind3(region,cscstation,cscring,cscchamber);
203  std::set<RPCDetId> myrolls3;
204  if (rollstoreCSC.find(ind3)!=rollstoreCSC.end()) {myrolls3=rollstoreCSC[ind3]; }
205  myrolls3.insert(rpcId);
206  rollstoreCSC[ind3]=myrolls3;
207 
208  }
209 
210  }
211 
212  }
213  }
214  }
215 }//beginRun
216 
217 
219 
221 
222 
226 
227  iSetup.get<MuonGeometryRecord>().get(rpcGeo);
228  iSetup.get<MuonGeometryRecord>().get(dtGeo);
229  iSetup.get<MuonGeometryRecord>().get(cscGeo);
230 
231  statistics->Fill(1);
232 
233  std::stringstream meIdRPC, meIdDT, meIdCSC;
234 
235  LogDebug("rpcefficiency") <<"\t Getting the RPC RecHits";
237  iEvent.getByToken(RPCRecHitLabel_,rpcHits);
238 
239  if(!rpcHits.isValid()) return;
240 
241  if(incldt){
242  LogDebug("rpcefficiency")<<"\t Getting the DT Segments";
244 
245  iEvent.getByToken(dt4DSegments, all4DSegments);
246 
247  if(all4DSegments.isValid()){
248 
249  if(all4DSegments->size()>0){
250 
251  if(all4DSegments->size()<=16) {statistics->Fill(2);}
252 
253  LogDebug("rpcefficiency")<<"\t Number of DT Segments in this event = "<<all4DSegments->size();
254 
255  std::map<DTChamberId,int> DTSegmentCounter;
257 
258  for (segment = all4DSegments->begin();segment!=all4DSegments->end(); ++segment){
259  DTSegmentCounter[segment->chamberId()]++;
260  }
261 
262  statistics->Fill(all4DSegments->size()+2);
263 
264  LogDebug("rpcefficiency")<<"\t Loop over all the 4D Segments";
265  for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){
266 
267  DTChamberId DTId = segment->chamberId();
268 
269 
270  if(DTSegmentCounter[DTId]==1 && DTId.station()!=4){
271 
272  int dtWheel = DTId.wheel();
273  int dtStation = DTId.station();
274  int dtSector = DTId.sector();
275 
276  LocalPoint segmentPosition= segment->localPosition();
277  LocalVector segmentDirection=segment->localDirection();
278 
279  const GeomDet* gdet=dtGeo->idToDet(segment->geographicalId());
280  const BoundPlane & DTSurface = gdet->surface();
281 
282  //check if the dimension of the segment is 4
283 
284  if(segment->dimension()==4){
285 
286  float Xo=segmentPosition.x();
287  float Yo=segmentPosition.y();
288  float Zo=segmentPosition.z();
289  float dx=segmentDirection.x();
290  float dy=segmentDirection.y();
291  float dz=segmentDirection.z();
292 
293  std::set<RPCDetId> rollsForThisDT = rollstoreDT[DTStationIndex(0,dtWheel,dtSector,dtStation)];
294 
295  LogDebug("rpcefficiency")<<"DT \t \t Loop over all the rolls asociated to this DT";
296  for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
297  const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
298  RPCDetId rpcId = rollasociated->id();
299  const BoundPlane & RPCSurface = rollasociated->surface();
300 
301  GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
302 
303  LocalPoint CenterRollinDTFrame = DTSurface.toLocal(CenterPointRollGlobal);
304 
305  float D=CenterRollinDTFrame.z();
306 
307  float X=Xo+dx*D/dz;
308  float Y=Yo+dy*D/dz;
309  float Z=D;
310 
311  const RectangularStripTopology* top_= dynamic_cast<const RectangularStripTopology*> (&(rollasociated->topology()));
312  LocalPoint xmin = top_->localPosition(0.);
313  LogDebug("rpcefficiency")<<"DT \t \t \t xmin of this Roll "<<xmin<<"cm";
314  LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
315  LogDebug("rpcefficiency")<<"DT \t \t \t xmax of this Roll "<<xmax<<"cm";
316  float rsize = fabs( xmax.x()-xmin.x() );
317  LogDebug("rpcefficiency")<<"DT \t \t \t Roll Size "<<rsize<<"cm";
318  float stripl = top_->stripLength();
319  float stripw = top_->pitch();
320 
321  float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
322 
323  if(extrapolatedDistance<=MaxD){
324 
325  GlobalPoint GlobalPointExtrapolated = DTSurface.toGlobal(LocalPoint(X,Y,Z));
326  LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
327 
328  if(fabs(PointExtrapolatedRPCFrame.z()) < 10. &&
329  fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
330  fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
331 
332  RPCDetId rollId = rollasociated->id();
333  RPCGeomServ rpcsrv(rollId);
334  std::string nameRoll = rpcsrv.name();
335  LogDebug("rpcefficiency")<<"DT \t \t \t \t The RPCName is "<<nameRoll;
336  const float stripPredicted =
337  rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
338 
339  LogDebug("rpcefficiency")<<"DT \t \t \t \t Candidate (from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
340  //---- HISTOGRAM STRIP PREDICTED FROM DT ----
341 
342  std::map<std::string, MonitorElement*> meMap=meCollection[rpcId.rawId()];
343  meIdDT.str("");
344  meIdDT<<"ExpectedOccupancyFromDT_"<<rollId.rawId();
345  meMap[meIdDT.str()]->Fill(stripPredicted);
346  //-----------------------------------------------------
347 
348 
349  //-------RecHitPart Just For Residual--------
350  int countRecHits = 0;
351  int cluSize = 0;
352  float minres = 3000.;
353 
354  typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
355  rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
357 
358  for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
359  countRecHits++;
360 
361  LocalPoint recHitPos=recHit->localPosition();
362  float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
363  LogDebug("rpcefficiency")<<"DT \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction.";
364  if(fabs(res)<fabs(minres)){
365  minres=res;
366  cluSize = recHit->clusterSize();
367  LogDebug("rpcefficiency")<<"DT \t \t \t \t \t \t New Min Res "<<res<<"cm.";
368  }
369  }
370 
371  if(countRecHits==0){
372  LogDebug("rpcefficiency") <<"DT \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT";
373  }else{
374  assert(minres!=3000);
375 
376  if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
377  LogDebug("rpcefficiency")<<"DT \t \t \t \t \t \t True!";
378 
379  if(rollId.station()==1&&rollId.layer()==1) {
380  if(cluSize==1*dupli) {hGlobalResClu1La[0]->Fill(minres);}
381  else if(cluSize==2*dupli){ hGlobalResClu2La[0]->Fill(minres);}
382  else if(cluSize==3*dupli){ hGlobalResClu3La[0]->Fill(minres);}
383  }else if(rollId.station()==1&&rollId.layer()==2){
384  if(cluSize==1*dupli) {hGlobalResClu1La[1]->Fill(minres);}
385  else if(cluSize==2*dupli){ hGlobalResClu2La[1]->Fill(minres);}
386  else if(cluSize==3*dupli){ hGlobalResClu3La[1]->Fill(minres);}
387  }else if(rollId.station()==2&&rollId.layer()==1){
388  if(cluSize==1*dupli) {hGlobalResClu1La[2]->Fill(minres);}
389  else if(cluSize==2*dupli){ hGlobalResClu2La[2]->Fill(minres);}
390  else if(cluSize==3*dupli){ hGlobalResClu3La[2]->Fill(minres);}
391  }else if(rollId.station()==2&&rollId.layer()==2){
392  if(cluSize==1*dupli) {hGlobalResClu1La[3]->Fill(minres);}
393  if(cluSize==2*dupli){ hGlobalResClu2La[3]->Fill(minres);}
394  else if(cluSize==3*dupli){ hGlobalResClu3La[3]->Fill(minres);}
395  }else if(rollId.station()==3){
396  if(cluSize==1*dupli) {hGlobalResClu1La[4]->Fill(minres);}
397  else if(cluSize==2*dupli){ hGlobalResClu2La[4]->Fill(minres);}
398  else if(cluSize==3*dupli){ hGlobalResClu3La[4]->Fill(minres);}
399  }
400  meIdRPC.str("");
401  meIdRPC<<"RPCDataOccupancyFromDT_"<<rollId.rawId();
402  meMap[meIdRPC.str()]->Fill(stripPredicted);
403  }
404  }
405  }else{
406  LogDebug("rpcefficiency")<<"DT \t \t \t \t No the prediction is outside of this roll";
407  }//Condition for the right match
408  }else{
409  LogDebug("rpcefficiency")<<"DT \t \t \t No, Exrtrapolation too long!, canceled";
410  }//D so big
411  }//loop over all the rolls asociated
412  }//Is the segment 4D?
413  }else {
414  LogDebug("rpcefficiency")<<"DT \t \t More than one segment in this chamber, or we are in Station 4";
415  }
416  }
417  } else {
418  LogDebug("rpcefficiency")<<"DT This Event doesn't have any DT4DDSegment"; //is ther more than 1 segment in this event?
419  }
420  }
421  }
422 
423  if(incldtMB4){
425  iEvent.getByToken(dt4DSegments, all4DSegments);
426 
427  if(all4DSegments.isValid() && all4DSegments->size()>0){
428 
429  std::map<DTChamberId,int> DTSegmentCounter;
431 
432  for (segment = all4DSegments->begin();segment!=all4DSegments->end(); ++segment){
433  DTSegmentCounter[segment->chamberId()]++;
434  }
435 
436  LogDebug("rpcefficiency")<<"MB4 \t \t Loop Over all4DSegments";
437  for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){
438 
439  DTChamberId DTId = segment->chamberId();
440 
441  if(DTSegmentCounter[DTId] == 1 && DTId.station()==4){
442  int dtWheel = DTId.wheel();
443  int dtStation = DTId.station();
444  int dtSector = DTId.sector();
445 
446  LocalPoint segmentPosition= segment->localPosition();
447  LocalVector segmentDirection=segment->localDirection();
448 
449  //check if the dimension of the segment is 2 and the station is 4
450 
451  if(segment->dimension()==2){
452  LocalVector segmentDirectionMB4=segmentDirection;
453  LocalPoint segmentPositionMB4=segmentPosition;
454 
455 
456  const BoundPlane& DTSurface4 = dtGeo->idToDet(DTId)->surface();
457 
459 
460  for(segMB3=all4DSegments->begin();segMB3!=all4DSegments->end();++segMB3){
461 
462  DTChamberId dtid3 = segMB3->chamberId();
463 
464  if(distsector_tmp(dtid3.sector(),DTId.sector())<=1
465  && dtid3.station()==3
466  && dtid3.wheel()==DTId.wheel()
467  && DTSegmentCounter[dtid3] == 1
468  && segMB3->dimension()==4){
469 
470  const GeomDet* gdet3=dtGeo->idToDet(segMB3->geographicalId());
471  const BoundPlane & DTSurface3 = gdet3->surface();
472 
473  LocalVector segmentDirectionMB3 = segMB3->localDirection();
474  GlobalPoint segmentPositionMB3inGlobal = DTSurface3.toGlobal(segMB3->localPosition());
475 
476  GlobalVector segDirMB4inGlobalFrame=DTSurface4.toGlobal(segmentDirectionMB4);
477  GlobalVector segDirMB3inGlobalFrame=DTSurface3.toGlobal(segmentDirectionMB3);
478 
479  float dx=segDirMB4inGlobalFrame.x();
480  float dy=segDirMB4inGlobalFrame.y();
481  // float dz=segDirMB4inGlobalFrame.z();
482 
483  float dx3=segDirMB3inGlobalFrame.x();
484  float dy3=segDirMB3inGlobalFrame.y();
485  // float dz3=segDirMB3inGlobalFrame.z();
486 
487  double cosAng=fabs(dx*dx3+dy*dy3/sqrt((dx3*dx3+dy3*dy3)*(dx*dx+dy*dy)));
488 
489  if(cosAng>MinCosAng){
490  if(dtSector==13){
491  dtSector=4;
492  }
493  if(dtSector==14){
494  dtSector=10;
495  }
496 
497  std::set<RPCDetId> rollsForThisDT = rollstoreDT[DTStationIndex(0,dtWheel,dtSector,dtStation)]; //It should be always 4
498 
499  assert(rollsForThisDT.size()>=1);
500 
501  for (std::set<RPCDetId>::iterator iteraRoll=rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
502  const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll); //roll asociado a MB4
503  const BoundPlane & RPCSurfaceRB4 = rollasociated->surface(); //surface MB4
504 
505  // RPCGeomServ rpcsrv(rpcId);
506  // std::string nameRoll = rpcsrv.name();
507 
508  GlobalPoint CenterPointRollGlobal=RPCSurfaceRB4.toGlobal(LocalPoint(0,0,0));
509  LocalPoint CenterRollinMB4Frame = DTSurface4.toLocal(CenterPointRollGlobal); //In MB4
510  LocalPoint segmentPositionMB3inMB4Frame = DTSurface4.toLocal(segmentPositionMB3inGlobal); //In MB4
511  LocalVector segmentDirectionMB3inMB4Frame = DTSurface4.toLocal(segDirMB3inGlobalFrame); //In MB4
512 
513  //The exptrapolation is done in MB4 frame. for local x and z is done from MB4,
514  float Dxz=CenterRollinMB4Frame.z();
515  float Xo4=segmentPositionMB4.x();
516  float dxl=segmentDirectionMB4.x(); //dx local for MB4 segment in MB4 Frame
517  float dzl=segmentDirectionMB4.z(); //dx local for MB4 segment in MB4 Frame
518 
519  float X=Xo4+dxl*Dxz/dzl; //In MB4 frame
520  float Z=Dxz;//In MB4 frame
521 
522  //for local y is done from MB3
523  float Yo34=segmentPositionMB3inMB4Frame.y();
524  float dy34 = segmentDirectionMB3inMB4Frame.y();
525  float dz34 = segmentDirectionMB3inMB4Frame.z();
526  float Dy=Dxz-(segmentPositionMB3inMB4Frame.z()); //Distance beetween the segment in MB3 and the RB4 surface
527 
528  float Y=Yo34+dy34*Dy/dz34;//In MB4 Frame
529 
530  const RectangularStripTopology* top_
531  =dynamic_cast<const RectangularStripTopology*>(&(rollasociated->topology())); //Topology roll asociated MB4
532  LocalPoint xmin = top_->localPosition(0.);
533  LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
534  float rsize = fabs( xmax.x()-xmin.x() );
535  float stripl = top_->stripLength();
536  float stripw = top_->pitch();
537 
538  float extrapolatedDistance = sqrt((Y-Yo34)*(Y-Yo34)+Dy*Dy);
539 
540  if(extrapolatedDistance<=MaxDrb4){
541 
542  GlobalPoint GlobalPointExtrapolated = DTSurface4.toGlobal(LocalPoint(X,Y,Z));
543  LocalPoint PointExtrapolatedRPCFrame = RPCSurfaceRB4.toLocal(GlobalPointExtrapolated);
544 
545  if(fabs(PointExtrapolatedRPCFrame.z()) < 5. &&
546  fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
547  fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
548 
549  RPCDetId rollId = rollasociated->id();
550  const float stripPredicted=rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
551 
552  LogDebug("rpcefficiency")<<"MB4 \t \t \t \t Candidate (from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
553  //--------- HISTOGRAM STRIP PREDICTED FROM DT MB4 -------------------
554 
555  std::map<std::string, MonitorElement*> meMap=meCollection[rollId.rawId()];
556  meIdDT.str("");
557  meIdDT<<"ExpectedOccupancyFromDT_"<<rollId.rawId();
558  meMap[meIdDT.str()]->Fill(stripPredicted);
559 
560  //-------RecHitPart Just For Residual--------
561  int countRecHits = 0;
562  int cluSize = 0;
563  float minres = 3000.;
564 
565  LogDebug("rpcefficiency")<<"MB4 \t \t \t \t Getting RecHits in Roll Asociated";
566  typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
567  rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
569 
570  for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
571  countRecHits++;
572  LocalPoint recHitPos=recHit->localPosition();
573  float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
574  LogDebug("rpcefficiency")<<"DT \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction.";
575  if(fabs(res)<fabs(minres)){
576  minres=res;
577  cluSize = recHit->clusterSize();
578  }
579  }
580 
581  if(countRecHits==0){
582  LogDebug("rpcefficiency") <<"MB4 \t \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT";
583  }else{
584  assert(minres!=3000);
585 
586  if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
587  assert(rollId.station()==4);
588  if(cluSize==1*dupli){ hGlobalResClu1La[5]->Fill(minres);}
589  else if(cluSize==2*dupli){ hGlobalResClu2La[5]->Fill(minres);}
590  else if(cluSize==3*dupli){ hGlobalResClu3La[5]->Fill(minres);}
591 
592  meIdRPC.str("");
593  meIdRPC<<"RPCDataOccupancyFromDT_"<<rollId.rawId();
594  meMap[meIdRPC.str()]->Fill(stripPredicted);
595  }
596  }
597  }else{
598  LogDebug("rpcefficiency")<<"MB4 \t \t \t \t No the prediction is outside of this roll";
599  }
600  }//Condition for the right match
601  else{
602  LogDebug("rpcefficiency")<<"MB4 \t \t \t No, Exrtrapolation too long!, canceled";
603  }
604  }//loop over all the rollsasociated
605  }else{
606  LogDebug("rpcefficiency")<<"MB4 \t \t \t \t I found segments in MB4 and MB3 adjacent or same wheel and sector but not compatibles Diferent Directions";
607  }
608  }else{//if dtid3.station()==3&&dtid3.sector()==DTId.sector()&&dtid3.wheel()==DTId.wheel()&&segMB3->dim()==4
609  LogDebug("rpcefficiency")<<"MB4 \t \t \t No the same station or same wheel or segment dim in mb3 not 4D";
610  }
611  }//loop over all the segments looking for one in MB3
612  }else{
613  LogDebug("rpcefficiency")<<"MB4 \t \t \t Is NOT a 2D Segment";
614  }
615  }else{
616  LogDebug("rpcefficiency")<<"MB4 \t \t \t \t There is not just one segment or is not in station 4";
617  }//De aca para abajo esta en dtpart.inl
618  }
619  }else{
620  LogDebug("rpcefficiency")<<"MB4 \t This event doesn't have 4D Segment";
621  }
622  }
623 
624 
625  if(inclcsc){
626  LogDebug("rpcefficiency") <<"\t Getting the CSC Segments";
627  edm::Handle<CSCSegmentCollection> allCSCSegments;
628 
629  iEvent.getByToken(cscSegments, allCSCSegments);
630 
631  if(allCSCSegments.isValid()){
632  if(allCSCSegments->size()>0){
633  statistics->Fill(18);
634 
635  LogDebug("rpcefficiency")<<"CSC \t Number of CSC Segments in this event = "<<allCSCSegments->size();
636 
637  std::map<CSCDetId,int> CSCSegmentsCounter;
639 
640  int segmentsInThisEventInTheEndcap=0;
641 
642  for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
643  CSCSegmentsCounter[segment->cscDetId()]++;
644  segmentsInThisEventInTheEndcap++;
645  }
646 
647  statistics->Fill(allCSCSegments->size()+18);
648 
649  LogDebug("rpcefficiency")<<"CSC \t loop over all the CSCSegments ";
650  for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
651  CSCDetId CSCId = segment->cscDetId();
652 
653  if(CSCSegmentsCounter[CSCId]==1 && CSCId.ring()!=1 && allCSCSegments->size()>=2){
654  LogDebug("rpcefficiency")<<"CSC \t \t yes";
655  int cscEndCap = CSCId.endcap();
656  int cscStation = CSCId.station();
657  int cscRing = CSCId.ring();
658  // int cscChamber = CSCId.chamber();
659  int rpcRegion = 1; if(cscEndCap==2) rpcRegion= -1;//Relation among the endcaps
660  int rpcRing = cscRing;
661  if(cscRing==4)rpcRing =1;
662  int rpcStation = cscStation;
663  int rpcSegment = CSCId.chamber();
664 
665  LocalPoint segmentPosition= segment->localPosition();
666  LocalVector segmentDirection=segment->localDirection();
667  float dz=segmentDirection.z();
668 
669  LogDebug("rpcefficiency")<<"CSC \t \t Is a good Segment? dim = 4, 4 <= nRecHits <= 10 Incident angle int range 45 < "<<acos(dz)*180/3.1415926<<" < 135? ";
670 
671  if(segment->dimension()==4 && (segment->nRecHits()<=10 && segment->nRecHits()>=4)&& acos(dz)*180/3.1415926 > 45. && acos(dz)*180/3.1415926 < 160. ){
672 
673  float Xo=segmentPosition.x();
674  float Yo=segmentPosition.y();
675  float Zo=segmentPosition.z();
676  float dx=segmentDirection.x();
677  float dy=segmentDirection.y();
678  float dz=segmentDirection.z();
679 
680 
681  LogDebug("rpcefficiency")<<"CSC \t \t Getting chamber from Geometry";
682  const CSCChamber* TheChamber=cscGeo->chamber(CSCId);
683  LogDebug("rpcefficiency")<<"CSC \t \t Getting ID from Chamber";
684  const CSCDetId TheId=TheChamber->id();
685  LogDebug("rpcefficiency")<<"CSC \t \t Printing The Id"<<TheId;
686  std::set<RPCDetId> rollsForThisCSC = rollstoreCSC[CSCStationIndex(rpcRegion,rpcStation,rpcRing,rpcSegment)];
687  LogDebug("rpcefficiency")<<"CSC \t \t Number of rolls for this CSC = "<<rollsForThisCSC.size();
688 
689  if(rpcRing!=1){
690 
691  //Loop over all the rolls
692  for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisCSC.begin();iteraRoll != rollsForThisCSC.end(); iteraRoll++){
693 
694  const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
695  RPCDetId rpcId = rollasociated->id();
696 
697  const BoundPlane & RPCSurface = rollasociated->surface();
698 
699  GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
700  LocalPoint CenterRollinCSCFrame = TheChamber->toLocal(CenterPointRollGlobal);
701 
702  float D=CenterRollinCSCFrame.z();
703 
704  float X=Xo+dx*D/dz;
705  float Y=Yo+dy*D/dz;
706  float Z=D;
707 
708  const TrapezoidalStripTopology* top_=dynamic_cast<const TrapezoidalStripTopology*>(&(rollasociated->topology()));
709  LocalPoint xmin = top_->localPosition(0.);
710  LogDebug("rpcefficiency")<<"CSC \t \t \t xmin of this Roll "<<xmin<<"cm";
711  LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
712  LogDebug("rpcefficiency")<<"CSC \t \t \t xmax of this Roll "<<xmax<<"cm";
713  float rsize = fabs( xmax.x()-xmin.x() );
714  LogDebug("rpcefficiency")<<"CSC \t \t \t Roll Size "<<rsize<<"cm";
715  float stripl = top_->stripLength();
716  float stripw = top_->pitch();
717 
718 
719  float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
720 
721 
722  if(extrapolatedDistance<=MaxD){
723 
724  GlobalPoint GlobalPointExtrapolated=TheChamber->toGlobal(LocalPoint(X,Y,Z));
725  LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
726 
727 
728  if(fabs(PointExtrapolatedRPCFrame.z()) < 10. &&
729  fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
730  fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
731 
732  RPCDetId rollId = rollasociated->id();
733  RPCGeomServ rpcsrv(rollId);
734  std::string nameRoll = rpcsrv.name();
735 
736  LogDebug("rpcefficiency")<<"CSC \t \t \t \t The RPCName is "<<nameRoll;
737 
738  const float stripPredicted =
739  rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
740 
741  LogDebug("rpcefficiency")<<"CSC \t \t \t \t \t Candidate"<<rollId<<" "<<"(from CSC Segment) STRIP---> "<<stripPredicted<< std::endl;
742  //--------- HISTOGRAM STRIP PREDICTED FROM CSC -------------------
743 
744  std::map<std::string, MonitorElement*> meMap=meCollection[rpcId.rawId()];
745  meIdCSC.str("");
746  meIdCSC<<"ExpectedOccupancyFromCSC_"<<rollId.rawId();
747  meMap[meIdCSC.str()]->Fill(stripPredicted);
748  //--------------------------------------------------------------------
749 
750 
751  //-------RecHitPart Just For Residual--------
752  int cluSize = 0;
753  int countRecHits = 0;
754  float minres = 3000.;
755 
756  LogDebug("rpcefficiency")<<"CSC \t \t \t \t \t Getting RecHits in Roll Asociated";
757  typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
758  rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
760 
761  for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
762 
763  countRecHits++;
764  LocalPoint recHitPos=recHit->localPosition();
765  float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
766  LogDebug("rpcefficiency")<<"CSC \t \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction.";
767  if(fabs(res)<fabs(minres)){
768  minres=res;
769  cluSize = recHit->clusterSize();
770  LogDebug("rpcefficiency")<<"CSC \t \t \t \t \t \t \t New Min Res "<<res<<"cm.";
771  }
772  }
773 
774  if(countRecHits==0){
775  LogDebug("rpcefficiency") <<"CSC \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT";
776  }else{
777  assert(minres!=3000);
778 
779  if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
780  LogDebug("rpcefficiency")<<"CSC \t \t \t \t \t \t True!";
781 
782  if(rollId.ring()==2&&rollId.roll()==1){
783  if(cluSize==1*dupli) { hGlobalResClu1R2A->Fill(minres);}
784  else if(cluSize==2*dupli) { hGlobalResClu2R2A->Fill(minres); }
785  else if(cluSize==3*dupli) {hGlobalResClu3R2A->Fill(minres);}
786  }
787  else if(rollId.ring()==2&&rollId.roll()==2){
788  if(cluSize==1*dupli) { hGlobalResClu1R2B->Fill(minres);}
789  else if(cluSize==2*dupli) { hGlobalResClu2R2B->Fill(minres); }
790  else if(cluSize==3*dupli) {hGlobalResClu3R2B->Fill(minres);}
791  }
792  else if(rollId.ring()==2&&rollId.roll()==3){
793  if(cluSize==1*dupli) { hGlobalResClu1R2C->Fill(minres);}
794  else if(cluSize==2*dupli) {hGlobalResClu2R2C->Fill(minres); }
795  else if(cluSize==3*dupli) { hGlobalResClu3R2C->Fill(minres);}
796  }
797  else if(rollId.ring()==3&&rollId.roll()==1){
798  if(cluSize==1*dupli) {hGlobalResClu1R3A->Fill(minres); }
799  else if(cluSize==2*dupli) { hGlobalResClu2R3A->Fill(minres); }
800  else if(cluSize==3*dupli) {hGlobalResClu3R3A->Fill(minres);}
801  }
802  else if(rollId.ring()==3&&rollId.roll()==2){
803  if(cluSize==1*dupli) { hGlobalResClu1R3B->Fill(minres); }
804  else if(cluSize==2*dupli) { hGlobalResClu2R3B->Fill(minres); }
805  else if(cluSize==3*dupli) { hGlobalResClu3R3B->Fill(minres);}
806  }
807  else if(rollId.ring()==3&&rollId.roll()==3){
808  if(cluSize==1*dupli) {hGlobalResClu1R3C->Fill(minres); }
809  else if(cluSize==2*dupli) {hGlobalResClu2R3C->Fill(minres); }
810  else if(cluSize==3*dupli) {hGlobalResClu3R3C->Fill(minres);}
811  }
812  meIdRPC.str("");
813  meIdRPC<<"RPCDataOccupancyFromCSC_"<<rollId.rawId();
814  meMap[meIdRPC.str()]->Fill(stripPredicted);
815  }
816  }
817 
818  }else{
819  LogDebug("rpcefficiency")<<"CSC \t \t \t \t No the prediction is outside of this roll";
820  }//Condition for the right match
821  }else{//if extrapolation distance D is not too long
822  LogDebug("rpcefficiency")<<"CSC \t \t \t No, Exrtrapolation too long!, canceled";
823  }//D so big
824  }//loop over the rolls asociated
825  }//Condition over the LS1 geometry!!!!
826  }//Is the segment 4D?
827  }else{
828  LogDebug("rpcefficiency")<<"CSC \t \t More than one segment in this chamber, or we are in Station Ring 1 or in Station 4";
829  }
830  }
831  }else{
832  LogDebug("rpcefficiency")<<"CSC This Event doesn't have any CSCSegment";
833  }
834  }
835  }
836 
837 }
838 
839 
840 
#define LogDebug(id)
const double Z[kNumberCalorimeter]
int chamber() const
Definition: CSCDetId.h:68
T getParameter(std::string const &) const
float strip(const LocalPoint &lp) const
Definition: RPCRoll.cc:71
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
const Topology & topology() const
Definition: RPCRoll.cc:30
int distsector_tmp(int sector1, int sector2)
MonitorElement * hGlobalResClu1R2A
virtual float stripLength() const
RPCEfficiency(const edm::ParameterSet &)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
CSCDetId id() const
Get the (concrete) DetId.
Definition: CSCChamber.h:37
void cd(void)
Definition: DQMStore.cc:265
MonitorElement * hGlobalResClu1R3A
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:52
assert(m_qm.get())
int nstrips() const
Definition: RPCRoll.cc:46
MonitorElement * hGlobalResClu3R2C
MonitorElement * hGlobalResClu2La[6]
T y() const
Definition: PV3DBase.h:63
#define X(str)
Definition: MuonsGrabber.cc:48
void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
MonitorElement * hGlobalResClu3R3C
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:67
std::map< int, std::map< std::string, MonitorElement * > > meCollection
void bookDetUnitSeg(DQMStore::IBooker &, RPCDetId &detId, int nstrips, std::string folder, std::map< std::string, MonitorElement * > &)
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
void Fill(long long x)
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
int endcap() const
Definition: CSCDetId.h:93
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
std::string EffRootFileName
RPCDetId id() const
Definition: RPCRoll.cc:24
int iEvent
Definition: GenABIO.cc:230
virtual std::string name()
Definition: RPCGeomServ.cc:20
int roll() const
Definition: RPCDetId.h:120
edm::EDGetTokenT< RPCRecHitCollection > RPCRecHitLabel_
MonitorElement * hGlobalResClu3La[6]
int ring() const
Definition: RPCDetId.h:72
T sqrt(T t)
Definition: SSEVec.h:18
std::string folderPath
T z() const
Definition: PV3DBase.h:64
MonitorElement * hGlobalResClu3R2A
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
const std::vector< const RPCRoll * > & rolls() const
Return the Rolls.
Definition: RPCChamber.cc:68
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double rangestripsRB4
MonitorElement * hGlobalResClu1R3B
MonitorElement * hGlobalResClu1R2C
MonitorElement * hGlobalResClu2R3C
MonitorElement * hGlobalResClu3R3A
bool isValid() const
Definition: HandleBase.h:75
std::map< CSCStationIndex, std::set< RPCDetId > > rollstoreCSC
MonitorElement * statistics
edm::EDGetTokenT< CSCSegmentCollection > cscSegments
std::map< DTStationIndex, std::set< RPCDetId > > rollstoreDT
edm::EDGetTokenT< DTRecSegment4DCollection > dt4DSegments
int ring() const
Definition: CSCDetId.h:75
int layer() const
Definition: RPCDetId.h:108
virtual int segment()
Definition: RPCGeomServ.cc:469
MonitorElement * hGlobalResClu2R2A
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:273
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:150
MonitorElement * hGlobalResClu2R3A
const T & get() const
Definition: EventSetup.h:56
MonitorElement * hGlobalResClu1La[6]
virtual LocalPoint localPosition(float strip) const
int sector() const
Sector id: the group of chambers at same phi (and increasing r)
Definition: RPCDetId.h:102
double rangestrips
virtual float pitch() const
int sector() const
Definition: DTChamberId.h:61
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
virtual void analyze(const edm::Event &, const edm::EventSetup &)
int station() const
Definition: CSCDetId.h:86
MonitorElement * hGlobalResClu1R3C
MonitorElement * hGlobalResClu2R2B
virtual LocalPoint localPosition(float strip) const
int station() const
Return the station number.
Definition: DTChamberId.h:51
MonitorElement * hGlobalResClu3R3B
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:45
T x() const
Definition: PV3DBase.h:62
virtual float stripLength() const
det heigth (strip length in the middle)
MonitorElement * hGlobalResClu3R2B
MonitorElement * hGlobalResClu2R3B
MonitorElement * hGlobalResClu1R2B
MonitorElement * hGlobalResClu2R2C
Definition: Run.h:43
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.
Definition: RPCDetId.h:63
int station() const
Definition: RPCDetId.h:96