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               bool compatiblesegments=false;
00309                 
00310               const BoundPlane& DTSurface4 = dtGeo->idToDet(DTId)->surface();
00311                 
00312               DTRecSegment4DCollection::const_iterator segMB3;  
00313                 
00314               if(debug) std::cout<<"MB4 \t \t \t \t Loop on segments in =sector && MB3 && adjacent sectors && y dim=4"<<std::endl;
00315               for(segMB3=all4DSegments->begin();segMB3!=all4DSegments->end();++segMB3){
00316                     
00317                 DTChamberId dtid3 = segMB3->chamberId();  
00318                 
00319                 if(debug) std::cout<<"MB4  \t \t \t \t Segment in Chamber ="<<dtid3<<std::endl;
00320 
00321                     
00322                 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.
00323                    && distwheel(dtid3.wheel(),DTId.wheel())<=1 //The we could have segments in neighbohr wheels in pp collisions 
00324                    && dtid3.station()==3
00325                    && DTSegmentCounter[dtid3] == 1
00326                    && segMB3->dimension()==4){
00327                   
00328                   if(debug) std::cout<<"MB4  \t \t \t \t distsector ="<<distsector(dtid3.sector(),DTId.sector())<<std::endl;
00329                   if(debug) std::cout<<"MB4  \t \t \t \t distwheel ="<<distwheel(dtid3.wheel(),DTId.wheel())<<std::endl;
00330 
00331                   const GeomDet* gdet3=dtGeo->idToDet(segMB3->geographicalId());
00332                   const BoundPlane & DTSurface3 = gdet3->surface();
00333 
00334                   LocalVector segmentDirectionMB3 =  segMB3->localDirection();
00335                   GlobalPoint segmentPositionMB3inGlobal = DTSurface3.toGlobal(segMB3->localPosition());
00336                   GlobalPoint segmentPositionMB4inGlobal = DTSurface4.toGlobal(segmentPosition);
00337                 
00338                   LocalVector segDirMB4inMB3Frame=DTSurface3.toLocal(DTSurface4.toGlobal(segmentDirectionMB4));
00339                   LocalVector segDirMB3inMB4Frame=DTSurface4.toLocal(DTSurface3.toGlobal(segmentDirectionMB3));
00340                 
00341                   GlobalVector segDirMB4inGlobalFrame=DTSurface4.toGlobal(segmentDirectionMB4);
00342                   GlobalVector segDirMB3inGlobalFrame=DTSurface3.toGlobal(segmentDirectionMB3);
00343                 
00344                   float dx=segDirMB4inGlobalFrame.x();
00345                   float dy=segDirMB4inGlobalFrame.y();
00346                                 
00347                   float dx3=segDirMB3inGlobalFrame.x();
00348                   float dy3=segDirMB3inGlobalFrame.y();
00349                                 
00350                   double cosAng=fabs(dx*dx3+dy*dy3/sqrt((dx3*dx3+dy3*dy3)*(dx*dx+dy*dy)));
00351 
00352                   if(debug) std::cout<<"MB4 \t \t \t \t cosAng"<<cosAng<<"Beetween "<<dtid3<<" and "<<DTId<<std::endl;
00353                 
00354                   if(debug){
00355                     std::cout<<"MB4 \t \t \t \t dx="<<dx<<" dy="<<dy<<std::endl;
00356                     std::cout<<"MB4 \t \t \t \t dx3="<<dx3<<" dy3="<<dy<<std::endl;
00357                     std::cout<<"MB4 \t \t \t \t cosAng="<<cosAng<<std::endl;
00358                   }
00359 
00360                   float DistanceBetweenSegments = ((segmentPositionMB3inGlobal) - (segmentPositionMB4inGlobal)).mag();
00361                   
00362                   if(cosAng>MinCosAng && DistanceBetweenSegments < MaxDistanceBetweenSegments){
00363                     
00364                     if(debug) std::cout<<"MB4 \t \t \t \t Distance between segments="<<DistanceBetweenSegments<<std::endl;
00365                   
00366                     compatiblesegments=true;
00367 
00368                     if(debug) std::cout<<"MB4 \t \t We found compatible Segments (similar direction and close enough) in "<<dtid3<<" and "<<DTId<<std::endl;
00369 
00370                     if(dtSector==13){
00371                       dtSector=4;
00372                     }
00373                     if(dtSector==14){
00374                       dtSector=10;
00375                     }
00376                    
00377                     if(debug)  std::cout<<"Calling to Object Map class"<<std::endl;
00378                     ObjectMap* TheObject = ObjectMap::GetInstance(iSetup);
00379                     if(debug) std::cout<<"Creating the DTIndex"<<std::endl;
00380                     DTStationIndex theindex(0,dtWheel,dtSector,dtStation);
00381                     if(debug) std::cout<<"Getting the Rolls for the given index"<<std::endl;
00382 
00383                     std::set<RPCDetId> rollsForThisDT = TheObject->GetInstance(iSetup)->GetRolls(theindex);
00384 
00385                     if(debug) std::cout<<"MB4 \t \t Number of rolls for this DT = "<<rollsForThisDT.size()<<std::endl;
00386                     
00387                     assert(rollsForThisDT.size()>=1);
00388                         
00389                     if(debug) std::cout<<"MB4  \t \t Loop over all the rolls asociated to this DT"<<std::endl;
00390                     for (std::set<RPCDetId>::iterator iteraRoll=rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++){
00391                       const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll); //roll asociado a MB4
00392                       RPCDetId rpcId = rollasociated->id();
00393                       const BoundPlane & RPCSurfaceRB4 = rollasociated->surface(); //surface MB4
00394 
00395                       RPCGeomServ rpcsrv(rpcId);
00396                       std::string nameRoll = rpcsrv.name();
00397 
00398                       if(debug) std::cout<<"MB4  \t \t \t RollName: "<<nameRoll<<std::endl;
00399                       if(debug) std::cout<<"MB4  \t \t \t Doing the extrapolation to this roll"<<std::endl;
00400                         
00401                       GlobalPoint CenterPointRollGlobal=RPCSurfaceRB4.toGlobal(LocalPoint(0,0,0));
00402                       LocalPoint CenterRollinMB4Frame = DTSurface4.toLocal(CenterPointRollGlobal); //In MB4
00403                       LocalPoint segmentPositionMB3inMB4Frame = DTSurface4.toLocal(segmentPositionMB3inGlobal); //In MB4
00404                       LocalPoint segmentPositionMB3inRB4Frame = RPCSurfaceRB4.toLocal(segmentPositionMB3inGlobal); //In MB4
00405                       LocalVector segmentDirectionMB3inMB4Frame = DTSurface4.toLocal(segDirMB3inGlobalFrame); //In MB4
00406                         
00407                       //The exptrapolation is done in MB4 frame. for local x and z is done from MB4,
00408                       float Dxz=CenterRollinMB4Frame.z();
00409                       float Xo4=segmentPositionMB4.x();
00410                       float dxl=segmentDirectionMB4.x(); //dx local for MB4 segment in MB4 Frame
00411                       float dzl=segmentDirectionMB4.z(); //dx local for MB4 segment in MB4 Frame
00412                         
00413                       float X=Xo4+dxl*Dxz/dzl; //In MB4 frame
00414                       float Z=Dxz;//In MB4 frame
00415                         
00416                       //for local y is done from MB3 
00417                       float Yo34 = segmentPositionMB3inMB4Frame.y();
00418                       float dy34 = segmentDirectionMB3inMB4Frame.y();
00419                       float dz34 = segmentDirectionMB3inMB4Frame.z();
00420                       float Dy=Dxz-(segmentPositionMB3inMB4Frame.z()); //Distance beetween the segment in MB3 and the RB4 surface
00421 
00422                       if(debug) std::cout<<"MB4 \t \t \t The distance to extrapolate in Y from MB3 is "<<Dy<<"cm"<<std::endl;
00423                         
00424                       float Y=Yo34+dy34*Dy/dz34;//In MB4 Frame
00425                           
00426                       const RectangularStripTopology* top_
00427                         =dynamic_cast<const RectangularStripTopology*>(&(rollasociated->topology())); //Topology roll asociated MB4
00428                       LocalPoint xmin = top_->localPosition(0.);
00429                       LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
00430                       float rsize = fabs( xmax.x()-xmin.x() );
00431                       float stripl = top_->stripLength();
00432                       float stripw = top_->pitch();
00433 
00434                             
00435                       if(debug) std::cout<<"MB4 \t \t \t Strip Lenght "<<stripl<<"cm"<<std::endl;
00436                       if(debug) std::cout<<"MB4 \t \t \t Strip Width "<<stripw<<"cm"<<std::endl;
00437 
00438                       if(debug) std::cout<<"MB4 \t \t \t X Predicted in MB4DTLocal= "<<X<<"cm"<<std::endl;
00439                       if(debug) std::cout<<"MB4 \t \t \t Y Predicted in MB4DTLocal= "<<Y<<"cm"<<std::endl;
00440                       if(debug) std::cout<<"MB4 \t \t \t Z Predicted in MB4DTLocal= "<<Z<<"cm"<<std::endl;
00441 
00442                       float extrapolatedDistance = sqrt((Y-Yo34)*(Y-Yo34)+Dy*Dy);
00443 
00444                       if(debug) std::cout<<"MB4 \t \t \t segmentPositionMB3inMB4Frame"<<segmentPositionMB3inMB4Frame<<std::endl;
00445                       if(debug) std::cout<<"MB4 \t \t \t segmentPositionMB4inMB4Frame"<<segmentPosition<<std::endl;
00446 
00447                       if(debug) std::cout<<"MB4 \t \t \t segmentDirMB3inMB4Frame"<<segDirMB3inMB4Frame<<std::endl;
00448                       if(debug) std::cout<<"MB4 \t \t \t segmentDirMB4inMB4Frame"<<segmentDirectionMB4<<std::endl;
00449                             
00450                       if(debug) std::cout<<"MB4 \t \t \t CenterRB4PositioninMB4Frame"<<CenterRollinMB4Frame<<std::endl;
00451                             
00452                       if(debug) std::cout<<"MB4 \t \t \t Is the extrapolation distance ="<<extrapolatedDistance<<"less than "<<MaxDrb4<<std::endl;
00453     
00454 
00455                       if(extrapolatedDistance<=MaxDrb4){ 
00456                         if(debug) std::cout<<"MB4 \t \t \t yes"<<std::endl;
00457 
00458                         GlobalPoint GlobalPointExtrapolated = DTSurface4.toGlobal(LocalPoint(X,Y,Z));
00459                                 
00460                         if(debug) std::cout<<"MB4 \t \t \t Point ExtraPolated in Global"<<GlobalPointExtrapolated<< std::endl;
00461                                 
00462                         LocalPoint PointExtrapolatedRPCFrame = RPCSurfaceRB4.toLocal(GlobalPointExtrapolated);
00463 
00464                         if(debug) std::cout<<"MB4 \t \t \t Point Extrapolated in RPCLocal"<<PointExtrapolatedRPCFrame<< std::endl;
00465                         if(debug) std::cout<<"MB4 \t \t \t Corner of the Roll = ("<<rsize*eyr<<","<<stripl*eyr<<")"<<std::endl;
00466                         if(debug) std::cout<<"MB4 \t \t \t Info About the Point Extrapolated in X Abs ("<<fabs(PointExtrapolatedRPCFrame.x())<<","
00467                                            <<fabs(PointExtrapolatedRPCFrame.y())<<","<<fabs(PointExtrapolatedRPCFrame.z())<<")"<<std::endl;
00468                           
00469                         if(debug) std::cout<<"MB4 \t \t \t Does the extrapolation go inside this roll?"<<std::endl;
00470                           
00471                         if(fabs(PointExtrapolatedRPCFrame.z()) < 5.  &&
00472                            fabs(PointExtrapolatedRPCFrame.x()) < rsize*eyr &&
00473                            fabs(PointExtrapolatedRPCFrame.y()) < stripl*eyr){
00474                           if(debug) std::cout<<"MB4 \t \t \t \t yes"<<std::endl;
00475                           if(debug) std::cout<<"MB4 \t \t \t \t Creating the RecHit"<<std::endl;
00476                           RPCRecHit RPCPointMB4(rpcId,0,PointExtrapolatedRPCFrame);
00477                           if(debug) std::cout<<"MB4 \t \t \t \t Clearing the RPCPointVector"<<std::endl;
00478                           RPCPointVector.clear();
00479                           if(debug) std::cout<<"MB4 \t \t \t \t Pushing Back"<<std::endl;
00480                           RPCPointVector.push_back(RPCPointMB4);
00481                           if(debug) std::cout<<"MB4 \t \t \t \t Putting for "<<rpcId<<std::endl;
00482                           if(debug) std::cout<<"MB4 \t \t \t \t Filling container with "<<nameRoll
00483                                              <<" Point.x="<<PointExtrapolatedRPCFrame.x()<<" Point.y="<<PointExtrapolatedRPCFrame.y()<<" size="<<RPCPointVector.size()<<std::endl;
00484                           if(debug) std::cout<<"MB4 \t \t \t \t Number of rolls already extrapolated in RB4 = "<<extrapolatedRolls.size()<<std::endl;
00485                           if(find (extrapolatedRolls.begin(),extrapolatedRolls.end(),rpcId.rawId()) == extrapolatedRolls.end()){
00486                             extrapolatedRolls.push_back(rpcId.rawId());
00487                             _ThePoints->put(rpcId,RPCPointVector.begin(),RPCPointVector.end());
00488                           }else{
00489                             if(debug) std::cout<<"MB4 \t \t \t \t roll already extrapolated "<<rpcId<<std::endl;
00490                           }
00491                           if(debug) std::cout<<"MB4 \t \t \t \t Extrapolations done after this point = "<<extrapolatedRolls.size()<<std::endl;
00492                           if(debug) for(uint32_t m=0;m<extrapolatedRolls.size();m++) std::cout<<"MB4 \t \t \t \t"<< extrapolatedRolls.at(m)<<std::endl;
00493                         }else{
00494                           if(debug) std::cout<<"MB4 \t \t \t \t No the prediction is outside of this roll"<<std::endl;
00495                         }
00496                       }//Condition for the right match
00497                       else{
00498                         if(debug) std::cout<<"MB4 \t \t \t No, Exrtrapolation too long!, canceled"<<std::endl;
00499                       }
00500                     }//loop over all the rollsasociated
00501                   }else{
00502                     compatiblesegments=false;
00503                     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;
00504                   }
00505                 }else{//if dtid3.station()==3&&dtid3.sector()==DTId.sector()&&dtid3.wheel()==DTId.wheel()&&segMB3->dim()==4
00506                   if(debug) std::cout<<"MB4 \t \t \t No the same station or same wheel or segment dim in mb3 not 4D"<<std::endl;
00507                 }
00508               }//loop over all the segments looking for one in MB3 
00509             }else{
00510               if(debug) std::cout<<"MB4 \t \t \t Is NOT a 2D Segment"<<std::endl;
00511             }
00512           }else{
00513             if(debug) std::cout<<"MB4 \t \t \t \t There is not just one segment or is not in station 4"<<std::endl;
00514           }//De aca para abajo esta en dtpart.inl
00515         }
00516       }else{
00517         if(debug) std::cout<<"MB4 \t This event doesn't have 4D Segment"<<std::endl;
00518       }
00519     }
00520   }
00521 
00522   
00523   
00524   
00525   /*
00526   clock_gettime(CLOCK_REALTIME, &stop_time);
00527   fs=start_time.tv_sec;
00528   fn=start_time.tv_nsec;
00529   ls=stop_time.tv_sec;
00530   ln=stop_time.tv_nsec;
00531   std::cout <<" =================|||| "<<ls-fs<<" sec "<<ln-fn<<" us"<<std::endl;
00532   */
00533 }
00534 
00535 
00536 
00537 DTSegtoRPC::~DTSegtoRPC(){
00538 }
/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/RecoLocalMuon/RPCRecHit/src/DTSegtoRPC.cc
