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 "RecoLocalMuon/RPCRecHit/interface/CSCSegtoRPC.h"
00014 #include <DataFormats/DetId/interface/DetId.h>
00015 #include <RecoLocalMuon/RPCRecHit/interface/TracktoRPC.h>
00016 #include <ctime>
00017 #include <TMath.h>
00018 
00019 ObjectMap2* ObjectMap2::mapInstance = NULL;
00020 
00021 ObjectMap2* ObjectMap2::GetInstance(const edm::EventSetup& iSetup){
00022   if (mapInstance == NULL){
00023     mapInstance = new ObjectMap2(iSetup);
00024   }
00025   return mapInstance;
00026 }
00027 
00028 ObjectMap2::ObjectMap2(const edm::EventSetup& iSetup){
00029   edm::ESHandle<RPCGeometry> rpcGeo;
00030   edm::ESHandle<DTGeometry> dtGeo;
00031   
00032   iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00033   iSetup.get<MuonGeometryRecord>().get(dtGeo);
00034   
00035   for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00036     if(dynamic_cast< RPCChamber* >( *it ) != 0 ){
00037       RPCChamber* ch = dynamic_cast< RPCChamber* >( *it ); 
00038       std::vector< const RPCRoll*> roles = (ch->rolls());
00039       for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){
00040         RPCDetId rpcId = (*r)->id();
00041         int region=rpcId.region();
00042         if(region==0){
00043           int wheel=rpcId.ring();
00044           int sector=rpcId.sector();
00045           int station=rpcId.station();
00046           DTStationIndex2 ind(region,wheel,sector,station);
00047           std::set<RPCDetId> myrolls;
00048           if (rollstoreDT.find(ind)!=rollstoreDT.end()) myrolls=rollstoreDT[ind];
00049           myrolls.insert(rpcId);
00050           rollstoreDT[ind]=myrolls;
00051         }
00052       }
00053     }
00054   }
00055 }
00056 
00057 int distsector2(int sector1,int sector2){
00058   if(sector1==13) sector1=4;
00059   if(sector1==14) sector1=10;
00060   
00061   if(sector2==13) sector2=4;
00062   if(sector2==14) sector2=10;
00063   
00064   int distance = std::abs(sector1 - sector2);
00065   if(distance>6) distance = 12-distance;
00066   return distance;
00067 }
00068 
00069 int distwheel2(int wheel1,int wheel2){
00070   int distance = std::abs(wheel1 - wheel2);
00071   return distance;
00072 }
00073 ObjectMap2CSC* ObjectMap2CSC::mapInstance = NULL;
00074 
00075 ObjectMap2CSC* ObjectMap2CSC::GetInstance(const edm::EventSetup& iSetup){
00076   if (mapInstance == NULL){
00077     mapInstance = new ObjectMap2CSC(iSetup);
00078   }
00079   return mapInstance;
00080 }
00081 
00082 ObjectMap2CSC::ObjectMap2CSC(const edm::EventSetup& iSetup){
00083   edm::ESHandle<RPCGeometry> rpcGeo;
00084   edm::ESHandle<CSCGeometry> cscGeo;
00085   
00086   iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00087   iSetup.get<MuonGeometryRecord>().get(cscGeo);
00088   
00089   for (TrackingGeometry::DetContainer::const_iterator it=rpcGeo->dets().begin();it<rpcGeo->dets().end();it++){
00090     if(dynamic_cast< RPCChamber* >( *it ) != 0 ){
00091       RPCChamber* ch = dynamic_cast< RPCChamber* >( *it ); 
00092       std::vector< const RPCRoll*> roles = (ch->rolls());
00093       for(std::vector<const RPCRoll*>::const_iterator r = roles.begin();r != roles.end(); ++r){
00094         RPCDetId rpcId = (*r)->id();
00095         int region=rpcId.region();
00096         if(region!=0){
00097           int station=rpcId.station();
00098           int ring=rpcId.ring();
00099           int cscring=ring;
00100           int cscstation=station;
00101           RPCGeomServ rpcsrv(rpcId);
00102           int rpcsegment = rpcsrv.segment();
00103           int cscchamber = rpcsegment; //FIX THIS ACCORDING TO RPCGeomServ::segment()Definition
00104           if((station==2||station==3)&&ring==3){//Adding Ring 3 of RPC to the CSC Ring 2
00105             cscring = 2;
00106           }
00107           CSCStationIndex2 ind(region,cscstation,cscring,cscchamber);
00108           std::set<RPCDetId> myrolls;
00109           if (rollstoreCSC.find(ind)!=rollstoreCSC.end()) myrolls=rollstoreCSC[ind];
00110           myrolls.insert(rpcId);
00111           rollstoreCSC[ind]=myrolls;
00112         }
00113       }
00114     }
00115   }
00116 }
00117 
00118 bool TracktoRPC::ValidRPCSurface(RPCDetId rpcid, LocalPoint LocalP, const edm::EventSetup& iSetup)
00119 {
00120   edm::ESHandle<RPCGeometry> rpcGeo;
00121   iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00122 
00123   const GeomDet *whichdet3 = rpcGeo->idToDet(rpcid.rawId());
00124   const RPCRoll *aroll = dynamic_cast<const RPCRoll *>(whichdet3);
00125   float locx=LocalP.x(), locy=LocalP.y();//, locz=LocalP.z();
00126   if(aroll->isBarrel())
00127   {
00128      const Bounds &rollbound = rpcGeo->idToDet((rpcid))->surface().bounds();
00129      float boundlength = rollbound.length();
00130      float boundwidth = rollbound.width();
00131 
00132      if(fabs(locx) < boundwidth/2 && fabs(locy) < boundlength/2 && locy > -boundlength/2) return true;
00133      else return false;
00134 
00135    }
00136    else if(aroll->isForward())
00137    {
00138      const Bounds &rollbound = rpcGeo->idToDet((rpcid))->surface().bounds();
00139      float boundlength = rollbound.length();
00140      float boundwidth = rollbound.width();
00141 
00142      float nminx = TMath::Pi()*(18*boundwidth/ TMath::Pi() - boundlength)/18;
00143      float ylimit = ((boundlength)/(boundwidth/2 - nminx/2))*fabs(locx) + boundlength/2 - ((boundlength)/(boundwidth/2 - nminx/2))*(boundwidth/2);
00144      if(ylimit < -boundlength/2 ) ylimit = -boundlength/2;
00145 
00146      if(fabs(locx) < boundwidth/2 && fabs(locy) < boundlength/2 && locy > ylimit) return true;
00147      else return false;
00148    } else return false;
00149 }
00150 
00151 TracktoRPC::TracktoRPC(edm::Handle<reco::TrackCollection> alltracks, const edm::EventSetup& iSetup,const edm::Event& iEvent,bool debug,const edm::ParameterSet& iConfig,edm::InputTag& tracklabel){ 
00152 
00153  _ThePoints = new RPCRecHitCollection();
00154 // if(alltracks->empty()) return;
00155 
00156  if(tracklabel.label().find("cosmic")==0) theTrackTransformer = new TrackTransformerForCosmicMuons(iConfig);
00157  else if(tracklabel.label().find("globalCosmic")==0) theTrackTransformer = new TrackTransformerForCosmicMuons(iConfig);
00158  else theTrackTransformer = new TrackTransformer(iConfig);
00159  theTrackTransformer->setServices(iSetup);  
00160 
00161  edm::ESHandle<RPCGeometry> rpcGeo;
00162  edm::ESHandle<DTGeometry> dtGeo;
00163  edm::ESHandle<CSCGeometry> cscGeo;
00164  
00165  iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny",thePropagator); 
00166  iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00167  iSetup.get<MuonGeometryRecord>().get(dtGeo);
00168  iSetup.get<MuonGeometryRecord>().get(cscGeo);
00169 
00170 std::vector<uint32_t> rpcput;
00171 double MaxD=999.;
00172 
00173 for (TrackCollection::const_iterator track = alltracks->begin(); track !=alltracks->end(); track++)
00174 {
00175  Trajectories trajectories = theTrackTransformer->transform(*track);
00176  if(debug) std::cout << "Building Trajectory from Track. " << std::endl;
00177 
00178  std::vector<uint32_t> rpcrolls;
00179  std::vector<uint32_t> rpcrolls2; 
00180  std::map<uint32_t, int> rpcNdtcsc;
00181  std::map<uint32_t, int> rpcrollCounter;
00182 
00183  float tInX = track->innerPosition().X(), tInY = track->innerPosition().Y(), tInZ = track->innerPosition().Z();
00184  float tOuX = track->outerPosition().X(), tOuY = track->outerPosition().Y(), tOuZ = track->outerPosition().Z();
00185  if(tInX > tOuX) { float temp=tOuX; tOuX=tInX; tInX=temp; }
00186  if(tInY > tOuY) { float temp=tOuY; tOuY=tInY; tInY=temp; }
00187  if(tInZ > tOuZ) { float temp=tOuZ; tOuZ=tInZ; tInZ=temp; }
00188 
00189  if(debug) std::cout << "in (x,y,z): ("<< tInX <<", "<< tInY <<", "<< tInZ << ")" << std::endl;
00190  if(debug) std::cout << "out (x,y,z): ("<< tOuX <<", "<< tOuY <<", "<< tOuZ << ")" << std::endl;
00191 
00192 if(debug) std::cout << "1. Search expeted RPC roll detid !!" << std::endl;
00193 for(trackingRecHit_iterator hit=track->recHitsBegin(); hit != track->recHitsEnd(); hit++)
00194  {
00195     if((*hit)->isValid())
00196     {
00197       DetId id = (*hit)->geographicalId();
00198 
00199        if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::DT)
00200        {
00201           const GeomDet *geomDet =  dtGeo->idToDet((*hit)->geographicalId());
00202           const DTLayer *dtlayer = dynamic_cast<const DTLayer *>(geomDet);
00203          if(dtlayer) for(Trajectories::const_iterator trajectory = trajectories.begin(); trajectory != trajectories.end(); ++trajectory)
00204           {
00205              const BoundPlane & DTSurface = dtlayer->surface();
00206              const GlobalPoint dcPoint = DTSurface.toGlobal(LocalPoint(0.,0.,0.));
00207 
00208              TrajectoryMeasurement tMt = trajectory->closestMeasurement(dcPoint);
00209              TrajectoryStateOnSurface upd2 = (tMt).updatedState();
00210              if(upd2.isValid())
00211              {
00212                 LocalPoint trajLP = upd2.localPosition();
00213                 LocalPoint trackLP = (*hit)->localPosition();
00214                 float dx = trajLP.x()-trackLP.x(), dy=trajLP.y()-trackLP.y();//, dz=trajLP.z()-trackLP.z();
00215                 if( dx>10. && dy>10.) continue;
00216 
00217                 DTChamberId dtid(geomDet->geographicalId().rawId());
00218                 int dtW=dtid.wheel(), dtS=dtid.sector(), dtT=dtid.station();
00219                 if(dtS==13) dtS=4; if(dtS==14) dtS=10;
00220                 ObjectMap2* TheObject = ObjectMap2::GetInstance(iSetup);
00221                 DTStationIndex2 theindex(0,dtW,dtS,dtT);
00222                 std::set<RPCDetId> rollsForThisDT = TheObject->GetInstance(iSetup)->GetRolls(theindex);
00223                 for(std::set<RPCDetId>::iterator iteraRoll = rollsForThisDT.begin();iteraRoll != rollsForThisDT.end(); iteraRoll++)
00224                 {                                 
00225                     const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00226                     RPCDetId rpcId = rollasociated->id();
00227 
00228                     TrajectoryStateOnSurface ptss =  thePropagator->propagate(upd2, rpcGeo->idToDet(rollasociated->id())->surface());
00229                     if(ptss.isValid()) if(ValidRPCSurface(rollasociated->id().rawId(), ptss.localPosition(), iSetup))
00230                     {
00231                       rpcrollCounter[rollasociated->id().rawId()]++;
00232                       bool check = true;
00233                       std::vector<uint32_t>::iterator rpcroll;
00234                       for( rpcroll=rpcrolls.begin() ; rpcroll < rpcrolls.end(); rpcroll++ )
00235                       if(rollasociated->id().rawId()== *rpcroll) check=false; 
00236                       if(check==true)
00237                       {
00238                         rpcrolls.push_back(rollasociated->id().rawId());
00239                         RPCGeomServ servId(rollasociated->id().rawId());
00240                         if(debug) std::cout << "1\t Barrel RPC roll" << rollasociated->id().rawId() << " "<< servId.name().c_str() <<std::endl;
00241                       }
00242                     }
00243                 }                                 
00244              }
00245           }
00246        }
00247        else if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::CSC) 
00248        {
00249           const GeomDet *geomDet =  cscGeo->idToDet((*hit)->geographicalId());
00250           const CSCLayer *csclayer = dynamic_cast<const CSCLayer *>(geomDet);
00251 
00252           CSCDetId cscid(geomDet->geographicalId().rawId());
00253           if(csclayer) for(Trajectories::const_iterator trajectory = trajectories.begin(); trajectory != trajectories.end(); ++trajectory)
00254           {
00255              const BoundPlane & CSCSurface = csclayer->surface();
00256              const GlobalPoint dcPoint = CSCSurface.toGlobal(LocalPoint(0.,0.,0.));
00257 
00258              TrajectoryMeasurement tMt = trajectory->closestMeasurement(dcPoint);
00259              TrajectoryStateOnSurface upd2 = (tMt).updatedState();
00260 
00261              if(upd2.isValid() && cscid.station()!=4 && cscid.ring()!=1 )
00262              {
00263                 LocalPoint trajLP = upd2.localPosition();
00264                 LocalPoint trackLP = (*hit)->localPosition();
00265                 float dx = trajLP.x()-trackLP.x(), dy=trajLP.y()-trackLP.y();//, dz=trajLP.z()-trackLP.z();
00266                 if( dx>10. && dy>10.) continue;
00267 
00268                 ObjectMap2CSC* TheObjectCSC = ObjectMap2CSC::GetInstance(iSetup);
00269                 int En = cscid.endcap(), St = cscid.station(), Ri = cscid.ring();
00270                 int rpcSegment = cscid.chamber();
00271                 if(En==2) En= -1; if(Ri==4) Ri =1; 
00272 
00273                 CSCStationIndex2 theindex(En,St,Ri,rpcSegment);
00274                 std::set<RPCDetId> rollsForThisCSC = TheObjectCSC->GetInstance(iSetup)->GetRolls(theindex);
00275                 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisCSC.begin();iteraRoll != rollsForThisCSC.end(); iteraRoll++)
00276                 {
00277                     const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
00278                     RPCDetId rpcId = rollasociated->id();
00279 
00280                     TrajectoryStateOnSurface ptss =  thePropagator->propagate(upd2, rpcGeo->idToDet(rollasociated->id())->surface());
00281                     if(ptss.isValid()) if(ValidRPCSurface(rollasociated->id().rawId(), ptss.localPosition(), iSetup))
00282                     {
00283                       rpcrollCounter[rollasociated->id().rawId()]++;
00284                       bool check = true;
00285                       std::vector<uint32_t>::iterator rpcroll;
00286                       for( rpcroll=rpcrolls.begin() ; rpcroll < rpcrolls.end(); rpcroll++ )
00287                       if(rollasociated->id().rawId()==*rpcroll) check=false; 
00288                       if(check==true)
00289                       {
00290                         rpcrolls.push_back(rollasociated->id().rawId());
00291                         RPCGeomServ servId(rollasociated->id().rawId());
00292                         if(debug) std::cout << "1\t Forward RPC roll" << rollasociated->id().rawId() << " "<< servId.name().c_str() <<std::endl;
00293                       }
00294                     }
00295                 }
00296              }
00297           }
00298        } else { if(debug) std::cout << "1\t The hit is not DT/CSC's.   " << std::endl;} 
00299     }
00300  }
00301  if(debug) std::cout << "First step OK!!\n2. Search nearest DT/CSC sufrace!!" << std::endl;
00302  std::vector<uint32_t>::iterator rpcroll;
00303  for( rpcroll=rpcrolls.begin() ; rpcroll < rpcrolls.end(); rpcroll++ )
00304  {
00305     RPCDetId rpcid(*rpcroll);
00306     const GlobalPoint &rGP = rpcGeo->idToDet(*rpcroll)->surface().toGlobal(LocalPoint(0,0,0));
00307     RPCGeomServ servId(rpcid);
00308     int rEn=rpcid.region(), rSe=rpcid.sector(), rWr=rpcid.ring(), rSt=rpcid.station(), rCh=servId.segment();
00309 
00310     if(rpcrollCounter[*rpcroll]<3) continue ;
00311     
00312     uint32_t dtcscid=0; double distance= MaxD;
00313 
00314    // if(rSt ==2 && rEn==0) MaxD=100;
00315    // else if(rSt ==3 && rEn==0) MaxD=100;
00316    // else if(rSt ==4 && rEn==0) MaxD =150;
00317     for(trackingRecHit_iterator hit=track->recHitsBegin(); hit != track->recHitsEnd(); hit++)
00318     {
00319         if((*hit)->isValid())
00320         {
00321             DetId id = (*hit)->geographicalId(); 
00322             if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::DT)
00323             {
00324                const GeomDet *geomDet =  dtGeo->idToDet((*hit)->geographicalId());
00325                //const DTLayer *dtlayer = dynamic_cast<const DTLayer *>(geomDet);
00326                const GlobalPoint &dtGP = dtGeo->idToDet((*hit)->geographicalId())->surface().toGlobal(LocalPoint(0,0,0));
00327                double dx = rGP.x()-dtGP.x(), dy = rGP.y()-dtGP.y(), dz = rGP.z()-dtGP.z();
00328                double distanceN = sqrt(dx*dx+dy*dy+dz*dz);
00329 
00330                DTChamberId dtid(geomDet->geographicalId().rawId());
00331                int Se = dtid.sector(), Wh = dtid.wheel(), St = dtid.station();
00332                if(Se == 13) Se=4; if(Se ==14) Se=10;  
00333 
00334                if( rEn==0&& (rSe-Se)==0 && (rWr-Wh) ==0 && (rSt-St)==0 && distanceN < distance)
00335                {
00336                    dtcscid=geomDet->geographicalId().rawId();
00337                    distance = distanceN;
00338                    if(debug) std::cout << "2\t DT "<< dtcscid << " Wheel : " << Wh << " station : " << St << " sector : " << Se << std::endl; 
00339                }
00340             }
00341             else if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::CSC)
00342             {
00343                const GeomDet *geomDet =  cscGeo->idToDet((*hit)->geographicalId());
00344                //const CSCLayer *csclayer = dynamic_cast<const CSCLayer *>(geomDet);
00345                const GlobalPoint &cscGP = cscGeo->idToDet((*hit)->geographicalId())->surface().toGlobal(LocalPoint(0,0,0));
00346                double dx = rGP.x()-cscGP.x(), dy = rGP.y()-cscGP.y(), dz = rGP.z()-cscGP.z();
00347                double distanceN = sqrt(dx*dx+dy*dy+dz*dz);
00348 
00349                CSCDetId cscid(geomDet->geographicalId().rawId());
00350                int En =cscid.endcap(), Ri=cscid.ring(), St=cscid.station(), Ch=cscid.chamber();
00351                if(En==2) En=-1; if(Ri==4) Ri=1;
00352 
00353                if((rEn-En)==0 && (rSt-St)==0 && (Ch-rCh) ==0 && rWr!=1 && rSt!=4 && distanceN < distance)
00354                {
00355                    dtcscid=geomDet->geographicalId().rawId();
00356                    distance = distanceN;
00357                    if(debug) std::cout << "2\t CSC " <<dtcscid <<" region : " << En << " station : " << St << " Ring : " << Ri << " chamber : " << Ch <<std::endl;
00358                }
00359             } 
00360          }
00361     }
00362     if(dtcscid != 0 && distance < MaxD)
00363     {
00364        rpcrolls2.push_back(*rpcroll);
00365        rpcNdtcsc[*rpcroll] = dtcscid;
00366     }
00367  } 
00368  if(debug) std::cout << "Second step OK!! \n3. Propagate to RPC from DT/CSC!!" << std::endl;  
00369  //std::map<uint32_t, int> rpcput;
00370  std::vector<uint32_t>::iterator rpcroll2;
00371  for( rpcroll2=rpcrolls2.begin() ; rpcroll2 < rpcrolls2.end(); rpcroll2++ )
00372  {
00373     bool check = true;
00374     std::vector<uint32_t>::iterator rpcput_;
00375     for( rpcput_=rpcput.begin() ; rpcput_ < rpcput.end(); rpcput_++ )
00376     if(*rpcroll2==*rpcput_) check = false;
00377 
00378     if(check == true)
00379     {
00380         uint32_t dtcscid = rpcNdtcsc[*rpcroll2];
00381         DetId id(dtcscid);
00382         if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::DT)
00383         {
00384            const GeomDet *geomDet =  dtGeo->idToDet(dtcscid);
00385            const DTLayer *dtlayer = dynamic_cast<const DTLayer *>(geomDet);
00386         
00387            if(dtlayer) for(Trajectories::const_iterator trajectory = trajectories.begin(); trajectory != trajectories.end(); ++trajectory)
00388            {
00389               const BoundPlane & DTSurface = dtlayer->surface();
00390               const GlobalPoint dcPoint = DTSurface.toGlobal(LocalPoint(0.,0.,0.));
00391 
00392               TrajectoryMeasurement tMt = trajectory->closestMeasurement(dcPoint);
00393               TrajectoryStateOnSurface upd2 = (tMt).updatedState();
00394               if(upd2.isValid())
00395               {
00396                  TrajectoryStateOnSurface ptss =  thePropagator->propagate(upd2, rpcGeo->idToDet(*rpcroll2)->surface());
00397                  if(ptss.isValid())  if(ValidRPCSurface(*rpcroll2, ptss.localPosition(), iSetup))
00398                  {
00399                    float rpcGPX = ptss.globalPosition().x();
00400                    float rpcGPY = ptss.globalPosition().y();
00401                    float rpcGPZ = ptss.globalPosition().z();
00402 
00403                    if(tInX > rpcGPX || tOuX < rpcGPX ) continue;
00404                    if(tInY > rpcGPY || tOuY < rpcGPY ) continue;
00405                    if(tInZ > rpcGPZ || tOuZ < rpcGPZ ) continue;
00406                  
00407                    const GeomDet *geomDet2 = rpcGeo->idToDet(*rpcroll2);
00408                    const RPCRoll *aroll = dynamic_cast<const RPCRoll *>(geomDet2);
00409                    const RectangularStripTopology* top_= dynamic_cast<const RectangularStripTopology*> (&(aroll->topology()));
00410                    LocalPoint xmin = top_->localPosition(0.);
00411                    LocalPoint xmax = top_->localPosition((float)aroll->nstrips());
00412                    float rsize = fabs( xmax.x()-xmin.x() );
00413                    float stripl = top_->stripLength();
00414                    //float stripw = top_->pitch();
00415                    float eyr=1;
00416 
00417                    float locx = ptss.localPosition().x(), locy = ptss.localPosition().y(), locz= ptss.localPosition().z();
00418                    if( locx < rsize*eyr && locy < stripl*eyr && locz < 1. )
00419                    {
00420                       RPCRecHit RPCPoint(*rpcroll2,0,LocalPoint(locx,locy,locz));
00421 
00422                       RPCGeomServ servId(*rpcroll2);
00423                       if(debug) std::cout << "3\t Barrel Expected RPC " << servId.name().c_str() <<
00424                         " \tLocalposition X: " << locx << ", Y: "<< locy << " GlobalPosition(x,y,z) (" << rpcGPX <<", "<< rpcGPY <<", " << rpcGPZ << ")"<< std::endl;        
00425                       RPCPointVector.clear();
00426                       RPCPointVector.push_back(RPCPoint);
00427                       _ThePoints->put(*rpcroll2,RPCPointVector.begin(),RPCPointVector.end());
00428                       rpcput.push_back(*rpcroll2);
00429                    }
00430                  }
00431               }
00432            }
00433         }
00434         else if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::CSC)
00435         {
00436            const GeomDet *geomDet4 =  cscGeo->idToDet(dtcscid);
00437            const CSCLayer *csclayer = dynamic_cast<const CSCLayer *>(geomDet4);
00438         
00439            if(csclayer) for(Trajectories::const_iterator trajectory = trajectories.begin(); trajectory != trajectories.end(); ++trajectory)
00440            {
00441               const BoundPlane & CSCSurface = csclayer->surface();
00442               const GlobalPoint dcPoint = CSCSurface.toGlobal(LocalPoint(0.,0.,0.));
00443 
00444               TrajectoryMeasurement tMt = trajectory->closestMeasurement(dcPoint);
00445               TrajectoryStateOnSurface upd2 = (tMt).updatedState();
00446               if(upd2.isValid())  
00447               {
00448                  TrajectoryStateOnSurface ptss =  thePropagator->propagate(upd2, rpcGeo->idToDet(*rpcroll2)->surface());
00449                  if(ptss.isValid()) if( ValidRPCSurface(*rpcroll2, ptss.localPosition(), iSetup))
00450                  {
00451                    float rpcGPX = ptss.globalPosition().x();
00452                    float rpcGPY = ptss.globalPosition().y();
00453                    float rpcGPZ = ptss.globalPosition().z();
00454 
00455                    if(tInX > rpcGPX || tOuX < rpcGPX ) continue;
00456                    if(tInY > rpcGPY || tOuY < rpcGPY ) continue;
00457                    if(tInZ > rpcGPZ || tOuZ < rpcGPZ ) continue;
00458 
00459                    RPCDetId rpcid(*rpcroll2); 
00460                    const GeomDet *geomDet3 = rpcGeo->idToDet(*rpcroll2);
00461                    const RPCRoll *aroll = dynamic_cast<const RPCRoll *>(geomDet3);
00462                    const TrapezoidalStripTopology* top_=dynamic_cast<const TrapezoidalStripTopology*>(&(aroll->topology()));
00463                    LocalPoint xmin = top_->localPosition(0.);
00464                    LocalPoint xmax = top_->localPosition((float)aroll->nstrips());
00465                    float rsize = fabs( xmax.x()-xmin.x() );
00466                    float stripl = top_->stripLength();
00467                    //float stripw = top_->pitch();
00468                   
00469                    float eyr=1;
00471                    float locx = ptss.localPosition().x(), locy = ptss.localPosition().y(), locz= ptss.localPosition().z();
00472                    if( locx < rsize*eyr && locy < stripl*eyr && locz < 1. )
00473                    {
00474                       RPCRecHit RPCPoint(*rpcroll2,0,LocalPoint(locx,locy,locz));
00475                       RPCGeomServ servId(*rpcroll2);
00476                       if(debug) std::cout << "3\t Forward Expected RPC " << servId.name().c_str() <<
00477                         " \tLocalposition X: " << locx << ", Y: "<< locy << " GlobalPosition(x,y,z) (" << rpcGPX <<", "<< rpcGPY <<", " << rpcGPZ << ")"<< std::endl;
00478                       RPCPointVector.clear();
00479                       RPCPointVector.push_back(RPCPoint);
00480                       _ThePoints->put(*rpcroll2,RPCPointVector.begin(),RPCPointVector.end());
00481                       rpcput.push_back(*rpcroll2);
00482                   
00483                    }
00484                  }
00485               }
00486            }       
00487         }
00488     }
00489  }
00490  if(debug) std::cout << "last steps OK!! " << std::endl; 
00491 }
00492 }
00493 
00494 TracktoRPC::~TracktoRPC(){
00495 }
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/RecoLocalMuon/RPCRecHit/src/TracktoRPC.cc
