CMSSW Reference Manual

00001 
00011 #include <DQM/RPCMonitorDigi/interface/RPCMonitorEfficiency.h>
00012 
00014 #include "FWCore/Framework/interface/ESHandle.h"
00015 
00016 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00017 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
00018 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
00019 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00020 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00021 
00023 #include "DataFormats/RPCDigi/interface/RPCDigi.h"
00024 #include "DataFormats/RPCDigi/interface/RPCDigiCollection.h"
00025 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
00026 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00027 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00028 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
00029 
00030 
00031 
00032 #include <cmath>
00033 
00034 
00035 class DTStationIndex{
00036 public: 
00037   DTStationIndex():_region(0),_wheel(0),_sector(0),_station(0){}
00038   DTStationIndex(int region, int wheel, int sector, int station) : 
00039     _region(region),
00040     _wheel(wheel),
00041     _sector(sector),
00042     _station(station){}
00043   ~DTStationIndex(){}
00044   int region() const {return _region;}
00045   int wheel() const {return _wheel;}
00046   int sector() const {return _sector;}
00047   int station() const {return _station;}
00048   bool operator<(const DTStationIndex& dtind) const{
00049     if(dtind.region()!=this->region())
00050       return dtind.region()<this->region();
00051     else if(dtind.wheel()!=this->wheel())
00052       return dtind.wheel()<this->wheel();
00053     else if(dtind.sector()!=this->sector())
00054       return dtind.sector()<this->sector();
00055     else if(dtind.station()!=this->station())
00056       return dtind.station()<this->station();
00057     return false;
00058   }
00059 private:
00060   int _region;
00061   int _wheel;
00062   int _sector;
00063   int _station; 
00064 };
00065 
00066 
00067 //**********************************************************************************************************
00068 
00069 
00070 
00071 RPCMonitorEfficiency::RPCMonitorEfficiency( const edm::ParameterSet& pset ){
00072   std::map<RPCDetId, int> buff;
00073   counter.clear();
00074   counter.reserve(3);
00075   //std::cout <<" Buff 1"<<std::endl;
00076   counter.push_back(buff);
00077   //std::cout <<" Buff 2"<<std::endl;
00078   counter.push_back(buff);
00079   //std::cout <<" Buff 3"<<std::endl;
00080   counter.push_back(buff);
00081   totalcounter.clear();
00082   totalcounter.reserve(3);
00083   totalcounter[0]=0;
00084   totalcounter[1]=0;
00085   totalcounter[2]=0;
00086   theRecHits4DLabel = pset.getParameter<std::string>("recHits4DLabel");
00087   digiLabel=pset.getParameter<std::string>("digiLabel");
00088   EffSaveRootFile  = pset.getUntrackedParameter<bool>("EffSaveRootFile", false); 
00089   EffSaveRootFileEventsInterval  = pset.getUntrackedParameter<int>("EffEventsInterval", 10000); 
00090   EffRootFileName  = pset.getUntrackedParameter<std::string>("EffRootFileName", "RPCEfficiency.root"); 
00091 
00093   dbe = edm::Service<DQMStore>().operator->();
00094   
00095 
00096   dbe->showDirStructure();
00097 
00098   _idList.clear(); 
00099   //ofrej.open("rejected.txt");
00100 }
00101 
00102 void RPCMonitorEfficiency::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup ){
00103   std::map<RPCDetId, int> buff;
00104   char layerLabel[128];
00105   char meIdRPC [128];
00106   char meIdDT [128];
00107   float dx=0.,dy=0.,dz=0.,Xo=0.,Yo=0.,X=0.,Y=0.,Z=0.,p1x=0.,p2x=0.,p3x=0.,p4x=0.,p1z=0.,p2z=0.,p3z=0.,p4z=0.,dx3=0.,dy3=0.,dz3=0.,Xo3=0.,Yo3=0.,x3=0.,x4=0.,z3=0.,z4=0.,m3=0.,m4=0.,xc=0.,zc=0.,b3=0.,b4=0.,w3=0.,w4=0.;
00108 
00109   float widestrip=5.;
00110   float widestripsRB4=8.;
00111   float circError=3.;
00112   float angle=0.01;
00113 
00114   edm::ESHandle<DTGeometry> dtGeo;
00115   iSetup.get<MuonGeometryRecord>().get(dtGeo);
00116   
00117   edm::ESHandle<RPCGeometry> rpcGeo;
00118   iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00119   
00120   edm::Handle<DTRecSegment4DCollection> all4DSegments;
00121   iEvent.getByLabel(theRecHits4DLabel, all4DSegments);
00122   
00123   edm::Handle<RPCDigiCollection> rpcDigis;
00124   iEvent.getByLabel(digiLabel, rpcDigis);
00125 
00126   std::map<DTStationIndex,std::set<RPCDetId> > rollstore;
00127   for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00128     RPCRoll* ir = dynamic_cast<RPCRoll*>(*it);
00129     RPCDetId rpcId = ir->id();
00130     int region=rpcId.region();
00131     int wheel=rpcId.ring();
00132     int sector=rpcId.sector();
00133     int station=rpcId.station();
00134     
00135     DTStationIndex ind(region,wheel,sector,station);
00136     
00137     std::set<RPCDetId> myrolls;
00138     if (rollstore.find(ind)!=rollstore.end()){
00139       myrolls=rollstore[ind];
00140     }
00141     myrolls.insert(rpcId);
00142     rollstore[ind]=myrolls;
00143   }
00144 
00145   
00146   if(all4DSegments->size()>0){
00147     //std::cout<<"Number of Segments in this event = "<<all4DSegments->size()<<std::endl;
00148     
00149 
00150     std::map<DTChamberId,int> scounter;
00151     DTRecSegment4DCollection::const_iterator segment;  
00152 
00153     
00154     for (segment = all4DSegments->begin();segment!=all4DSegments->end(); ++segment){
00155       scounter[segment->chamberId()]++;
00156     }    
00157     
00158     //std::cout<<"Loop over all the 4D Segments"<<std::endl;
00159     //loop over all the 4D Segments
00160 
00161     for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){ 
00162       DTChamberId DTId = segment->chamberId();
00163       //std::cout<<"\t This Segment is in Chamber id: "<<DTId<<std::endl;
00164       //std::cout<<"\t Number of segments in this DT = "<<scounter[DTId]<<std::endl;
00165       //std::cout<<"\t DT Segment Dimension "<<segment->dimension()<<std::endl; 
00166       //std::cout<<"\t Is the only in this DT?"<<std::endl;
00167 
00168       //there must be only one segment per Chamber
00169       if(scounter[DTId] == 1){  
00170         //std::cout<<"\t \t yes"<<std::endl;
00171         int dtWheel = DTId.wheel();
00172         int dtStation = DTId.station();
00173         int dtSector = DTId.sector();
00174         
00175         LocalPoint segmentPosition= segment->localPosition();
00176         LocalVector segmentDirection=segment->localDirection();
00177 
00178         const GeomDet* gdet=dtGeo->idToDet(segment->geographicalId());
00179         const BoundPlane & DTSurface = gdet->surface();
00180 
00181         //check if the dimension of the segment is 4 
00182 
00183         if(segment->dimension()==4){
00184           Xo=segmentPosition.x();
00185           Yo=segmentPosition.y();
00186           dx=segmentDirection.x();
00187           dy=segmentDirection.y();
00188           dz=segmentDirection.z();
00189           //std::cout<<"\t \t Loop over all the rolls asociated to this DT"<<std::endl;
00190           std::set<RPCDetId> rollsForThisDT = 
00191             rollstore[DTStationIndex(0,dtWheel,dtSector,dtStation)];
00192           //Loop over all the rolls
00193           for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00194             const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00195             //To get the roll's surface
00196             const BoundPlane & RPCSurface = rollasociated->surface(); 
00197             //std::cout<<"\t \t RollID: "<<rollasociated->id()<<std::endl;
00198             //std::cout<<"\t \t Doing the extrapolation"<<std::endl;
00199             //std::cout<<"\t \t DT Segment Direction in DTLocal "<<segmentDirection<<std::endl;
00200             //std::cout<<"\t \t DT Segment Point in DTLocal "<<segmentPosition<<std::endl;
00201             
00202             GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
00203             //std::cout<<"\t \t Center (0,0,0) of the Roll in Global"<<CenterPointRollGlobal<<std::endl;
00204             
00205             LocalPoint CenterRollinDTFrame = DTSurface.toLocal(CenterPointRollGlobal);
00206             //std::cout<<"\t \t Center (0,0,0) Roll In DTLocal"<<CenterRollinDTFrame<<std::endl;
00207             
00208             float D=CenterRollinDTFrame.z();
00209             //std::cout<<"\t \t D="<<D<<"cm"<<std::endl;
00210             
00211             X=Xo+dx*D/dz;
00212             Y=Yo+dy*D/dz;
00213             Z=D;
00214             
00215             const RectangularStripTopology* top_= dynamic_cast<const RectangularStripTopology*> (&(rollasociated->topology()));
00216             LocalPoint xmin = top_->localPosition(0.);
00217             LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00218             float rsize = fabs( xmax.x()-xmin.x() )*0.5;
00219             float stripl = top_->stripLength();
00220                   
00221             //std::cout<<"\t \t X Predicted in DTLocal= "<<X<<"cm"<<std::endl;
00222             //std::cout<<"\t \t Y Predicted in DTLocal= "<<Y<<"cm"<<std::endl;
00223             //std::cout<<"\t \t Z Predicted in DTLocal= "<<Z<<"cm"<<std::endl;
00224             
00225             GlobalPoint GlobalPointExtrapolated = 
00226               DTSurface.toGlobal(LocalPoint(X,Y,Z));
00227             //std::cout<<"\t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
00228             
00229             LocalPoint PointExtrapolatedRPCFrame =
00230               RPCSurface.toLocal(GlobalPointExtrapolated);
00231             //std::cout<<"\t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
00232             //std::cout<<"\t \t Does the extrapolation go inside this roll?"<<std::endl;
00233             //conditions to find the right roll to extrapolate
00234             
00235             if(fabs(PointExtrapolatedRPCFrame.z()) < 0.01 && fabs(PointExtrapolatedRPCFrame.x()) < rsize && fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){ 
00236               //std::cout<<"\t \t \t yes"<<std::endl;   
00237               //getting the number of the strip
00238               const float stripPredicted = 
00239                 rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.)); 
00240               
00241               //std::cout<<"\t \t \t Candidate"<<rollasociated->id()<<" "<<"(from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
00242               
00243               //--------- HISTOGRAM STRIP PREDICTED FROM DT  -------------------
00244               
00245               RPCDetId  rollId = rollasociated->id();
00246               uint32_t id = rollId.rawId();
00247 
00248               RPCDetId otherRollId1,otherRollId2;
00249 
00250               if(rollId.roll() == 1){
00251                 RPCDetId tempRollId1(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),2);
00252                 RPCDetId tempRollId2(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),3);
00253                 otherRollId1 = tempRollId1;
00254                 otherRollId2 = tempRollId2;
00255               }
00256               else if(rollId.roll() == 2){
00257 
00258                 RPCDetId tempRollId1(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),1);
00259                 RPCDetId tempRollId2(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),3);
00260                 otherRollId1 = tempRollId1;
00261                 otherRollId2 = tempRollId2;
00262               }
00263               else if(rollId.roll() == 3){
00264                 RPCDetId tempRollId1(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),1);
00265                 RPCDetId tempRollId2(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),2);
00266                 otherRollId1 = tempRollId1;
00267                 otherRollId2 = tempRollId2;
00268               }
00269 
00270               RPCDigiCollection::Range rpcRangeDigi1=rpcDigis->get(otherRollId1);
00271               RPCDigiCollection::Range rpcRangeDigi2=rpcDigis->get(otherRollId2);
00272 
00273               _idList.push_back(id);
00274               
00275               char detUnitLabel[128];
00276               sprintf(detUnitLabel ,"%d",id);
00277               sprintf(layerLabel ,"layer%d_subsector%d_roll%d",rollId.layer(),rollId.subsector(),rollId.roll());
00278               
00279               std::map<uint32_t, std::map<std::string,MonitorElement*> >::iterator meItr = meCollection.find(id);
00280               if (meItr == meCollection.end()){
00281                 meCollection[id] = bookDetUnitMEEff(rollId);
00282                 //std::cout << "\t \t \t Create new histograms  for "<<layerLabel<<std::endl;
00283               }
00284               
00285               std::map<std::string, MonitorElement*> meMap=meCollection[id];
00286               sprintf(meIdDT,"ExpectedOccupancyFromDT_%s",detUnitLabel);
00287               meMap[meIdDT]->Fill(stripPredicted);
00288 
00289               sprintf(meIdDT,"ExpectedOccupancy2DFromDT_%s",detUnitLabel);
00290               meMap[meIdDT]->Fill(stripPredicted,Y);
00291 
00292               //std::cout << "\t \t \t One for counterPREDICT"<<std::endl;
00293               totalcounter[0]++;
00294               buff=counter[0];
00295               buff[rollId]++;
00296               counter[0]=buff;
00297               //-------------------------------------------------------------------
00298               
00299               //std::cout<<"\t \t \t We have a Candidate let's see in the digis!"<<std::endl;
00300               
00301               bool anycoincidence=false;
00302               int stripDetected = 0;
00303               RPCDigiCollection::Range rpcRangeDigi=rpcDigis->get(rollasociated->id());
00304 
00305               int stripCounter = 0;
00306 
00307               for (RPCDigiCollection::const_iterator digiIt = rpcRangeDigi.first;digiIt!=rpcRangeDigi.second;++digiIt){//loop over the digis in the event
00308                 stripCounter++;
00309                 //std::cout<<"\t \t \t \t Digi "<<*digiIt<<std::endl;//print the digis in the event
00310                 stripDetected=digiIt->strip();
00311 
00312                 float res = (float)(stripDetected) - stripPredicted;
00313                 sprintf(meIdRPC,"RPCResiduals_%s",detUnitLabel);
00314                 meMap[meIdRPC]->Fill(res);
00315                 
00316                 sprintf(meIdRPC,"RPCResiduals2D_%s",detUnitLabel);
00317                 meMap[meIdRPC]->Fill(res,Y);
00318                 
00319                 if(res > 7)
00320                   std::cout<<"STRANGE     "<<"EVENTO NUM = "<<iEvent.id().event()<<"   Residuo = "<<res<<"   Strip Num = "<<stripDetected<<"   Strip totali = "<<stripCounter<<std::endl;
00321                 //compare the strip Detected with the predicted
00322                 if(fabs((float)(stripDetected) - stripPredicted) < widestrip){
00323                   //std::cout <<"\t \t \t \t COINCEDENCE Predict "
00324                   //    <<stripPredicted<<" Detect "
00325                   //    <<stripDetected<<std::endl;
00326                   anycoincidence=true;
00327                   //break;//funciona solo para hacerlo mas rapido
00328                   //We can not divide two diferents things
00329                 }
00330               }
00331               if (anycoincidence) {
00332                 sprintf(meIdRPC,"ExpectedOccupancyFromDT_forCrT_%s",detUnitLabel);
00333                 meMap[meIdRPC]->Fill(stripPredicted);
00334                 sprintf(meIdRPC,"RealDetectedOccupancy_%s",detUnitLabel);
00335                 meMap[meIdRPC]->Fill(stripDetected);
00336                 sprintf(meIdRPC,"YExpectedOccupancyFromDT_%s",detUnitLabel);
00337                 meMap[meIdRPC]->Fill(Y);
00338 
00339                 for (RPCDigiCollection::const_iterator digiIt1 = rpcRangeDigi1.first;digiIt1!=rpcRangeDigi1.second;++digiIt1){
00340 
00341                   if(fabs(stripDetected - digiIt1->strip()) <= 1){
00342                     
00343                     sprintf(meIdRPC,"XCrossTalk_1_%s",detUnitLabel);
00344                     meMap[meIdRPC]->Fill(stripPredicted);
00345 
00346                     sprintf(meIdRPC,"XDetectCrossTalk_1_%s",detUnitLabel);
00347                     meMap[meIdRPC]->Fill(stripDetected);
00348 
00349                     sprintf(meIdRPC,"YCrossTalk_1_%s",detUnitLabel);
00350                     meMap[meIdRPC]->Fill(Y);
00351                     break;
00352                   }
00353                 }
00354 
00355                 for (RPCDigiCollection::const_iterator digiIt2 = rpcRangeDigi2.first;digiIt2!=rpcRangeDigi2.second;++digiIt2){
00356 
00357                   if(fabs(stripDetected - digiIt2->strip()) <= 1){
00358                     
00359                     sprintf(meIdRPC,"XCrossTalk_2_%s",detUnitLabel);
00360                     meMap[meIdRPC]->Fill(stripPredicted);
00361 
00362                     sprintf(meIdRPC,"XDetectCrossTalk_2_%s",detUnitLabel);
00363                     meMap[meIdRPC]->Fill(stripDetected);
00364 
00365                     sprintf(meIdRPC,"YCrossTalk_2_%s",detUnitLabel);
00366                     meMap[meIdRPC]->Fill(Y);
00367                     break;
00368                   }
00369                 }
00370 
00371                 sprintf(meIdRPC,"RPCDataOccupancy_%s",detUnitLabel);
00372                 meMap[meIdRPC]->Fill(stripPredicted);
00373 
00374                 sprintf(meIdRPC,"RPCDataOccupancy2D_%s",detUnitLabel);
00375                 meMap[meIdRPC]->Fill(stripPredicted,Y);
00376 
00377                 totalcounter[1]++;
00378                 buff=counter[1];
00379                 buff[rollId]++;
00380                 counter[1]=buff;                
00381               }
00382               else {
00383                 //std::cout <<"\t \t \t \t XXXXX THIS PREDICTION DOESN'T HAVE ANY CORRESPONDENCE WITH THE DATA"<<std::endl;
00384                 totalcounter[2]++;
00385                 buff=counter[2];
00386                 buff[rollId]++;
00387                 counter[2]=buff;                
00388                 //std::cout << "\t \t \t \t One for counterFAIL"<<std::endl;
00389                 //ofrej<<"Wh "<<dtWheel<<"\t  St "<<dtStation
00390                 //     <<"\t Se "<<dtSector<<"\t Event "
00391                 //    <<iEvent.id().event()<<std::endl;
00392               }
00393             }
00394             else {
00395               //std::cout<<"\t \t \t no"<<std::endl;
00396             }//Condition for the right match
00397           }//loop over all the rolls
00398 
00399 
00400           // dedicated RB4 analysis part, that misses DT 4D segments
00401 
00402 
00403         }else if(segment->dimension()==2&&dtStation==4){
00404           
00405           LocalVector segmentDirectionMB4=segmentDirection;
00406           LocalPoint segmentPositionMB4=segmentPosition;
00407           
00408 
00409           //std::cout<<"\t \t 2D in RB4"<<DTId<<" with D="<<segment->dimension()<<localDirection<<segmentPositionMB4<<std::endl;          
00410           bool compatiblesegments=false;
00411           Xo=segmentPositionMB4.x();
00412           dx=segmentDirectionMB4.x();
00413           dz=segmentDirectionMB4.z();
00414           //std::cout<<"\t \t Loop over all the segments"<<std::endl;     
00415           DTRecSegment4DCollection::const_iterator segMB3;  
00416 
00417           const BoundPlane& DTSurface4 = dtGeo->idToDet(DTId)->surface();
00418           w4 = DTSurface4.bounds().thickness()*0.5; // along local Z
00419 
00420           for(segMB3=all4DSegments->begin();segMB3!=all4DSegments->end();++segMB3){
00421             DTChamberId dtid = segMB3->chamberId();
00422             if(dtid.station()==3){
00423               const GeomDet* gdet3=dtGeo->idToDet(segMB3->geographicalId());
00424               const BoundPlane & DTSurface3 = gdet3->surface();
00425               w3 = DTSurface3.bounds().thickness()*0.5; // along local Z
00426               
00427               dx3=segMB3->localDirection().x();
00428               dy3=segMB3->localDirection().y();
00429               dz3=segMB3->localDirection().z();
00430               
00431               //LocalVector segDirMB4inMB3Frame=DTSurface3.toLocal(DTSurface4.toGlobal(segmentDirectionMB4));
00432 
00433               if(fabs(dx-dx3)<=angle&&fabs(dz-dz3)<=angle){//same direction?
00434                 compatiblesegments=true;
00435               }else{
00436                 //They don't have the same local dir, and the segments in diferent sectors,compatibles with circle?
00437                 Xo3=segMB3->localPosition().x();
00438                 Yo3=segMB3->localPosition().y();
00439 
00440                 //Do we have segments compatibles with a muon?
00441                 
00442                 // pos=seg.position()
00443                 // dir=seg.direction()
00444                 // dest =pos + dir*height*cos(dir.theta())
00445                 
00446                 //std::cout<<"Ancho de MB3 w3="<<w3<<" |Ancho de MB4 w4= "<<w4<<std::endl;
00447                 //std::cout<<"Informacion Inicial MB4"<<segment->localDirection()<<" "<<segment->localPosition()<<std::endl;
00448                 //std::cout<<"Informacion Inicial MB3"<<segMB3->localDirection()<<" "<<segMB3->localPosition()<<std::endl;
00449                 compatiblesegments=false;
00450                 
00451                 //This in MB3 Local
00452                 
00453                 p1x=Xo3+dx3*w3/dz3;
00454                 p1z=w3;
00455                 
00456                 p2x=Xo3-dx3*w3/dz3;
00457                 p2z=-w3;
00458                 
00459                 //This in MB4 Local
00460                 
00461                 p3x=Xo+dx*w4/dz;
00462                 p3z=w4;
00463                 
00464                 p4x=Xo-dx*w4/dz;
00465                 p4z=-w4;
00466                 
00467                 LocalPoint P1=LocalPoint(p1x,0,p1z);
00468                 LocalPoint P2=LocalPoint(p2x,0,p2z);
00469                 LocalPoint P3=LocalPoint(p3x,0,p3z);
00470                 LocalPoint P4=LocalPoint(p4x,0,p4z);
00471                 
00472                 //std::cout<<"Points in MB3="<<P1<<" --- "<<P2<<std::endl;
00473                 //std::cout<<"Points in MB4="<<P3<<" --- "<<P4<<std::endl;
00474                 
00475                 //Now we have to convert to global
00476                 
00477                 LocalPoint P1g=P1;
00478                 LocalPoint P2g=P2;
00479                 
00480                 LocalPoint P3g=DTSurface3.toLocal(DTSurface.toGlobal(P3));
00481                 LocalPoint P4g=DTSurface3.toLocal(DTSurface.toGlobal(P4));
00482                 
00483                 //std::cout<<"Points in MB3 ="<<P1g<<" --- "<<P2g<<std::endl;
00484                 //std::cout<<"Points in MB4 in MB3Frame="<<P3g<<" --- "<<P4g<<std::endl;
00485                 
00486                 
00487                 float dx3g=P1g.x()-P2g.x();
00488                 float dz3g=P1g.z()-P2g.z();
00489               
00490                 float dxg=P3g.x()-P4g.x();
00491                 float dzg=P3g.z()-P4g.z();
00492                 
00493                 //std::cout<<"dx3g="<<dx3g<<" dz3g="<<dz3g<<std::endl;
00494                 //std::cout<<"dxg="<<dxg<<" dzg="<<dzg<<std::endl;
00495                 
00496                 m3=-dx3g/dz3g;
00497                 m4=-dxg/dzg;
00498                 
00499                 x3=(P1g.x()+P2g.x())*0.5;
00500                 z3=(P1g.z()+P2g.z())*0.5;
00501                 
00502                 x4=(P3g.x()+P4g.x())*0.5;
00503                 z4=(P3g.z()+P4g.z())*0.5;
00504                 
00505                 b3=z3-m3*x3;
00506                 b4=z4-m4*x4;
00507                 
00508                 if(m3!=m4){
00509                   xc=(b4-b3)/(m3-m4);
00510                   zc=m3*xc+b3;
00511                 
00512                   GlobalPoint Pc=GlobalPoint(xc,0,zc);
00513                   
00514                   //std::cout<<Pc<<std::endl;
00515                   
00516                   float distance=fabs((GlobalPoint(P2g.x(),0,P2g.z())-GlobalPoint(Pc.x(),0,Pc.z())).mag()-(GlobalPoint(P3g.x(),0,P3g.z())-GlobalPoint(Pc.x(),0,Pc.z())).mag());
00517                   
00518                   if(distance<circError){
00519                     compatiblesegments=true;
00520                     //std::cout<<"YES in the same circle... p2 a C="<<P1g<<P2g<<P3g<<P4g<<"Distancia "<<distance<<std::endl;
00521                   }
00522                   else{
00523                     //std::cout<<"NOT in the same circle... p2 a C="<<P1g<<P2g<<P3g<<P4g<<"Distancia "<<distance<<std::endl;
00524                     compatiblesegments=false;
00525                   }
00526                 }
00527                 else{
00528                   //std::cout<<"We have the same slope m3="<<m3<<" m4="<<m4<<std::endl;
00529                   if(fabs(b3-b4)<=circError){
00530                     compatiblesegments=true;
00531                     //std::cout<<"and the segments are in a line"<<std::endl;
00532                   }else{std::cout<<"But we don't have the same intercept b3="<<b3<<" b4="<<b4<<std::endl;}
00533                   compatiblesegments=false;
00534                   system("sleep 30");
00535                 }
00536               }
00537                 
00538               //conditions in MB3
00539               if(scounter[dtid]==1 && compatiblesegments){
00540                 //std::cout<<"********\t \t \t In the same event there is a segment in RB3 "<<dtid<<" with D="<<segMB3->dimension()<<segMB3->localDirection()<<segMB3->localPosition()<<"scounter "<<scounter[dtid]<<std::endl;
00541                 
00542                 std::set<RPCDetId> rollsForThisDT = 
00543                   rollstore[DTStationIndex(0,dtWheel,dtSector,dtStation)];
00544                 //Loop over all the rolls asociated to RB4
00545                 for (std::set<RPCDetId>::iterator iteraRoll=rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00546                   const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00547                   //To get the roll's surface
00548                   const BoundPlane & RPCSurfaceRB4 = rollasociated->surface(); 
00549                   
00550                   const GeomDet* gdet=dtGeo->idToDet(segMB3->geographicalId());
00551                   const BoundPlane & DTSurfaceMB3 = gdet->surface();
00552                   
00553                   //std::cout<<"\t \t \t RollID: should be RB4"<<rollasociated->id()<<std::endl;
00554                   //std::cout<<"\t \t \t Making the extrapolation"<<std::endl;
00555                   //std::cout<<"\t \t \t DT Segment Direction in MB3 DTLocal "<<segMB3->localDirection()<<std::endl;
00556                   //std::cout<<"\t \t \t DT Segment Point in MB3 DTLocal "<<segMB3->localPosition()<<std::endl;
00557                   
00558                   GlobalPoint CenterPointRollGlobal=RPCSurfaceRB4.toGlobal(LocalPoint(0,0,0));
00559                   //std::cout<<"\t \t \t Center (0,0,0) of the RB4 Roll in Global"<<CenterPointRollGlobal<<std::endl;
00560                   
00561                   LocalPoint CenterRollinDTFrame = DTSurfaceMB3.toLocal(CenterPointRollGlobal);
00562                   //std::cout<<"\t \t \t Center (0,0,0) Roll In DT MB3 Local"<<CenterRollinDTFrame<<std::endl;
00563                   
00564                   float D=CenterRollinDTFrame.z();
00565                   //std::cout<<"\t \t \t D="<<D<<"cm"<<std::endl;
00566                   
00567                   X=Xo3+dx3*D/dz3;
00568                   Y=Yo3+dy3*D/dz3;
00569                   Z=D;
00570                   
00571                   const RectangularStripTopology* top_=dynamic_cast<const RectangularStripTopology*>(&(rollasociated->topology()));
00572                   LocalPoint xmin = top_->localPosition(0.);
00573                   LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00574                   float rsize = fabs( xmax.x()-xmin.x() )*0.5;
00575                   float stripl = top_->stripLength();
00576                   
00577                   //std::cout<<"\t \t \t X Predicted in DT MB3 Local= "<<X<<"cm"<<std::endl;
00578                   //std::cout<<"\t \t \t Y Predicted in DT MB3 Local= "<<Y<<"cm"<<std::endl;
00579                   //std::cout<<"\t \t \t Z Predicted in DT MB3 Local= "<<Z<<"cm"<<std::endl;
00580                   
00581                   GlobalPoint GlobalPointExtrapolated = DTSurfaceMB3.toGlobal(LocalPoint(X,Y,Z));
00582                   //std::cout<<"\t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
00583                   
00584                   LocalPoint PointExtrapolatedRPCFrame = RPCSurfaceRB4.toLocal(GlobalPointExtrapolated);
00585                   //std::cout<<"\t \t \t Point Extrapolated in RPC RB4 Local"<<PointExtrapolatedRPCFrame<< std::endl;
00586                   
00587                   //std::cout<<"\t \t \t Does the extrapolation go inside this roll?"<<std::endl;
00588                   //conditions to find the right roll to extrapolate
00589                   if(fabs(PointExtrapolatedRPCFrame.z()) < 0.01  &&
00590                      fabs(PointExtrapolatedRPCFrame.x()) < rsize &&
00591                      fabs(PointExtrapolatedRPCFrame.y()) < stripl*0.5){ 
00592                     
00593                     //std::cout<<"\t \t \t \t yes"<<std::endl;
00594                     //getting the number of the strip
00595                     const float stripPredicted=
00596                       rollasociated->strip(LocalPoint(PointExtrapolatedRPCFrame.x(),PointExtrapolatedRPCFrame.y(),0.)); 
00597                     
00598                     //std::cout<<"\t \t \t \t Candidate"<<rollasociated->id()<<" "<<"(from DT Segment) STRIP---> "<<stripPredicted<< std::endl;
00599                     
00600                     //--------- HISTOGRAM STRIP PREDICTED FROM DT  -------------------
00601                     
00602                     RPCDetId  rollId = rollasociated->id();
00603                     uint32_t id = rollId.rawId();
00604                     
00605                     RPCDetId otherRollId1,otherRollId2;
00606                     
00607                     if(rollId.roll() == 1){
00608                       RPCDetId tempRollId1(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),2);
00609                       RPCDetId tempRollId2(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),3);
00610                       otherRollId1 = tempRollId1;
00611                       otherRollId2 = tempRollId2;
00612                     }
00613                     else if(rollId.roll() == 2){
00614                       
00615                       RPCDetId tempRollId1(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),1);
00616                       RPCDetId tempRollId2(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),3);
00617                       otherRollId1 = tempRollId1;
00618                       otherRollId2 = tempRollId2;
00619                     }
00620                     else if(rollId.roll() == 3){
00621                       RPCDetId tempRollId1(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),1);
00622                       RPCDetId tempRollId2(rollId.region(),rollId.ring(),rollId.station(),rollId.sector(),rollId.layer(),rollId.subsector(),2);
00623                       otherRollId1 = tempRollId1;
00624                       otherRollId2 = tempRollId2;
00625                     }
00626                     RPCDigiCollection::Range rpcRangeDigi1=rpcDigis->get(otherRollId1);
00627                     RPCDigiCollection::Range rpcRangeDigi2=rpcDigis->get(otherRollId2);
00628 
00629                     _idList.push_back(id);
00630                     
00631                     char detUnitLabel[128];
00632                     sprintf(detUnitLabel ,"%d",id);
00633                     sprintf(layerLabel ,"layer%d_subsector%d_roll%d",rollId.layer(),rollId.subsector(),rollId.roll());
00634                     
00635                     std::map<uint32_t, std::map<std::string,MonitorElement*> >::iterator meItr = meCollection.find(id);
00636                     if (meItr == meCollection.end()){
00637                       meCollection[id] = bookDetUnitMEEff(rollId);
00638                       std::cout << "\t \t \t \t Create new histograms  for "<<layerLabel<<std::endl;
00639                     }
00640                     
00641                     std::map<std::string, MonitorElement*> meMap=meCollection[id];
00642                     sprintf(meIdDT,"ExpectedOccupancyFromDT_%s",detUnitLabel);
00643                     meMap[meIdDT]->Fill(stripPredicted);
00644 
00645                     sprintf(meIdDT,"ExpectedOccupancy2DFromDT_%s",detUnitLabel);
00646                     meMap[meIdDT]->Fill(stripPredicted,Y);
00647 
00648                     //std::cout << "\t \t \t \t One for counterPREDICT"<<std::endl;
00649                     totalcounter[0]++;
00650                     buff=counter[0];
00651                     buff[rollId]++;
00652                     counter[0]=buff;            
00653                     
00654                     bool anycoincidence=false;
00655                     int stripDetected = 0;
00656                     RPCDigiCollection::Range rpcRangeDigi = 
00657                       rpcDigis->get(rollasociated->id());
00658                                    
00659                     int stripCounter = 0;
00660 
00661                     //loop over the digis in the event
00662                     for (RPCDigiCollection::const_iterator digiIt = rpcRangeDigi.first;digiIt!=rpcRangeDigi.second;++digiIt){
00663                       stripCounter++;
00664                       //std::cout<<"\t \t \t \t \t Digi "<<*digiIt<<std::endl;//print the digis in the event
00665                       stripDetected=digiIt->strip();
00666 
00667                       stripCounter++;
00668                       float res = (float)(stripDetected) - stripPredicted;
00669                       sprintf(meIdRPC,"RPCResiduals_%s",detUnitLabel);
00670                       meMap[meIdRPC]->Fill(res);
00671                 
00672                       sprintf(meIdRPC,"RPCResiduals2D_%s",detUnitLabel);
00673                       meMap[meIdRPC]->Fill(res,Y);
00674                       
00675                       if(res > 7)
00676                         std::cout<<"STRANGE     "<<"EVENTO NUM = "<<iEvent.id().event()<<"   Residuo = "<<res<<"   Strip Num = "<<stripDetected<<"   Strip totali = "<<stripCounter<<std::endl;
00677                       if(fabs((float)(stripDetected) - stripPredicted)<widestripsRB4){//Detected Vs Predicted
00678                         //std::cout <<"\t \t \t \t \t COINCEDENCE Predict "<<stripPredicted<<" Detect "<<stripDetected<<std::endl;
00679                         anycoincidence=true;
00680                         //break;//funciona solo para hacerlo mas rapido
00681                       }
00682                     }
00683                     if (anycoincidence){
00684 
00685                       sprintf(meIdRPC,"ExpectedOccupancyFromDT_forCrT_%s",detUnitLabel);
00686                       meMap[meIdRPC]->Fill(stripPredicted);
00687                       sprintf(meIdRPC,"RealDetectedOccupancy_%s",detUnitLabel);
00688                       meMap[meIdRPC]->Fill(stripDetected);
00689                       sprintf(meIdRPC,"YExpectedOccupancyFromDT_%s",detUnitLabel);
00690                       meMap[meIdRPC]->Fill(Y);
00691 
00692                       for (RPCDigiCollection::const_iterator digiIt1 = rpcRangeDigi1.first;digiIt1!=rpcRangeDigi1.second;++digiIt1){
00693 
00694                         if(fabs(stripDetected - digiIt1->strip()) <= 1){
00695                           
00696                           sprintf(meIdRPC,"XCrossTalk_1_%s",detUnitLabel);
00697                           meMap[meIdRPC]->Fill(stripPredicted);
00698                           
00699                           sprintf(meIdRPC,"XDetectCrossTalk_1_%s",detUnitLabel);
00700                           meMap[meIdRPC]->Fill(stripDetected);
00701                           
00702                           sprintf(meIdRPC,"YCrossTalk_1_%s",detUnitLabel);
00703                           meMap[meIdRPC]->Fill(Y);
00704                           break;
00705                         }
00706                       }
00707 
00708                       for (RPCDigiCollection::const_iterator digiIt2 = rpcRangeDigi2.first;digiIt2!=rpcRangeDigi2.second;++digiIt2){
00709                         
00710                         if(fabs(stripDetected - digiIt2->strip()) <= 1){
00711                           
00712                           sprintf(meIdRPC,"XCrossTalk_2_%s",detUnitLabel);
00713                           meMap[meIdRPC]->Fill(stripPredicted);
00714                           
00715                           sprintf(meIdRPC,"XDetectCrossTalk_2_%s",detUnitLabel);
00716                           meMap[meIdRPC]->Fill(stripDetected);
00717                           
00718                           sprintf(meIdRPC,"YCrossTalk_2_%s",detUnitLabel);
00719                           meMap[meIdRPC]->Fill(Y);
00720                           break;
00721                         }
00722                       }
00723 
00724                       sprintf(meIdRPC,"RPCDataOccupancy_%s",detUnitLabel);
00725                       meMap[meIdRPC]->Fill(stripPredicted);
00726 
00727                       sprintf(meIdRPC,"RPCDataOccupancy2D_%s",detUnitLabel);
00728                       meMap[meIdRPC]->Fill(stripPredicted,Y);
00729 
00730                       totalcounter[1]++;
00731                       buff=counter[1];
00732                       buff[rollId]++;
00733                       counter[1]=buff;          
00734                     }
00735                     else{
00736                       totalcounter[2]++;
00737                       buff=counter[2];
00738                       buff[rollId]++;
00739                       counter[2]=buff;          
00740                       //std::cout <<"\t \t \t \t \t XXXXX THIS PREDICTION DOESN'T HAVE ANY CORRESPONDENCE WITH THE DATA"<<std::endl;
00741                       //std::cout << "\t \t \t \t \t One for counterFAIL"<<std::endl;
00742                       //ofrej<<"Wh "<<dtWheel<<"\t St "<<dtStation
00743                       //         <<"\t Se "<<dtSector<<"\t Event "
00744                       //         <<iEvent.id().event()<<std::endl;
00745                       
00746                     }
00747                   }
00748                   else{
00749                     //std::cout<<"\t \t \t \t no"<<std::endl;
00750                   }//Condition for the right match
00751                   
00752                 }//loop over all the rolls FOR RB3 -------------------------------
00753               }//Si tenemos solo un segmento en la chamber
00754             }//Condition Segment in MB3
00755           }//Loop avoer all the segments to see if it is the right MB3
00756         }
00757         else{
00758           //std::cout<<"\t \t Strange Segment"<<std::endl;
00759         }//Is not a 4D Segment neither a 2D in MB4
00760       }
00761       else {
00762         //std::cout<<"\t \t no"<<std::endl;
00763       }//There is one segment in the chamber?
00764     }//loop over the segments
00765   }else {
00766     //std::cout<<"This Event doesn't have any DT4DSegment"<<std::endl;
00767   }//is ther more than 1 segment in this event?
00768   
00769 }
00770 
00771 void RPCMonitorEfficiency::endJob(void){
00772   std::map<RPCDetId, int> pred = counter[0];
00773   std::map<RPCDetId, int> obse = counter[1];
00774   std::map<RPCDetId, int> reje = counter[2];
00775 
00776   std::map<RPCDetId, int>::iterator irpc;
00777   for (irpc=pred.begin(); irpc!=pred.end();irpc++){
00778     RPCDetId id=irpc->first;
00779     int p=pred[id]; 
00780     int o=obse[id]; 
00781     int r=reje[id]; 
00782     assert(p==o+r);
00783     float ef = float(o)/float(p);
00784     float er = sqrt(ef*(1.-ef)/float(p));
00785     std::cout <<"\n "<<id<<"\t Predicted "<<p<<"\t Observed "<<o<<"\t Eff = "<<ef*100.<<" % +/- "<<er*100.<<" %";
00786     if(ef<0.8){
00787       std::cout<<"\t \t Warning!";
00788     } 
00789   }
00790   
00791   float tote = float(totalcounter[1])/float(totalcounter[0]);
00792   float totr = sqrt(tote*(1.-tote)/float(totalcounter[0]));
00793   std::cout <<"\n\n \t \t TOTAL EFFICIENCY \t Predicted "<<totalcounter[1]<<"\t Observed "<<totalcounter[0]<<"\t Eff = "<<tote*100.<<"\t +/- \t"<<totr*100.<<"%"<<std::endl;
00794   std::cout <<totalcounter[1]<<" "<<totalcounter[0]<<" flagcode"<<std::endl;
00795 
00796   
00797   std::vector<uint32_t>::iterator meIt;
00798   for(meIt = _idList.begin(); meIt != _idList.end(); ++meIt){
00799 
00800     char detUnitLabel[128];
00801     char meIdRPC [128];
00802     char meIdDT [128];
00803     char meIdDTCrT [128];
00804     char effIdRPC [128];
00805     char meIdRPCCrT [128];
00806 
00807     char XCrTalkId_1 [128];
00808     char XCrTalkId_2 [128];
00809     char effXCrTalkId_1 [128];
00810     char effXCrTalkId_2 [128];
00811 
00812     char XCrTalkDetId_1 [128];
00813     char XCrTalkDetId_2 [128];
00814     char effXCrTalkDetId_1 [128];
00815     char effXCrTalkDetId_2 [128];
00816 
00817     char YCrTalkId_1 [128];
00818     char YCrTalkId_2 [128];
00819     char effYCrTalkId_1 [128];
00820     char effYCrTalkId_2 [128];
00821 
00822     char YCrTalkDTId [128];
00823 
00824     char meIdRPC_2D [128];
00825     char meIdDT_2D [128];
00826     char effIdRPC_2D [128];
00827 
00828     sprintf(detUnitLabel ,"%d",*meIt);
00829     sprintf(meIdRPC,"RPCDataOccupancy_%s",detUnitLabel);
00830     sprintf(meIdDT,"ExpectedOccupancyFromDT_%s",detUnitLabel);
00831     sprintf(meIdDTCrT,"ExpectedOccupancyFromDT_forCrT_%s",detUnitLabel);
00832     sprintf(YCrTalkDTId,"YExpectedOccupancyFromDT_%s",detUnitLabel);
00833     sprintf(meIdRPCCrT,"RealDetectedOccupancy_%s",detUnitLabel);
00834 
00835     sprintf(effIdRPC,"EfficienyFromDTExtrapolation_%s",detUnitLabel);
00836 
00837     sprintf(meIdRPC_2D,"RPCDataOccupancy2D_%s",detUnitLabel);
00838     sprintf(meIdDT_2D,"ExpectedOccupancy2DFromDT_%s",detUnitLabel);
00839     sprintf(effIdRPC_2D,"EfficienyFromDT2DExtrapolation_%s",detUnitLabel);
00840 
00841     sprintf(XCrTalkId_1,"XCrossTalk_1_%s",detUnitLabel);
00842     sprintf(XCrTalkId_2,"XCrossTalk_2_%s",detUnitLabel);
00843     sprintf(YCrTalkId_1,"YCrossTalk_1_%s",detUnitLabel);
00844     sprintf(YCrTalkId_2,"YCrossTalk_2_%s",detUnitLabel);
00845     sprintf(XCrTalkDetId_1,"XDetectCrossTalk_1_%s",detUnitLabel);
00846     sprintf(XCrTalkDetId_2,"XDetectCrossTalk_2_%s",detUnitLabel);
00847 
00848     sprintf(effXCrTalkId_1,"XCrossTalkFromDTExtrapolation_1_%s",detUnitLabel);
00849     sprintf(effXCrTalkId_2,"XCrossTalkFromDTExtrapolation_2_%s",detUnitLabel);
00850     sprintf(effXCrTalkDetId_1,"XCrossTalkFromDetectedStrip_1_%s",detUnitLabel);
00851     sprintf(effXCrTalkDetId_2,"XCrossTalkFromDetectedStrip_2_%s",detUnitLabel);
00852     sprintf(effYCrTalkId_1,"YCrossTalkFromDTExtrapolation_1_%s",detUnitLabel);
00853     sprintf(effYCrTalkId_2,"YCrossTalkFromDTExtrapolation_2_%s",detUnitLabel);
00854 
00855     std::map<std::string, MonitorElement*> meMap=meCollection[*meIt];
00856 
00857     for(unsigned int i=1;i<=100;++i){
00858       
00859       if(meMap[meIdDT]->getBinContent(i) != 0){
00860         float eff = meMap[meIdRPC]->getBinContent(i)/meMap[meIdDT]->getBinContent(i);
00861         float erreff = sqrt(eff*(1-eff)/meMap[meIdDT]->getBinContent(i));
00862         meMap[effIdRPC]->setBinContent(i,eff*100.);
00863         meMap[effIdRPC]->setBinError(i,erreff*100.);
00864       }
00865     }
00866     for(unsigned int i=1;i<=100;++i){
00867       for(unsigned int j=1;j<=200;++j){
00868         if(meMap[meIdDT_2D]->getBinContent(i,j) != 0){
00869           float eff = meMap[meIdRPC_2D]->getBinContent(i,j)/meMap[meIdDT_2D]->getBinContent(i,j);
00870           float erreff = sqrt(eff*(1-eff)/meMap[meIdDT_2D]->getBinContent(i,j));
00871           meMap[effIdRPC_2D]->setBinContent(i,j,eff*100.);
00872           meMap[effIdRPC_2D]->setBinError(i,j,erreff*100.);
00873         }
00874       }
00875     }
00876 
00877     //--------------------  CROSS TALK PLOT ---------------------------------------------------
00878 
00879     //-------------------- With predicted STRIP -----------------------------------------------
00880 
00881     for(unsigned int i=1;i<=100;++i){
00882       
00883       if(meMap[meIdDTCrT]->getBinContent(i) != 0){
00884         float crt = meMap[XCrTalkId_1]->getBinContent(i)/meMap[meIdDTCrT]->getBinContent(i);
00885         float errcrt = sqrt(crt*(1-crt)/meMap[meIdDTCrT]->getBinContent(i));
00886         meMap[effXCrTalkId_1]->setBinContent(i,crt*100.);
00887         meMap[effXCrTalkId_1]->setBinError(i,errcrt*100.);
00888       }
00889     }
00890 
00891     for(unsigned int i=1;i<=100;++i){
00892       
00893       if(meMap[meIdDTCrT]->getBinContent(i) != 0){
00894         float crt = meMap[XCrTalkId_2]->getBinContent(i)/meMap[meIdDTCrT]->getBinContent(i);
00895         float errcrt = sqrt(crt*(1-crt)/meMap[meIdDTCrT]->getBinContent(i));
00896         meMap[effXCrTalkId_2]->setBinContent(i,crt*100.);
00897         meMap[effXCrTalkId_2]->setBinError(i,errcrt*100.);
00898       }
00899     }
00900 
00901     //-------------------- With detected STRIP -----------------------------------------------
00902 
00903     for(unsigned int i=1;i<=100;++i){
00904       if(meMap[meIdRPCCrT]->getBinContent(i) != 0){
00905         float crt = meMap[XCrTalkDetId_1]->getBinContent(i)/meMap[meIdRPCCrT]->getBinContent(i);
00906         float errcrt = sqrt(crt*(1-crt)/meMap[meIdRPCCrT]->getBinContent(i));
00907         meMap[effXCrTalkDetId_1]->setBinContent(i,crt*100.);
00908         meMap[effXCrTalkDetId_1]->setBinError(i,errcrt*100.);
00909       }
00910     }
00911 
00912     for(unsigned int i=1;i<=100;++i){
00913       
00914       if(meMap[meIdRPCCrT]->getBinContent(i) != 0){
00915         float crt = meMap[XCrTalkDetId_2]->getBinContent(i)/meMap[meIdRPCCrT]->getBinContent(i);
00916         float errcrt = sqrt(crt*(1-crt)/meMap[meIdRPCCrT]->getBinContent(i));
00917         meMap[effXCrTalkDetId_2]->setBinContent(i,crt*100.);
00918         meMap[effXCrTalkDetId_2]->setBinError(i,errcrt*100.);
00919       }
00920     }
00921 
00922     //-------------------- With Y coordinate -------------------------------------------------
00923 
00924     for(unsigned int i=1;i<=200;++i){
00925       
00926       if(meMap[YCrTalkDTId]->getBinContent(i) != 0){
00927         float crt = meMap[YCrTalkId_1]->getBinContent(i)/meMap[YCrTalkDTId]->getBinContent(i);
00928         float errcrt = sqrt(crt*(1-crt)/meMap[YCrTalkDTId]->getBinContent(i));
00929         meMap[effYCrTalkId_1]->setBinContent(i,crt*100.);
00930         meMap[effYCrTalkId_1]->setBinError(i,errcrt*100.);
00931       }
00932     }
00933 
00934     for(unsigned int i=1;i<=200;++i){
00935       
00936       if(meMap[YCrTalkDTId]->getBinContent(i) != 0){
00937         float crt = meMap[YCrTalkId_2]->getBinContent(i)/meMap[YCrTalkDTId]->getBinContent(i);
00938         float errcrt = sqrt(crt*(1-crt)/meMap[YCrTalkDTId]->getBinContent(i));
00939         meMap[effYCrTalkId_2]->setBinContent(i,crt*100.);
00940         meMap[effYCrTalkId_2]->setBinError(i,errcrt*100.);
00941       }
00942     }
00943   }
00944 
00945   if(EffSaveRootFile) dbe->save(EffRootFileName);
00946   //  theFile->Write();
00947   //  theile->Close();
00948   std::cout<<"End Job"<<std::endl;
00949 }
00950 
00951 RPCMonitorEfficiency::~RPCMonitorEfficiency(){}
00952 
00953 
00954 
00955 
00956 
00957
RPCMonitorEfficiency.cc
