CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoLocalMuon/RPCRecHit/src/CSCSegtoRPC.cc

Go to the documentation of this file.
00001 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
00002 #include <Geometry/CSCGeometry/interface/CSCGeometry.h>
00003 #include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
00004 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
00005 #include <Geometry/Records/interface/MuonGeometryRecord.h>
00006 #include <Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h>
00007 #include <FWCore/Framework/interface/EDAnalyzer.h>
00008 #include "FWCore/Framework/interface/ESHandle.h"
00009 #include <Geometry/RPCGeometry/interface/RPCGeomServ.h>
00010 #include <DataFormats/RPCRecHit/interface/RPCRecHit.h>
00011 #include <DataFormats/RPCRecHit/interface/RPCRecHitCollection.h>
00012 #include <RecoLocalMuon/RPCRecHit/interface/CSCSegtoRPC.h>
00013 
00014 ObjectMapCSC* ObjectMapCSC::mapInstance = NULL;
00015 
00016 ObjectMapCSC* ObjectMapCSC::GetInstance(const edm::EventSetup& iSetup){
00017   if (mapInstance == NULL){
00018     mapInstance = new ObjectMapCSC(iSetup);
00019   }
00020   return mapInstance;
00021 }
00022 
00023 ObjectMapCSC::ObjectMapCSC(const edm::EventSetup& iSetup){
00024   edm::ESHandle<RPCGeometry> rpcGeo;
00025   edm::ESHandle<CSCGeometry> cscGeo;
00026   
00027   iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00028   iSetup.get<MuonGeometryRecord>().get(cscGeo);
00029   
00030   for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00031     if(dynamic_cast< RPCChamber* >( *it ) != 0 ){
00032       RPCChamber* ch = dynamic_cast< RPCChamber* >( *it ); 
00033       std::vector< const RPCRoll*> roles = (ch->rolls());
00034       for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){
00035         RPCDetId rpcId = (*r)->id();
00036         int region=rpcId.region();
00037         if(region!=0){
00038           int station=rpcId.station();
00039           int ring=rpcId.ring();
00040           int cscring=ring;
00041           int cscstation=station;
00042           RPCGeomServ rpcsrv(rpcId);
00043           int rpcsegment = rpcsrv.segment();
00044           int cscchamber = rpcsegment; //FIX THIS ACCORDING TO RPCGeomServ::segment()Definition
00045           if((station==2||station==3)&&ring==3){//Adding Ring 3 of RPC to the CSC Ring 2
00046             cscring = 2;
00047           }
00048           CSCStationIndex ind(region,cscstation,cscring,cscchamber);
00049           std::set<RPCDetId> myrolls;
00050           if (rollstoreCSC.find(ind)!=rollstoreCSC.end()) myrolls=rollstoreCSC[ind];
00051           myrolls.insert(rpcId);
00052           rollstoreCSC[ind]=myrolls;
00053         }
00054       }
00055     }
00056   }
00057 }
00058   
00059 CSCSegtoRPC::CSCSegtoRPC(edm::Handle<CSCSegmentCollection> allCSCSegments, const edm::EventSetup& iSetup,const edm::Event& iEvent, bool debug, double eyr){
00060   
00061   edm::ESHandle<RPCGeometry> rpcGeo;
00062   edm::ESHandle<CSCGeometry> cscGeo;
00063   
00064   iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00065   iSetup.get<MuonGeometryRecord>().get(cscGeo);
00066   
00067   MaxD=80.;
00068 
00069   if(debug) std::cout<<"CSC \t Number of CSC Segments in this event = "<<allCSCSegments->size()<<std::endl;
00070 
00071   _ThePoints = new RPCRecHitCollection();
00072 
00073   if(allCSCSegments->size()==0){
00074     if(debug) std::cout<<"CSC 0 segments skiping event"<<std::endl;
00075   }else {
00076     std::map<CSCDetId,int> CSCSegmentsCounter;
00077     CSCSegmentCollection::const_iterator segment;
00078       
00079     int segmentsInThisEventInTheEndcap=0;
00080       
00081     for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
00082       CSCSegmentsCounter[segment->cscDetId()]++;
00083       segmentsInThisEventInTheEndcap++;
00084     }    
00085       
00086     if(debug) std::cout<<"CSC \t loop over all the CSCSegments "<<std::endl;
00087     for (segment = allCSCSegments->begin();segment!=allCSCSegments->end(); ++segment){
00088       CSCDetId CSCId = segment->cscDetId();
00089         
00090       if(debug) std::cout<<"CSC \t \t This Segment is in Chamber id: "<<CSCId<<std::endl;
00091       if(debug) std::cout<<"CSC \t \t Number of segments in this CSC = "<<CSCSegmentsCounter[CSCId]<<std::endl;
00092       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;
00093 
00094       if(CSCSegmentsCounter[CSCId]==1 && CSCId.station()!=4 && CSCId.ring()!=1 && allCSCSegments->size()>=2){
00095         if(debug) std::cout<<"CSC \t \t yes"<<std::endl;
00096         int cscEndCap = CSCId.endcap();
00097         int cscStation = CSCId.station();
00098         int cscRing = CSCId.ring();
00099         int rpcRegion = 1; if(cscEndCap==2) rpcRegion= -1;//Relacion entre las endcaps
00100         int rpcRing = cscRing;
00101         if(cscRing==4)rpcRing =1;
00102         int rpcStation = cscStation;
00103         int rpcSegment = CSCId.chamber();
00104         
00105         LocalPoint segmentPosition= segment->localPosition();
00106         LocalVector segmentDirection=segment->localDirection();
00107         float dz=segmentDirection.z();
00108 
00109         if(debug) std::cout<<"CSC \t \t \t Information about the segment" 
00110                            <<"RecHits ="<<segment->nRecHits()
00111                            <<"Angle ="<<acos(dz)*180/3.1415926<<std::endl;
00112                       
00113         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;
00114 
00115         if((segment->dimension()==4) && (segment->nRecHits()<=10 && segment->nRecHits()>=4)){
00116           //&& acos(dz)*180/3.1415926 > 45. && acos(dz)*180/3.1415926 < 135.){ 
00117           //&& segment->chi2()< ??)Add 3 segmentes in the endcaps???
00118 
00119 
00120           if(debug) std::cout<<"CSC \t \t yes"<<std::endl;
00121           if(debug) std::cout<<"CSC \t \t CSC Segment Dimension "<<segment->dimension()<<std::endl; 
00122             
00123           float Xo=segmentPosition.x();
00124           float Yo=segmentPosition.y();
00125           float Zo=segmentPosition.z();
00126           float dx=segmentDirection.x();
00127           float dy=segmentDirection.y();
00128           float dz=segmentDirection.z();
00129 
00130           if(debug)  std::cout<<"Calling to Object Map class"<<std::endl;
00131           ObjectMapCSC* TheObjectCSC = ObjectMapCSC::GetInstance(iSetup);
00132           if(debug) std::cout<<"Creating the CSCIndex"<<std::endl;
00133           CSCStationIndex theindex(rpcRegion,rpcStation,rpcRing,rpcSegment);
00134           if(debug) std::cout<<"Getting the Rolls for the given index"<<std::endl;
00135         
00136           std::set<RPCDetId> rollsForThisCSC = TheObjectCSC->GetInstance(iSetup)->GetRolls(theindex);
00137                 
00138            
00139           if(debug) std::cout<<"CSC \t \t Getting chamber from Geometry"<<std::endl;
00140           const CSCChamber* TheChamber=cscGeo->chamber(CSCId); 
00141           if(debug) std::cout<<"CSC \t \t Getting ID from Chamber"<<std::endl;
00142           const CSCDetId TheId=TheChamber->id();
00143 
00144           if(debug) std::cout<<"CSC \t \t Number of rolls for this CSC = "<<rollsForThisCSC.size()<<std::endl;
00145 
00146           if(debug) std::cout<<"CSC \t \t Printing The Id"<<TheId<<std::endl;
00147 
00148           if(rpcRing!=1&&rpcStation!=4){//They don't exist!
00149           
00150             assert(rollsForThisCSC.size()>=1);
00151 
00152             if(debug) std::cout<<"CSC \t \t Loop over all the rolls asociated to this CSC"<<std::endl;      
00153             for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisCSC.begin();iteraRoll != rollsForThisCSC.end(); iteraRoll++){
00154               const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00155               RPCDetId rpcId = rollasociated->id();
00156                 
00157               if(debug) std::cout<<"CSC \t \t \t We are in the roll getting the surface"<<rpcId<<std::endl;
00158               const BoundPlane & RPCSurface = rollasociated->surface(); 
00159 
00160               if(debug) std::cout<<"CSC \t \t \t RollID: "<<rpcId<<std::endl;
00161                 
00162               if(debug) std::cout<<"CSC \t \t \t Doing the extrapolation to this roll"<<std::endl;
00163               if(debug) std::cout<<"CSC \t \t \t CSC Segment Direction in CSCLocal "<<segmentDirection<<std::endl;
00164               if(debug) std::cout<<"CSC \t \t \t CSC Segment Point in CSCLocal "<<segmentPosition<<std::endl;  
00165                 
00166               GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
00167               if(debug) std::cout<<"CSC \t \t \t Center (0,0,0) of the Roll in Global"<<CenterPointRollGlobal<<std::endl;
00168               GlobalPoint CenterPointCSCGlobal = TheChamber->toGlobal(LocalPoint(0,0,0));
00169               if(debug) std::cout<<"CSC \t \t \t Center (0,0,0) of the CSC in Global"<<CenterPointCSCGlobal<<std::endl;
00170               GlobalPoint segmentPositionInGlobal=TheChamber->toGlobal(segmentPosition); //new way to convert to global
00171               if(debug) std::cout<<"CSC \t \t \t Segment Position in Global"<<segmentPositionInGlobal<<std::endl;
00172               LocalPoint CenterRollinCSCFrame = TheChamber->toLocal(CenterPointRollGlobal);
00173 
00174               if(debug){//to check CSC RPC phi relation!
00175                 float rpcphi=0;
00176                 float cscphi=0;
00177                   
00178                 (CenterPointRollGlobal.barePhi()<0)? 
00179                   rpcphi = 2*3.141592+CenterPointRollGlobal.barePhi():rpcphi=CenterPointRollGlobal.barePhi();
00180                   
00181                 (CenterPointCSCGlobal.barePhi()<0)? 
00182                   cscphi = 2*3.1415926536+CenterPointCSCGlobal.barePhi():cscphi=CenterPointCSCGlobal.barePhi();
00183 
00184                 float df=fabs(cscphi-rpcphi); 
00185                 float dr=fabs(CenterPointRollGlobal.perp()-CenterPointCSCGlobal.perp());
00186                 float diffz=CenterPointRollGlobal.z()-CenterPointCSCGlobal.z();
00187                 float dfg=df*180./3.14159265;
00188 
00189                 if(debug) std::cout<<"CSC \t \t \t z of RPC="<<CenterPointRollGlobal.z()<<"z of CSC"<<CenterPointCSCGlobal.z()<<" dfg="<<dfg<<std::endl;
00190                   
00191                 RPCGeomServ rpcsrv(rpcId);
00192                   
00193                 if(dr>200.||fabs(dz)>55.||dfg>1.){ 
00194                   //if(rpcRegion==1&&dfg>1.&&dr>100.){  
00195                   if (debug) std::cout
00196                     <<"\t \t \t CSC Station= "<<CSCId.station()
00197                     <<" Ring= "<<CSCId.ring()
00198                     <<" Chamber= "<<CSCId.chamber()
00199                     <<" cscphi="<<cscphi*180/3.14159265
00200                     <<"\t RPC Station= "<<rpcId.station()
00201                     <<" ring= "<<rpcId.ring()
00202                     <<" segment =-> "<<rpcsrv.segment()
00203                     <<" rollphi="<<rpcphi*180/3.14159265
00204                     <<"\t dfg="<<dfg
00205                     <<" dz="<<diffz
00206                     <<" dr="<<dr
00207                     <<std::endl;
00208                     
00209                 }
00210               }
00211               
00212               float D=CenterRollinCSCFrame.z();
00213                   
00214               float X=Xo+dx*D/dz;
00215               float Y=Yo+dy*D/dz;
00216               float Z=D;
00217 
00218               const TrapezoidalStripTopology* top_=dynamic_cast<const TrapezoidalStripTopology*>(&(rollasociated->topology()));
00219               LocalPoint xmin = top_->localPosition(0.);
00220               if(debug) std::cout<<"CSC \t \t \t xmin of this  Roll "<<xmin<<"cm"<<std::endl;
00221               LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00222               if(debug) std::cout<<"CSC \t \t \t xmax of this  Roll "<<xmax<<"cm"<<std::endl;
00223               float rsize = fabs( xmax.x()-xmin.x() );
00224               if(debug) std::cout<<"CSC \t \t \t Roll Size "<<rsize<<"cm"<<std::endl;
00225               float stripl = top_->stripLength();
00226               float stripw = top_->pitch();
00227 
00228               if(debug) std::cout<<"CSC \t \t \t Strip Lenght "<<stripl<<"cm"<<std::endl;
00229               if(debug) std::cout<<"CSC \t \t \t Strip Width "<<stripw<<"cm"<<std::endl;
00230 
00231               if(debug) std::cout<<"CSC \t \t \t X Predicted in CSCLocal= "<<X<<"cm"<<std::endl;
00232               if(debug) std::cout<<"CSC \t \t \t Y Predicted in CSCLocal= "<<Y<<"cm"<<std::endl;
00233               if(debug) std::cout<<"CSC \t \t \t Z Predicted in CSCLocal= "<<Z<<"cm"<<std::endl;
00234           
00235               float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
00236 
00237               if(debug) std::cout<<"CSC \t \t \t Is the distance of extrapolation less than MaxD? ="<<extrapolatedDistance<<"cm"<<" MaxD="<<MaxD<<"cm"<<std::endl;
00238           
00239               if(extrapolatedDistance<=MaxD){ 
00240 
00241                 if(debug) std::cout<<"CSC \t \t \t yes"<<std::endl;
00242 
00243                 GlobalPoint GlobalPointExtrapolated=TheChamber->toGlobal(LocalPoint(X,Y,Z));
00244                 if(debug) std::cout<<"CSC \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
00245 
00246               
00247                 LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
00248 
00249                 if(debug) std::cout<<"CSC \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
00250                 if(debug) std::cout<<"CSC \t \t \t Corner of the Roll = ("<<rsize*eyr<<","<<stripl*eyr<<")"<<std::endl;
00251                 if(debug) std::cout<<"CSC \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.x())<<","
00252                                    <<fabs(PointExtrapolatedRPCFrame.y())<<","<<fabs(PointExtrapolatedRPCFrame.z())<<")"<<std::endl;
00253                 if(debug) std::cout<<"CSC \t \t \t dz="
00254                                    <<fabs(PointExtrapolatedRPCFrame.z())<<" dx="
00255                                    <<fabs(PointExtrapolatedRPCFrame.x())<<" dy="
00256                                    <<fabs(PointExtrapolatedRPCFrame.y())<<std::endl;
00257                   
00258                 if(debug) std::cout<<"CSC \t \t \t Does the extrapolation go inside this roll????"<<std::endl;
00259 
00260                 if(fabs(PointExtrapolatedRPCFrame.z()) < 1. && 
00261                    fabs(PointExtrapolatedRPCFrame.x()) < rsize*eyr && 
00262                    fabs(PointExtrapolatedRPCFrame.y()) < stripl*eyr){ 
00263                   if(debug) std::cout<<"CSC \t \t \t \t yes"<<std::endl;
00264                   if(debug) std::cout<<"CSC \t \t \t \t Creating the RecHit"<<std::endl;
00265                   RPCRecHit RPCPoint(rpcId,0,PointExtrapolatedRPCFrame);
00266                   if(debug) std::cout<<"CSC \t \t \t \t Clearing the vector"<<std::endl;        
00267                   RPCPointVector.clear();
00268                   if(debug) std::cout<<"CSC \t \t \t \t Pushing back"<<std::endl;       
00269                   RPCPointVector.push_back(RPCPoint); 
00270                   if(debug) std::cout<<"CSC \t \t \t \t Putting the vector"<<std::endl; 
00271                   _ThePoints->put(rpcId,RPCPointVector.begin(),RPCPointVector.end());
00272                 }
00273               }
00274             }
00275           }
00276         }
00277       }
00278     }
00279   }
00280 }
00281 
00282 CSCSegtoRPC::~CSCSegtoRPC(){
00283 
00284 }