Reference Manual

Go to the documentation of this file.
00001 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
00002 #include <Geometry/DTGeometry/interface/DTGeometry.h>
00003 #include <DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h>
00004 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
00005 #include <Geometry/Records/interface/MuonGeometryRecord.h>
00006 #include <Geometry/CommonTopologies/interface/RectangularStripTopology.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/DTSegtoRPC.h>
00013 #include <ctime>
00014 
00015 ObjectMap* ObjectMap::mapInstance = NULL;
00016 
00017 ObjectMap* ObjectMap::GetInstance(const edm::EventSetup& iSetup){
00018   if (mapInstance == NULL){
00019     mapInstance = new ObjectMap(iSetup);
00020   }
00021   return mapInstance;
00022 }
00023 
00024 ObjectMap::ObjectMap(const edm::EventSetup& iSetup){
00025   edm::ESHandle<RPCGeometry> rpcGeo;
00026   edm::ESHandle<DTGeometry> dtGeo;
00027   
00028   iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00029   iSetup.get<MuonGeometryRecord>().get(dtGeo);
00030   
00031   for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00032     if(dynamic_cast< RPCChamber* >( *it ) != 0 ){
00033       RPCChamber* ch = dynamic_cast< RPCChamber* >( *it ); 
00034       std::vector< const RPCRoll*> roles = (ch->rolls());
00035       for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){
00036         RPCDetId rpcId = (*r)->id();
00037         int region=rpcId.region();
00038         if(region==0){
00039           int wheel=rpcId.ring();
00040           int sector=rpcId.sector();
00041           int station=rpcId.station();
00042           DTStationIndex ind(region,wheel,sector,station);
00043           std::set<RPCDetId> myrolls;
00044           if (rollstoreDT.find(ind)!=rollstoreDT.end()) myrolls=rollstoreDT[ind];
00045           myrolls.insert(rpcId);
00046           rollstoreDT[ind]=myrolls;
00047         }
00048       }
00049     }
00050   }
00051 }
00052 
00053 int distsector(int sector1,int sector2){
00054   if(sector1==13) sector1=4;
00055   if(sector1==14) sector1=10;
00056   
00057   if(sector2==13) sector2=4;
00058   if(sector2==14) sector2=10;
00059   
00060   int distance = std::abs(sector1 - sector2);
00061   if(distance>6) distance = 12-distance;
00062   return distance;
00063 }
00064 
00065 int distwheel(int wheel1,int wheel2){
00066   int distance = std::abs(wheel1 - wheel2);
00067   return distance;
00068 }
00069 
00070 DTSegtoRPC::DTSegtoRPC(edm::Handle<DTRecSegment4DCollection> all4DSegments, const edm::EventSetup& iSetup,const edm::Event& iEvent,bool debug,double eyr){
00071 
00072   /*
00073   MinCosAng=iConfig.getUntrackedParameter<double>("MinCosAng",0.95);
00074   MaxD=iConfig.getUntrackedParameter<double>("MaxD",80.);
00075   MaxDrb4=iConfig.getUntrackedParameter<double>("MaxDrb4",150.);
00076   */
00077   incldt=true;
00078   incldtMB4=true;
00079  
00080   //By now hard coded parameters
00081   MinCosAng=0.85;
00082   MaxD=80.;
00083   MaxDrb4=150.;
00084   MaxDistanceBetweenSegments = 150;
00085   /*
00086 
00087   //These should be always true expect for debuggin porpouses
00088   incldt=true;
00089   incldtMB4=true;
00090 
00091   
00092     struct timespec start_time, stop_time;
00093     time_t fs;
00094     time_t fn;
00095     time_t ls;
00096     time_t ln;
00097     clock_gettime(CLOCK_REALTIME, &start_time);
00098   */
00099 
00100   _ThePoints = new RPCRecHitCollection();
00101 
00102   if(all4DSegments->size()>8){
00103     if(debug) std::cout<<"Too many segments in this event we are not doing the extrapolation"<<std::endl;
00104   }else{ 
00105     edm::ESHandle<RPCGeometry> rpcGeo;
00106     edm::ESHandle<DTGeometry> dtGeo;
00107   
00108     iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00109     iSetup.get<MuonGeometryRecord>().get(dtGeo);
00110 
00111     /*
00112       clock_gettime(CLOCK_REALTIME, &stop_time);
00113       fs=start_time.tv_sec;
00114       fn=start_time.tv_nsec;
00115       ls=stop_time.tv_sec;
00116       ln=stop_time.tv_nsec;
00117       std::cout <<" =================|| "<<ls-fs<<" sec "<<ln-fn<<" us"<<std::endl;
00118       clock_gettime(CLOCK_REALTIME, &start_time);
00119     */
00120   
00121     std::map<DTChamberId,int> DTSegmentCounter;
00122     DTRecSegment4DCollection::const_iterator segment;  
00123   
00124     for (segment = all4DSegments->begin();segment!=all4DSegments->end(); ++segment){
00125       DTSegmentCounter[segment->chamberId()]++;
00126     }
00127   
00128 
00129     /*
00130       clock_gettime(CLOCK_REALTIME, &stop_time);
00131       fs=start_time.tv_sec;
00132       fn=start_time.tv_nsec;
00133       ls=stop_time.tv_sec;
00134       ln=stop_time.tv_nsec;
00135       if(debug) std::cout <<" =================||| "<<ls-fs<<" sec "<<ln-fn<<" us"<<std::endl;
00136       clock_gettime(CLOCK_REALTIME, &start_time);
00137     */
00138 
00139     if(incldt){
00140       for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){
00141       
00142         DTChamberId DTId = segment->chamberId();
00143       
00144         if(debug) std::cout<<"DT  \t \t This Segment is in Chamber id: "<<DTId<<std::endl;
00145         if(debug) std::cout<<"DT  \t \t Number of segments in this DT = "<<DTSegmentCounter[DTId]<<std::endl;
00146         if(debug) std::cout<<"DT  \t \t Is the only one in this DT? and is not in the 4th Station?"<<std::endl;
00147       
00148         if(DTSegmentCounter[DTId]!=1 || DTId.station()==4){
00149           if(debug) std::cout<<"DT \t \t More than one segment in this chamber, or we are in Station 4"<<std::endl;
00150           continue;
00151         }
00152       
00153         int dtWheel = DTId.wheel();
00154         int dtStation = DTId.station();
00155         int dtSector = DTId.sector();      
00156       
00157         LocalPoint segmentPosition= segment->localPosition();
00158         LocalVector segmentDirection=segment->localDirection();
00159       
00160         const GeomDet* gdet=dtGeo->idToDet(segment->geographicalId());
00161         const BoundPlane & DTSurface = gdet->surface();
00162       
00163         //check if the dimension of the segment is 4 
00164       
00165         if(debug) std::cout<<"DT  \t \t Is the segment 4D?"<<std::endl;
00166       
00167         if(segment->dimension()!=4){
00168           if(debug) std::cout<<"DT  \t \t no"<<std::endl;
00169           continue; 
00170         }
00171       
00172         if(debug) std::cout<<"DT  \t \t yes"<<std::endl;
00173         if(debug) std::cout<<"DT  \t \t DT Segment Dimension "<<segment->dimension()<<std::endl; 
00174       
00175         float Xo=segmentPosition.x();
00176         float Yo=segmentPosition.y();
00177         float Zo=segmentPosition.z();
00178         float dx=segmentDirection.x();
00179         float dy=segmentDirection.y();
00180         float dz=segmentDirection.z();
00181       
00182         if(debug)  std::cout<<"Calling to Object Map class"<<std::endl;
00183         ObjectMap* TheObject = ObjectMap::GetInstance(iSetup);
00184         if(debug) std::cout<<"Creating the DTIndex"<<std::endl;
00185         DTStationIndex theindex(0,dtWheel,dtSector,dtStation);
00186         if(debug) std::cout<<"Getting the Rolls for the given index"<<std::endl;
00187       
00188         std::set<RPCDetId> rollsForThisDT = TheObject->GetInstance(iSetup)->GetRolls(theindex);
00189       
00190         if(debug) std::cout<<"DT  \t \t Number of rolls for this DT = "<<rollsForThisDT.size()<<std::endl;
00191       
00192         assert(rollsForThisDT.size()>=1);
00193       
00194         if(debug) std::cout<<"DT  \t \t Loop over all the rolls asociated to this DT"<<std::endl;
00195         for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00196           const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00197           RPCDetId rpcId = rollasociated->id();
00198           const BoundPlane & RPCSurface = rollasociated->surface(); 
00199         
00200           RPCGeomServ rpcsrv(rpcId);
00201           std::string nameRoll = rpcsrv.name();
00202         
00203           if(debug) std::cout<<"DT  \t \t \t RollName: "<<nameRoll<<std::endl;
00204           if(debug) std::cout<<"DT  \t \t \t Doing the extrapolation to this roll"<<std::endl;
00205           if(debug) std::cout<<"DT  \t \t \t DT Segment Direction in DTLocal "<<segmentDirection<<std::endl;
00206           if(debug) std::cout<<"DT  \t \t \t DT Segment Point in DTLocal "<<segmentPosition<<std::endl;
00207         
00208           GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0,0,0));
00209         
00210           LocalPoint CenterRollinDTFrame = DTSurface.toLocal(CenterPointRollGlobal);
00211         
00212           if(debug) std::cout<<"DT  \t \t \t Center (0,0,0) Roll In DTLocal"<<CenterRollinDTFrame<<std::endl;
00213           if(debug) std::cout<<"DT  \t \t \t Center (0,0,0) of the Roll in Global"<<CenterPointRollGlobal<<std::endl;
00214         
00215           float D=CenterRollinDTFrame.z();
00216         
00217           float X=Xo+dx*D/dz;
00218           float Y=Yo+dy*D/dz;
00219           float Z=D;
00220         
00221           const RectangularStripTopology* top_= dynamic_cast<const RectangularStripTopology*> (&(rollasociated->topology()));
00222           LocalPoint xmin = top_->localPosition(0.);
00223           if(debug) std::cout<<"DT  \t \t \t xmin of this  Roll "<<xmin<<"cm"<<std::endl;
00224           LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00225           if(debug) std::cout<<"DT  \t \t \t xmax of this  Roll "<<xmax<<"cm"<<std::endl;
00226           float rsize = fabs( xmax.x()-xmin.x() );
00227           if(debug) std::cout<<"DT  \t \t \t Roll Size "<<rsize<<"cm"<<std::endl;
00228           float stripl = top_->stripLength();
00229         
00230           float stripw = top_->pitch();
00231         
00232           if(debug) std::cout<<"DT  \t \t \t Strip Lenght "<<stripl<<"cm"<<std::endl;
00233           if(debug) std::cout<<"DT  \t \t \t Strip Width "<<stripw<<"cm"<<std::endl;
00234           if(debug) std::cout<<"DT  \t \t \t X Predicted in DTLocal= "<<X<<"cm"<<std::endl;
00235           if(debug) std::cout<<"DT  \t \t \t Y Predicted in DTLocal= "<<Y<<"cm"<<std::endl;
00236           if(debug) std::cout<<"DT  \t \t \t Z Predicted in DTLocal= "<<Z<<"cm"<<std::endl;
00237         
00238           float extrapolatedDistance = sqrt((X-Xo)*(X-Xo)+(Y-Yo)*(Y-Yo)+(Z-Zo)*(Z-Zo));
00239         
00240           if(debug) std::cout<<"DT  \t \t \t Is the distance of extrapolation less than MaxD? ="<<extrapolatedDistance<<"cm"<<"MaxD="<<MaxD<<"cm"<<std::endl;
00241         
00242           if(extrapolatedDistance<=MaxD){ 
00243             if(debug) std::cout<<"DT  \t \t \t yes"<<std::endl;   
00244             GlobalPoint GlobalPointExtrapolated = DTSurface.toGlobal(LocalPoint(X,Y,Z));
00245             if(debug) std::cout<<"DT  \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
00246             LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
00247           
00248             if(debug) std::cout<<"DT  \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
00249             if(debug) std::cout<<"DT  \t \t \t Corner of the Roll = ("<<rsize*eyr<<","<<stripl*eyr<<")"<<std::endl;
00250             if(debug) std::cout<<"DT \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.x())<<","
00251                                <<fabs(PointExtrapolatedRPCFrame.y())<<","<<fabs(PointExtrapolatedRPCFrame.z())<<")"<<std::endl;
00252             if(debug) std::cout<<"DT  \t \t \t Does the extrapolation go inside this roll?"<<std::endl;
00253           
00254             if(fabs(PointExtrapolatedRPCFrame.z()) < 1. && 
00255                fabs(PointExtrapolatedRPCFrame.x()) < rsize*eyr && 
00256                fabs(PointExtrapolatedRPCFrame.y()) < stripl*eyr){
00257               if(debug) std::cout<<"DT  \t \t \t \t yes"<<std::endl;    
00258               if(debug) std::cout<<"DT  \t \t \t \t Creating the RecHit"<<std::endl;    
00259               RPCRecHit RPCPoint(rpcId,0,PointExtrapolatedRPCFrame);
00260               if(debug) std::cout<<"DT  \t \t \t \t Clearing the vector"<<std::endl;    
00261               RPCPointVector.clear();
00262               if(debug) std::cout<<"DT  \t \t \t \t Pushing back"<<std::endl;   
00263               RPCPointVector.push_back(RPCPoint); 
00264               if(debug) std::cout<<"DT  \t \t \t \t Putting the vector"<<std::endl;     
00265               _ThePoints->put(rpcId,RPCPointVector.begin(),RPCPointVector.end());
00266 
00267               if(debug) std::cout<<"DT \t \t \t \t Filling container with "<<nameRoll
00268                                  <<" Point.x="<<PointExtrapolatedRPCFrame.x()<<" Point.y="<<PointExtrapolatedRPCFrame.y()<<" size="<<RPCPointVector.size()<<std::endl;
00269             
00270             }else {
00271               if(debug) std::cout<<"DT \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00272             }//Condition for the right match
00273           }else{
00274             if(debug) std::cout<<"DT \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00275           }//D so big
00276         }//loop over all the rolls asociated
00277       }
00278     }
00279   
00280     if(incldtMB4){
00281       if(all4DSegments->size()>0){     
00282         if(debug) std::cout<<"MB4 \t \t Loop Over all4DSegments "<<all4DSegments->size()<<std::endl;
00283         extrapolatedRolls.clear();
00284         for (segment = all4DSegments->begin(); segment != all4DSegments->end(); ++segment){ 
00285     
00286           DTChamberId DTId = segment->chamberId();
00287         
00288           if(debug) std::cout<<"MB4 \t \t This Segment is in Chamber id: "<<DTId<<std::endl;
00289           if(debug) std::cout<<"MB4 \t \t Number of segments in this DT = "<<DTSegmentCounter[DTId]<<std::endl;
00290           if(debug) std::cout<<"MB4 \t \t \t Is the only one in this DT? and is in the Station 4?"<<std::endl;
00291 
00292           if(DTSegmentCounter[DTId] == 1 && DTId.station()==4){
00293 
00294             if(debug) std::cout<<"MB4 \t \t \t yes"<<std::endl;
00295             int dtWheel = DTId.wheel();
00296             int dtStation = DTId.station();
00297             int dtSector = DTId.sector();
00298       
00299             LocalPoint segmentPosition= segment->localPosition();
00300             LocalVector segmentDirection=segment->localDirection();
00301             
00302             if(debug) std::cout<<"MB4 \t \t \t \t The Segment in MB4 is 2D?"<<std::endl;
00303             if(segment->dimension()==2){
00304               if(debug) std::cout<<"MB4 \t \t \t \t yes"<<std::endl;
00305               LocalVector segmentDirectionMB4=segmentDirection;
00306               LocalPoint segmentPositionMB4=segmentPosition;
00307             
00308               const BoundPlane& DTSurface4 = dtGeo->idToDet(DTId)->surface();
00309                 
00310               DTRecSegment4DCollection::const_iterator segMB3;  
00311                 
00312               if(debug) std::cout<<"MB4 \t \t \t \t Loop on segments in =sector && MB3 && adjacent sectors && y dim=4"<<std::endl;
00313               for(segMB3=all4DSegments->begin();segMB3!=all4DSegments->end();++segMB3){
00314                     
00315                 DTChamberId dtid3 = segMB3->chamberId();  
00316                 
00317                 if(debug) std::cout<<"MB4  \t \t \t \t Segment in Chamber ="<<dtid3<<std::endl;
00318 
00319                     
00320                 if(distsector(dtid3.sector(),DTId.sector())<=1 //The DT sector could be 13 or 14 and because is corrected in the calculation of the distance.
00321                    && distwheel(dtid3.wheel(),DTId.wheel())<=1 //The we could have segments in neighbohr wheels in pp collisions 
00322                    && dtid3.station()==3
00323                    && DTSegmentCounter[dtid3] == 1
00324                    && segMB3->dimension()==4){
00325                   
00326                   if(debug) std::cout<<"MB4  \t \t \t \t distsector ="<<distsector(dtid3.sector(),DTId.sector())<<std::endl;
00327                   if(debug) std::cout<<"MB4  \t \t \t \t distwheel ="<<distwheel(dtid3.wheel(),DTId.wheel())<<std::endl;
00328 
00329                   const GeomDet* gdet3=dtGeo->idToDet(segMB3->geographicalId());
00330                   const BoundPlane & DTSurface3 = gdet3->surface();
00331 
00332                   LocalVector segmentDirectionMB3 =  segMB3->localDirection();
00333                   GlobalPoint segmentPositionMB3inGlobal = DTSurface3.toGlobal(segMB3->localPosition());
00334                   GlobalPoint segmentPositionMB4inGlobal = DTSurface4.toGlobal(segmentPosition);
00335                 
00336                   //LocalVector segDirMB4inMB3Frame=DTSurface3.toLocal(DTSurface4.toGlobal(segmentDirectionMB4));
00337                   LocalVector segDirMB3inMB4Frame=DTSurface4.toLocal(DTSurface3.toGlobal(segmentDirectionMB3));
00338                 
00339                   GlobalVector segDirMB4inGlobalFrame=DTSurface4.toGlobal(segmentDirectionMB4);
00340                   GlobalVector segDirMB3inGlobalFrame=DTSurface3.toGlobal(segmentDirectionMB3);
00341                 
00342                   float dx=segDirMB4inGlobalFrame.x();
00343                   float dy=segDirMB4inGlobalFrame.y();
00344                                 
00345                   float dx3=segDirMB3inGlobalFrame.x();
00346                   float dy3=segDirMB3inGlobalFrame.y();
00347                                 
00348                   double cosAng=fabs(dx*dx3+dy*dy3/sqrt((dx3*dx3+dy3*dy3)*(dx*dx+dy*dy)));
00349 
00350                   if(debug) std::cout<<"MB4 \t \t \t \t cosAng"<<cosAng<<"Beetween "<<dtid3<<" and "<<DTId<<std::endl;
00351                 
00352                   if(debug){
00353                     std::cout<<"MB4 \t \t \t \t dx="<<dx<<" dy="<<dy<<std::endl;
00354                     std::cout<<"MB4 \t \t \t \t dx3="<<dx3<<" dy3="<<dy<<std::endl;
00355                     std::cout<<"MB4 \t \t \t \t cosAng="<<cosAng<<std::endl;
00356                   }
00357 
00358                   float DistanceBetweenSegments = ((segmentPositionMB3inGlobal) - (segmentPositionMB4inGlobal)).mag();
00359                   
00360                   if(cosAng>MinCosAng && DistanceBetweenSegments < MaxDistanceBetweenSegments){
00361                     
00362                     if(debug) std::cout<<"MB4 \t \t \t \t Distance between segments="<<DistanceBetweenSegments<<std::endl;
00363                   
00364                     if(debug) std::cout<<"MB4 \t \t We found compatible Segments (similar direction and close enough) in "<<dtid3<<" and "<<DTId<<std::endl;
00365 
00366                     if(dtSector==13){
00367                       dtSector=4;
00368                     }
00369                     if(dtSector==14){
00370                       dtSector=10;
00371                     }
00372                    
00373                     if(debug)  std::cout<<"Calling to Object Map class"<<std::endl;
00374                     ObjectMap* TheObject = ObjectMap::GetInstance(iSetup);
00375                     if(debug) std::cout<<"Creating the DTIndex"<<std::endl;
00376                     DTStationIndex theindex(0,dtWheel,dtSector,dtStation);
00377                     if(debug) std::cout<<"Getting the Rolls for the given index"<<std::endl;
00378 
00379                     std::set<RPCDetId> rollsForThisDT = TheObject->GetInstance(iSetup)->GetRolls(theindex);
00380 
00381                     if(debug) std::cout<<"MB4 \t \t Number of rolls for this DT = "<<rollsForThisDT.size()<<std::endl;
00382                     
00383                     assert(rollsForThisDT.size()>=1);
00384                         
00385                     if(debug) std::cout<<"MB4  \t \t Loop over all the rolls asociated to this DT"<<std::endl;
00386                     for (std::set<RPCDetId>::iterator iteraRoll=rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00387                       const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll); //roll asociado a MB4
00388                       RPCDetId rpcId = rollasociated->id();
00389                       const BoundPlane & RPCSurfaceRB4 = rollasociated->surface(); //surface MB4
00390 
00391                       RPCGeomServ rpcsrv(rpcId);
00392                       std::string nameRoll = rpcsrv.name();
00393 
00394                       if(debug) std::cout<<"MB4  \t \t \t RollName: "<<nameRoll<<std::endl;
00395                       if(debug) std::cout<<"MB4  \t \t \t Doing the extrapolation to this roll"<<std::endl;
00396                         
00397                       GlobalPoint CenterPointRollGlobal=RPCSurfaceRB4.toGlobal(LocalPoint(0,0,0));
00398                       LocalPoint CenterRollinMB4Frame = DTSurface4.toLocal(CenterPointRollGlobal); //In MB4
00399                       LocalPoint segmentPositionMB3inMB4Frame = DTSurface4.toLocal(segmentPositionMB3inGlobal); //In MB4
00400                       //LocalPoint segmentPositionMB3inRB4Frame = RPCSurfaceRB4.toLocal(segmentPositionMB3inGlobal); //In MB4
00401                       LocalVector segmentDirectionMB3inMB4Frame = DTSurface4.toLocal(segDirMB3inGlobalFrame); //In MB4
00402                         
00403                       //The exptrapolation is done in MB4 frame. for local x and z is done from MB4,
00404                       float Dxz=CenterRollinMB4Frame.z();
00405                       float Xo4=segmentPositionMB4.x();
00406                       float dxl=segmentDirectionMB4.x(); //dx local for MB4 segment in MB4 Frame
00407                       float dzl=segmentDirectionMB4.z(); //dx local for MB4 segment in MB4 Frame
00408                         
00409                       float X=Xo4+dxl*Dxz/dzl; //In MB4 frame
00410                       float Z=Dxz;//In MB4 frame
00411                         
00412                       //for local y is done from MB3 
00413                       float Yo34 = segmentPositionMB3inMB4Frame.y();
00414                       float dy34 = segmentDirectionMB3inMB4Frame.y();
00415                       float dz34 = segmentDirectionMB3inMB4Frame.z();
00416                       float Dy=Dxz-(segmentPositionMB3inMB4Frame.z()); //Distance beetween the segment in MB3 and the RB4 surface
00417 
00418                       if(debug) std::cout<<"MB4 \t \t \t The distance to extrapolate in Y from MB3 is "<<Dy<<"cm"<<std::endl;
00419                         
00420                       float Y=Yo34+dy34*Dy/dz34;//In MB4 Frame
00421                           
00422                       const RectangularStripTopology* top_
00423                         =dynamic_cast<const RectangularStripTopology*>(&(rollasociated->topology())); //Topology roll asociated MB4
00424                       LocalPoint xmin = top_->localPosition(0.);
00425                       LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00426                       float rsize = fabs( xmax.x()-xmin.x() );
00427                       float stripl = top_->stripLength();
00428                       float stripw = top_->pitch();
00429 
00430                             
00431                       if(debug) std::cout<<"MB4 \t \t \t Strip Lenght "<<stripl<<"cm"<<std::endl;
00432                       if(debug) std::cout<<"MB4 \t \t \t Strip Width "<<stripw<<"cm"<<std::endl;
00433 
00434                       if(debug) std::cout<<"MB4 \t \t \t X Predicted in MB4DTLocal= "<<X<<"cm"<<std::endl;
00435                       if(debug) std::cout<<"MB4 \t \t \t Y Predicted in MB4DTLocal= "<<Y<<"cm"<<std::endl;
00436                       if(debug) std::cout<<"MB4 \t \t \t Z Predicted in MB4DTLocal= "<<Z<<"cm"<<std::endl;
00437 
00438                       float extrapolatedDistance = sqrt((Y-Yo34)*(Y-Yo34)+Dy*Dy);
00439 
00440                       if(debug) std::cout<<"MB4 \t \t \t segmentPositionMB3inMB4Frame"<<segmentPositionMB3inMB4Frame<<std::endl;
00441                       if(debug) std::cout<<"MB4 \t \t \t segmentPositionMB4inMB4Frame"<<segmentPosition<<std::endl;
00442 
00443                       if(debug) std::cout<<"MB4 \t \t \t segmentDirMB3inMB4Frame"<<segDirMB3inMB4Frame<<std::endl;
00444                       if(debug) std::cout<<"MB4 \t \t \t segmentDirMB4inMB4Frame"<<segmentDirectionMB4<<std::endl;
00445                             
00446                       if(debug) std::cout<<"MB4 \t \t \t CenterRB4PositioninMB4Frame"<<CenterRollinMB4Frame<<std::endl;
00447                             
00448                       if(debug) std::cout<<"MB4 \t \t \t Is the extrapolation distance ="<<extrapolatedDistance<<"less than "<<MaxDrb4<<std::endl;
00449     
00450 
00451                       if(extrapolatedDistance<=MaxDrb4){ 
00452                         if(debug) std::cout<<"MB4 \t \t \t yes"<<std::endl;
00453 
00454                         GlobalPoint GlobalPointExtrapolated = DTSurface4.toGlobal(LocalPoint(X,Y,Z));
00455                                 
00456                         if(debug) std::cout<<"MB4 \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
00457                                 
00458                         LocalPoint PointExtrapolatedRPCFrame = RPCSurfaceRB4.toLocal(GlobalPointExtrapolated);
00459 
00460                         if(debug) std::cout<<"MB4 \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
00461                         if(debug) std::cout<<"MB4 \t \t \t Corner of the Roll = ("<<rsize*eyr<<","<<stripl*eyr<<")"<<std::endl;
00462                         if(debug) std::cout<<"MB4 \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.x())<<","
00463                                            <<fabs(PointExtrapolatedRPCFrame.y())<<","<<fabs(PointExtrapolatedRPCFrame.z())<<")"<<std::endl;
00464                           
00465                         if(debug) std::cout<<"MB4 \t \t \t Does the extrapolation go inside this roll?"<<std::endl;
00466                           
00467                         if(fabs(PointExtrapolatedRPCFrame.z()) < 5.  &&
00468                            fabs(PointExtrapolatedRPCFrame.x()) < rsize*eyr &&
00469                            fabs(PointExtrapolatedRPCFrame.y()) < stripl*eyr){
00470                           if(debug) std::cout<<"MB4 \t \t \t \t yes"<<std::endl;
00471                           if(debug) std::cout<<"MB4 \t \t \t \t Creating the RecHit"<<std::endl;
00472                           RPCRecHit RPCPointMB4(rpcId,0,PointExtrapolatedRPCFrame);
00473                           if(debug) std::cout<<"MB4 \t \t \t \t Clearing the RPCPointVector"<<std::endl;
00474                           RPCPointVector.clear();
00475                           if(debug) std::cout<<"MB4 \t \t \t \t Pushing Back"<<std::endl;
00476                           RPCPointVector.push_back(RPCPointMB4);
00477                           if(debug) std::cout<<"MB4 \t \t \t \t Putting for "<<rpcId<<std::endl;
00478                           if(debug) std::cout<<"MB4 \t \t \t \t Filling container with "<<nameRoll
00479                                              <<" Point.x="<<PointExtrapolatedRPCFrame.x()<<" Point.y="<<PointExtrapolatedRPCFrame.y()<<" size="<<RPCPointVector.size()<<std::endl;
00480                           if(debug) std::cout<<"MB4 \t \t \t \t Number of rolls already extrapolated in RB4 = "<<extrapolatedRolls.size()<<std::endl;
00481                           if(find (extrapolatedRolls.begin(),extrapolatedRolls.end(),rpcId.rawId()) == extrapolatedRolls.end()){
00482                             extrapolatedRolls.push_back(rpcId.rawId());
00483                             _ThePoints->put(rpcId,RPCPointVector.begin(),RPCPointVector.end());
00484                           }else{
00485                             if(debug) std::cout<<"MB4 \t \t \t \t roll already extrapolated "<<rpcId<<std::endl;
00486                           }
00487                           if(debug) std::cout<<"MB4 \t \t \t \t Extrapolations done after this point = "<<extrapolatedRolls.size()<<std::endl;
00488                           if(debug) for(uint32_t m=0;m<extrapolatedRolls.size();m++) std::cout<<"MB4 \t \t \t \t"<< extrapolatedRolls.at(m)<<std::endl;
00489                         }else{
00490                           if(debug) std::cout<<"MB4 \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00491                         }
00492                       }//Condition for the right match
00493                       else{
00494                         if(debug) std::cout<<"MB4 \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00495                       }
00496                     }//loop over all the rollsasociated
00497                   }else{
00498                     if(debug) std::cout<<"MB4 \t \t \t \t I found segments in MB4 and MB3 adjacent wheel and/or sector but not compatibles, Diferent Directions"<<std::endl;
00499                   }
00500                 }else{//if dtid3.station()==3&&dtid3.sector()==DTId.sector()&&dtid3.wheel()==DTId.wheel()&&segMB3->dim()==4
00501                   if(debug) std::cout<<"MB4 \t \t \t No the same station or same wheel or segment dim in mb3 not 4D"<<std::endl;
00502                 }
00503               }//loop over all the segments looking for one in MB3 
00504             }else{
00505               if(debug) std::cout<<"MB4 \t \t \t Is NOT a 2D Segment"<<std::endl;
00506             }
00507           }else{
00508             if(debug) std::cout<<"MB4 \t \t \t \t There is not just one segment or is not in station 4"<<std::endl;
00509           }//De aca para abajo esta en dtpart.inl
00510         }
00511       }else{
00512         if(debug) std::cout<<"MB4 \t This event doesn't have 4D Segment"<<std::endl;
00513       }
00514     }
00515   }
00516 
00517   
00518   
00519   
00520   /*
00521   clock_gettime(CLOCK_REALTIME, &stop_time);
00522   fs=start_time.tv_sec;
00523   fn=start_time.tv_nsec;
00524   ls=stop_time.tv_sec;
00525   ln=stop_time.tv_nsec;
00526   std::cout <<" =================|||| "<<ls-fs<<" sec "<<ln-fn<<" us"<<std::endl;
00527   */
00528 }
00529 
00530 
00531 
00532 DTSegtoRPC::~DTSegtoRPC(){
00533 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/RecoLocalMuon/RPCRecHit/src/DTSegtoRPC.cc
