#include <MuonTrackResidualAnalyzer.h>

Inheritance diagram for MuonTrackResidualAnalyzer:

Classes
struct	RadiusComparatorInOut

Public Types
enum	EtaRange { all, barrel, endcap }

Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

typedef WorkerT< EDAnalyzer >	WorkerType

Public Member Functions
void	analyze (const edm::Event &event, const edm::EventSetup &eventSetup)

virtual void	beginJob ()

virtual void	endJob ()

	MuonTrackResidualAnalyzer (const edm::ParameterSet &pset)
	Constructor. More...

virtual	~MuonTrackResidualAnalyzer ()
	Destructor. More...

Public Member Functions inherited from edm::EDAnalyzer
	EDAnalyzer ()

std::string	workerType () const

virtual	~EDAnalyzer ()

Private Member Functions
void	computeResolution (Trajectory &trajectory, std::map< DetId, const PSimHit * > &hitIdMap, HResolution1DRecHit *histos)

bool	isInTheAcceptance (double eta)

std::map< DetId, const PSimHit * >	mapMuSimHitsPerId (edm::Handle< edm::PSimHitContainer > dtSimhits, edm::Handle< edm::PSimHitContainer > cscSimhits, edm::Handle< edm::PSimHitContainer > rpcSimhits)

void	mapMuSimHitsPerId (edm::Handle< edm::PSimHitContainer > simhits, std::map< DetId, const PSimHit * > &hitIdMap)

Private Attributes
edm::InputTag	cscSimHitLabel

DQMStore *	dbe_

std::string	dirName_

edm::InputTag	dtSimHitLabel

HResolution1DRecHit *	h1DRecHitRes

HResolution1DRecHit *	h1DSimHitRes

MonitorElement *	hDeltaPtVsEtaSim

MonitorElement *	hDeltaPtVsEtaSim2

MonitorElement *	hDPtRef

MonitorElement *	hSimHitsPerTrack

MonitorElement *	hSimHitsPerTrackVsEta

std::string	out

edm::InputTag	rpcSimHitLabel

edm::InputTag	theDataType

MeasurementEstimator *	theEstimator

EtaRange	theEtaRange

int	theMuonSimHitNumberPerEvent

edm::InputTag	theMuonTrackLabel

edm::InputTag	theSeedCollectionLabel

MuonServiceProxy *	theService

std::vector< const PSimHit * >	theSimHitContainer

unsigned int	theSimTkId

KFUpdator *	theUpdator

Additional Inherited Members
Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

Protected Member Functions inherited from edm::EDAnalyzer
CurrentProcessingContext const *	currentContext () const

Detailed Description

No description available.

Date:: 2010/02/20 21:02:36

Revision:: 1.7

Author: R. Bellan - INFN Torino ricca.nosp@m.rdo..nosp@m.bella.nosp@m.n@ce.nosp@m.rn.ch

Definition at line 44 of file MuonTrackResidualAnalyzer.h.

Member Enumeration Documentation

enum MuonTrackResidualAnalyzer::EtaRange

Enumerator
all
barrel
endcap

Definition at line 47 of file MuonTrackResidualAnalyzer.h.

47 {all,barrel,endcap};

MuonTrackResidualAnalyzer::all

Definition: MuonTrackResidualAnalyzer.h:47

MuonTrackResidualAnalyzer::endcap

Definition: MuonTrackResidualAnalyzer.h:47

MuonTrackResidualAnalyzer::barrel

Definition: MuonTrackResidualAnalyzer.h:47

Constructor & Destructor Documentation

MuonTrackResidualAnalyzer::MuonTrackResidualAnalyzer ( const edm::ParameterSet & pset )

Constructor.

Definition at line 38 of file MuonTrackResidualAnalyzer.cc.

