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
00510 const BoundPlane& DTSurface4 = dtGeo->idToDet(DTId)->surface();
00511
00512 DTRecSegment4DCollection::const_iterator segMB3;
00513
00514 for(segMB3=all4DSegments->begin();segMB3!=all4DSegments->end();++segMB3){
00515
00516 DTChamberId dtid3 = segMB3->chamberId();
00517
00518 if(distsector_tmp(dtid3.sector(),DTId.sector())<=1
00519 && dtid3.station()==3
00520 && dtid3.wheel()==DTId.wheel()
00521 && DTSegmentCounter[dtid3] == 1
00522 && segMB3->dimension()==4){
00523
00524 const GeomDet* gdet3=dtGeo->idToDet(segMB3->geographicalId());
00525 const BoundPlane & DTSurface3 = gdet3->surface();
00526
00527 LocalVector segmentDirectionMB3 = segMB3->localDirection();
00528 GlobalPoint segmentPositionMB3inGlobal = DTSurface3.toGlobal(segMB3->localPosition());
00529
00530 GlobalVector segDirMB4inGlobalFrame=DTSurface4.toGlobal(segmentDirectionMB4);
00531 GlobalVector segDirMB3inGlobalFrame=DTSurface3.toGlobal(segmentDirectionMB3);
00532
00533 float dx=segDirMB4inGlobalFrame.x();
00534 float dy=segDirMB4inGlobalFrame.y();
00535
00536
00537 float dx3=segDirMB3inGlobalFrame.x();
00538 float dy3=segDirMB3inGlobalFrame.y();
00539
00540
00541 double cosAng=fabs(dx*dx3+dy*dy3/sqrt((dx3*dx3+dy3*dy3)*(dx*dx+dy*dy)));
00542
00543 if(cosAng>MinCosAng){
00544 if(dtSector==13){
00545 dtSector=4;
00546 }
00547 if(dtSector==14){
00548 dtSector=10;
00549 }
00550
00551 std::set<RPCDetId> rollsForThisDT = rollstoreDT[DTStationIndex(0,dtWheel,dtSector,dtStation)];
00552
00553 assert(rollsForThisDT.size()>=1);
00554
00555 for (std::set<RPCDetId>::iterator iteraRoll=rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00556 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00557 const BoundPlane & RPCSurfaceRB4 = rollasociated->surface();
00558
00559
00560
00561
00562 GlobalPoint CenterPointRollGlobal=RPCSurfaceRB4.toGlobal(LocalPoint(0,0,0));
00563 LocalPoint CenterRollinMB4Frame = DTSurface4.toLocal(CenterPointRollGlobal);
00564 LocalPoint segmentPositionMB3inMB4Frame = DTSurface4.toLocal(segmentPositionMB3inGlobal);
00565 LocalVector segmentDirectionMB3inMB4Frame = DTSurface4.toLocal(segDirMB3inGlobalFrame);
00566
00567
00568 float Dxz=CenterRollinMB4Frame.z();
00569 float Xo4=segmentPositionMB4.x();
00570 float dxl=segmentDirectionMB4.x();
00571 float dzl=segmentDirectionMB4.z();
00572
00573 float X=Xo4+dxl*Dxz/dzl;
00574 float Z=Dxz;
00575
00576
00577 float Yo34=segmentPositionMB3inMB4Frame.y();
00578 float dy34 = segmentDirectionMB3inMB4Frame.y();
00579 float dz34 = segmentDirectionMB3inMB4Frame.z();
00580 float Dy=Dxz-(segmentPositionMB3inMB4Frame.z());
00581
00582 float Y=Yo34+dy34*Dy/dz34;
00583
00584 const RectangularStripTopology* top_
00585 =dynamic_cast<const RectangularStripTopology*>(&(rollasociated->topology()));
00586 LocalPoint xmin = top_->localPosition(0.);
00587 LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00588 float rsize = fabs( xmax.x()-xmin.x() );
00589 float stripl = top_->stripLength();
00590 float stripw = top_->pitch();
00591
00592 float extrapolatedDistance = sqrt((Y-Yo34)*(Y-Yo34)+Dy*Dy);
00593
00594 if(extrapolatedDistance<=MaxDrb4){
00595
00596 GlobalPoint GlobalPointExtrapolated = DTSurface4.toGlobal(LocalPoint(X,Y,Z));
00597 LocalPoint PointExtrapolatedRPCFrame = RPCSurfaceRB4.toLocal(GlobalPointExtrapolated);
00598
00599 if(fabs(PointExtrapolatedRPCFrame.z()) < 5. &&
00600 fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
00601 fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
00602
00603 RPCDetId rollId = rollasociated->id();
00604
00605
00606
00607
00608 const float stripPredicted=
00609 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
00610
00611 if(debug) std::cout<<"MB4 \t \t \t \t Candidate (from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
00612
00613
00614 std::map<std::string, MonitorElement*> meMap=meCollection[rollId.rawId()];
00615
00616 meIdDT.str("");
00617 meIdDT<<"ExpectedOccupancyFromDT_"<<rollId.rawId();
00618 meMap[meIdDT.str()]->Fill(stripPredicted);
00619
00620
00621
00622
00623 int countRecHits = 0;
00624 int cluSize = 0;
00625 float minres = 3000.;
00626
00627 if(debug) std::cout<<"MB4 \t \t \t \t Getting RecHits in Roll Asociated"<<std::endl;
00628 typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
00629 rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
00630 RPCRecHitCollection::const_iterator recHit;
00631
00632 for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
00633 countRecHits++;
00634 LocalPoint recHitPos=recHit->localPosition();
00635 float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
00636 if(debug) std::cout<<"DT \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction."<<std::endl;
00637 if(fabs(res)<fabs(minres)){
00638 minres=res;
00639 cluSize = recHit->clusterSize();
00640 }
00641 }
00642
00643 if(countRecHits==0){
00644 if(debug) std::cout <<"MB4 \t \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT"<<std::endl;
00645 }else{
00646 assert(minres!=3000);
00647
00648 if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
00649 assert(rollId.station()==4);
00650 if(cluSize==1*dupli){ hGlobalResClu1La[5]->Fill(minres);}
00651 else if(cluSize==2*dupli){ hGlobalResClu2La[5]->Fill(minres);}
00652 else if(cluSize==3*dupli){ hGlobalResClu3La[5]->Fill(minres);}
00653
00654 meIdRPC.str("");
00655 meIdRPC<<"RPCDataOccupancyFromDT_"<<rollId.rawId();
00656 meMap[meIdRPC.str()]->Fill(stripPredicted);
00657 }
00658 }
00659 }else{
00660 if(debug) std::cout<<"MB4 \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00661 }
00662 }
00663 else{
00664 if(debug) std::cout<<"MB4 \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00665 }
00666 }
00667 }else{
00668 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;
00669 }
00670 }else{
00671 if(debug) std::cout<<"MB4 \t \t \t No the same station or same wheel or segment dim in mb3 not 4D"<<std::endl;
00672 }
00673 }
00674 }else{
00675 if(debug) std::cout<<"MB4 \t \t \t Is NOT a 2D Segment"<<std::endl;
00676 }
00677 }else{
00678 if(debug) std::cout<<"MB4 \t \t \t \t There is not just one segment or is not in station 4"<<std::endl;
00679 }
00680 }
00681 }else{
00682 if(debug) std::cout<<"MB4 \t This event doesn't have 4D Segment"<<std::endl;
00683 }
00684
00685 }
00686
00687
00688 if(inclcsc){
00689 if(debug) std::cout <<"\t Getting the CSC Segments"<<std::endl;
00690 edm::Handle<CSCSegmentCollection> allCSCSegments;
00691
00692 iEvent.getByLabel(cscSegments, allCSCSegments);
00693
00694 if(allCSCSegments.isValid()){
00695 if(allCSCSegments->size()>0){
00696 statistics->Fill(18);
00697
00698 if(debug) std::cout<<"CSC \t Number of CSC Segments in this event = "<<allCSCSegments->size()<<std::endl;
00699
00700 std::map<CSCDetId,int> CSCSegmentsCounter;
00701 CSCSegmentCollection::const_iterator segment;
00702
00703 int segmentsInThisEventInTheEndcap=0;
00704
00705 for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
00706 CSCSegmentsCounter[segment->cscDetId()]++;
00707 segmentsInThisEventInTheEndcap++;
00708 }
00709
00710 statistics->Fill(allCSCSegments->size()+18);
00711
00712 if(debug) std::cout<<"CSC \t loop over all the CSCSegments "<<std::endl;
00713 for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
00714 CSCDetId CSCId = segment->cscDetId();
00715
00716 if(CSCSegmentsCounter[CSCId]==1 && CSCId.station()!=4 && CSCId.ring()!=1 && allCSCSegments->size()>=2){
00717 if(debug) std::cout<<"CSC \t \t yes"<<std::endl;
00718 int cscEndCap = CSCId.endcap();
00719 int cscStation = CSCId.station();
00720 int cscRing = CSCId.ring();
00721
00722 int rpcRegion = 1; if(cscEndCap==2) rpcRegion= -1;
00723 int rpcRing = cscRing;
00724 if(cscRing==4)rpcRing =1;
00725 int rpcStation = cscStation;
00726 int rpcSegment = CSCId.chamber();
00727
00728 LocalPoint segmentPosition= segment->localPosition();
00729 LocalVector segmentDirection=segment->localDirection();
00730 float dz=segmentDirection.z();
00731
00732 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;
00733
00734 if(segment->dimension()==4 && (segment->nRecHits()<=10 && segment->nRecHits()>=4)&& acos(dz)*180/3.1415926 > 45. && acos(dz)*180/3.1415926 < 160. ){
00735
00736 float Xo=segmentPosition.x();
00737 float Yo=segmentPosition.y();
00738 float Zo=segmentPosition.z();
00739 float dx=segmentDirection.x();
00740 float dy=segmentDirection.y();
00741 float dz=segmentDirection.z();
00742
00743
00744 if(debug) std::cout<<"CSC \t \t Getting chamber from Geometry"<<std::endl;
00745 const CSCChamber* TheChamber=cscGeo->chamber(CSCId);
00746 if(debug) std::cout<<"CSC \t \t Getting ID from Chamber"<<std::endl;
00747 const CSCDetId TheId=TheChamber->id();
00748 if(debug) std::cout<<"CSC \t \t Printing The Id"<<TheId<<std::endl;
00749 std::set<RPCDetId> rollsForThisCSC = rollstoreCSC[CSCStationIndex(rpcRegion,rpcStation,rpcRing,rpcSegment)];
00750 if(debug) std::cout<<"CSC \t \t Number of rolls for this CSC = "<<rollsForThisCSC.size()<<std::endl;
00751
00752 if(rpcRing!=1&&rpcStation!=4){
00753
00754
00755 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisCSC.begin();iteraRoll != rollsForThisCSC.end(); iteraRoll++){
00756
00757 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00758 RPCDetId rpcId = rollasociated->id();
00759
00760 const BoundPlane & RPCSurface = rollasociated->surface();
00761
00762 GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
00763 LocalPoint CenterRollinCSCFrame = TheChamber->toLocal(CenterPointRollGlobal);
00764
00765 float D=CenterRollinCSCFrame.z();
00766
00767 float X=Xo+dx*D/dz;
00768 float Y=Yo+dy*D/dz;
00769 float Z=D;
00770
00771 const TrapezoidalStripTopology* top_=dynamic_cast<const TrapezoidalStripTopology*>(&(rollasociated->topology()));
00772 LocalPoint xmin = top_->localPosition(0.);
00773 if(debug) std::cout<<"CSC \t \t \t xmin of this Roll "<<xmin<<"cm"<<std::endl;
00774 LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00775 if(debug) std::cout<<"CSC \t \t \t xmax of this Roll "<<xmax<<"cm"<<std::endl;
00776 float rsize = fabs( xmax.x()-xmin.x() );
00777 if(debug) std::cout<<"CSC \t \t \t Roll Size "<<rsize<<"cm"<<std::endl;
00778 float stripl = top_->stripLength();
00779 float stripw = top_->pitch();
00780
00781
00782 float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
00783
00784
00785 if(extrapolatedDistance<=MaxD){
00786
00787 GlobalPoint GlobalPointExtrapolated=TheChamber->toGlobal(LocalPoint(X,Y,Z));
00788 LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
00789
00790
00791 if(fabs(PointExtrapolatedRPCFrame.z()) < 10. &&
00792 fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
00793 fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
00794
00795 RPCDetId rollId = rollasociated->id();
00796 RPCGeomServ rpcsrv(rollId);
00797 std::string nameRoll = rpcsrv.name();
00798
00799 if(debug) std::cout<<"CSC \t \t \t \t The RPCName is "<<nameRoll<<std::endl;
00800
00801 const float stripPredicted =
00802 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
00803
00804 if(debug) std::cout<<"CSC \t \t \t \t \t Candidate"<<rollId<<" "<<"(from CSC Segment) STRIP---> "<<stripPredicted<< std::endl;
00805
00806
00807 std::map<std::string, MonitorElement*> meMap=meCollection[rpcId.rawId()];
00808 meIdCSC.str("");
00809 meIdCSC<<"ExpectedOccupancyFromCSC_"<<rollId.rawId();
00810 meMap[meIdCSC.str()]->Fill(stripPredicted);
00811
00812
00813
00814
00815 int cluSize = 0;
00816 int countRecHits = 0;
00817 float minres = 3000.;
00818
00819 if(debug) std::cout<<"CSC \t \t \t \t \t Getting RecHits in Roll Asociated"<<std::endl;
00820 typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
00821 rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
00822 RPCRecHitCollection::const_iterator recHit;
00823
00824 for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
00825
00826 countRecHits++;
00827 LocalPoint recHitPos=recHit->localPosition();
00828 float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
00829 if(debug) std::cout<<"CSC \t \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction."<<std::endl;
00830 if(fabs(res)<fabs(minres)){
00831 minres=res;
00832 cluSize = recHit->clusterSize();
00833 if(debug) std::cout<<"CSC \t \t \t \t \t \t \t New Min Res "<<res<<"cm."<<std::endl;
00834 }
00835 }
00836
00837 if(countRecHits==0){
00838 if(debug) std::cout <<"CSC \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT"<<std::endl;
00839 }else{
00840 assert(minres!=3000);
00841
00842 if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
00843 if(debug) std::cout<<"CSC \t \t \t \t \t \t True!"<<std::endl;
00844
00845 if(rollId.ring()==2&&rollId.roll()==1){
00846 if(cluSize==1*dupli) hGlobalResClu1R2A->Fill(minres);
00847 else if(cluSize==2*dupli) hGlobalResClu2R2A->Fill(minres);
00848 else if(cluSize==3*dupli) hGlobalResClu3R2A->Fill(minres);
00849 }
00850 else if(rollId.ring()==2&&rollId.roll()==2){
00851 if(cluSize==1*dupli) hGlobalResClu1R2B->Fill(minres);
00852 else if(cluSize==2*dupli) hGlobalResClu2R2B->Fill(minres);
00853 else if(cluSize==3*dupli) hGlobalResClu3R2B->Fill(minres);
00854 }
00855 else if(rollId.ring()==2&&rollId.roll()==3){
00856 if(cluSize==1*dupli) hGlobalResClu1R2C->Fill(minres);
00857 else if(cluSize==2*dupli) hGlobalResClu2R2C->Fill(minres);
00858 else if(cluSize==3*dupli) hGlobalResClu3R2C->Fill(minres);
00859 }
00860 else if(rollId.ring()==3&&rollId.roll()==1){
00861 if(cluSize==1*dupli) hGlobalResClu1R3A->Fill(minres);
00862 else if(cluSize==2*dupli) hGlobalResClu2R3A->Fill(minres);
00863 else if(cluSize==3*dupli) hGlobalResClu3R3A->Fill(minres);
00864 }
00865 else if(rollId.ring()==3&&rollId.roll()==2){
00866 if(cluSize==1*dupli) hGlobalResClu1R3B->Fill(minres);
00867 else if(cluSize==2*dupli) hGlobalResClu2R3B->Fill(minres);
00868 else if(cluSize==3*dupli) hGlobalResClu3R3B->Fill(minres);
00869 }
00870 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);
00871 }
00872 meIdRPC.str("");
00873 meIdRPC<<"RPCDataOccupancyFromCSC_"<<rollId.rawId();
00874 meMap[meIdRPC.str()]->Fill(stripPredicted);
00875 }
00876 }
00877
00878 }else{
00879 if(debug) std::cout<<"CSC \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00880 }
00881 }else{
00882 if(debug) std::cout<<"CSC \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00883 }
00884 }
00885 }
00886 }
00887 }else{
00888 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;
00889 }
00890 }
00891 }else{
00892 if(debug) std::cout<<"CSC This Event doesn't have any CSCSegment"<<std::endl;
00893 }
00894 }
00895 }
00896
00897 }
00898
00899
00900 void RPCEfficiency::endRun(const edm::Run& r, const edm::EventSetup& iSetup){
00901 if (EffSaveRootFile){
00902 dbe->save(EffRootFileName);
00903 }
00904 }
00905
00906
00907 void RPCEfficiency::endJob(){
00908 dbe =0;
00909 }
1.7.3