00001
00002
00003
00004
00005
00006
00007
00008 #include "DQM/RPCMonitorDigi/interface/RPCEfficiency.h"
00009
00010 #include <memory>
00011 #include "FWCore/Framework/interface/MakerMacros.h"
00012 #include <DataFormats/RPCDigi/interface/RPCDigiCollection.h>
00013 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
00014 #include <DataFormats/MuonDetId/interface/RPCDetId.h>
00015 #include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
00016 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
00017 #include <Geometry/RPCGeometry/interface/RPCGeomServ.h>
00018 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
00019 #include <Geometry/Records/interface/MuonGeometryRecord.h>
00020 #include <Geometry/CommonTopologies/interface/RectangularStripTopology.h>
00021 #include <Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h>
00022
00023 #include <cmath>
00024 #include "TFile.h"
00025 #include "TH1F.h"
00026 #include "TH2F.h"
00027 #include "TCanvas.h"
00028 #include "TAxis.h"
00029 #include "TString.h"
00030
00031 void RPCEfficiency::beginJob(){
00032
00033 }
00034
00035 int distsector(int sector1,int sector2){
00036
00037 if(sector1==13) sector1=4;
00038 if(sector1==14) sector1=10;
00039
00040 if(sector2==13) sector2=4;
00041 if(sector2==14) sector2=10;
00042
00043 int distance = abs(sector1 - sector2);
00044 if(distance>6) distance = 12-distance;
00045 return distance;
00046 }
00047
00048
00049 RPCEfficiency::RPCEfficiency(const edm::ParameterSet& iConfig){
00050 incldt=iConfig.getUntrackedParameter<bool>("incldt",true);
00051 incldtMB4=iConfig.getUntrackedParameter<bool>("incldtMB4",true);
00052 inclcsc=iConfig.getUntrackedParameter<bool>("inclcsc",true);
00053 debug=iConfig.getUntrackedParameter<bool>("debug",false);
00054
00055 rangestrips = iConfig.getUntrackedParameter<double>("rangestrips",4.);
00056 rangestripsRB4=iConfig.getUntrackedParameter<double>("rangestripsRB4",4.);
00057 dupli = iConfig.getUntrackedParameter<int>("DuplicationCorrection",2);
00058 MinCosAng=iConfig.getUntrackedParameter<double>("MinCosAng",0.96);
00059 MaxD=iConfig.getUntrackedParameter<double>("MaxD",80.);
00060 MaxDrb4=iConfig.getUntrackedParameter<double>("MaxDrb4",150.);
00061 muonRPCDigis=iConfig.getUntrackedParameter<std::string>("muonRPCDigis","muonRPCDigis");
00062 RPCRecHitLabel_ = iConfig.getParameter<edm::InputTag>("RecHitLabel");
00063 cscSegments=iConfig.getUntrackedParameter<std::string>("cscSegments","cscSegments");
00064 dt4DSegments=iConfig.getUntrackedParameter<std::string>("dt4DSegments","dt4DSegments");
00065
00066 folderPath=iConfig.getUntrackedParameter<std::string>("folderPath","RPC/RPCEfficiency/");
00067
00068 nameInLog = iConfig.getUntrackedParameter<std::string>("moduleLogName", "RPC_Eff");
00069 EffSaveRootFile = iConfig.getUntrackedParameter<bool>("EffSaveRootFile", false);
00070 EffRootFileName = iConfig.getUntrackedParameter<std::string>("EffRootFileName", "RPCEfficiency.root");
00071
00072
00073
00074 dbe = edm::Service<DQMStore>().operator->();
00075
00076 std::string folder = folderPath+"MuonSegEff";
00077 dbe->setCurrentFolder(folder);
00078 statistics = dbe->book1D("Statistics","All Statistics",33,0.5,33.5);
00079
00080 if(debug) std::cout<<"booking Global histograms with "<<folderPath<<std::endl;
00081
00082 folder = folderPath+"MuonSegEff/"+"Residuals/Barrel";
00083 dbe->setCurrentFolder(folder);
00084
00085
00086 hGlobalResClu1La1 = dbe->book1D("GlobalResidualsClu1La1","RPC Residuals Layer 1 Cluster Size 1",101,-10.,10.);
00087 hGlobalResClu1La2 = dbe->book1D("GlobalResidualsClu1La2","RPC Residuals Layer 2 Cluster Size 1",101,-10.,10.);
00088 hGlobalResClu1La3 = dbe->book1D("GlobalResidualsClu1La3","RPC Residuals Layer 3 Cluster Size 1",101,-10.,10.);
00089 hGlobalResClu1La4 = dbe->book1D("GlobalResidualsClu1La4","RPC Residuals Layer 4 Cluster Size 1",101,-10.,10.);
00090 hGlobalResClu1La5 = dbe->book1D("GlobalResidualsClu1La5","RPC Residuals Layer 5 Cluster Size 1",101,-10.,10.);
00091 hGlobalResClu1La6 = dbe->book1D("GlobalResidualsClu1La6","RPC Residuals Layer 6 Cluster Size 1",101,-10.,10.);
00092
00093 hGlobalResClu2La1 = dbe->book1D("GlobalResidualsClu2La1","RPC Residuals Layer 1 Cluster Size 2",101,-10.,10.);
00094 hGlobalResClu2La2 = dbe->book1D("GlobalResidualsClu2La2","RPC Residuals Layer 2 Cluster Size 2",101,-10.,10.);
00095 hGlobalResClu2La3 = dbe->book1D("GlobalResidualsClu2La3","RPC Residuals Layer 3 Cluster Size 2",101,-10.,10.);
00096 hGlobalResClu2La4 = dbe->book1D("GlobalResidualsClu2La4","RPC Residuals Layer 4 Cluster Size 2",101,-10.,10.);
00097 hGlobalResClu2La5 = dbe->book1D("GlobalResidualsClu2La5","RPC Residuals Layer 5 Cluster Size 2",101,-10.,10.);
00098 hGlobalResClu2La6 = dbe->book1D("GlobalResidualsClu2La6","RPC Residuals Layer 6 Cluster Size 2",101,-10.,10.);
00099
00100 hGlobalResClu3La1 = dbe->book1D("GlobalResidualsClu3La1","RPC Residuals Layer 1 Cluster Size 3",101,-10.,10.);
00101 hGlobalResClu3La2 = dbe->book1D("GlobalResidualsClu3La2","RPC Residuals Layer 2 Cluster Size 3",101,-10.,10.);
00102 hGlobalResClu3La3 = dbe->book1D("GlobalResidualsClu3La3","RPC Residuals Layer 3 Cluster Size 3",101,-10.,10.);
00103 hGlobalResClu3La4 = dbe->book1D("GlobalResidualsClu3La4","RPC Residuals Layer 4 Cluster Size 3",101,-10.,10.);
00104 hGlobalResClu3La5 = dbe->book1D("GlobalResidualsClu3La5","RPC Residuals Layer 5 Cluster Size 3",101,-10.,10.);
00105 hGlobalResClu3La6 = dbe->book1D("GlobalResidualsClu3La6","RPC Residuals Layer 6 Cluster Size 3",101,-10.,10.);
00106
00107 if(debug) std::cout<<"Booking Residuals for EndCap"<<std::endl;
00108 folder = folderPath+"MuonSegEff/Residuals/EndCap";
00109 dbe->setCurrentFolder(folder);
00110
00111
00112 hGlobalResClu1R3C = dbe->book1D("GlobalResidualsClu1R3C","RPC Residuals Ring 3 Roll C Cluster Size 1",101,-10.,10.);
00113 hGlobalResClu1R3B = dbe->book1D("GlobalResidualsClu1R3B","RPC Residuals Ring 3 Roll B Cluster Size 1",101,-10.,10.);
00114 hGlobalResClu1R3A = dbe->book1D("GlobalResidualsClu1R3A","RPC Residuals Ring 3 Roll A Cluster Size 1",101,-10.,10.);
00115 hGlobalResClu1R2C = dbe->book1D("GlobalResidualsClu1R2C","RPC Residuals Ring 2 Roll C Cluster Size 1",101,-10.,10.);
00116 hGlobalResClu1R2B = dbe->book1D("GlobalResidualsClu1R2B","RPC Residuals Ring 2 Roll B Cluster Size 1",101,-10.,10.);
00117 hGlobalResClu1R2A = dbe->book1D("GlobalResidualsClu1R2A","RPC Residuals Ring 2 Roll A Cluster Size 1",101,-10.,10.);
00118
00119 hGlobalResClu2R3C = dbe->book1D("GlobalResidualsClu2R3C","RPC Residuals Ring 3 Roll C Cluster Size 2",101,-10.,10.);
00120 hGlobalResClu2R3B = dbe->book1D("GlobalResidualsClu2R3B","RPC Residuals Ring 3 Roll B Cluster Size 2",101,-10.,10.);
00121 hGlobalResClu2R3A = dbe->book1D("GlobalResidualsClu2R3A","RPC Residuals Ring 3 Roll A Cluster Size 2",101,-10.,10.);
00122 hGlobalResClu2R2C = dbe->book1D("GlobalResidualsClu2R2C","RPC Residuals Ring 2 Roll C Cluster Size 2",101,-10.,10.);
00123 hGlobalResClu2R2B = dbe->book1D("GlobalResidualsClu2R2B","RPC Residuals Ring 2 Roll B Cluster Size 2",101,-10.,10.);
00124 hGlobalResClu2R2A = dbe->book1D("GlobalResidualsClu2R2A","RPC Residuals Ring 2 Roll A Cluster Size 2",101,-10.,10.);
00125
00126 hGlobalResClu3R3C = dbe->book1D("GlobalResidualsClu3R3C","RPC Residuals Ring 3 Roll C Cluster Size 3",101,-10.,10.);
00127 hGlobalResClu3R3B = dbe->book1D("GlobalResidualsClu3R3B","RPC Residuals Ring 3 Roll B Cluster Size 3",101,-10.,10.);
00128 hGlobalResClu3R3A = dbe->book1D("GlobalResidualsClu3R3A","RPC Residuals Ring 3 Roll A Cluster Size 3",101,-10.,10.);
00129 hGlobalResClu3R2C = dbe->book1D("GlobalResidualsClu3R2C","RPC Residuals Ring 2 Roll C Cluster Size 3",101,-10.,10.);
00130 hGlobalResClu3R2B = dbe->book1D("GlobalResidualsClu3R2B","RPC Residuals Ring 2 Roll B Cluster Size 3",101,-10.,10.);
00131 hGlobalResClu3R2A = dbe->book1D("GlobalResidualsClu3R2A","RPC Residuals Ring 2 Roll A Cluster Size 3",101,-10.,10.);
00132
00133
00134 if(debug) ofrej.open("rejected.txt");
00135
00136 if(debug) std::cout<<"Rejected done"<<std::endl;
00137
00138 }
00139
00140 void RPCEfficiency::beginRun(const edm::Run& run, const edm::EventSetup& iSetup){
00141
00142 iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00143 iSetup.get<MuonGeometryRecord>().get(dtGeo);
00144 iSetup.get<MuonGeometryRecord>().get(cscGeo);
00145
00146 for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00147 if(dynamic_cast< RPCChamber* >( *it ) != 0 ){
00148 RPCChamber* ch = dynamic_cast< RPCChamber* >( *it );
00149 std::vector< const RPCRoll*> roles = (ch->rolls());
00150 for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){
00151 RPCDetId rpcId = (*r)->id();
00152 int region=rpcId.region();
00153
00154 RPCGeomServ rpcsrv(rpcId);
00155 std::string nameRoll = rpcsrv.name();
00156 if(debug) std::cout<<"Booking for "<<rpcId.rawId()<<std::endl;
00157 meCollection[rpcId.rawId()] = bookDetUnitSeg(rpcId,(*r)->nstrips(),folderPath+"MuonSegEff/");
00158
00159 if(region==0&&(incldt||incldtMB4)){
00160
00161 int wheel=rpcId.ring();
00162 int sector=rpcId.sector();
00163 int station=rpcId.station();
00164 DTStationIndex ind(region,wheel,sector,station);
00165 std::set<RPCDetId> myrolls;
00166 if (rollstoreDT.find(ind)!=rollstoreDT.end()) myrolls=rollstoreDT[ind];
00167 myrolls.insert(rpcId);
00168 rollstoreDT[ind]=myrolls;
00169
00170 }
00171 if(region!=0 && inclcsc){
00172 int region=rpcId.region();
00173 int station=rpcId.station();
00174 int ring=rpcId.ring();
00175 int cscring=ring;
00176 int cscstation=station;
00177 RPCGeomServ rpcsrv(rpcId);
00178 int rpcsegment = rpcsrv.segment();
00179 int cscchamber = rpcsegment;
00180 if((station==2||station==3)&&ring==3){
00181 cscring = 2;
00182 }
00183
00184 CSCStationIndex ind(region,cscstation,cscring,cscchamber);
00185 std::set<RPCDetId> myrolls;
00186 if (rollstoreCSC.find(ind)!=rollstoreCSC.end()){
00187 myrolls=rollstoreCSC[ind];
00188 }
00189 myrolls.insert(rpcId);
00190 rollstoreCSC[ind]=myrolls;
00191 }
00192 }
00193 }
00194 }
00195 for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00196 if( dynamic_cast< RPCChamber* >( *it ) != 0 ){
00197
00198 RPCChamber* ch = dynamic_cast< RPCChamber* >( *it );
00199 std::vector< const RPCRoll*> roles = (ch->rolls());
00200 for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){
00201 RPCDetId rpcId = (*r)->id();
00202
00203 int region=rpcId.region();
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232 if(region!=0 && inclcsc && (rpcId.ring()==2 || rpcId.ring()==3)){
00233 int region=rpcId.region();
00234 int station=rpcId.station();
00235 int ring=rpcId.ring();
00236 int cscring = ring;
00237
00238 if((station==2||station==3)&&ring==3) cscring = 2;
00239
00240
00241 int cscstation=station;
00242 RPCGeomServ rpcsrv(rpcId);
00243 int rpcsegment = rpcsrv.segment();
00244
00245
00246 int cscchamber = rpcsegment+1;
00247 if(cscchamber==37)cscchamber=1;
00248 CSCStationIndex ind(region,cscstation,cscring,cscchamber);
00249 std::set<RPCDetId> myrolls;
00250 if (rollstoreCSC.find(ind)!=rollstoreCSC.end())myrolls=rollstoreCSC[ind];
00251 myrolls.insert(rpcId);
00252 rollstoreCSC[ind]=myrolls;
00253
00254 cscchamber = rpcsegment-1;
00255 if(cscchamber==0)cscchamber=36;
00256 CSCStationIndex indDos(region,cscstation,cscring,cscchamber);
00257 std::set<RPCDetId> myrollsDos;
00258 if (rollstoreCSC.find(indDos)!=rollstoreCSC.end()) myrollsDos=rollstoreCSC[indDos];
00259 myrollsDos.insert(rpcId);
00260 rollstoreCSC[indDos]=myrollsDos;
00261 }
00262 }
00263 }
00264 }
00265 }
00266
00267 RPCEfficiency::~RPCEfficiency()
00268 {
00269
00270 }
00271
00272 void RPCEfficiency::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00273 {
00274
00275 statistics->Fill(1);
00276
00277 char meIdRPC [128];
00278 char meIdDT [128];
00279 char meIdCSC [128];
00280
00281 if(debug) std::cout <<"\t Getting the RPC RecHits"<<std::endl;
00282 edm::Handle<RPCRecHitCollection> rpcHits;
00283 iEvent.getByLabel(RPCRecHitLabel_,rpcHits);
00284
00285
00286 if(rpcHits.isValid()){
00287 if(incldt){
00288 if(debug) std::cout<<"\t Getting the DT Segments"<<std::endl;
00289 edm::Handle<DTRecSegment4DCollection> all4DSegments;
00290
00291 iEvent.getByLabel(dt4DSegments, all4DSegments);
00292
00293 if(all4DSegments.isValid()){
00294 if(all4DSegments->size()>0){
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(debug) std::cout<<"DT \t \t This Segment is in Chamber id: "<<DTId<<std::endl;
00315 if(debug) std::cout<<"DT \t \t Number of segments in this DT = "<<DTSegmentCounter[DTId]<<std::endl;
00316 if(debug) std::cout<<"DT \t \t Is the only one in this DT? and is not in the 4th Station?"<<std::endl;
00317
00318
00319 if(DTSegmentCounter[DTId]==1 && DTId.station()!=4){
00320
00321 int dtWheel = DTId.wheel();
00322 int dtStation = DTId.station();
00323 int dtSector = DTId.sector();
00324
00325 LocalPoint segmentPosition= segment->localPosition();
00326 LocalVector segmentDirection=segment->localDirection();
00327
00328 const GeomDet* gdet=dtGeo->idToDet(segment->geographicalId());
00329 const BoundPlane & DTSurface = gdet->surface();
00330
00331
00332
00333 if(debug) std::cout<<"DT \t \t Is the segment 4D?"<<std::endl;
00334
00335 if(segment->dimension()==4){
00336
00337 if(debug) std::cout<<"DT \t \t yes"<<std::endl;
00338 if(debug) std::cout<<"DT \t \t DT Segment Dimension "<<segment->dimension()<<std::endl;
00339
00340 float Xo=segmentPosition.x();
00341 float Yo=segmentPosition.y();
00342 float Zo=segmentPosition.z();
00343 float dx=segmentDirection.x();
00344 float dy=segmentDirection.y();
00345 float dz=segmentDirection.z();
00346
00347 std::set<RPCDetId> rollsForThisDT = rollstoreDT[DTStationIndex(0,dtWheel,dtSector,dtStation)];
00348
00349 if(debug) std::cout<<"DT \t \t Number of rolls for this DT = "<<rollsForThisDT.size()<<std::endl;
00350
00351 assert(rollsForThisDT.size()>=1);
00352
00353 if(debug) std::cout<<"DT \t \t Loop over all the rolls asociated to this DT"<<std::endl;
00354 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00355 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00356 RPCDetId rpcId = rollasociated->id();
00357 const BoundPlane & RPCSurface = rollasociated->surface();
00358
00359 RPCGeomServ rpcsrv(rpcId);
00360 std::string nameRoll = rpcsrv.name();
00361
00362 if(debug) std::cout<<"DT \t \t \t RollName: "<<nameRoll<<std::endl;
00363 if(debug) std::cout<<"DT \t \t \t Doing the extrapolation to this roll"<<std::endl;
00364 if(debug) std::cout<<"DT \t \t \t DT Segment Direction in DTLocal "<<segmentDirection<<std::endl;
00365 if(debug) std::cout<<"DT \t \t \t DT Segment Point in DTLocal "<<segmentPosition<<std::endl;
00366
00367 GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
00368
00369 LocalPoint CenterRollinDTFrame = DTSurface.toLocal(CenterPointRollGlobal);
00370
00371 if(debug) std::cout<<"DT \t \t \t Center (0,0,0) Roll In DTLocal"<<CenterRollinDTFrame<<std::endl;
00372 if(debug) std::cout<<"DT \t \t \t Center (0,0,0) of the Roll in Global"<<CenterPointRollGlobal<<std::endl;
00373
00374 float D=CenterRollinDTFrame.z();
00375
00376 float X=Xo+dx*D/dz;
00377 float Y=Yo+dy*D/dz;
00378 float Z=D;
00379
00380 const RectangularStripTopology* top_= dynamic_cast<const RectangularStripTopology*> (&(rollasociated->topology()));
00381 LocalPoint xmin = top_->localPosition(0.);
00382 if(debug) std::cout<<"DT \t \t \t xmin of this Roll "<<xmin<<"cm"<<std::endl;
00383 LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00384 if(debug) std::cout<<"DT \t \t \t xmax of this Roll "<<xmax<<"cm"<<std::endl;
00385 float rsize = fabs( xmax.x()-xmin.x() );
00386 if(debug) std::cout<<"DT \t \t \t Roll Size "<<rsize<<"cm"<<std::endl;
00387 float stripl = top_->stripLength();
00388 float stripw = top_->pitch();
00389
00390 if(debug) std::cout<<"DT \t \t \t Strip Lenght "<<stripl<<"cm"<<std::endl;
00391 if(debug) std::cout<<"DT \t \t \t Strip Width "<<stripw<<"cm"<<std::endl;
00392 if(debug) std::cout<<"DT \t \t \t X Predicted in DTLocal= "<<X<<"cm"<<std::endl;
00393 if(debug) std::cout<<"DT \t \t \t Y Predicted in DTLocal= "<<Y<<"cm"<<std::endl;
00394 if(debug) std::cout<<"DT \t \t \t Z Predicted in DTLocal= "<<Z<<"cm"<<std::endl;
00395
00396 float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
00397
00398 if(debug) std::cout<<"DT \t \t \t Is the distance of extrapolation less than MaxD? ="<<extrapolatedDistance<<"cm"<<"MaxD="<<MaxD<<"cm"<<std::endl;
00399
00400 if(extrapolatedDistance<=MaxD){
00401 if(debug) std::cout<<"DT \t \t \t yes"<<std::endl;
00402 GlobalPoint GlobalPointExtrapolated = DTSurface.toGlobal(LocalPoint(X,Y,Z));
00403 if(debug) std::cout<<"DT \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
00404 LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
00405
00406 if(debug) std::cout<<"DT \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
00407 if(debug) std::cout<<"DT \t \t \t Corner of the Roll = ("<<rsize*0.5<<","<<stripl*0.5<<")"<<std::endl;
00408 if(debug) std::cout<<"DT \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.x())<<","
00409 <<fabs(PointExtrapolatedRPCFrame.y())<<","<<fabs(PointExtrapolatedRPCFrame.z())<<")"<<std::endl;
00410 if(debug) std::cout<<"DT \t \t \t Does the extrapolation go inside this roll?"<<std::endl;
00411
00412 if(fabs(PointExtrapolatedRPCFrame.z()) < 10. &&
00413 fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
00414 fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
00415
00416 if(debug) std::cout<<"DT \t \t \t \t yes"<<std::endl;
00417
00418 RPCDetId rollId = rollasociated->id();
00419
00420 RPCGeomServ rpcsrv(rollId);
00421 std::string nameRoll = rpcsrv.name();
00422 if(debug) std::cout<<"DT \t \t \t \t The RPCName is "<<nameRoll<<std::endl;
00423 const float stripPredicted =
00424 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
00425
00426 if(debug) std::cout<<"DT \t \t \t \t Candidate (from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
00427
00428
00429 std::map<std::string, MonitorElement*> meMap=meCollection[rpcId.rawId()];
00430
00431 sprintf(meIdDT,"ExpectedOccupancyFromDT_%d",rollId.rawId());
00432 if(debug) std::cout<<"DT \t \t \t \t Filling Expected for "<<meIdDT<<" with "<<stripPredicted<<std::endl;
00433 if(fabs(stripPredicted-rollasociated->nstrips())<1.) if(debug) std::cout<<"DT \t \t \t \t Extrapolating near last strip, Event"<<iEvent.id()<<" stripPredicted="<<stripPredicted<<" Number of strips="<<rollasociated->nstrips()<<std::endl;
00434 if(fabs(stripPredicted)<1.) if(debug) std::cout<<"DT \t \t \t \t Extrapolating near first strip, Event"<<iEvent.id()<<" stripPredicted="<<stripPredicted<<" Number of strips="<<rollasociated->nstrips()<<std::endl;
00435 meMap[meIdDT]->Fill(stripPredicted);
00436
00437
00438
00439
00440 int countRecHits = 0;
00441 int cluSize = 0;
00442 float minres = 3000.;
00443
00444 if(debug) std::cout<<"DT \t \t \t \t Getting RecHits in Roll Asociated"<<std::endl;
00445 typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
00446 rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
00447 RPCRecHitCollection::const_iterator recHit;
00448
00449 for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
00450 countRecHits++;
00451
00452 sprintf(meIdRPC,"BXDistribution_%d",rollasociated->id().rawId());
00453 meMap[meIdRPC]->Fill(recHit->BunchX());
00454
00455 LocalPoint recHitPos=recHit->localPosition();
00456 float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
00457 if(debug) std::cout<<"DT \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction."<<std::endl;
00458 if(fabs(res)<fabs(minres)){
00459 minres=res;
00460 cluSize = recHit->clusterSize();
00461 if(debug) std::cout<<"DT \t \t \t \t \t \t New Min Res "<<res<<"cm."<<std::endl;
00462 }
00463 }
00464
00465 bool anycoincidence=false;
00466
00467 if(countRecHits==0){
00468 if(debug) std::cout <<"DT \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT"<<std::endl;
00469 }else{
00470 assert(minres!=3000);
00471
00472 if(debug) std::cout<<"DT \t \t \t \t \t PointExtrapolatedRPCFrame.x="<<PointExtrapolatedRPCFrame.x()<<" Minimal Residual="<<minres<<std::endl;
00473 if(debug) std::cout<<"DT \t \t \t \t \t Minimal Residual less than stripw*rangestrips? minres="<<minres<<" range="<<rangestrips<<" stripw="<<stripw<<" cluSize"<<cluSize<<" <=compare minres with"<<(rangestrips+cluSize*0.5)*stripw<<std::endl;
00474 if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
00475 if(debug) std::cout<<"DT \t \t \t \t \t \t True!"<<std::endl;
00476 anycoincidence=true;
00477 }
00478 }
00479 if(anycoincidence){
00480 if(debug) std::cout<<"DT \t \t \t \t \t At least one RecHit inside the range, Predicted="<<stripPredicted<<" minres="<<minres<<"cm range="<<rangestrips<<"strips stripw="<<stripw<<"cm"<<std::endl;
00481 if(debug) std::cout<<"DT \t \t \t \t \t Norm of Cosine Directors="<<dx*dx+dy*dy+dz*dz<<"~1?"<<std::endl;
00482
00483 float cosal = dx/sqrt(dx*dx+dz*dz);
00484 if(debug) std::cout<<"DT \t \t \t \t \t Angle="<<acos(cosal)*180/3.1415926<<" degree"<<std::endl;
00485 if(debug) std::cout<<"DT \t \t \t \t \t Filling the Residuals Histogram for globals with "<<minres<<"And the angular incidence with Cos Alpha="<<cosal<<std::endl;
00486 if(rollId.station()==1&&rollId.layer()==1) { if(cluSize==1*dupli) {hGlobalResClu1La1->Fill(minres);}if(cluSize==2*dupli){ hGlobalResClu2La1->Fill(minres);} else if(cluSize==3*dupli){ hGlobalResClu3La1->Fill(minres);}}
00487 else if(rollId.station()==1&&rollId.layer()==2){ if(cluSize==1*dupli) {hGlobalResClu1La2->Fill(minres);}if(cluSize==2*dupli){ hGlobalResClu2La2->Fill(minres);} else if(cluSize==3*dupli){ hGlobalResClu3La2->Fill(minres);}}
00488 else if(rollId.station()==2&&rollId.layer()==1){ if(cluSize==1*dupli) {hGlobalResClu1La3->Fill(minres);}if(cluSize==2*dupli){ hGlobalResClu2La3->Fill(minres);} else if(cluSize==3*dupli){ hGlobalResClu3La3->Fill(minres);}}
00489 else if(rollId.station()==2&&rollId.layer()==2){ if(cluSize==1*dupli) {hGlobalResClu1La4->Fill(minres);}if(cluSize==2*dupli){ hGlobalResClu2La4->Fill(minres);} else if(cluSize==3*dupli){ hGlobalResClu3La4->Fill(minres);}}
00490 else if(rollId.station()==3) { if(cluSize==1*dupli) {hGlobalResClu1La5->Fill(minres);}if(cluSize==2*dupli){ hGlobalResClu2La5->Fill(minres);} else if(cluSize==3*dupli){ hGlobalResClu3La5->Fill(minres);}}
00491
00492 sprintf(meIdRPC,"RPCDataOccupancyFromDT_%d",rollId.rawId());
00493 if(debug) std::cout<<"DT \t \t \t \t \t COINCIDENCE!!! Event="<<iEvent.id()<<" Filling RPC Data Occupancy for "<<meIdRPC<<" with "<<stripPredicted<<std::endl;
00494 meMap[meIdRPC]->Fill(stripPredicted);
00495 }
00496 else{
00497 RPCGeomServ rpcsrv(rollasociated->id());
00498 std::string nameRoll = rpcsrv.name();
00499 if(debug) std::cout<<"DT \t \t \t \t \t A roll was ineficient in event "<<iEvent.id().event()<<std::endl;
00500 if(debug) ofrej<<"DTs \t Wh "<<dtWheel
00501 <<"\t St "<<dtStation
00502 <<"\t Se "<<dtSector
00503 <<"\t Roll "<<nameRoll
00504 <<"\t Event "
00505 <<iEvent.id().event()
00506 <<"\t Run "
00507 <<iEvent.id().run()
00508 <<std::endl;
00509 }
00510 }else{
00511 if(debug) std::cout<<"DT \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00512 }
00513 }else{
00514 if(debug) std::cout<<"DT \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00515 }
00516 }
00517 }
00518 }else {
00519 if(debug) std::cout<<"DT \t \t More than one segment in this chamber, or we are in Station 4"<<std::endl;
00520 }
00521 }
00522 } else {
00523 if(debug) std::cout<<"DT This Event doesn't have any DT4DDSegment"<<std::endl;
00524 }
00525 }
00526 }
00527
00528 if(incldtMB4){
00529 if(debug) std::cout <<"MB4 \t Getting ALL the DT Segments"<<std::endl;
00530 edm::Handle<DTRecSegment4DCollection> all4DSegments;
00531 iEvent.getByLabel(dt4DSegments, all4DSegments);
00532
00533 iEvent.getByLabel(dt4DSegments, all4DSegments);
00534
00535 if(all4DSegments.isValid()){
00536 if(all4DSegments->size()>0){
00537 std::map<DTChamberId,int> DTSegmentCounter;
00538 DTRecSegment4DCollection::const_iterator segment;
00539
00540 for (segment = all4DSegments->begin();segment!=all4DSegments->end(); ++segment){
00541 DTSegmentCounter[segment->chamberId()]++;
00542 }
00543
00544 if(debug) std::cout<<"MB4 \t \t Loop Over all4DSegments"<<std::endl;
00545 for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){
00546
00547 DTChamberId DTId = segment->chamberId();
00548
00549 if(debug) std::cout<<"MB4 \t \t This Segment is in Chamber id: "<<DTId<<std::endl;
00550 if(debug) std::cout<<"MB4 \t \t Number of segments in this DT = "<<DTSegmentCounter[DTId]<<std::endl;
00551 if(debug) std::cout<<"MB4 \t \t \t Is the only one in this DT? and is in the Station 4?"<<std::endl;
00552
00553 if(DTSegmentCounter[DTId] == 1 && DTId.station()==4){
00554
00555 if(debug) std::cout<<"MB4 \t \t \t yes"<<std::endl;
00556 int dtWheel = DTId.wheel();
00557 int dtStation = DTId.station();
00558 int dtSector = DTId.sector();
00559
00560 LocalPoint segmentPosition= segment->localPosition();
00561 LocalVector segmentDirection=segment->localDirection();
00562
00563
00564
00565
00566 if(debug) std::cout<<"MB4 \t \t \t \t The Segment in MB4 is 2D?"<<std::endl;
00567 if(segment->dimension()==2){
00568 if(debug) std::cout<<"MB4 \t \t \t \t yes"<<std::endl;
00569 LocalVector segmentDirectionMB4=segmentDirection;
00570 LocalPoint segmentPositionMB4=segmentPosition;
00571
00572 bool compatiblesegments=false;
00573
00574 const BoundPlane& DTSurface4 = dtGeo->idToDet(DTId)->surface();
00575
00576 DTRecSegment4DCollection::const_iterator segMB3;
00577
00578 if(debug) std::cout<<"MB4 \t \t \t \t Loop on segments in =sector && MB3 && adjacent sectors && y dim=4"<<std::endl;
00579 for(segMB3=all4DSegments->begin();segMB3!=all4DSegments->end();++segMB3){
00580
00581 DTChamberId dtid3 = segMB3->chamberId();
00582
00583 if(distsector(dtid3.sector(),DTId.sector())<=1
00584 && dtid3.station()==3
00585 && dtid3.wheel()==DTId.wheel()
00586 && DTSegmentCounter[dtid3] == 1
00587 && segMB3->dimension()==4){
00588
00589 if(debug) std::cout<<"MB4 \t \t \t \t distsector ="<<distsector(dtid3.sector(),DTId.sector())<<std::endl;
00590
00591 const GeomDet* gdet3=dtGeo->idToDet(segMB3->geographicalId());
00592 const BoundPlane & DTSurface3 = gdet3->surface();
00593
00594 LocalVector segmentDirectionMB3 = segMB3->localDirection();
00595 GlobalPoint segmentPositionMB3inGlobal = DTSurface3.toGlobal(segMB3->localPosition());
00596
00597
00598 LocalVector segDirMB4inMB3Frame=DTSurface3.toLocal(DTSurface4.toGlobal(segmentDirectionMB4));
00599 LocalVector segDirMB3inMB4Frame=DTSurface4.toLocal(DTSurface3.toGlobal(segmentDirectionMB3));
00600
00601 GlobalVector segDirMB4inGlobalFrame=DTSurface4.toGlobal(segmentDirectionMB4);
00602 GlobalVector segDirMB3inGlobalFrame=DTSurface3.toGlobal(segmentDirectionMB3);
00603
00604 float dx=segDirMB4inGlobalFrame.x();
00605 float dy=segDirMB4inGlobalFrame.y();
00606 float dz=segDirMB4inGlobalFrame.z();
00607
00608 float dx3=segDirMB3inGlobalFrame.x();
00609 float dy3=segDirMB3inGlobalFrame.y();
00610 float dz3=segDirMB3inGlobalFrame.z();
00611
00612 double cosAng=fabs(dx*dx3+dy*dy3/sqrt((dx3*dx3+dy3*dy3)*(dx*dx+dy*dy)));
00613
00614 if(debug) std::cout<<"MB4 \t \t \t \t cosAng"<<cosAng<<"Beetween "<<dtid3<<" and "<<DTId<<std::endl;
00615
00616 if(fabs(cosAng)>1.){
00617 if(debug) std::cout<<"dx="<<dx<<" dz="<<dz<<std::endl;
00618 if(debug) std::cout<<"dx3="<<dx3<<" dz3="<<dz<<std::endl;
00619 if(debug) std::cout<<cosAng<<std::endl;
00620 }
00621
00622 if(cosAng>MinCosAng){
00623 compatiblesegments=true;
00624 if(dtSector==13){
00625 dtSector=4;
00626 }
00627 if(dtSector==14){
00628 dtSector=10;
00629 }
00630
00631 std::set<RPCDetId> rollsForThisDT = rollstoreDT[DTStationIndex(0,dtWheel,dtSector,dtStation)];
00632
00633 if(debug) std::cout<<"MB4 \t \t Number of rolls for this DT = "<<rollsForThisDT.size()<<std::endl;
00634
00635 assert(rollsForThisDT.size()>=1);
00636
00637 if(debug) std::cout<<"MB4 \t \t Loop over all the rolls asociated to this DT"<<std::endl;
00638
00639 for (std::set<RPCDetId>::iterator iteraRoll=rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00640 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00641 RPCDetId rpcId = rollasociated->id();
00642 const BoundPlane & RPCSurfaceRB4 = rollasociated->surface();
00643
00644 RPCGeomServ rpcsrv(rpcId);
00645 std::string nameRoll = rpcsrv.name();
00646
00647 if(debug) std::cout<<"MB4 \t \t \t RollName: "<<nameRoll<<std::endl;
00648 if(debug) std::cout<<"MB4 \t \t \t Doing the extrapolation to this roll"<<std::endl;
00649
00650 GlobalPoint CenterPointRollGlobal=RPCSurfaceRB4.toGlobal(LocalPoint(0,0,0));
00651 LocalPoint CenterRollinMB4Frame = DTSurface4.toLocal(CenterPointRollGlobal);
00652 LocalPoint segmentPositionMB3inMB4Frame = DTSurface4.toLocal(segmentPositionMB3inGlobal);
00653 LocalPoint segmentPositionMB3inRB4Frame = RPCSurfaceRB4.toLocal(segmentPositionMB3inGlobal);
00654 LocalVector segmentDirectionMB3inMB4Frame = DTSurface4.toLocal(segDirMB3inGlobalFrame);
00655
00656
00657 float Dxz=CenterRollinMB4Frame.z();
00658 float Xo4=segmentPositionMB4.x();
00659 float dxl=segmentDirectionMB4.x();
00660 float dzl=segmentDirectionMB4.z();
00661
00662 float X=Xo4+dxl*Dxz/dzl;
00663 float Z=Dxz;
00664
00665
00666 float Yo34=segmentPositionMB3inMB4Frame.y();
00667 float dy34 = segmentDirectionMB3inMB4Frame.y();
00668 float dz34 = segmentDirectionMB3inMB4Frame.z();
00669 float Dy=Dxz-(segmentPositionMB3inMB4Frame.z());
00670
00671 if(debug) std::cout<<"MB4 \t \t \t The distance to extrapolate in Y from MB3 is "<<Dy<<"cm"<<std::endl;
00672
00673 float Y=Yo34+dy34*Dy/dz34;
00674
00675 const RectangularStripTopology* top_
00676 =dynamic_cast<const RectangularStripTopology*>(&(rollasociated->topology()));
00677 LocalPoint xmin = top_->localPosition(0.);
00678 LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00679 float rsize = fabs( xmax.x()-xmin.x() );
00680 float stripl = top_->stripLength();
00681 float stripw = top_->pitch();
00682
00683
00684 if(debug) std::cout<<"MB4 \t \t \t Strip Lenght "<<stripl<<"cm"<<std::endl;
00685 if(debug) std::cout<<"MB4 \t \t \t Strip Width "<<stripw<<"cm"<<std::endl;
00686
00687 if(debug) std::cout<<"MB4 \t \t \t X Predicted in MB4DTLocal= "<<X<<"cm"<<std::endl;
00688 if(debug) std::cout<<"MB4 \t \t \t Y Predicted in MB4DTLocal= "<<Y<<"cm"<<std::endl;
00689 if(debug) std::cout<<"MB4 \t \t \t Z Predicted in MB4DTLocal= "<<Z<<"cm"<<std::endl;
00690
00691 float extrapolatedDistance = sqrt((Y-Yo34)*(Y-Yo34)+Dy*Dy);
00692
00693 if(debug) std::cout<<"MB4 \t \t \t segmentPositionMB3inMB4Frame"<<segmentPositionMB3inMB4Frame<<std::endl;
00694 if(debug) std::cout<<"MB4 \t \t \t segmentPositionMB4inMB4Frame"<<segmentPosition<<std::endl;
00695
00696 if(debug) std::cout<<"MB4 \t \t \t segmentDirMB3inMB4Frame"<<segDirMB3inMB4Frame<<std::endl;
00697 if(debug) std::cout<<"MB4 \t \t \t segmentDirMB4inMB4Frame"<<segmentDirectionMB4<<std::endl;
00698
00699 if(debug) std::cout<<"MB4 \t \t \t CenterRB4PositioninMB4Frame"<<CenterRollinMB4Frame<<std::endl;
00700
00701 if(debug) std::cout<<"MB4 \t \t \t Is the extrapolation distance ="<<extrapolatedDistance<<"less than "<<MaxDrb4<<std::endl;
00702
00703
00704 if(extrapolatedDistance<=MaxDrb4){
00705 if(debug) std::cout<<"MB4 \t \t \t yes"<<std::endl;
00706
00707 GlobalPoint GlobalPointExtrapolated = DTSurface4.toGlobal(LocalPoint(X,Y,Z));
00708
00709 if(debug) std::cout<<"MB4 \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
00710
00711 LocalPoint PointExtrapolatedRPCFrame = RPCSurfaceRB4.toLocal(GlobalPointExtrapolated);
00712
00713 if(debug) std::cout<<"MB4 \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
00714 if(debug) std::cout<<"MB4 \t \t \t Corner of the Roll = ("<<rsize*0.5<<","<<stripl*0.5<<")"<<std::endl;
00715 if(debug) std::cout<<"MB4 \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.x())<<","
00716 <<fabs(PointExtrapolatedRPCFrame.y())<<","<<fabs(PointExtrapolatedRPCFrame.z())<<")"<<std::endl;
00717
00718 if(debug) std::cout<<"MB4 \t \t \t Does the extrapolation go inside this roll?"<<std::endl;
00719
00720 if(fabs(PointExtrapolatedRPCFrame.z()) < 5. &&
00721 fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
00722 fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
00723
00724 if(debug) std::cout<<"MB4 \t \t \t \t yes"<<std::endl;
00725
00726 RPCDetId rollId = rollasociated->id();
00727
00728 RPCGeomServ rpcsrv(rollId);
00729 std::string nameRoll = rpcsrv.name();
00730 if(debug) std::cout<<"MB4 \t \t \t \t \t The RPCName is "<<nameRoll<<std::endl;
00731 const float stripPredicted=
00732 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
00733
00734 if(debug) std::cout<<"MB4 \t \t \t \t Candidate (from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
00735
00736
00737 std::map<std::string, MonitorElement*> meMap=meCollection[rollId.rawId()];
00738
00739 if(debug) std::cout<<"MB4 \t \t \t \t \t Filling Expected"<<std::endl;
00740 sprintf(meIdDT,"ExpectedOccupancyFromDT_%d",rollId.rawId());
00741 meMap[meIdDT]->Fill(stripPredicted);
00742
00743
00744
00745
00746 int countRecHits = 0;
00747 int cluSize = 0;
00748 float minres = 3000.;
00749
00750 if(debug) std::cout<<"MB4 \t \t \t \t Getting RecHits in Roll Asociated"<<std::endl;
00751 typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
00752 rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
00753 RPCRecHitCollection::const_iterator recHit;
00754
00755 for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
00756 countRecHits++;
00757 LocalPoint recHitPos=recHit->localPosition();
00758 float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
00759 if(debug) std::cout<<"DT \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction."<<std::endl;
00760 if(fabs(res)<fabs(minres)){
00761 minres=res;
00762 cluSize = recHit->clusterSize();
00763 }
00764 }
00765
00766 bool anycoincidence=false;
00767
00768 if(countRecHits==0){
00769 if(debug) std::cout <<"MB4 \t \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT"<<std::endl;
00770 }else{
00771 assert(minres!=3000);
00772
00773 if(debug) std::cout<<"MB4 \t \t \t \t \t \t PointExtrapolatedRPCFrame.x="<<PointExtrapolatedRPCFrame.x()<<" Minimal Residual ="<<minres<<std::endl;
00774 if(debug) std::cout<<"MB4 \t \t \t \t \t \t Minimal Residual less than stripw*rangestrips? minres="<<minres<<" range="<<rangestrips<<" stripw="<<stripw<<" cluSize"<<cluSize<<" <=compare minres with"<<(rangestrips+cluSize*0.5)*stripw<<std::endl;
00775 if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
00776 if(debug) std::cout<<"MB4 \t \t \t \t \t \t \t True!"<<std::endl;
00777 anycoincidence=true;
00778 }
00779 }
00780 if(anycoincidence){
00781 if(debug) std::cout<<"MB4 \t \t \t \t \t At least one RecHit inside the range, Predicted="<<stripPredicted<<" minres="<<minres<<"cm range="<<rangestrips<<"strips stripw="<<stripw<<"cm"<<std::endl;
00782 if(debug) std::cout<<"MB4 \t \t \t \t \t Norm of Cosine Directors="<<dx3*dx3+dy3*dy3+dz3*dz3<<"~1?"<<std::endl;
00783
00784 float cosal = dx/sqrt(dx*dx+dz*dz);
00785 if(debug) std::cout<<"MB4 \t \t \t \t \t Angle="<<acos(cosal)*180/3.1415926<<" degree"<<std::endl;
00786 if(debug) std::cout<<"MB4 \t \t \t \t \t Filling the Residuals Histogram for globals with "<<minres<<"And the angular incidence with Cos Theta="<<-1*dz<<std::endl;
00787 assert(rollId.station()==4);
00788 if(cluSize==1*dupli){ hGlobalResClu1La6->Fill(minres);}
00789 else if(cluSize==2*dupli){ hGlobalResClu2La6->Fill(minres);}
00790 else if(cluSize==3*dupli){ hGlobalResClu3La6->Fill(minres);}
00791
00792 sprintf(meIdRPC,"RPCDataOccupancyFromDT_%d",rollId.rawId());
00793 if(debug) std::cout<<"MB4 \t \t \t \t \t \t COINCIDENCE!!! Event="<<iEvent.id()<<"Filling RPC Data Occupancy for "<<meIdRPC<<" with "<<stripPredicted<<std::endl;
00794 meMap[meIdRPC]->Fill(stripPredicted);
00795 }
00796 else{
00797 RPCGeomServ rpcsrv(rollasociated->id());
00798 std::string nameRoll = rpcsrv.name();
00799 if(debug) std::cout<<"MB4 \t \t \t \t \t \t A roll was ineficient in event"<<iEvent.id().event()<<std::endl;
00800 if(debug) ofrej<<"MB4 \t Wh "<<dtWheel
00801 <<"\t St "<<dtStation
00802 <<"\t Se "<<dtSector
00803 <<"\t Roll "<<nameRoll
00804 <<"\t Event "
00805 <<iEvent.id().event()
00806 <<"\t Run "
00807 <<iEvent.id().run()
00808 <<std::endl;
00809 }
00810 }else{
00811 if(debug) std::cout<<"MB4 \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00812 }
00813 }
00814 else{
00815 if(debug) std::cout<<"MB4 \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00816 }
00817 }
00818 }else{
00819 compatiblesegments=false;
00820 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;
00821 }
00822 }else{
00823 if(debug) std::cout<<"MB4 \t \t \t No the same station or same wheel or segment dim in mb3 not 4D"<<std::endl;
00824 }
00825 }
00826 }else{
00827 if(debug) std::cout<<"MB4 \t \t \t Is NOT a 2D Segment"<<std::endl;
00828 }
00829 }else{
00830 if(debug) std::cout<<"MB4 \t \t \t \t There is not just one segment or is not in station 4"<<std::endl;
00831 }
00832 }
00833 }else{
00834 if(debug) std::cout<<"MB4 \t This event doesn't have 4D Segment"<<std::endl;
00835 }
00836 }
00837 }
00838
00839
00840 if(inclcsc){
00841 if(debug) std::cout <<"\t Getting the CSC Segments"<<std::endl;
00842 edm::Handle<CSCSegmentCollection> allCSCSegments;
00843
00844 iEvent.getByLabel(cscSegments, allCSCSegments);
00845
00846 if(allCSCSegments.isValid()){
00847 if(allCSCSegments->size()>0){
00848 statistics->Fill(18);
00849
00850 if(debug) std::cout<<"CSC \t Number of CSC Segments in this event = "<<allCSCSegments->size()<<std::endl;
00851
00852 std::map<CSCDetId,int> CSCSegmentsCounter;
00853 CSCSegmentCollection::const_iterator segment;
00854
00855 int segmentsInThisEventInTheEndcap=0;
00856
00857 for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
00858 CSCSegmentsCounter[segment->cscDetId()]++;
00859 segmentsInThisEventInTheEndcap++;
00860 }
00861
00862 statistics->Fill(allCSCSegments->size()+18);
00863
00864 if(debug) std::cout<<"CSC \t loop over all the CSCSegments "<<std::endl;
00865 for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
00866 CSCDetId CSCId = segment->cscDetId();
00867
00868 if(debug) std::cout<<"CSC \t \t This Segment is in Chamber id: "<<CSCId<<std::endl;
00869 if(debug) std::cout<<"CSC \t \t Number of segments in this CSC = "<<CSCSegmentsCounter[CSCId]<<std::endl;
00870 if(debug) std::cout<<"CSC \t \t Is the only one in this CSC? is not ind the ring 1 or station 4? Are there more than 2 segments in the event?"<<std::endl;
00871
00872 if(CSCSegmentsCounter[CSCId]==1 && CSCId.station()!=4 && CSCId.ring()!=1 && allCSCSegments->size()>=2){
00873 if(debug) std::cout<<"CSC \t \t yes"<<std::endl;
00874 int cscEndCap = CSCId.endcap();
00875 int cscStation = CSCId.station();
00876 int cscRing = CSCId.ring();
00877 int cscChamber = CSCId.chamber();
00878 int rpcRegion = 1; if(cscEndCap==2) rpcRegion= -1;
00879 int rpcRing = cscRing;
00880 if(cscRing==4)rpcRing =1;
00881 int rpcStation = cscStation;
00882 int rpcSegment = CSCId.chamber();
00883
00884 LocalPoint segmentPosition= segment->localPosition();
00885 LocalVector segmentDirection=segment->localDirection();
00886 float dz=segmentDirection.z();
00887
00888 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;
00889
00890 if(segment->dimension()==4 && (segment->nRecHits()<=10 && segment->nRecHits()>=4)&& acos(dz)*180/3.1415926 > 45. && acos(dz)*180/3.1415926 < 160. ){
00891
00892
00893
00894 if(debug) std::cout<<"CSC \t \t yes"<<std::endl;
00895 if(debug) std::cout<<"CSC \t \t CSC Segment Dimension "<<segment->dimension()<<std::endl;
00896
00897 float Xo=segmentPosition.x();
00898 float Yo=segmentPosition.y();
00899 float Zo=segmentPosition.z();
00900 float dx=segmentDirection.x();
00901 float dy=segmentDirection.y();
00902 float dz=segmentDirection.z();
00903
00904
00905 if(debug) std::cout<<"CSC \t \t Getting chamber from Geometry"<<std::endl;
00906 const CSCChamber* TheChamber=cscGeo->chamber(CSCId);
00907 if(debug) std::cout<<"CSC \t \t Getting ID from Chamber"<<std::endl;
00908 const CSCDetId TheId=TheChamber->id();
00909 if(debug) std::cout<<"CSC \t \t Printing The Id"<<TheId<<std::endl;
00910 std::set<RPCDetId> rollsForThisCSC = rollstoreCSC[CSCStationIndex(rpcRegion,rpcStation,rpcRing,rpcSegment)];
00911 if(debug) std::cout<<"CSC \t \t Number of rolls for this CSC = "<<rollsForThisCSC.size()<<std::endl;
00912
00913 if(debug) std::cout<<"CSC \t \t Loop over all the rolls asociated to this CSC"<<std::endl;
00914
00915 if(rpcRing!=1&&rpcStation!=4){
00916
00917 if(rollsForThisCSC.size()==0){
00918 if(debug) std::cout<<"CSC Fail for CSCId="<<TheId<<" rpcRegion="<<rpcRegion<<" rpcStation="<<rpcStation<<" rpcRing="<<rpcRing<<" rpcSegment="<<rpcSegment<<std::endl;
00919 }
00920
00921 assert(rollsForThisCSC.size()>=1);
00922
00923
00924 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisCSC.begin();iteraRoll != rollsForThisCSC.end(); iteraRoll++){
00925 const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00926 RPCDetId rpcId = rollasociated->id();
00927
00928 if(debug) std::cout<<"CSC \t \t \t We are in the roll getting the surface"<<rpcId<<std::endl;
00929 const BoundPlane & RPCSurface = rollasociated->surface();
00930
00931 if(debug) std::cout<<"CSC \t \t \t RollID: "<<rpcId<<std::endl;
00932
00933 if(debug) std::cout<<"CSC \t \t \t Doing the extrapolation to this roll"<<std::endl;
00934 if(debug) std::cout<<"CSC \t \t \t CSC Segment Direction in CSCLocal "<<segmentDirection<<std::endl;
00935 if(debug) std::cout<<"CSC \t \t \t CSC Segment Point in CSCLocal "<<segmentPosition<<std::endl;
00936
00937 GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
00938 if(debug) std::cout<<"CSC \t \t \t Center (0,0,0) of the Roll in Global"<<CenterPointRollGlobal<<std::endl;
00939 GlobalPoint CenterPointCSCGlobal = TheChamber->toGlobal(LocalPoint(0,0,0));
00940 if(debug) std::cout<<"CSC \t \t \t Center (0,0,0) of the CSC in Global"<<CenterPointCSCGlobal<<std::endl;
00941 GlobalPoint segmentPositionInGlobal=TheChamber->toGlobal(segmentPosition);
00942 if(debug) std::cout<<"CSC \t \t \t Segment Position in Global"<<segmentPositionInGlobal<<std::endl;
00943 LocalPoint CenterRollinCSCFrame = TheChamber->toLocal(CenterPointRollGlobal);
00944
00945 if(debug){
00946 float rpcphi=0;
00947 float cscphi=0;
00948
00949 (CenterPointRollGlobal.barePhi()<0)?
00950 rpcphi = 2*3.141592+CenterPointRollGlobal.barePhi():rpcphi=CenterPointRollGlobal.barePhi();
00951
00952 (CenterPointCSCGlobal.barePhi()<0)?
00953 cscphi = 2*3.1415926536+CenterPointCSCGlobal.barePhi():cscphi=CenterPointCSCGlobal.barePhi();
00954
00955 float df=fabs(cscphi-rpcphi);
00956 float dr=fabs(CenterPointRollGlobal.perp()-CenterPointCSCGlobal.perp());
00957 float diffz=CenterPointRollGlobal.z()-CenterPointCSCGlobal.z();
00958 float dfg=df*180./3.14159265;
00959
00960 if(debug) std::cout<<"CSC \t \t \t z of RPC="<<CenterPointRollGlobal.z()<<"z of CSC"<<CenterPointCSCGlobal.z()<<" dfg="<<dfg<<std::endl;
00961
00962
00963 RPCGeomServ rpcsrv(rpcId);
00964
00965
00966 if(dr>200.||fabs(dz)>55.||dfg>1.){
00967
00968 if (debug) std::cout
00969 <<"\t \t \t CSC Station= "<<CSCId.station()
00970 <<" Ring= "<<CSCId.ring()
00971 <<" Chamber= "<<CSCId.chamber()
00972 <<" cscphi="<<cscphi*180/3.14159265
00973 <<"\t RPC Station= "<<rpcId.station()
00974 <<" ring= "<<rpcId.ring()
00975 <<" segment =-> "<<rpcsrv.name()
00976 <<" rollphi="<<rpcphi*180/3.14159265
00977 <<"\t dfg="<<dfg
00978 <<" dz="<<diffz
00979 <<" dr="<<dr
00980 <<std::endl;
00981
00982 }
00983 }
00984
00985
00986
00987
00988 float D=CenterRollinCSCFrame.z();
00989
00990 float X=Xo+dx*D/dz;
00991 float Y=Yo+dy*D/dz;
00992 float Z=D;
00993
00994 const TrapezoidalStripTopology* top_=dynamic_cast<const TrapezoidalStripTopology*>(&(rollasociated->topology()));
00995 LocalPoint xmin = top_->localPosition(0.);
00996 if(debug) std::cout<<"CSC \t \t \t xmin of this Roll "<<xmin<<"cm"<<std::endl;
00997 LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00998 if(debug) std::cout<<"CSC \t \t \t xmax of this Roll "<<xmax<<"cm"<<std::endl;
00999 float rsize = fabs( xmax.x()-xmin.x() );
01000 if(debug) std::cout<<"CSC \t \t \t Roll Size "<<rsize<<"cm"<<std::endl;
01001 float stripl = top_->stripLength();
01002 float stripw = top_->pitch();
01003
01004 if(debug) std::cout<<"CSC \t \t \t Strip Lenght "<<stripl<<"cm"<<std::endl;
01005 if(debug) std::cout<<"CSC \t \t \t Strip Width "<<stripw<<"cm"<<std::endl;
01006
01007 if(debug) std::cout<<"CSC \t \t \t X Predicted in CSCLocal= "<<X<<"cm"<<std::endl;
01008 if(debug) std::cout<<"CSC \t \t \t Y Predicted in CSCLocal= "<<Y<<"cm"<<std::endl;
01009 if(debug) std::cout<<"CSC \t \t \t Z Predicted in CSCLocal= "<<Z<<"cm"<<std::endl;
01010
01011 float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
01012
01013 if(debug) std::cout<<"CSC \t \t \t Is the distance of extrapolation less than MaxD? ="<<extrapolatedDistance<<"cm"<<"MaxD="<<MaxD<<"cm"<<std::endl;
01014
01015 if(extrapolatedDistance<=MaxD){
01016
01017 if(debug) std::cout<<"CSC \t \t \t yes"<<std::endl;
01018
01019 GlobalPoint GlobalPointExtrapolated=TheChamber->toGlobal(LocalPoint(X,Y,Z));
01020 if(debug) std::cout<<"CSC \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
01021
01022
01023 LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
01024
01025 if(debug) std::cout<<"CSC \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
01026 if(debug) std::cout<<"CSC \t \t \t Corner of the Roll = ("<<rsize*0.5<<","<<stripl*0.5<<")"<<std::endl;
01027 if(debug) std::cout<<"CSC \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.x())<<","
01028 <<fabs(PointExtrapolatedRPCFrame.y())<<","<<fabs(PointExtrapolatedRPCFrame.z())<<")"<<std::endl;
01029 if(debug) std::cout<<"CSC \t \t \t dz="
01030 <<fabs(PointExtrapolatedRPCFrame.z())<<" dx="
01031 <<fabs(PointExtrapolatedRPCFrame.x())<<" dy="
01032 <<fabs(PointExtrapolatedRPCFrame.y())<<std::endl;
01033
01034 if(debug) std::cout<<"CSC \t \t \t Does the extrapolation go inside this roll????"<<std::endl;
01035
01036 if(fabs(PointExtrapolatedRPCFrame.z()) < 10. &&
01037 fabs(PointExtrapolatedRPCFrame.x()) < rsize*0.5 &&
01038 fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){
01039
01040 if(debug) std::cout<<"CSC \t \t \t \t yes"<<std::endl;
01041
01042 RPCDetId rollId = rollasociated->id();
01043
01044 RPCGeomServ rpcsrv(rollId);
01045 std::string nameRoll = rpcsrv.name();
01046 if(debug) std::cout<<"CSC \t \t \t \t The RPCName is "<<nameRoll<<std::endl;
01047
01048 const float stripPredicted =
01049 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.));
01050
01051 if(debug) std::cout<<"CSC \t \t \t \t \t Candidate"<<rollId<<" "<<"(from CSC Segment) STRIP---> "<<stripPredicted<< std::endl;
01052
01053
01054 std::map<std::string, MonitorElement*> meMap=meCollection[rpcId.rawId()];
01055
01056 if(debug) std::cout<<"CSC \t \t \t \t Filling Expected"<<std::endl;
01057 sprintf(meIdCSC,"ExpectedOccupancyFromCSC_%d",rollId.rawId());
01058 meMap[meIdCSC]->Fill(stripPredicted);
01059
01060
01061
01062
01063 int cluSize = 0;
01064 int countRecHits = 0;
01065 float minres = 3000.;
01066
01067 if(debug) std::cout<<"CSC \t \t \t \t \t Getting RecHits in Roll Asociated"<<std::endl;
01068 typedef std::pair<RPCRecHitCollection::const_iterator, RPCRecHitCollection::const_iterator> rangeRecHits;
01069 rangeRecHits recHitCollection = rpcHits->get(rollasociated->id());
01070 RPCRecHitCollection::const_iterator recHit;
01071
01072 for (recHit = recHitCollection.first; recHit != recHitCollection.second ; recHit++) {
01073
01074 sprintf(meIdRPC,"BXDistribution_%d",rollasociated->id().rawId());
01075 meMap[meIdRPC]->Fill(recHit->BunchX());
01076
01077 countRecHits++;
01078 LocalPoint recHitPos=recHit->localPosition();
01079 float res=PointExtrapolatedRPCFrame.x()- recHitPos.x();
01080 if(debug) std::cout<<"CSC \t \t \t \t \t \t Found Rec Hit at "<<res<<"cm of the prediction."<<std::endl;
01081 if(fabs(res)<fabs(minres)){
01082 minres=res;
01083 cluSize = recHit->clusterSize();
01084 if(debug) std::cout<<"CSC \t \t \t \t \t \t \t New Min Res "<<res<<"cm."<<std::endl;
01085 }
01086 }
01087
01088 bool anycoincidence = false;
01089
01090 if(countRecHits==0){
01091 if(debug) std::cout <<"CSC \t \t \t \t \t THIS ROLL DOESN'T HAVE ANY RECHIT"<<std::endl;
01092 }else{
01093 assert(minres!=3000);
01094
01095 if(debug) std::cout<<"CSC \t \t \t \t \t PointExtrapolatedRPCFrame.x="<<PointExtrapolatedRPCFrame.x()<<" Minimal Residual"<<minres<<std::endl;
01096 if(debug) std::cout<<"CSC \t \t \t \t \t Minimal Residual less than stripw*rangestrips? minres="<<minres<<" range="<<rangestrips<<" stripw="<<stripw<<" cluSize"<<cluSize<<" <=compare minres with"<<(rangestrips+cluSize*0.5)*stripw<<std::endl;
01097 if(fabs(minres)<=(rangestrips+cluSize*0.5)*stripw){
01098 if(debug) std::cout<<"CSC \t \t \t \t \t \t True!"<<std::endl;
01099 anycoincidence=true;
01100 }
01101 }
01102 if(anycoincidence){
01103 if(debug) std::cout<<"CSC \t \t \t \t \t At least one RecHit inside the range, Predicted="<<stripPredicted<<" minres="<<minres<<"cm range="<<rangestrips<<"strips stripw="<<stripw<<"cm"<<std::endl;
01104 if(debug) std::cout<<"CSC \t \t \t \t \t Norm of Cosine Directors="<<dx*dx+dy*dy+dz*dz<<"~1?"<<std::endl;
01105
01106 float cosal = dx/sqrt(dx*dx+dz*dz);
01107 if(debug) std::cout<<"CSC \t \t \t \t \t Angle="<<acos(cosal)*180/3.1415926<<" degree"<<std::endl;
01108 if(debug) std::cout<<"CSC \t \t \t \t \t Filling the Residuals Histogram for globals with "<<minres<<"And the angular incidence with Cos Theta="<<-1*dz<<std::endl;
01109 if(rollId.ring()==2&&rollId.roll()==1){if(cluSize==1*dupli) hGlobalResClu1R2A->Fill(minres); if(cluSize==2*dupli) hGlobalResClu2R2A->Fill(minres); if(cluSize==3*dupli) hGlobalResClu3R2A->Fill(minres);}
01110 if(rollId.ring()==2&&rollId.roll()==2){if(cluSize==1*dupli) hGlobalResClu1R2B->Fill(minres); if(cluSize==2*dupli) hGlobalResClu2R2B->Fill(minres); if(cluSize==3*dupli) hGlobalResClu3R2B->Fill(minres);}
01111 if(rollId.ring()==2&&rollId.roll()==3){if(cluSize==1*dupli) hGlobalResClu1R2C->Fill(minres); if(cluSize==2*dupli) hGlobalResClu2R2C->Fill(minres); if(cluSize==3*dupli) hGlobalResClu3R2C->Fill(minres);}
01112 if(rollId.ring()==3&&rollId.roll()==1){if(cluSize==1*dupli) hGlobalResClu1R3A->Fill(minres); if(cluSize==2*dupli) hGlobalResClu2R3A->Fill(minres); if(cluSize==3*dupli) hGlobalResClu3R3A->Fill(minres);}
01113 if(rollId.ring()==3&&rollId.roll()==2){if(cluSize==1*dupli) hGlobalResClu1R3B->Fill(minres); if(cluSize==2*dupli) hGlobalResClu2R3B->Fill(minres); if(cluSize==3*dupli) hGlobalResClu3R3B->Fill(minres);}
01114 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);}
01115
01116 sprintf(meIdRPC,"RPCDataOccupancyFromCSC_%d",rollId.rawId());
01117 if(debug) std::cout <<"CSC \t \t \t \t \t \t COINCEDENCE!!! Event="<<iEvent.id()<<"Filling RPC Data Occupancy for "<<meIdRPC<<" with "<<stripPredicted<<std::endl;
01118 meMap[meIdRPC]->Fill(stripPredicted);
01119 }
01120 else{
01121 RPCGeomServ rpcsrv(rollasociated->id());
01122 std::string nameRoll = rpcsrv.name();
01123 if(debug) std::cout<<"CSC \t \t \t \t \t \t A roll was ineficient in event"<<iEvent.id().event()<<std::endl;
01124 if(debug) ofrej<<"CSC \t EndCap "<<rpcRegion
01125 <<"\t cscStation "<<cscStation
01126 <<"\t cscRing "<<cscRing
01127 <<"\t cscChamber "<<cscChamber
01128 <<"\t Roll "<<nameRoll
01129 <<"\t Event "<<iEvent.id().event()
01130 <<"\t CSCId "<<CSCId
01131 <<"\t Event "
01132 <<iEvent.id().event()
01133 <<"\t Run "
01134 <<iEvent.id().run()
01135 <<std::endl;
01136 }
01137 }else{
01138 if(debug) std::cout<<"CSC \t \t \t \t No the prediction is outside of this roll"<<std::endl;
01139 }
01140 }else{
01141 if(debug) std::cout<<"CSC \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
01142 }
01143 }
01144 }
01145 }
01146 }else{
01147 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;
01148 }
01149 }
01150 }else{
01151 if(debug) std::cout<<"CSC This Event doesn't have any CSCSegment"<<std::endl;
01152 }
01153 }
01154 }
01155 }
01156 }
01157
01158
01159 void RPCEfficiency::endRun(const edm::Run& r, const edm::EventSetup& iSetup){
01160 if (EffSaveRootFile){
01161 dbe->save(EffRootFileName);
01162 }
01163 }
01164
01165
01166 void RPCEfficiency::endJob(){
01167 dbe =0;
01168 }
1.7.3