References Chi2MeasurementEstimatorESProducer_cfi::Chi2MeasurementEstimator, dbe_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), LogDebug, MuonServiceProxy_cff::MuonServiceProxy, cppFunctionSkipper::operator, and dbtoconf::out.

                                                                                {
   
     // service parameters
   ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
   // the services
   theService = new MuonServiceProxy(serviceParameters);
   
   theMuonTrackLabel = pset.getParameter<InputTag>("MuonTrack");
   theSeedCollectionLabel = pset.getParameter<InputTag>("MuonSeed");
 
   cscSimHitLabel = pset.getParameter<InputTag>("CSCSimHit");
   dtSimHitLabel = pset.getParameter<InputTag>("DTSimHit");
   rpcSimHitLabel = pset.getParameter<InputTag>("RPCSimHit");
 
   dbe_ = edm::Service<DQMStore>().operator->();
   out = pset.getUntrackedParameter<string>("rootFileName");
   dirName_ = pset.getUntrackedParameter<std::string>("dirName");
   
   // Sim or Real
   theDataType = pset.getParameter<InputTag>("DataType"); 
   if(theDataType.label() != "RealData" && theDataType.label() != "SimData")
     LogDebug("MuonTrackResidualAnalyzer")<<"Error in Data Type!!";
 
   theEtaRange = (EtaRange) pset.getParameter<int>("EtaRange");
 
   theUpdator = new KFUpdator();
   theEstimator = new Chi2MeasurementEstimator(100000.);
 
   theMuonSimHitNumberPerEvent = 0;
 }

MuonTrackResidualAnalyzer::~MuonTrackResidualAnalyzer ( )

virtual

Destructor.

Definition at line 70 of file MuonTrackResidualAnalyzer.cc.

                                                      {
   delete theUpdator;
   delete theEstimator;
   delete theService;
 }

Member Function Documentation

void MuonTrackResidualAnalyzer::analyze	(	const edm::Event &	event,
		const edm::EventSetup &	eventSetup
	)

virtual

Implements edm::EDAnalyzer.

Definition at line 117 of file MuonTrackResidualAnalyzer.cc.

