#include <MuonTrackResidualAnalyzer.h>

Inheritance diagram for MuonTrackResidualAnalyzer:

Classes
struct	RadiusComparatorInOut
Public Types
enum	EtaRange { all, barrel, endcap }
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.
virtual	~MuonTrackResidualAnalyzer ()
	Destructor.
Private Member Functions
void	computeResolution (Trajectory &trajectory, std::map< DetId, const PSimHit * > &hitIdMap, HResolution1DRecHit *histos)
bool	isInTheAcceptance (double eta)
void	mapMuSimHitsPerId (edm::Handle< edm::PSimHitContainer > simhits, std::map< DetId, const PSimHit * > &hitIdMap)
std::map< DetId, const PSimHit * >	mapMuSimHitsPerId (edm::Handle< edm::PSimHitContainer > dtSimhits, edm::Handle< edm::PSimHitContainer > cscSimhits, edm::Handle< edm::PSimHitContainer > rpcSimhits)
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

Detailed Description

No description available.

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

Revision:: 1.7

Author:: R. Bellan - INFN Torino <riccardo.bellan@cern.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.

{all,barrel,endcap};

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, dbtoconf::out, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                               {
  
    // 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 DQMStore::book1D(), DQMStore::book2D(), DQMStore::cd(), dbe_, TrackerOfflineValidation_Dqm_cff::dirName, edm::InputTag::instance(), edm::InputTag::label(), LogDebug, edm::InputTag::process(), linker::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 551 of file DQMGenericClient.cc.

References FitSlicesYTool::getFittedMeanWithError(), FitSlicesYTool::getFittedSigmaWithError(), and j.

{
  if ( ! theDQM->dirExists(startDir) ) {
    if ( verbose_ >= 2 || (verbose_ == 1 && !isWildcardUsed_) ) {
      LogInfo("DQMGenericClient") << "computeResolution() : "
                                  << "Cannot find sub-directory " << startDir << endl;
    }
    return;
  }

  theDQM->cd();

  ME* srcME = theDQM->get(startDir+"/"+srcName);
  if ( !srcME ) {
    if ( verbose_ >= 2 || (verbose_ == 1 && !isWildcardUsed_) ) {
      LogInfo("DQMGenericClient") << "computeResolution() : "
                                  << "No source ME '" << srcName << "' found\n";
    }
    return;
  }

  TH2F* hSrc = srcME->getTH2F();
  if ( !hSrc ) {
    if ( verbose_ >= 2 || (verbose_ == 1 && !isWildcardUsed_) ) {
      LogInfo("DQMGenericClient") << "computeResolution() : "
                                  << "Cannot create TH2F from source-ME\n";
    }
    return;
  }

  const int nBin = hSrc->GetNbinsX();

  string newDir = startDir;
  string newPrefix = namePrefix;
  string::size_type shiftPos;
  if ( string::npos != (shiftPos = namePrefix.rfind('/')) ) {
    newDir += "/"+namePrefix.substr(0, shiftPos);
    newPrefix.erase(0, shiftPos+1);
  }

  theDQM->setCurrentFolder(newDir);

  float * lowedgesfloats = new float[nBin+1];
  ME* meanME;
  ME* sigmaME;
  if (hSrc->GetXaxis()->GetXbins()->GetSize())
  {
    for (int j=0; j<nBin+1; ++j)
      lowedgesfloats[j] = (float)hSrc->GetXaxis()->GetXbins()->GetAt(j);
    meanME = theDQM->book1D(newPrefix+"_Mean", titlePrefix+" Mean", nBin, lowedgesfloats);
    sigmaME = theDQM->book1D(newPrefix+"_Sigma", titlePrefix+" Sigma", nBin, lowedgesfloats);
  }
  else
  {
    meanME = theDQM->book1D(newPrefix+"_Mean", titlePrefix+" Mean", nBin,
                            hSrc->GetXaxis()->GetXmin(),
                            hSrc->GetXaxis()->GetXmax());
    sigmaME = theDQM->book1D(newPrefix+"_Sigma", titlePrefix+" Sigma", nBin,
                            hSrc->GetXaxis()->GetXmin(),
                            hSrc->GetXaxis()->GetXmax());                            
  }
  
  if (meanME && sigmaME)
  {
    meanME->setEfficiencyFlag();
    sigmaME->setEfficiencyFlag();

    if (! resLimitedFit_ ) {
      FitSlicesYTool fitTool(srcME);
      fitTool.getFittedMeanWithError(meanME);
      fitTool.getFittedSigmaWithError(sigmaME);
    } else {
      limitedFit(srcME,meanME,sigmaME);
    }
  }
  delete[] lowedgesfloats;
}

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;}
  }
}

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

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;
}