Reference Manual

Go to the documentation of this file.
00001 //-------------------------------------------------
00002 //
00017 //
00018 //--------------------------------------------------
00019 
00020 // Class Header
00021 #include "RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGenerator.h"
00022 
00023 // Data Formats 
00024 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeed.h"
00025 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
00026 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
00027 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00028 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00029 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTExtendedCand.h"
00030 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"
00031 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h"
00032 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
00033 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
00034 #include "DataFormats/Common/interface/Handle.h"
00035 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
00036 #include "DataFormats/Math/interface/deltaR.h"
00037 
00038 #include "CLHEP/Vector/ThreeVector.h"
00039 
00040 #include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"
00041 
00042 // Framework
00043 #include "FWCore/Framework/interface/EventSetup.h"
00044 #include "FWCore/Framework/interface/Event.h"
00045 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00046 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00047 
00048 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
00049 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
00050 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00051 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00052 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00053 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
00054 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00055 
00056 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
00057 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00058 
00059 using namespace std;
00060 using namespace edm;
00061 using namespace l1extra;
00062 
00063 // constructors
00064 L2MuonSeedGenerator::L2MuonSeedGenerator(const edm::ParameterSet& iConfig) : 
00065   theSource(iConfig.getParameter<InputTag>("InputObjects")),
00066   theL1GMTReadoutCollection(iConfig.getParameter<InputTag>("GMTReadoutCollection")),
00067   thePropagatorName(iConfig.getParameter<string>("Propagator")),
00068   theL1MinPt(iConfig.getParameter<double>("L1MinPt")),
00069   theL1MaxEta(iConfig.getParameter<double>("L1MaxEta")),
00070   theL1MinQuality(iConfig.getParameter<unsigned int>("L1MinQuality")),
00071   useOfflineSeed(iConfig.getUntrackedParameter<bool>("UseOfflineSeed", false)){
00072   
00073   // service parameters
00074   ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
00075   
00076   // the services
00077   theService = new MuonServiceProxy(serviceParameters);
00078 
00079   // the estimator
00080   theEstimator = new Chi2MeasurementEstimator(10000.);
00081 
00082   if(useOfflineSeed)
00083     theOfflineSeedLabel = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
00084 
00085   produces<L2MuonTrajectorySeedCollection>(); 
00086 }
00087 
00088 // destructor
00089 L2MuonSeedGenerator::~L2MuonSeedGenerator(){
00090   if (theService) delete theService;
00091   if (theEstimator) delete theEstimator;
00092 }
00093 
00094 void L2MuonSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00095 {
00096   const std::string metname = "Muon|RecoMuon|L2MuonSeedGenerator";
00097   MuonPatternRecoDumper debug;
00098 
00099   auto_ptr<L2MuonTrajectorySeedCollection> output(new L2MuonTrajectorySeedCollection());
00100   
00101   // Muon particles and GMT readout collection
00102   edm::Handle<L1MuGMTReadoutCollection> gmtrc_handle;
00103   iEvent.getByLabel(theL1GMTReadoutCollection,gmtrc_handle);
00104   L1MuGMTReadoutRecord const& gmtrr = gmtrc_handle.product()->getRecord(0);
00105 
00106   edm::Handle<L1MuonParticleCollection> muColl;
00107   iEvent.getByLabel(theSource, muColl);
00108   LogTrace(metname) << "Number of muons " << muColl->size() << endl;
00109 
00110   edm::Handle<edm::View<TrajectorySeed> > offlineSeedHandle;
00111   vector<int> offlineSeedMap;
00112   if(useOfflineSeed) {
00113     iEvent.getByLabel(theOfflineSeedLabel, offlineSeedHandle);
00114     LogTrace(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl;
00115     offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
00116   }
00117 
00118   L1MuonParticleCollection::const_iterator it;
00119   L1MuonParticleRef::key_type l1ParticleIndex = 0;
00120 
00121   for(it = muColl->begin(); it != muColl->end(); ++it,++l1ParticleIndex) {
00122     
00123     const L1MuGMTExtendedCand muonCand = (*it).gmtMuonCand();
00124     unsigned int quality = 0;
00125     bool valid_charge = false;;
00126 
00127     if ( muonCand.empty() ) {
00128       LogWarning(metname) << "L2MuonSeedGenerator: WARNING, no L1MuGMTCand! " << endl;
00129       LogWarning(metname) << "L2MuonSeedGenerator:   this should make sense only within MC tests" << endl;
00130       // FIXME! Temporary to handle the MC input
00131       quality = 7;
00132       valid_charge = true;
00133     }
00134     else {
00135       quality =  muonCand.quality();
00136       valid_charge = muonCand.charge_valid();
00137     }
00138     
00139     float pt    =  (*it).pt();
00140     float eta   =  (*it).eta();
00141     float theta =  2*atan(exp(-eta));
00142     float phi   =  (*it).phi();      
00143     int charge  =  (*it).charge();
00144     // Set charge=0 for the time being if the valid charge bit is zero
00145     if (!valid_charge) charge = 0;
00146     bool barrel = !(*it).isForward();
00147 
00148     // Get a better eta and charge from regional information
00149     // Phi has the same resolution in GMT than regionally, is not it?
00150     if ( !(muonCand.empty()) ) {
00151       int idx = -1;
00152       vector<L1MuRegionalCand> rmc;
00153       if ( !muonCand.isRPC() ) {
00154             idx = muonCand.getDTCSCIndex();
00155             if (muonCand.isFwd()) rmc = gmtrr.getCSCCands();
00156             else rmc = gmtrr.getDTBXCands();
00157       } else {
00158             idx = muonCand.getRPCIndex();
00159             if (muonCand.isFwd()) rmc = gmtrr.getFwdRPCCands();
00160             else rmc = gmtrr.getBrlRPCCands();
00161       }
00162       if (idx>=0) {
00163             eta = rmc[idx].etaValue();
00164             //phi = rmc[idx].phiValue();
00165             // Use this charge if the valid charge bit is zero
00166             if (!valid_charge) charge = rmc[idx].chargeValue();
00167       }
00168     }
00169 
00170     if ( pt < theL1MinPt || fabs(eta) > theL1MaxEta ) continue;
00171     
00172     LogTrace(metname) << "New L2 Muon Seed";
00173     LogTrace(metname) << "Pt = " << pt << " GeV/c";
00174     LogTrace(metname) << "eta = " << eta;
00175     LogTrace(metname) << "theta = " << theta << " rad";
00176     LogTrace(metname) << "phi = " << phi << " rad";
00177     LogTrace(metname) << "charge = "<< charge;
00178     LogTrace(metname) << "In Barrel? = "<< barrel;
00179     
00180     if ( quality <= theL1MinQuality ) continue;
00181     LogTrace(metname) << "quality = "<< quality; 
00182     
00183     // Update the services
00184     theService->update(iSetup);
00185 
00186     const DetLayer *detLayer = 0;
00187     float radius = 0.;
00188   
00189     CLHEP::Hep3Vector vec(0.,1.,0.);
00190     vec.setTheta(theta);
00191     vec.setPhi(phi);
00192         
00193     // Get the det layer on which the state should be put
00194     if ( barrel ){
00195       LogTrace(metname) << "The seed is in the barrel";
00196       
00197       // MB2
00198       DetId id = DTChamberId(0,2,0);
00199       detLayer = theService->detLayerGeometry()->idToLayer(id);
00200       LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
00201       
00202       const BoundSurface* sur = &(detLayer->surface());
00203       const BoundCylinder* bc = dynamic_cast<const BoundCylinder*>(sur);
00204 
00205       radius = fabs(bc->radius()/sin(theta));
00206 
00207       LogTrace(metname) << "radius "<<radius;
00208 
00209       if ( pt < 3.5 ) pt = 3.5;
00210     }
00211     else { 
00212       LogTrace(metname) << "The seed is in the endcap";
00213       
00214       DetId id;
00215       // ME2
00216       if ( theta < Geom::pi()/2. )
00217         id = CSCDetId(1,2,0,0,0); 
00218       else
00219         id = CSCDetId(2,2,0,0,0); 
00220       
00221       detLayer = theService->detLayerGeometry()->idToLayer(id);
00222       LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
00223 
00224       radius = fabs(detLayer->position().z()/cos(theta));      
00225       
00226       if( pt < 1.0) pt = 1.0;
00227     }
00228         
00229     vec.setMag(radius);
00230     
00231     GlobalPoint pos(vec.x(),vec.y(),vec.z());
00232       
00233     GlobalVector mom(pt*cos(phi), pt*sin(phi), pt*cos(theta)/sin(theta));
00234 
00235     GlobalTrajectoryParameters param(pos,mom,charge,&*theService->magneticField());
00236     AlgebraicSymMatrix55 mat;
00237     
00238     mat[0][0] = (0.25/pt)*(0.25/pt);  // sigma^2(charge/abs_momentum)
00239     if ( !barrel ) mat[0][0] = (0.4/pt)*(0.4/pt);
00240 
00241     //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
00242     if (!valid_charge) mat[0][0] = (1./pt)*(1./pt);
00243     
00244     mat[1][1] = 0.05*0.05;        // sigma^2(lambda)
00245     mat[2][2] = 0.2*0.2;          // sigma^2(phi)
00246     mat[3][3] = 20.*20.;          // sigma^2(x_transverse))
00247     mat[4][4] = 20.*20.;          // sigma^2(y_transverse))
00248     
00249     CurvilinearTrajectoryError error(mat);
00250 
00251     const FreeTrajectoryState state(param,error);
00252    
00253     LogTrace(metname) << "Free trajectory State from the parameters";
00254     LogTrace(metname) << debug.dumpFTS(state);
00255 
00256     // Propagate the state on the MB2/ME2 surface
00257     TrajectoryStateOnSurface tsos = theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
00258    
00259     LogTrace(metname) << "State after the propagation on the layer";
00260     LogTrace(metname) << debug.dumpLayer(detLayer);
00261     LogTrace(metname) << debug.dumpFTS(state);
00262 
00263     if (tsos.isValid()) {
00264       // Get the compatible dets on the layer
00265       std::vector< pair<const GeomDet*,TrajectoryStateOnSurface> > 
00266         detsWithStates = detLayer->compatibleDets(tsos, 
00267                                                   *theService->propagator(thePropagatorName), 
00268                                                   *theEstimator);   
00269       if (detsWithStates.size()){
00270 
00271         TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
00272         const GeomDet *newTSOSDet = detsWithStates.front().first;
00273         
00274         LogTrace(metname) << "Most compatible det";
00275         LogTrace(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
00276 
00277         LogDebug(metname) << "L1 info: Det and State:";
00278         LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
00279 
00280         if (newTSOS.isValid()){
00281 
00282           //LogDebug(metname) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", " 
00283           //                  << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")";
00284           LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag() << ", phi=" 
00285                             << newTSOS.globalPosition().phi() << ", eta=" << newTSOS.globalPosition().eta() << ")";
00286           LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge()*newTSOS.globalMomentum().perp() << ", phi=" 
00287                             << newTSOS.globalMomentum().phi() << ", eta=" << newTSOS.globalMomentum().eta() << ")";
00288 
00289           //LogDebug(metname) << "State on it";
00290           //LogDebug(metname) << debug.dumpTSOS(newTSOS);
00291 
00292           //PTrajectoryStateOnDet seedTSOS;
00293           edm::OwnVector<TrackingRecHit> container;
00294 
00295           if(useOfflineSeed) {
00296             const TrajectorySeed *assoOffseed = 
00297               associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS);
00298 
00299             if(assoOffseed!=0) {
00300               PTrajectoryStateOnDet const & seedTSOS = assoOffseed->startingState();
00301               TrajectorySeed::const_iterator 
00302                 tsci  = assoOffseed->recHits().first,
00303                 tscie = assoOffseed->recHits().second;
00304               for(; tsci!=tscie; ++tsci) {
00305                 container.push_back(*tsci);
00306               }
00307               output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
00308                                                      L1MuonParticleRef(muColl,l1ParticleIndex)));
00309             }
00310             else {
00311               // convert the TSOS into a PTSOD
00312               PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( newTSOS,newTSOSDet->geographicalId().rawId());
00313               output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
00314                                                      L1MuonParticleRef(muColl,l1ParticleIndex)));
00315             }
00316           }
00317           else {
00318             // convert the TSOS into a PTSOD
00319             PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( newTSOS,newTSOSDet->geographicalId().rawId());
00320             output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
00321                                                    L1MuonParticleRef(muColl,l1ParticleIndex)));
00322           }
00323         }
00324       }
00325     }
00326   }
00327   
00328   iEvent.put(output);
00329 }
00330 
00331 
00332 // FIXME: does not resolve ambiguities yet! 
00333 const TrajectorySeed* L2MuonSeedGenerator::associateOfflineSeedToL1( edm::Handle<edm::View<TrajectorySeed> > & offseeds, 
00334                                                                      std::vector<int> & offseedMap, 
00335                                                                      TrajectoryStateOnSurface & newTsos) {
00336 
00337   const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGenerator";
00338   MuonPatternRecoDumper debugtmp;
00339 
00340   edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
00341   const TrajectorySeed *selOffseed = 0;
00342   double bestDr = 99999.;
00343   unsigned int nOffseed(0);
00344   int lastOffseed(-1);
00345 
00346   for(offseed=offseeds->begin(); offseed!=endOffseed; ++offseed, ++nOffseed) {
00347     if(offseedMap[nOffseed]!=0) continue;
00348     GlobalPoint glbPos = theService->trackingGeometry()->idToDet(offseed->startingState().detId())->surface().toGlobal(offseed->startingState().parameters().position());
00349     GlobalVector glbMom = theService->trackingGeometry()->idToDet(offseed->startingState().detId())->surface().toGlobal(offseed->startingState().parameters().momentum());
00350 
00351     // Preliminary check
00352     double preDr = deltaR( newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi() );
00353     if(preDr > 1.0) continue;
00354 
00355     const FreeTrajectoryState offseedFTS(glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField()); 
00356     TrajectoryStateOnSurface offseedTsos = theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
00357     LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
00358     LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
00359     //LogDebug(metlabel) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", " 
00360     //                   << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")" << std::endl;
00361     LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" 
00362                        << offseedFTS.position().phi() << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
00363     LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge()*offseedFTS.momentum().perp() << ", phi=" 
00364                        << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")" << std::endl << std::endl;
00365     //LogDebug(metlabel) << debugtmp.dumpFTS(offseedFTS) << std::endl;
00366 
00367     if(offseedTsos.isValid()) {
00368       LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
00369       //LogDebug(metlabel) << "(x, y, z) = (" << offseedTsos.globalPosition().x() << ", " 
00370       //                   << offseedTsos.globalPosition().y() << ", " << offseedTsos.globalPosition().z() << ")" << std::endl;
00371       LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag() << ", phi=" 
00372                          << offseedTsos.globalPosition().phi() << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
00373       LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge()*offseedTsos.globalMomentum().perp() << ", phi=" 
00374                          << offseedTsos.globalMomentum().phi() << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl << std::endl;
00375       //LogDebug(metlabel) << debugtmp.dumpTSOS(offseedTsos) << std::endl;
00376       double newDr = deltaR( newTsos.globalPosition().eta(),     newTsos.globalPosition().phi(), 
00377                              offseedTsos.globalPosition().eta(), offseedTsos.globalPosition().phi() );
00378       LogDebug(metlabel) << "   -- DR = " << newDr << std::endl;
00379       if( newDr<0.3 && newDr<bestDr ) {
00380         LogDebug(metlabel) << "          --> OK! " << newDr << std::endl << std::endl;
00381         selOffseed = &*offseed;
00382         bestDr = newDr;
00383         offseedMap[nOffseed] = 1;
00384         if(lastOffseed>-1) offseedMap[lastOffseed] = 0;
00385         lastOffseed = nOffseed;
00386       }
00387       else {
00388         LogDebug(metlabel) << "          --> Rejected. " << newDr << std::endl << std::endl;
00389       }
00390     }
00391     else {
00392       LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
00393     }
00394   }
00395 
00396   return selOffseed;
00397 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGenerator.cc