References abs, debug, LogDebug, LogTrace, mag(), query::result, and GeomDet::toGlobal().

                                                                                               {
   LogDebug("MuonTrackResidualAnalyzer")<<"Analyze";
 
   // Update the services
   theService->update(eventSetup);
   MuonPatternRecoDumper debug;
   theMuonSimHitNumberPerEvent = 0;
 
   // Get the SimHit collection from the event
   Handle<PSimHitContainer> dtSimHits;
   event.getByLabel(dtSimHitLabel.instance(),dtSimHitLabel.label(), dtSimHits);
 
   Handle<PSimHitContainer> cscSimHits;
   event.getByLabel(cscSimHitLabel.instance(),cscSimHitLabel.label(), cscSimHits);
 
   Handle<PSimHitContainer> rpcSimHits;
   event.getByLabel(rpcSimHitLabel.instance(),rpcSimHitLabel.label(), rpcSimHits);
 
   Handle<SimTrackContainer> simTracks;
 
   // FIXME Add the tracker one??
 
   // Map simhits per DetId
   map<DetId, const PSimHit* > muonSimHitsPerId =
     mapMuSimHitsPerId(dtSimHits,cscSimHits,rpcSimHits);
   
   hSimHitsPerTrack->Fill(theMuonSimHitNumberPerEvent);
   
   double etaSim=0;
 
   if(theDataType.label() == "SimData"){
     
     // Get the SimTrack collection from the event
     event.getByLabel(theDataType.instance(),simTracks);
     
     // Loop over the Sim tracks
     SimTrackContainer::const_iterator simTrack;
     
     for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
       if (abs((*simTrack).type()) == 13){
     hSimHitsPerTrackVsEta->Fill((*simTrack).momentum().eta(),theMuonSimHitNumberPerEvent);
     etaSim = (*simTrack).momentum().eta();
     theSimTkId = (*simTrack).trackId();
 
       }
   }
 
   
   // Get the RecTrack collection from the event
   Handle<reco::TrackCollection> muonTracks;
   event.getByLabel(theMuonTrackLabel, muonTracks);
 
   reco::TrackCollection::const_iterator muonTrack;
   
   // Loop over the Rec tracks
   for (muonTrack = muonTracks->begin(); muonTrack != muonTracks->end(); ++muonTrack) {
     
     reco::TransientTrack track(*muonTrack,&*theService->magneticField(),theService->trackingGeometry());   
 
     TrajectoryStateOnSurface outerTSOS = track.outermostMeasurementState();
     TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
 
     TransientTrackingRecHit::ConstRecHitContainer result;
     
     // SimHit Energy loss analysis  
     double momAtEntry = -150., momAtExit = -150.;
 
     if(theSimHitContainer.size()>1){
 
       const GeomDet *geomDetA = theService->trackingGeometry()->idToDet(DetId(theSimHitContainer.front()->detUnitId()));
       double distA = geomDetA->toGlobal(theSimHitContainer.front()->localPosition()).mag();
 
       const GeomDet *geomDetB = theService->trackingGeometry()->idToDet(DetId(theSimHitContainer.back()->detUnitId()));
       double distB = geomDetB->toGlobal(theSimHitContainer.back()->localPosition()).mag();
 
       LogTrace("MuonTrackResidualAnalyzer")<<"Inner SimHit: "<<theSimHitContainer.front()->particleType()
                        <<" Pt: "<<theSimHitContainer.front()->momentumAtEntry().perp()
                        <<" E: "<<theSimHitContainer.front()->momentumAtEntry().perp()
                        <<" R: "<<distA<<endl;
       LogTrace("MuonTrackResidualAnalyzer")<<"Outer SimHit: "<<theSimHitContainer.back()->particleType()
                        <<" Pt: "<<theSimHitContainer.back()->momentumAtEntry().perp()
                        <<" E: "<<theSimHitContainer.front()->momentumAtEntry().perp()
                        <<" R: "<<distB<<endl;
       
       momAtEntry = theSimHitContainer.front()->momentumAtEntry().perp();
       momAtExit =  theSimHitContainer.back()->momentumAtEntry().perp();
     }
     
     trackingRecHit_iterator rhFirst = track.recHitsBegin();
     trackingRecHit_iterator rhLast = track.recHitsEnd()-1;
     map<DetId,const PSimHit*>::const_iterator itFirst = muonSimHitsPerId.find((*rhFirst)->geographicalId());
     map<DetId,const PSimHit*>::const_iterator itLast = muonSimHitsPerId.find((*rhLast)->geographicalId());
     
     double momAtEntry2 = -150, momAtExit2 = -150.;
     if (itFirst != muonSimHitsPerId.end() )
       momAtEntry2 = itFirst->second->momentumAtEntry().perp();
     else {
       LogDebug("MuonTrackResidualAnalyzer")<<"No first sim hit found";
       ++rhFirst;
       itFirst = muonSimHitsPerId.find((*rhFirst)->geographicalId());
       if (itFirst != muonSimHitsPerId.end() )
     momAtEntry2 = itFirst->second->momentumAtEntry().perp();
       else{
     LogDebug("MuonTrackResidualAnalyzer")<<"No second sim hit found";
     // continue;
       }
     }
     
     if (itLast != muonSimHitsPerId.end() )
       momAtExit2 = itLast->second->momentumAtEntry().perp();
     else {
       LogDebug("MuonTrackResidualAnalyzer")<<"No last sim hit found";
       --rhLast;
       itLast = muonSimHitsPerId.find((*rhLast)->geographicalId());
       if (itLast != muonSimHitsPerId.end() )
     momAtExit2 = itLast->second->momentumAtEntry().perp();
       else{
     LogDebug("MuonTrackResidualAnalyzer")<<"No last but one sim hit found";
     // continue;
       }
     }
     
     if(etaSim){
       if(momAtEntry >=0 && momAtExit >= 0)
     hDeltaPtVsEtaSim->Fill(etaSim,momAtEntry-momAtExit);
       if(momAtEntry2 >=0 && momAtExit2 >= 0)
     hDeltaPtVsEtaSim2->Fill(etaSim,momAtEntry2-momAtExit2);
     }
     else
       LogDebug("MuonTrackResidualAnalyzer")<<"NO SimTrack'eta";
     //
 
     // computeResolution(trajectoryBW,muonSimHitsPerId,h1DSimHitRes);
     // computeResolution(smoothed,muonSimHitsPerId,h1DSimHitRes);
     
   }
 }

void MuonTrackResidualAnalyzer::beginJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 77 of file MuonTrackResidualAnalyzer.cc.

