00001
00002
00003
00004
00005
00006
00007
00008 #include "DQM/RPCMonitorDigi/interface/RPCEfficiency.h"
00009 #include <sstream>
00010 #include "FWCore/Framework/interface/MakerMacros.h"
00011 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
00012 #include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
00013 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
00014 #include <Geometry/RPCGeometry/interface/RPCGeomServ.h>
00015 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
00016 #include <Geometry/Records/interface/MuonGeometryRecord.h>
00017 #include <Geometry/CommonTopologies/interface/RectangularStripTopology.h>
00018 #include <Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h>
00019
00020
00021 void RPCEfficiency::beginJob(){}
00022
00023 int distsector_tmp(int sector1,int sector2){
00024
00025 if(sector1==13) sector1=4;
00026 if(sector1==14) sector1=10;
00027
00028 if(sector2==13) sector2=4;
00029 if(sector2==14) sector2=10;
00030
00031 int distance = abs(sector1 - sector2);
00032 if(distance>6) distance = 12-distance;
00033 return distance;
00034 }
00035
00036
00037 RPCEfficiency::RPCEfficiency(const edm::ParameterSet& iConfig){
00038 incldt=iConfig.getUntrackedParameter<bool>("incldt",true);
00039 incldtMB4=iConfig.getUntrackedParameter<bool>("incldtMB4",true);
00040 inclcsc=iConfig.getUntrackedParameter<bool>("inclcsc",true);
00041 debug=iConfig.getUntrackedParameter<bool>("debug",false);
00042
00043 rangestrips = iConfig.getUntrackedParameter<double>("rangestrips",4.);
00044 rangestripsRB4=iConfig.getUntrackedParameter<double>("rangestripsRB4",4.);
00045 dupli = iConfig.getUntrackedParameter<int>("DuplicationCorrection",2);
00046 MinCosAng=iConfig.getUntrackedParameter<double>("MinCosAng",0.96);
00047 MaxD=iConfig.getUntrackedParameter<double>("MaxD",80.);
00048 MaxDrb4=iConfig.getUntrackedParameter<double>("MaxDrb4",150.);
00049
00050
00051
00052 cscSegments=iConfig.getParameter<edm::InputTag>("cscSegments");
00053 dt4DSegments=iConfig.getParameter<edm::InputTag>("dt4DSegments");
00054 RPCRecHitLabel_ = iConfig.getParameter<edm::InputTag>("RecHitLabel");
00055
00056
00057 folderPath=iConfig.getUntrackedParameter<std::string>("folderPath","RPC/RPCEfficiency/");
00058
00059 EffSaveRootFile = iConfig.getUntrackedParameter<bool>("EffSaveRootFile", false);
00060 EffRootFileName = iConfig.getUntrackedParameter<std::string>("EffRootFileName", "RPCEfficiency.root");
00061
00062
00063
00064 dbe = edm::Service<DQMStore>().operator->();
00065
00066 std::string folder;
00067 dbe->setCurrentFolder(folderPath);
00068 statistics = dbe->book1D("Statistics","All Statistics",33,0.5,33.5);
00069
00070 statistics->setBinLabel(1,"Events ",1);
00071 statistics->setBinLabel(2,"Events with DT seg",1);
00072 statistics->setBinLabel(3,"1 DT seg",1);
00073 statistics->setBinLabel(4,"2 DT seg",1);
00074 statistics->setBinLabel(5,"3 DT seg",1);
00075 statistics->setBinLabel(6,"4 DT seg",1);
00076 statistics->setBinLabel(7,"5 DT seg",1);
00077 statistics->setBinLabel(8,"6 DT seg",1);
00078 statistics->setBinLabel(9,"7 DT seg",1);
00079 statistics->setBinLabel(10,"8 DT seg",1);
00080 statistics->setBinLabel(11,"9 DT seg",1);
00081 statistics->setBinLabel(12,"10 DT seg",1);
00082 statistics->setBinLabel(13,"11 DT seg",1);
00083 statistics->setBinLabel(14,"12 DT seg",1);
00084 statistics->setBinLabel(15,"13 DT seg",1);
00085 statistics->setBinLabel(16,"14 DT seg",1);
00086 statistics->setBinLabel(17,"15 DT seg",1);
00087 statistics->setBinLabel(18,"Events with CSC seg",1);
00088 statistics->setBinLabel(16+3,"1 CSC seg",1);
00089 statistics->setBinLabel(16+4,"2 CSC seg",1);
00090 statistics->setBinLabel(16+5,"3 CSC seg",1);
00091 statistics->setBinLabel(16+6,"4 CSC seg",1);
00092 statistics->setBinLabel(16+7,"5 CSC seg",1);
00093 statistics->setBinLabel(16+8,"6 CSC seg",1);
00094 statistics->setBinLabel(16+9,"7 CSC seg",1);
00095 statistics->setBinLabel(16+10,"8 CSC seg",1);
00096 statistics->setBinLabel(16+11,"9 CSC seg",1);
00097 statistics->setBinLabel(16+12,"10 CSC seg",1);
00098 statistics->setBinLabel(16+13,"11 CSC seg",1);
00099 statistics->setBinLabel(16+14,"12 CSC seg",1);
00100 statistics->setBinLabel(16+15,"13 CSC seg",1);
00101 statistics->setBinLabel(16+16,"14 CSC seg",1);
00102 statistics->setBinLabel(16+17,"15 CSC seg",1);
00103
00104 if(debug) std::cout<<"booking Global histograms with "<<folderPath<<std::endl;
00105
00106 folder = folderPath+"MuonSegEff/"+"Residuals/Barrel";
00107 dbe->setCurrentFolder(folder);
00108
00109
00110 std::stringstream histoName, histoTitle;
00111
00112 for (int layer = 1 ; layer<= 6 ;layer++){
00113 histoName.str("");
00114 histoTitle.str("");
00115 histoName<<"GlobalResidualsClu1La"<<layer;
00116 histoTitle<<"RPC Residuals Layer "<<layer<<" Cluster Size 1";
00117 hGlobalResClu1La[layer-1] = dbe->book1D(histoName.str(), histoTitle.str(),101,-10.,10.);
00118
00119 histoName.str("");
00120 histoTitle.str("");
00121 histoName<<"GlobalResidualsClu2La"<<layer;
00122 histoTitle<<"RPC Residuals Layer "<<layer<<" Cluster Size 2";
00123 hGlobalResClu2La[layer-1] = dbe->book1D(histoName.str(), histoTitle.str(),101,-10.,10.);
00124
00125 histoName.str("");
00126 histoTitle.str("");
00127 histoName<<"GlobalResidualsClu3La"<<layer;
00128 histoTitle<<"RPC Residuals Layer "<<layer<<" Cluster Size 3";
00129 hGlobalResClu3La[layer-1] = dbe->book1D(histoName.str(), histoTitle.str(),101,-10.,10.);
00130
00131 }
00132
00133 if(debug) std::cout<<"Booking Residuals for EndCap"<<std::endl;
00134 folder = folderPath+"MuonSegEff/Residuals/EndCap";
00135 dbe->setCurrentFolder(folder);
00136
00137
00138
00139 hGlobalResClu1R3C = dbe->book1D("GlobalResidualsClu1R3C","RPC Residuals Ring 3 Roll C Cluster Size 1",101,-10.,10.);
00140 hGlobalResClu1R3B = dbe->book1D("GlobalResidualsClu1R3B","RPC Residuals Ring 3 Roll B Cluster Size 1",101,-10.,10.);
00141 hGlobalResClu1R3A = dbe->book1D("GlobalResidualsClu1R3A","RPC Residuals Ring 3 Roll A Cluster Size 1",101,-10.,10.);
00142 hGlobalResClu1R2C = dbe->book1D("GlobalResidualsClu1R2C","RPC Residuals Ring 2 Roll C Cluster Size 1",101,-10.,10.);
00143 hGlobalResClu1R2B = dbe->book1D("GlobalResidualsClu1R2B","RPC Residuals Ring 2 Roll B Cluster Size 1",101,-10.,10.);
00144 hGlobalResClu1R2A = dbe->book1D("GlobalResidualsClu1R2A","RPC Residuals Ring 2 Roll A Cluster Size 1",101,-10.,10.);
00145
00146 hGlobalResClu2R3C = dbe->book1D("GlobalResidualsClu2R3C","RPC Residuals Ring 3 Roll C Cluster Size 2",101,-10.,10.);
00147 hGlobalResClu2R3B = dbe->book1D("GlobalResidualsClu2R3B","RPC Residuals Ring 3 Roll B Cluster Size 2",101,-10.,10.);
00148 hGlobalResClu2R3A = dbe->book1D("GlobalResidualsClu2R3A","RPC Residuals Ring 3 Roll A Cluster Size 2",101,-10.,10.);
00149 hGlobalResClu2R2C = dbe->book1D("GlobalResidualsClu2R2C","RPC Residuals Ring 2 Roll C Cluster Size 2",101,-10.,10.);
00150 hGlobalResClu2R2B = dbe->book1D("GlobalResidualsClu2R2B","RPC Residuals Ring 2 Roll B Cluster Size 2",101,-10.,10.);
00151 hGlobalResClu2R2A = dbe->book1D("GlobalResidualsClu2R2A","RPC Residuals Ring 2 Roll A Cluster Size 2",101,-10.,10.);
00152
00153 hGlobalResClu3R3C = dbe->book1D("GlobalResidualsClu3R3C","RPC Residuals Ring 3 Roll C Cluster Size 3",101,-10.,10.);
00154 hGlobalResClu3R3B = dbe->book1D("GlobalResidualsClu3R3B","RPC Residuals Ring 3 Roll B Cluster Size 3",101,-10.,10.);
00155 hGlobalResClu3R3A = dbe->book1D("GlobalResidualsClu3R3A","RPC Residuals Ring 3 Roll A Cluster Size 3",101,-10.,10.);
00156 hGlobalResClu3R2C = dbe->book1D("GlobalResidualsClu3R2C","RPC Residuals Ring 2 Roll C Cluster Size 3",101,-10.,10.);
00157 hGlobalResClu3R2B = dbe->book1D("GlobalResidualsClu3R2B","RPC Residuals Ring 2 Roll B Cluster Size 3",101,-10.,10.);
00158 hGlobalResClu3R2A = dbe->book1D("GlobalResidualsClu3R2A","RPC Residuals Ring 2 Roll A Cluster Size 3",101,-10.,10.);
00159
00160 }
00161
00162 void RPCEfficiency::beginRun(const edm::Run& run, const edm::EventSetup& iSetup){
00163
00164 edm::ESHandle<RPCGeometry> rpcGeo;
00165 iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00166
00167
00168 for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00169 if(dynamic_cast< RPCChamber* >( *it ) != 0 ){
00170 RPCChamber* ch = dynamic_cast< RPCChamber* >( *it );
00171 std::vector< const RPCRoll*> roles = (ch->rolls());
00172 for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){
00173
00174 RPCDetId rpcId = (*r)->id();
00175 int region=rpcId.region();
00176
00177
00178
00179
00180 if(debug) std::cout<<"Booking for "<<rpcId.rawId()<<std::endl;
00181
00182 bookDetUnitSeg(rpcId,(*r)->nstrips(),folderPath+"MuonSegEff/", meCollection[rpcId.rawId()] );
00183
00184 if(region==0&&(incldt||incldtMB4)){
00185
00186 int wheel=rpcId.ring();
00187 int sector=rpcId.sector();
00188 int station=rpcId.station();
00189 DTStationIndex ind(region,wheel,sector,station);
00190 std::set<RPCDetId> myrolls;
00191 if (rollstoreDT.find(ind)!=rollstoreDT.end()) myrolls=rollstoreDT[ind];
00192 myrolls.insert(rpcId);
00193 rollstoreDT[ind]=myrolls;
00194
00195 }else if(region!=0 && inclcsc){
00196 int station=rpcId.station();
00197 int ring=rpcId.ring();
00198 int cscring=ring;
00199 int cscstation=station;
00200 RPCGeomServ rpcsrv(rpcId);
00201 int rpcsegment = rpcsrv.segment();
00202 int cscchamber = rpcsegment;
00203 if((station==2||station==3)&&ring==3){
00204 cscring = 2;
00205 }
00206
00207 CSCStationIndex ind(region,cscstation,cscring,cscchamber);
00208 std::set<RPCDetId> myrolls;
00209 if (rollstoreCSC.find(ind)!=rollstoreCSC.end()){
00210 myrolls=rollstoreCSC[ind];
00211 }
00212 myrolls.insert(rpcId);
00213 rollstoreCSC[ind]=myrolls;
00214 }
00215 }
00216 }
00217 }
00218 for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00219 if( dynamic_cast< RPCChamber* >( *it ) != 0 ){
00220
00221 RPCChamber* ch = dynamic_cast< RPCChamber* >( *it );
00222 std::vector< const RPCRoll*> roles = (ch->rolls());
00223 for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){
00224 RPCDetId rpcId = (*r)->id();
00225
00226 int region=rpcId.region();
00227
00228 if(region!=0 && inclcsc && (rpcId.ring()==2 || rpcId.ring()==3)){
00229 int region=rpcId.region();
00230 int station=rpcId.station();
00231 int ring=rpcId.ring();
00232 int cscring = ring;
00233
00234 if((station==2||station==3)&&ring==3) cscring = 2;
00235
00236
00237 int cscstation=station;
00238 RPCGeomServ rpcsrv(rpcId);
00239 int rpcsegment = rpcsrv.segment();
00240
00241 int cscchamber = rpcsegment+1;
00242 if(cscchamber==37)cscchamber=1;
00243 CSCStationIndex ind(region,cscstation,cscring,cscchamber);
00244 std::set<RPCDetId> myrolls;
00245 if (rollstoreCSC.find(ind)!=rollstoreCSC.end())myrolls=rollstoreCSC[ind];
00246 myrolls.insert(rpcId);
00247 rollstoreCSC[ind]=myrolls;
00248
00249 cscchamber = rpcsegment-1;
00250 if(cscchamber==0)cscchamber=36;
00251 CSCStationIndex indDos(region,cscstation,cscring,cscchamber);
00252 std::set<RPCDetId> myrollsDos;
00253 if (rollstoreCSC.find(indDos)!=rollstoreCSC.end()) myrollsDos=rollstoreCSC[indDos];
00254 myrollsDos.insert(rpcId);
00255 rollstoreCSC[indDos]=myrollsDos;
00256 }
00257 }
00258 }
00259 }
00260 }
00261
00262 RPCEfficiency::~RPCEfficiency(){}
00263
00264 void RPCEfficiency::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup){
00265
00266
00267 edm::ESHandle<RPCGeometry> rpcGeo;
00268 edm::ESHandle<DTGeometry> dtGeo;
00269 edm::ESHandle<CSCGeometry> cscGeo;
00270
00271 iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00272 iSetup.get<MuonGeometryRecord>().get(dtGeo);
00273 iSetup.get<MuonGeometryRecord>().get(cscGeo);
00274
00275 statistics->Fill(1);
00276
00277 std::stringstream meIdRPC, meIdDT, meIdCSC;
00278
00279 if(debug) std::cout <<"\t Getting the RPC RecHits"<<std::endl;
00280 edm::Handle<RPCRecHitCollection> rpcHits;
00281 iEvent.getByLabel(RPCRecHitLabel_,rpcHits);
00282
00283 if(!rpcHits.isValid()) return;
00284
00285 if(incldt){
00286 if(debug) std::cout<<"\t Getting the DT Segments"<<std::endl;
00287 edm::Handle<DTRecSegment4DCollection> all4DSegments;
00288
00289 iEvent.getByLabel(dt4DSegments, all4DSegments);
00290
00291 if(all4DSegments.isValid()){
00292
00293 if(all4DSegments->size()>0){
00294
00295 if(all4DSegments->size()<=16) statistics->Fill(2);
00296
00297 if(debug) std::cout<<"\t Number of DT Segments in this event = "<<all4DSegments->size()<<std::endl;
00298
00299 std::map<DTChamberId,int> DTSegmentCounter;
00300 DTRecSegment4DCollection::const_iterator segment;
00301
00302 for (segment = all4DSegments->begin();segment!=all4DSegments->end(); ++segment){
00303 DTSegmentCounter[segment->chamberId()]++;
00304 }
00305
00306 statistics->Fill(all4DSegments->size()+2);
00307
00308 if(debug) std::cout<<"\t Loop over all the 4D Segments"<<std::endl;
00309 for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){
00310
00311 DTChamberId DTId = segment->chamberId();
00312
00313
00314 if(DTSegmentCounter[DTId]==1 && DTId.station()!=4){
00315
00316 int dtWheel = DTId.wheel();
00317 int dtStation = DTId.station();
00318 int dtSector = DTId.sector();
00319
00320 LocalPoint segmentPosition= segment->localPosition();
00321 LocalVector segmentDirection=segment->localDirection();
00322
00323 const GeomDet* gdet=dtGeo->idToDet(segment->geographicalId());
00324 const BoundPlane & DTSurface = gdet->surface();
00325
00326
00327
00328 if(segment->dimension()==4){
00329
00330 float Xo=segmentPosition.x();
00331 float Yo=segmentPosition.y();
00332 float Zo=segmentPosition.z();
00333 float dx=segmentDirection.x();
00334 float dy=segmentDirection.y();
00335 float dz=segmentDirection.z();
00336
00337 std::set<RPCDetId> rollsForThisDT = rollstoreDT[DTStationIndex(0,dtWheel,dtSector,dtStation)];
00338
00339 if(debug) std::cout<<"DT \t \t Loop over all the rolls asociated to this DT"<<std::endl;
00340 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00341 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00342 RPCDetId rpcId = rollasociated->id();
00343 const BoundPlane & RPCSurface = rollasociated->surface();
00344
00345
00346
00347
00348 GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
00349
00350 LocalPoint CenterRollinDTFrame = DTSurface.toLocal(CenterPointRollGlobal);
00351
00352 float D=CenterRollinDTFrame.z();
00353
00354 float X=Xo+dx*D/dz;
00355 float Y=Yo+dy*D/dz;
00356 float Z=D;
00357
00358 const RectangularStripTopology* top_= dynamic_cast<const RectangularStripTopology*> (&(rollasociated->topology()));
00359 LocalPoint xmin = top_->localPosition(0.);
00360 if(debug) std::cout<<"DT \t \t \t xmin of this Roll "<<xmin<<"cm"<<std::endl;
00361 LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00362 if(debug) std::cout<<"DT \t \t \t xmax of this Roll "<<xmax<<"cm"<<std::endl;
00363 float rsize = fabs( xmax.x()-xmin.x() );
00364 if(debug) std::cout<<"DT \t \t \t Roll Size "<<rsize<<"cm"<<std::endl;
00365 float stripl = top_->stripLength();
00366 float stripw = top_->pitch();
00367
00368 float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
00369
00370 if(extrapolatedDistance<=MaxD){
00371
00372 GlobalPoint GlobalPointExtrapolated = DTSurface.toGlobal(LocalPoint(X,Y,Z));
00373 LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
00374
00375 if(fabs(PointExtrapolatedRPCFrame.z()) < 10. &&
00376 fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
00377 fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
00378
00379 RPCDetId rollId = rollasociated->id();
00380 RPCGeomServ rpcsrv(rollId);
00381 std::string nameRoll = rpcsrv.name();
00382 if(debug) std::cout<<"DT \t \t \t \t The RPCName is "<<nameRoll<<std::endl;
00383 const float stripPredicted =
00384 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
00385
00386 if(debug) std::cout<<"DT \t \t \t \t Candidate (from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
00387
00388
00389 std::map<std::string, MonitorElement*> meMap=meCollection[rpcId.rawId()];
00390 meIdDT.str("");
00391 meIdDT<<"ExpectedOccupancyFromDT_"<<rollId.rawId();
00392 meMap[meIdDT.str()]->Fill(stripPredicted);
00393
00394
00395
00396
00397 int countRecHits = 0;
00398 int cluSize = 0;
00399 float minres = 3000.;
00400
00401 typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
00402 rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
00403 RPCRecHitCollection::const_iterator recHit;
00404
00405 for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
00406 countRecHits++;
00407
00408
00409
00410
00411 LocalPoint recHitPos=recHit->localPosition();
00412 float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
00413 if(debug) std::cout<<"DT \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction."<<std::endl;
00414 if(fabs(res)<fabs(minres)){
00415 minres=res;
00416 cluSize = recHit->clusterSize();
00417 if(debug) std::cout<<"DT \t \t \t \t \t \t New Min Res "<<res<<"cm."<<std::endl;
00418 }
00419 }
00420
00421 if(countRecHits==0){
00422 if(debug) std::cout <<"DT \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT"<<std::endl;
00423 }else{
00424 assert(minres!=3000);
00425
00426 if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
00427 if(debug) std::cout<<"DT \t \t \t \t \t \t True!"<<std::endl;
00428
00429
00430
00431 if(rollId.station()==1&&rollId.layer()==1) {
00432 if(cluSize==1*dupli) {hGlobalResClu1La[0]->Fill(minres);}
00433 else if(cluSize==2*dupli){ hGlobalResClu2La[0]->Fill(minres);}
00434 else if(cluSize==3*dupli){ hGlobalResClu3La[0]->Fill(minres);}}
00435 else if(rollId.station()==1&&rollId.layer()==2){
00436 if(cluSize==1*dupli) {hGlobalResClu1La[1]->Fill(minres);}
00437 else if(cluSize==2*dupli){ hGlobalResClu2La[1]->Fill(minres);}
00438 else if(cluSize==3*dupli){ hGlobalResClu3La[1]->Fill(minres);}}
00439 else if(rollId.station()==2&&rollId.layer()==1){
00440 if(cluSize==1*dupli) {hGlobalResClu1La[2]->Fill(minres);}
00441 else if(cluSize==2*dupli){ hGlobalResClu2La[2]->Fill(minres);}
00442 else if(cluSize==3*dupli){ hGlobalResClu3La[2]->Fill(minres);}
00443 }
00444 else if(rollId.station()==2&&rollId.layer()==2){
00445 if(cluSize==1*dupli) {hGlobalResClu1La[3]->Fill(minres);}
00446 if(cluSize==2*dupli){ hGlobalResClu2La[3]->Fill(minres);}
00447 else if(cluSize==3*dupli){ hGlobalResClu3La[3]->Fill(minres);}
00448 }
00449 else if(rollId.station()==3){
00450 if(cluSize==1*dupli) {hGlobalResClu1La[4]->Fill(minres);}
00451 else if(cluSize==2*dupli){ hGlobalResClu2La[4]->Fill(minres);}
00452 else if(cluSize==3*dupli){ hGlobalResClu3La[4]->Fill(minres);}
00453 }
00454 meIdRPC.str("");
00455 meIdRPC<<"RPCDataOccupancyFromDT_"<<rollId.rawId();
00456 meMap[meIdRPC.str()]->Fill(stripPredicted);
00457 }
00458 }
00459 }else{
00460 if(debug) std::cout<<"DT \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00461 }
00462 }else{
00463 if(debug) std::cout<<"DT \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00464 }
00465 }
00466 }
00467 }else {
00468 if(debug) std::cout<<"DT \t \t More than one segment in this chamber, or we are in Station 4"<<std::endl;
00469 }
00470 }
00471 } else {
00472 if(debug) std::cout<<"DT This Event doesn't have any DT4DDSegment"<<std::endl;
00473 }
00474 }
00475 }
00476
00477 if(incldtMB4){
00478 edm::Handle<DTRecSegment4DCollection> all4DSegments;
00479 iEvent.getByLabel(dt4DSegments, all4DSegments);
00480
00481 if(all4DSegments.isValid() && all4DSegments->size()>0){
00482
00483 std::map<DTChamberId,int> DTSegmentCounter;
00484 DTRecSegment4DCollection::const_iterator segment;
00485
00486 for (segment = all4DSegments->begin();segment!=all4DSegments->end(); ++segment){
00487 DTSegmentCounter[segment->chamberId()]++;
00488 }
00489
00490 if(debug) std::cout<<"MB4 \t \t Loop Over all4DSegments"<<std::endl;
00491 for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){
00492
00493 DTChamberId DTId = segment->chamberId();
00494
00495 if(DTSegmentCounter[DTId] == 1 && DTId.station()==4){
00496 int dtWheel = DTId.wheel();
00497 int dtStation = DTId.station();
00498 int dtSector = DTId.sector();
00499
00500 LocalPoint segmentPosition= segment->localPosition();
00501 LocalVector segmentDirection=segment->localDirection();
00502
00503
00504
00505 if(segment->dimension()==2){
00506 LocalVector segmentDirectionMB4=segmentDirection;
00507 LocalPoint segmentPositionMB4=segmentPosition;
00508
00509 bool compatiblesegments=false;
00510
00511 const BoundPlane& DTSurface4 = dtGeo->idToDet(DTId)->surface();
00512
00513 DTRecSegment4DCollection::const_iterator segMB3;
00514
00515 for(segMB3=all4DSegments->begin();segMB3!=all4DSegments->end();++segMB3){
00516
00517 DTChamberId dtid3 = segMB3->chamberId();
00518
00519 if(distsector_tmp(dtid3.sector(),DTId.sector())<=1
00520 && dtid3.station()==3
00521 && dtid3.wheel()==DTId.wheel()
00522 && DTSegmentCounter[dtid3] == 1
00523 && segMB3->dimension()==4){
00524
00525 const GeomDet* gdet3=dtGeo->idToDet(segMB3->geographicalId());
00526 const BoundPlane & DTSurface3 = gdet3->surface();
00527
00528 LocalVector segmentDirectionMB3 = segMB3->localDirection();
00529 GlobalPoint segmentPositionMB3inGlobal = DTSurface3.toGlobal(segMB3->localPosition());
00530
00531 LocalVector segDirMB4inMB3Frame=DTSurface3.toLocal(DTSurface4.toGlobal(segmentDirectionMB4));
00532 LocalVector segDirMB3inMB4Frame=DTSurface4.toLocal(DTSurface3.toGlobal(segmentDirectionMB3));
00533
00534 GlobalVector segDirMB4inGlobalFrame=DTSurface4.toGlobal(segmentDirectionMB4);
00535 GlobalVector segDirMB3inGlobalFrame=DTSurface3.toGlobal(segmentDirectionMB3);
00536
00537 float dx=segDirMB4inGlobalFrame.x();
00538 float dy=segDirMB4inGlobalFrame.y();
00539
00540
00541 float dx3=segDirMB3inGlobalFrame.x();
00542 float dy3=segDirMB3inGlobalFrame.y();
00543
00544
00545 double cosAng=fabs(dx*dx3+dy*dy3/sqrt((dx3*dx3+dy3*dy3)*(dx*dx+dy*dy)));
00546
00547 if(cosAng>MinCosAng){
00548 compatiblesegments=true;
00549 if(dtSector==13){
00550 dtSector=4;
00551 }
00552 if(dtSector==14){
00553 dtSector=10;
00554 }
00555
00556 std::set<RPCDetId> rollsForThisDT = rollstoreDT[DTStationIndex(0,dtWheel,dtSector,dtStation)];
00557
00558 assert(rollsForThisDT.size()>=1);
00559
00560 for (std::set<RPCDetId>::iterator iteraRoll=rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00561 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00562 RPCDetId rpcId = rollasociated->id();
00563 const BoundPlane & RPCSurfaceRB4 = rollasociated->surface();
00564
00565
00566
00567
00568 GlobalPoint CenterPointRollGlobal=RPCSurfaceRB4.toGlobal(LocalPoint(0,0,0));
00569 LocalPoint CenterRollinMB4Frame = DTSurface4.toLocal(CenterPointRollGlobal);
00570 LocalPoint segmentPositionMB3inMB4Frame = DTSurface4.toLocal(segmentPositionMB3inGlobal);
00571 LocalPoint segmentPositionMB3inRB4Frame = RPCSurfaceRB4.toLocal(segmentPositionMB3inGlobal);
00572 LocalVector segmentDirectionMB3inMB4Frame = DTSurface4.toLocal(segDirMB3inGlobalFrame);
00573
00574
00575 float Dxz=CenterRollinMB4Frame.z();
00576 float Xo4=segmentPositionMB4.x();
00577 float dxl=segmentDirectionMB4.x();
00578 float dzl=segmentDirectionMB4.z();
00579
00580 float X=Xo4+dxl*Dxz/dzl;
00581 float Z=Dxz;
00582
00583
00584 float Yo34=segmentPositionMB3inMB4Frame.y();
00585 float dy34 = segmentDirectionMB3inMB4Frame.y();
00586 float dz34 = segmentDirectionMB3inMB4Frame.z();
00587 float Dy=Dxz-(segmentPositionMB3inMB4Frame.z());
00588
00589 float Y=Yo34+dy34*Dy/dz34;
00590
00591 const RectangularStripTopology* top_
00592 =dynamic_cast<const RectangularStripTopology*>(&(rollasociated->topology()));
00593 LocalPoint xmin = top_->localPosition(0.);
00594 LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00595 float rsize = fabs( xmax.x()-xmin.x() );
00596 float stripl = top_->stripLength();
00597 float stripw = top_->pitch();
00598
00599 float extrapolatedDistance = sqrt((Y-Yo34)*(Y-Yo34)+Dy*Dy);
00600
00601 if(extrapolatedDistance<=MaxDrb4){
00602
00603 GlobalPoint GlobalPointExtrapolated = DTSurface4.toGlobal(LocalPoint(X,Y,Z));
00604 LocalPoint PointExtrapolatedRPCFrame = RPCSurfaceRB4.toLocal(GlobalPointExtrapolated);
00605
00606 if(fabs(PointExtrapolatedRPCFrame.z()) < 5. &&
00607 fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
00608 fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
00609
00610 RPCDetId rollId = rollasociated->id();
00611
00612
00613
00614
00615 const float stripPredicted=
00616 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
00617
00618 if(debug) std::cout<<"MB4 \t \t \t \t Candidate (from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
00619
00620
00621 std::map<std::string, MonitorElement*> meMap=meCollection[rollId.rawId()];
00622
00623 meIdDT.str("");
00624 meIdDT<<"ExpectedOccupancyFromDT_"<<rollId.rawId();
00625 meMap[meIdDT.str()]->Fill(stripPredicted);
00626
00627
00628
00629
00630 int countRecHits = 0;
00631 int cluSize = 0;
00632 float minres = 3000.;
00633
00634 if(debug) std::cout<<"MB4 \t \t \t \t Getting RecHits in Roll Asociated"<<std::endl;
00635 typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
00636 rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
00637 RPCRecHitCollection::const_iterator recHit;
00638
00639 for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
00640 countRecHits++;
00641 LocalPoint recHitPos=recHit->localPosition();
00642 float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
00643 if(debug) std::cout<<"DT \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction."<<std::endl;
00644 if(fabs(res)<fabs(minres)){
00645 minres=res;
00646 cluSize = recHit->clusterSize();
00647 }
00648 }
00649
00650 if(countRecHits==0){
00651 if(debug) std::cout <<"MB4 \t \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT"<<std::endl;
00652 }else{
00653 assert(minres!=3000);
00654
00655 if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
00656 assert(rollId.station()==4);
00657 if(cluSize==1*dupli){ hGlobalResClu1La[5]->Fill(minres);}
00658 else if(cluSize==2*dupli){ hGlobalResClu2La[5]->Fill(minres);}
00659 else if(cluSize==3*dupli){ hGlobalResClu3La[5]->Fill(minres);}
00660
00661 meIdRPC.str("");
00662 meIdRPC<<"RPCDataOccupancyFromDT_"<<rollId.rawId();
00663 meMap[meIdRPC.str()]->Fill(stripPredicted);
00664 }
00665 }
00666 }else{
00667 if(debug) std::cout<<"MB4 \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00668 }
00669 }
00670 else{
00671 if(debug) std::cout<<"MB4 \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00672 }
00673 }
00674 }else{
00675 compatiblesegments=false;
00676 if(debug) std::cout<<"MB4 \t \t \t \t I found segments in MB4 and MB3 adjacent or same wheel and sector but not compatibles Diferent Directions"<<std::endl;
00677 }
00678 }else{
00679 if(debug) std::cout<<"MB4 \t \t \t No the same station or same wheel or segment dim in mb3 not 4D"<<std::endl;
00680 }
00681 }
00682 }else{
00683 if(debug) std::cout<<"MB4 \t \t \t Is NOT a 2D Segment"<<std::endl;
00684 }
00685 }else{
00686 if(debug) std::cout<<"MB4 \t \t \t \t There is not just one segment or is not in station 4"<<std::endl;
00687 }
00688 }
00689 }else{
00690 if(debug) std::cout<<"MB4 \t This event doesn't have 4D Segment"<<std::endl;
00691 }
00692
00693 }
00694
00695
00696 if(inclcsc){
00697 if(debug) std::cout <<"\t Getting the CSC Segments"<<std::endl;
00698 edm::Handle<CSCSegmentCollection> allCSCSegments;
00699
00700 iEvent.getByLabel(cscSegments, allCSCSegments);
00701
00702 if(allCSCSegments.isValid()){
00703 if(allCSCSegments->size()>0){
00704 statistics->Fill(18);
00705
00706 if(debug) std::cout<<"CSC \t Number of CSC Segments in this event = "<<allCSCSegments->size()<<std::endl;
00707
00708 std::map<CSCDetId,int> CSCSegmentsCounter;
00709 CSCSegmentCollection::const_iterator segment;
00710
00711 int segmentsInThisEventInTheEndcap=0;
00712
00713 for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
00714 CSCSegmentsCounter[segment->cscDetId()]++;
00715 segmentsInThisEventInTheEndcap++;
00716 }
00717
00718 statistics->Fill(allCSCSegments->size()+18);
00719
00720 if(debug) std::cout<<"CSC \t loop over all the CSCSegments "<<std::endl;
00721 for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
00722 CSCDetId CSCId = segment->cscDetId();
00723
00724 if(CSCSegmentsCounter[CSCId]==1 && CSCId.station()!=4 && CSCId.ring()!=1 && allCSCSegments->size()>=2){
00725 if(debug) std::cout<<"CSC \t \t yes"<<std::endl;
00726 int cscEndCap = CSCId.endcap();
00727 int cscStation = CSCId.station();
00728 int cscRing = CSCId.ring();
00729
00730 int rpcRegion = 1; if(cscEndCap==2) rpcRegion= -1;
00731 int rpcRing = cscRing;
00732 if(cscRing==4)rpcRing =1;
00733 int rpcStation = cscStation;
00734 int rpcSegment = CSCId.chamber();
00735
00736 LocalPoint segmentPosition= segment->localPosition();
00737 LocalVector segmentDirection=segment->localDirection();
00738 float dz=segmentDirection.z();
00739
00740 if(debug) std::cout<<"CSC \t \t Is a good Segment? dim = 4, 4 <= nRecHits <= 10 Incident angle int range 45 < "<<acos(dz)*180/3.1415926<<" < 135? "<<std::endl;
00741
00742 if(segment->dimension()==4 && (segment->nRecHits()<=10 && segment->nRecHits()>=4)&& acos(dz)*180/3.1415926 > 45. && acos(dz)*180/3.1415926 < 160. ){
00743
00744 float Xo=segmentPosition.x();
00745 float Yo=segmentPosition.y();
00746 float Zo=segmentPosition.z();
00747 float dx=segmentDirection.x();
00748 float dy=segmentDirection.y();
00749 float dz=segmentDirection.z();
00750
00751
00752 if(debug) std::cout<<"CSC \t \t Getting chamber from Geometry"<<std::endl;
00753 const CSCChamber* TheChamber=cscGeo->chamber(CSCId);
00754 if(debug) std::cout<<"CSC \t \t Getting ID from Chamber"<<std::endl;
00755 const CSCDetId TheId=TheChamber->id();
00756 if(debug) std::cout<<"CSC \t \t Printing The Id"<<TheId<<std::endl;
00757 std::set<RPCDetId> rollsForThisCSC = rollstoreCSC[CSCStationIndex(rpcRegion,rpcStation,rpcRing,rpcSegment)];
00758 if(debug) std::cout<<"CSC \t \t Number of rolls for this CSC = "<<rollsForThisCSC.size()<<std::endl;
00759
00760 if(rpcRing!=1&&rpcStation!=4){
00761
00762
00763 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisCSC.begin();iteraRoll != rollsForThisCSC.end(); iteraRoll++){
00764
00765 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00766 RPCDetId rpcId = rollasociated->id();
00767
00768 const BoundPlane & RPCSurface = rollasociated->surface();
00769
00770 GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
00771 GlobalPoint CenterPointCSCGlobal = TheChamber->toGlobal(LocalPoint(0,0,0));
00772 GlobalPoint segmentPositionInGlobal=TheChamber->toGlobal(segmentPosition);
00773 LocalPoint CenterRollinCSCFrame = TheChamber->toLocal(CenterPointRollGlobal);
00774
00775 float D=CenterRollinCSCFrame.z();
00776
00777 float X=Xo+dx*D/dz;
00778 float Y=Yo+dy*D/dz;
00779 float Z=D;
00780
00781 const TrapezoidalStripTopology* top_=dynamic_cast<const TrapezoidalStripTopology*>(&(rollasociated->topology()));
00782 LocalPoint xmin = top_->localPosition(0.);
00783 if(debug) std::cout<<"CSC \t \t \t xmin of this Roll "<<xmin<<"cm"<<std::endl;
00784 LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00785 if(debug) std::cout<<"CSC \t \t \t xmax of this Roll "<<xmax<<"cm"<<std::endl;
00786 float rsize = fabs( xmax.x()-xmin.x() );
00787 if(debug) std::cout<<"CSC \t \t \t Roll Size "<<rsize<<"cm"<<std::endl;
00788 float stripl = top_->stripLength();
00789 float stripw = top_->pitch();
00790
00791
00792 float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
00793
00794
00795 if(extrapolatedDistance<=MaxD){
00796
00797 GlobalPoint GlobalPointExtrapolated=TheChamber->toGlobal(LocalPoint(X,Y,Z));
00798 LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
00799
00800
00801 if(fabs(PointExtrapolatedRPCFrame.z()) < 10. &&
00802 fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
00803 fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
00804
00805 RPCDetId rollId = rollasociated->id();
00806 RPCGeomServ rpcsrv(rollId);
00807 std::string nameRoll = rpcsrv.name();
00808
00809 if(debug) std::cout<<"CSC \t \t \t \t The RPCName is "<<nameRoll<<std::endl;
00810
00811 const float stripPredicted =
00812 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
00813
00814 if(debug) std::cout<<"CSC \t \t \t \t \t Candidate"<<rollId<<" "<<"(from CSC Segment) STRIP---> "<<stripPredicted<< std::endl;
00815
00816
00817 std::map<std::string, MonitorElement*> meMap=meCollection[rpcId.rawId()];
00818 meIdCSC.str("");
00819 meIdCSC<<"ExpectedOccupancyFromCSC_"<<rollId.rawId();
00820 meMap[meIdCSC.str()]->Fill(stripPredicted);
00821
00822
00823
00824
00825 int cluSize = 0;
00826 int countRecHits = 0;
00827 float minres = 3000.;
00828
00829 if(debug) std::cout<<"CSC \t \t \t \t \t Getting RecHits in Roll Asociated"<<std::endl;
00830 typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
00831 rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
00832 RPCRecHitCollection::const_iterator recHit;
00833
00834 for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
00835
00836 countRecHits++;
00837 LocalPoint recHitPos=recHit->localPosition();
00838 float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
00839 if(debug) std::cout<<"CSC \t \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction."<<std::endl;
00840 if(fabs(res)<fabs(minres)){
00841 minres=res;
00842 cluSize = recHit->clusterSize();
00843 if(debug) std::cout<<"CSC \t \t \t \t \t \t \t New Min Res "<<res<<"cm."<<std::endl;
00844 }
00845 }
00846
00847 if(countRecHits==0){
00848 if(debug) std::cout <<"CSC \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT"<<std::endl;
00849 }else{
00850 assert(minres!=3000);
00851
00852 if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
00853 if(debug) std::cout<<"CSC \t \t \t \t \t \t True!"<<std::endl;
00854
00855 if(rollId.ring()==2&&rollId.roll()==1){
00856 if(cluSize==1*dupli) hGlobalResClu1R2A->Fill(minres);
00857 else if(cluSize==2*dupli) hGlobalResClu2R2A->Fill(minres);
00858 else if(cluSize==3*dupli) hGlobalResClu3R2A->Fill(minres);
00859 }
00860 else if(rollId.ring()==2&&rollId.roll()==2){
00861 if(cluSize==1*dupli) hGlobalResClu1R2B->Fill(minres);
00862 else if(cluSize==2*dupli) hGlobalResClu2R2B->Fill(minres);
00863 else if(cluSize==3*dupli) hGlobalResClu3R2B->Fill(minres);
00864 }
00865 else if(rollId.ring()==2&&rollId.roll()==3){
00866 if(cluSize==1*dupli) hGlobalResClu1R2C->Fill(minres);
00867 else if(cluSize==2*dupli) hGlobalResClu2R2C->Fill(minres);
00868 else if(cluSize==3*dupli) hGlobalResClu3R2C->Fill(minres);
00869 }
00870 else if(rollId.ring()==3&&rollId.roll()==1){
00871 if(cluSize==1*dupli) hGlobalResClu1R3A->Fill(minres);
00872 else if(cluSize==2*dupli) hGlobalResClu2R3A->Fill(minres);
00873 else if(cluSize==3*dupli) hGlobalResClu3R3A->Fill(minres);
00874 }
00875 else if(rollId.ring()==3&&rollId.roll()==2){
00876 if(cluSize==1*dupli) hGlobalResClu1R3B->Fill(minres);
00877 else if(cluSize==2*dupli) hGlobalResClu2R3B->Fill(minres);
00878 else if(cluSize==3*dupli) hGlobalResClu3R3B->Fill(minres);
00879 }
00880 else if(rollId.ring()==3&&rollId.roll()==3){if(cluSize==1*dupli) hGlobalResClu1R3C->Fill(minres); if(cluSize==2*dupli) hGlobalResClu2R3C->Fill(minres); if(cluSize==3*dupli) hGlobalResClu3R3C->Fill(minres);
00881 }
00882 meIdRPC.str("");
00883 meIdRPC<<"RPCDataOccupancyFromCSC_"<<rollId.rawId();
00884 meMap[meIdRPC.str()]->Fill(stripPredicted);
00885 }
00886 }
00887
00888 }else{
00889 if(debug) std::cout<<"CSC \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00890 }
00891 }else{
00892 if(debug) std::cout<<"CSC \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00893 }
00894 }
00895 }
00896 }
00897 }else{
00898 if(debug) std::cout<<"CSC \t \t More than one segment in this chamber, or we are in Station Ring 1 or in Station 4"<<std::endl;
00899 }
00900 }
00901 }else{
00902 if(debug) std::cout<<"CSC This Event doesn't have any CSCSegment"<<std::endl;
00903 }
00904 }
00905 }
00906
00907 }
00908
00909
00910 void RPCEfficiency::endRun(const edm::Run& r, const edm::EventSetup& iSetup){
00911 if (EffSaveRootFile){
00912 dbe->save(EffRootFileName);
00913 }
00914 }
00915
00916
00917 void RPCEfficiency::endJob(){
00918 dbe =0;
00919 }
1.7.3