References algo, DQMStore::book1D(), DQMStore::book2D(), DQMStore::cd(), dbe_, TrackerOfflineValidation_Dqm_cff::dirName, edm::InputTag::instance(), edm::InputTag::label(), LogDebug, edm::InputTag::process(), python.rootplot.root2matplotlib::replace(), DQMStore::setCurrentFolder(), and DQMStore::showDirStructure().

                                         {
   LogDebug("MuonTrackResidualAnalyzer")<<"Begin Job";
   
   dbe_->showDirStructure();
   
   dbe_->cd();
   InputTag algo = theMuonTrackLabel;
   string dirName=dirName_;
   if (algo.process()!="")
     dirName+=algo.process()+"_";
   if(algo.label()!="")
     dirName+=algo.label()+"_";
   if(algo.instance()!="")
     dirName+=algo.instance()+"";
   if (dirName.find("Tracks")<dirName.length()){
     dirName.replace(dirName.find("Tracks"),6,"");
   }
   std::replace(dirName.begin(), dirName.end(), ':', '_');
   dbe_->setCurrentFolder(dirName.c_str());
   
   
   hDPtRef = dbe_->book1D("DeltaPtRef","P^{in}_{t}-P^{in ref}",10000,-20,20);
   
   // Resolution wrt the 1D Rec Hits
   //  h1DRecHitRes = new HResolution1DRecHit("TotalRec");
   
   // Resolution wrt the 1d Sim Hits
   //  h1DSimHitRes = new HResolution1DRecHit("TotalSim");
   
   hSimHitsPerTrack  = dbe_->book1D("SimHitsPerTrack","Number of sim hits per track",55,0,55);
   hSimHitsPerTrackVsEta  = dbe_->book2D("SimHitsPerTrackVsEta","Number of sim hits per track VS #eta",120,-3.,3.,55,0,55);
   hDeltaPtVsEtaSim = dbe_->book2D("DeltaPtVsEtaSim","#Delta P_{t} vs #eta gen, sim quantity",120,-3.,3.,500,-250.,250.);
   hDeltaPtVsEtaSim2 = dbe_->book2D("DeltaPtVsEtaSim2","#Delta P_{t} vs #eta gen, sim quantity",120,-3.,3.,500,-250.,250.);
 }

void MuonTrackResidualAnalyzer::computeResolution	(	Trajectory &	trajectory,
		std::map< DetId, const PSimHit * > &	hitIdMap,
		HResolution1DRecHit *	histos
	)

private

Definition at line 329 of file MuonTrackResidualAnalyzer.cc.

References abs, gather_cfg::cout, data, PSimHit::detUnitId(), diffTreeTool::diff, PV3DBase< T, PVType, FrameType >::eta(), HResolution1DRecHit::Fill(), PSimHit::localPosition(), Trajectory::measurements(), PSimHit::momentumAtEntry(), PV3DBase< T, PVType, FrameType >::phi(), GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                       {
   
   Trajectory::DataContainer data = trajectory.measurements();
 
   for(Trajectory::DataContainer::const_iterator datum = data.begin();
       datum != data.end(); ++datum){
     
     GlobalPoint fitPoint = datum->updatedState().globalPosition();
 
     // FIXME!
     //     double errX = datum->updatedState().cartesianError().matrix()[0][0];
     //     double errY = datum->updatedState().cartesianError().matrix()[1][1];
     //     double errZ = datum->updatedState().cartesianError().matrix()[2][2];
     //
     double errX = datum->updatedState().localError().matrix()(3,3);
     double errY = datum->updatedState().localError().matrix()(4,4);
     double errZ = 1.;
 
     map<DetId,const PSimHit*>::const_iterator it = hitIdMap.find(datum->recHit()->geographicalId());
  
     
     if (it == hitIdMap.end() ) continue; // FIXME! Put a counter
     
     const PSimHit* simhit = it->second;
     
     LocalPoint simHitPoint = simhit->localPosition();
 
     const GeomDet* geomDet = theService->trackingGeometry()->idToDet(DetId(simhit->detUnitId()));
 
     LocalPoint fitLocalPoint = geomDet->toLocal(fitPoint);
     
     LocalVector diff = fitLocalPoint-simHitPoint;
 
     
     cout << "SimHit position "<< simHitPoint << endl;
     cout << "Fit position "<<  fitLocalPoint << endl;
     cout << "Fit position2 "<< datum->updatedState().localPosition() << endl;
     cout << "Errors on the fit position: (" << errX << "," << errY << "," << errZ << ")"<<endl;
     cout << "Resolution on x: " << diff.x()/abs(simHitPoint.x()) << endl;
     cout << "Resolution on y: " << diff.y()/abs(simHitPoint.y()) << endl;
     cout << "Resolution on z: " << diff.z()/abs(simHitPoint.z()) << endl;
 
     cout << "Eta direction: "<< simhit->momentumAtEntry().eta() <<" eta position: " << simHitPoint.eta() << endl;
     cout << "Phi direction: "<< simhit->momentumAtEntry().phi() <<" phi position: " << simHitPoint.phi() << endl;
 
     
     histos->Fill( simHitPoint.x(), simHitPoint.y(), simHitPoint.z(), 
           diff.x(), diff.y(), diff.z(),
           errX,  errY, errZ,
           simhit->momentumAtEntry().eta(), simhit->momentumAtEntry().phi());
           // simHitPoint.eta(), simHitPoint.phi() ); // FIXME!
   }
 
 
 }

void MuonTrackResidualAnalyzer::endJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 112 of file MuonTrackResidualAnalyzer.cc.

References dbe_, dbtoconf::out, and DQMStore::save().

                                       {
   if ( out.size() != 0 && dbe_ ) dbe_->save(out);
 }

bool MuonTrackResidualAnalyzer::isInTheAcceptance ( double eta )

private

Definition at line 255 of file MuonTrackResidualAnalyzer.cc.

References abs, cond::ecalcond::all, Reference_intrackfit_cff::barrel, Reference_intrackfit_cff::endcap, and LogDebug.

                                                            {
   switch(theEtaRange){
   case all:
     return ( abs(eta) <= 2.4 ) ? true : false;
   case barrel:
     return ( abs(eta) < 1.1 ) ? true : false;
   case endcap:
     return ( abs(eta) >= 1.1 && abs(eta) <= 2.4 ) ? true : false;  
   default:
     {LogDebug("MuonTrackResidualAnalyzer")<<"No correct Eta range selected!! "; return false;}
   }
 }

map< DetId, const PSimHit * > MuonTrackResidualAnalyzer::mapMuSimHitsPerId	(	edm::Handle< edm::PSimHitContainer >	dtSimhits,
		edm::Handle< edm::PSimHitContainer >	cscSimhits,
		edm::Handle< edm::PSimHitContainer >	rpcSimhits
	)

private

Definition at line 270 of file MuonTrackResidualAnalyzer.cc.

References prof2calltree::count, debug, MuonPatternRecoDumper::dumpMuonId(), LogDebug, and LogTrace.

                                                               {
   
   MuonPatternRecoDumper debug;
   
   map<DetId,const PSimHit*> hitIdMap;
   theSimHitContainer.clear();
 
   mapMuSimHitsPerId(dtSimhits,hitIdMap);
   mapMuSimHitsPerId(cscSimhits,hitIdMap);
   mapMuSimHitsPerId(rpcSimhits,hitIdMap);
 
   if(theSimHitContainer.size() >1)
     stable_sort(theSimHitContainer.begin(),theSimHitContainer.end(),RadiusComparatorInOut(theService->trackingGeometry()));
 
   LogDebug("MuonTrackResidualAnalyzer")<<"Sim Hit list";
   int count=1;
   for(vector<const PSimHit*>::const_iterator it = theSimHitContainer.begin(); 
       it != theSimHitContainer.end(); ++it){
     LogTrace("MuonTrackResidualAnalyzer")<<count 
     << " " 
     << " Process Type: " << (*it)->processType()
     << " " 
     << debug.dumpMuonId(DetId( (*it)->detUnitId() ))<<endl;
   }  
 
   return hitIdMap;
 }

void MuonTrackResidualAnalyzer::mapMuSimHitsPerId	(	edm::Handle< edm::PSimHitContainer >	simhits,
		std::map< DetId, const PSimHit * > &	hitIdMap
	)

private

Definition at line 301 of file MuonTrackResidualAnalyzer.cc.

References abs, MuonSubdetId::DT, LogDebug, and DetId::rawId().

                                                               {
   
   for(PSimHitContainer::const_iterator simhit = simhits->begin();
       simhit != simhits->end(); ++simhit) {
     
     if ( abs(simhit->particleType()) != 13 && theSimTkId != simhit->trackId()) continue; 
     
     theSimHitContainer.push_back(&*simhit);
     DetId id = DetId(simhit->detUnitId());
     
     if(id.subdetId() == MuonSubdetId::DT){
       DTLayerId lId(id.rawId());
       id = DetId(lId.rawId());
     }
 
     map<DetId,const PSimHit*>::const_iterator it = hitIdMap.find(id);
     
     if (it == hitIdMap.end() )
       hitIdMap[id] = &*simhit;
     else
       LogDebug("MuonTrackResidualAnalyzer")<<"TWO muons in the same sensible volume!!";
     
     ++theMuonSimHitNumberPerEvent;
   }
 }