|
|
|
#include <DTSegment4DQuality.h>
Public Member Functions | |
| void | analyze (const edm::Event &event, const edm::EventSetup &eventSetup) |
| Perform the real analysis. | |
| DTSegment4DQuality (const edm::ParameterSet &pset) | |
| Constructor. | |
| void | endJob () |
| void | endLuminosityBlock (edm::LuminosityBlock const &lumiSeg, edm::EventSetup const &c) |
| virtual | ~DTSegment4DQuality () |
| Destructor. | |
Private Attributes | |
| DQMStore * | dbe_ |
| bool | debug |
| bool | doall |
| HRes4DHit * | h4DHit |
| HRes4DHit * | h4DHit_W0 |
| HRes4DHit * | h4DHit_W1 |
| HRes4DHit * | h4DHit_W2 |
| HRes4DHit * | h4DHitWS [3][4] |
| HEff4DHit * | hEff_All |
| HEff4DHit * | hEff_W0 |
| HEff4DHit * | hEff_W1 |
| HEff4DHit * | hEff_W2 |
| HEff4DHit * | hEffWS [3][4] |
| MonitorElement * | hHitMult [3][4] |
| MonitorElement * | ht0 [3][4] |
| bool | local |
| std::string | rootFileName |
| edm::InputTag | segment4DLabel |
| double | sigmaResAlpha |
| double | sigmaResBeta |
| double | sigmaResX |
| double | sigmaResY |
| edm::InputTag | simHitLabel |
Basic analyzer class which accesses 4D DTSegments and plots resolution comparing reconstructed and simulated quantities
Only true 4D segments are considered. Station 4 segments are not looked at. FIXME: Add flag to consider also segments with only phi view? Possible bias?
Residual/pull plots are filled for the reco segment with alpha closest to the simulated muon direction (defined from muon simhits in the chamber).
Efficiencies are defined as reconstructed 4D segments with alpha, beta, x, y, within 5 sigma relative to the sim muon, with sigmas specified in the config. Note that loss of even only one of the two views is considered as inefficiency!
Definition at line 43 of file DTSegment4DQuality.h.
| DTSegment4DQuality::DTSegment4DQuality | ( | const edm::ParameterSet & | pset | ) |
Constructor.
Definition at line 39 of file DTSegment4DQuality.cc.
References dbe_, debug, DTHitQualityUtils::debug, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), mergeVDriftHistosByStation::name, cppFunctionSkipper::operator, dtTPAnalyzer_cfg::rootFileName, alignCSCRings::s, and w().
{
// Get the debug parameter for verbose output
debug = pset.getUntrackedParameter<bool>("debug");
DTHitQualityUtils::debug = debug;
rootFileName = pset.getUntrackedParameter<string>("rootFileName");
// the name of the simhit collection
simHitLabel = pset.getUntrackedParameter<InputTag>("simHitLabel");
// the name of the 4D rec hit collection
segment4DLabel = pset.getUntrackedParameter<InputTag>("segment4DLabel");
//sigma resolution on position
sigmaResX = pset.getParameter<double>("sigmaResX");
sigmaResY = pset.getParameter<double>("sigmaResY");
//sigma resolution on angle
sigmaResAlpha = pset.getParameter<double>("sigmaResAlpha");
sigmaResBeta = pset.getParameter<double>("sigmaResBeta");
doall = pset.getUntrackedParameter<bool>("doall", false);
local = pset.getUntrackedParameter<bool>("local", false);
// Create the root file
//theFile = new TFile(rootFileName.c_str(), "RECREATE");
//theFile->cd();
// ----------------------
// get hold of back-end interface
dbe_ = 0;
dbe_ = Service<DQMStore>().operator->();
if ( dbe_ ) {
if (debug) {
dbe_->setVerbose(1);
} else {
dbe_->setVerbose(0);
}
}
if ( dbe_ ) {
if ( debug ) dbe_->showDirStructure();
}
h4DHit= new HRes4DHit ("All",dbe_,doall,local);
h4DHit_W0= new HRes4DHit ("W0",dbe_,doall,local);
h4DHit_W1= new HRes4DHit ("W1",dbe_,doall,local);
h4DHit_W2= new HRes4DHit ("W2",dbe_,doall,local);
if(doall) {
hEff_All= new HEff4DHit ("All",dbe_);
hEff_W0= new HEff4DHit ("W0",dbe_);
hEff_W1= new HEff4DHit ("W1",dbe_);
hEff_W2= new HEff4DHit ("W2",dbe_);
}
if (local) {
// Plots with finer granularity, not to be included in DQM
TString name="W";
for (long w=0;w<=2;++w) {
for (long s=1;s<=4;++s){
// FIXME station 4 is not filled
TString nameWS=(name+w+"_St"+s);
h4DHitWS[w][s-1] = new HRes4DHit(nameWS.Data(),dbe_,doall,local);
hEffWS[w][s-1] = new HEff4DHit (nameWS.Data(),dbe_);
dbe_->setCurrentFolder("DT/4DSegments/");
hHitMult[w][s-1] = dbe_->book2D("4D_" +nameWS+ "_hNHits", "NHits",12,0,12, 6,0,6);
ht0[w][s-1] = dbe_->book2D("4D_" +nameWS+ "_ht0", "t0",200,-25,25,200,-25,25);
}
}
}
}
| DTSegment4DQuality::~DTSegment4DQuality | ( | ) | [virtual] |
| void DTSegment4DQuality::analyze | ( | const edm::Event & | event, |
| const edm::EventSetup & | eventSetup | ||
| ) | [virtual] |
Perform the real analysis.
Implements edm::EDAnalyzer.
Definition at line 142 of file DTSegment4DQuality.cc.
References abs, funct::cos(), gather_cfg::cout, debug, PV3DBase< T, PVType, FrameType >::eta(), HEff4DHit::Fill(), HRes4DHit::Fill(), DTHitQualityUtils::findMuSimSegment(), DTHitQualityUtils::findMuSimSegmentDirAndPos(), DTHitQualityUtils::findSegmentAlphaAndBeta(), edm::EventSetup::get(), interpolateCardsSimple::histo, edm::HandleBase::isValid(), DTRecSegment2D::localDirection(), DTRecSegment4D::localDirection(), DTRecSegment2D::localDirectionError(), DTRecSegment4D::localDirectionError(), DTRecSegment4D::localPosition(), DTRecSegment2D::localPosition(), DTRecSegment4D::localPositionError(), DTRecSegment2D::localPositionError(), DTHitQualityUtils::mapMuSimHitsPerWire(), DTHitQualityUtils::mapSimHitsPerWire(), PV3DBase< T, PVType, FrameType >::phi(), DTRecSegment2D::recHits(), trackerHits::simHits, mathSSE::sqrt(), DTSLRecSegment2D::superLayerId(), DTRecSegment2D::t0(), PV3DBase< T, PVType, FrameType >::theta(), GeomDet::toGlobal(), GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), PV3DBase< T, PVType, FrameType >::z(), and DTRecSegment4D::zSegment().
{
//theFile->cd();
// Get the DT Geometry
ESHandle<DTGeometry> dtGeom;
eventSetup.get<MuonGeometryRecord>().get(dtGeom);
// Get the SimHit collection from the event
edm::Handle<PSimHitContainer> simHits;
event.getByLabel(simHitLabel, simHits); //FIXME: second string to be removed
//Map simHits by chamber
map<DTChamberId, PSimHitContainer > simHitsPerCh;
for(PSimHitContainer::const_iterator simHit = simHits->begin();
simHit != simHits->end(); simHit++){
// Create the id of the chamber (the simHits in the DT known their wireId)
DTChamberId chamberId = (((DTWireId(simHit->detUnitId())).layerId()).superlayerId()).chamberId();
// Fill the map
simHitsPerCh[chamberId].push_back(*simHit);
}
// Get the 4D rechits from the event
Handle<DTRecSegment4DCollection> segment4Ds;
event.getByLabel(segment4DLabel, segment4Ds);
if(!segment4Ds.isValid()) {
if(debug) cout << "[DTSegment4DQuality]**Warning: no 4D Segments with label: " <<segment4DLabel
<< " in this event, skipping!" << endl;
return;
}
// Loop over all chambers containing a segment
DTRecSegment4DCollection::id_iterator chamberId;
for (chamberId = segment4Ds->id_begin();
chamberId != segment4Ds->id_end();
++chamberId){
if((*chamberId).station() == 4)
continue; //use DTSegment2DSLPhiQuality to analyze MB4 performaces
//------------------------- simHits ---------------------------//
//Get simHits of each chamber
PSimHitContainer simHits = simHitsPerCh[(*chamberId)];
// Map simhits per wire
map<DTWireId, PSimHitContainer > simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
map<DTWireId, const PSimHit*> muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
int nMuSimHit = muSimHitPerWire.size();
if(nMuSimHit == 0 || nMuSimHit == 1) {
if(debug && nMuSimHit == 1)
cout << "[DTSegment4DQuality] Only " << nMuSimHit << " mu SimHit in this chamber, skipping!" << endl;
continue; // If no or only one mu SimHit is found skip this chamber
}
if(debug)
cout << "=== Chamber " << (*chamberId) << " has " << nMuSimHit << " SimHits" << endl;
//Find outer and inner mu SimHit to build a segment
pair<const PSimHit*, const PSimHit*> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire);
//Find direction and position of the sim Segment in Chamber RF
pair<LocalVector, LocalPoint> dirAndPosSimSegm = DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit,
(*chamberId),&(*dtGeom));
LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
const DTChamber* chamber = dtGeom->chamber(*chamberId);
GlobalPoint simSegmGlobalPos = chamber->toGlobal(simSegmLocalPos);
GlobalVector simSegmGlobalDir = chamber->toGlobal(simSegmLocalDir);
//phi and theta angle of simulated segment in Chamber RF
float alphaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
float betaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).second;
//x,y position of simulated segment in Chamber RF
float xSimSeg = simSegmLocalPos.x();
float ySimSeg = simSegmLocalPos.y();
//Position (in eta,phi coordinates) in lobal RF
float etaSimSeg = simSegmGlobalPos.eta();
float phiSimSeg = simSegmGlobalPos.phi();
if(debug)
cout<<" Simulated segment: local direction "<<simSegmLocalDir<<endl
<<" local position "<<simSegmLocalPos<<endl
<<" alpha "<<alphaSimSeg<<endl
<<" beta "<<betaSimSeg<<endl;
//---------------------------- recHits --------------------------//
// Get the range of rechit for the corresponding chamberId
bool recHitFound = false;
DTRecSegment4DCollection::range range = segment4Ds->get(*chamberId);
int nsegm = distance(range.first, range.second);
if(debug)
cout << " Chamber: " << *chamberId << " has " << nsegm
<< " 4D segments" << endl;
if (nsegm!=0) {
// Find the best RecHit: look for the 4D RecHit with the phi angle closest
// to that of segment made of SimHits.
// RecHits must have delta alpha and delta position within 5 sigma of
// the residual distribution (we are looking for residuals of segments
// usefull to the track fit) for efficency purpose
const DTRecSegment4D* bestRecHit = 0;
bool bestRecHitFound = false;
double deltaAlpha = 99999;
// Loop over the recHits of this chamberId
for (DTRecSegment4DCollection::const_iterator segment4D = range.first;
segment4D!=range.second;
++segment4D){
if (local) {
const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment();
const DTSLRecSegment2D* zSeg = (*segment4D).zSegment();
float t0phi = -999;
float t0z = -999;
int nHitPhi=0;
int nHitZ=0;
if (phiSeg) {
t0phi = phiSeg->t0();
nHitPhi = phiSeg->recHits().size();
}
if (zSeg) {
t0z = zSeg->t0();
nHitZ = zSeg->recHits().size();
}
hHitMult[abs((*chamberId).wheel())][(*chamberId).station()-1]->Fill(nHitPhi,nHitZ);
ht0[abs((*chamberId).wheel())][(*chamberId).station()-1]->Fill(t0phi,t0z);
// Uncomment to skip segments w/o t0 computed
// if (!(t0phi>-998&&t0z>-998)) continue
}
// Check the dimension
if((*segment4D).dimension() != 4) {
if(debug)cout << "[DTSegment4DQuality]***Error: This is not 4D segment!!!" << endl;
continue;
}
// Segment Local Direction and position (in Chamber RF)
LocalVector recSegDirection = (*segment4D).localDirection();
//LocalPoint recSegPosition = (*segment4D).localPosition();
float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first;
if(debug)
cout << " RecSegment direction: " << recSegDirection << endl
<< " position : " << (*segment4D).localPosition() << endl
<< " alpha : " << recSegAlpha << endl
<< " beta : " << DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).second << endl;
if(fabs(recSegAlpha - alphaSimSeg) < deltaAlpha) {
deltaAlpha = fabs(recSegAlpha - alphaSimSeg);
bestRecHit = &(*segment4D);
bestRecHitFound = true;
}
} // End of Loop over all 4D RecHits
if(bestRecHitFound) {
// Best rechit direction and position in Chamber RF
LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
// Errors on x and y
LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
float alphaBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first;
float betaBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).second;
// Errors on alpha and beta
// Get position and direction using the rx projection (so in SL
// reference frame). Note that x (and y) are swapped wrt to Chamber
// frame
//if (bestRecHit->hasZed() ) {
const DTSLRecSegment2D * zedRecSeg = bestRecHit->zSegment();
LocalPoint bestRecHitLocalPosRZ = zedRecSeg->localPosition();
LocalVector bestRecHitLocalDirRZ = zedRecSeg->localDirection();
// Errors on x and y
LocalError bestRecHitLocalPosErrRZ = zedRecSeg->localPositionError();
LocalError bestRecHitLocalDirErrRZ = zedRecSeg->localDirectionError();
float alphaBestRHitRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDirRZ).first;
// Get SimSeg position and Direction in rZ SL frame
const DTSuperLayer* sl = dtGeom->superLayer(zedRecSeg->superLayerId());
LocalPoint simSegLocalPosRZTmp = sl->toLocal(simSegmGlobalPos);
LocalVector simSegLocalDirRZ = sl->toLocal(simSegmGlobalDir);
LocalPoint simSegLocalPosRZ =
simSegLocalPosRZTmp + simSegLocalDirRZ*(-simSegLocalPosRZTmp.z()/(cos(simSegLocalDirRZ.theta())));
float alphaSimSegRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegLocalDirRZ).first;
if (debug) cout <<
"RZ SL: recPos " << bestRecHitLocalPosRZ <<
"recDir " << bestRecHitLocalDirRZ <<
"recAlpha " << alphaBestRHitRZ << endl <<
"RZ SL: simPos " << simSegLocalPosRZ <<
"simDir " << simSegLocalDirRZ <<
"simAlpha " << alphaSimSegRZ << endl ;
//}
if(fabs(alphaBestRHit - alphaSimSeg) < 5*sigmaResAlpha &&
fabs(betaBestRHit - betaSimSeg) < 5*sigmaResBeta &&
fabs(bestRecHitLocalPos.x() - xSimSeg) < 5*sigmaResX &&
fabs(bestRecHitLocalPos.y() - ySimSeg) < 5*sigmaResY) {
recHitFound = true;
}
// Fill Residual histos
HRes4DHit *histo=0;
if((*chamberId).wheel() == 0)
histo = h4DHit_W0;
else if(abs((*chamberId).wheel()) == 1)
histo = h4DHit_W1;
else if(abs((*chamberId).wheel()) == 2)
histo = h4DHit_W2;
histo->Fill(alphaSimSeg,
alphaBestRHit,
betaSimSeg,
betaBestRHit,
xSimSeg,
bestRecHitLocalPos.x(),
ySimSeg,
bestRecHitLocalPos.y(),
etaSimSeg,
phiSimSeg,
bestRecHitLocalPosRZ.x(),
simSegLocalPosRZ.x(),
alphaBestRHitRZ,
alphaSimSegRZ,
sqrt(bestRecHitLocalDirErr.xx()),
sqrt(bestRecHitLocalDirErr.yy()),
sqrt(bestRecHitLocalPosErr.xx()),
sqrt(bestRecHitLocalPosErr.yy()),
sqrt(bestRecHitLocalDirErrRZ.xx()),
sqrt(bestRecHitLocalPosErrRZ.xx())
);
h4DHit->Fill(alphaSimSeg,
alphaBestRHit,
betaSimSeg,
betaBestRHit,
xSimSeg,
bestRecHitLocalPos.x(),
ySimSeg,
bestRecHitLocalPos.y(),
etaSimSeg,
phiSimSeg,
bestRecHitLocalPosRZ.x(),
simSegLocalPosRZ.x(),
alphaBestRHitRZ,
alphaSimSegRZ,
sqrt(bestRecHitLocalDirErr.xx()),
sqrt(bestRecHitLocalDirErr.yy()),
sqrt(bestRecHitLocalPosErr.xx()),
sqrt(bestRecHitLocalPosErr.yy()),
sqrt(bestRecHitLocalDirErrRZ.xx()),
sqrt(bestRecHitLocalPosErrRZ.xx())
);
if (local) h4DHitWS[abs((*chamberId).wheel())][(*chamberId).station()-1]->Fill(alphaSimSeg,
alphaBestRHit,
betaSimSeg,
betaBestRHit,
xSimSeg,
bestRecHitLocalPos.x(),
ySimSeg,
bestRecHitLocalPos.y(),
etaSimSeg,
phiSimSeg,
bestRecHitLocalPosRZ.x(),
simSegLocalPosRZ.x(),
alphaBestRHitRZ,
alphaSimSegRZ,
sqrt(bestRecHitLocalDirErr.xx()),
sqrt(bestRecHitLocalDirErr.yy()),
sqrt(bestRecHitLocalPosErr.xx()),
sqrt(bestRecHitLocalPosErr.yy()),
sqrt(bestRecHitLocalDirErrRZ.xx()),
sqrt(bestRecHitLocalPosErrRZ.xx())
);
} //end of if(bestRecHitFound)
} //end of if(nsegm!=0)
// Fill Efficiency plot
if(doall){
HEff4DHit *heff = 0;
if((*chamberId).wheel() == 0)
heff = hEff_W0;
else if(abs((*chamberId).wheel()) == 1)
heff = hEff_W1;
else if(abs((*chamberId).wheel()) == 2)
heff = hEff_W2;
heff->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound);
hEff_All->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound);
if (local) hEffWS[abs((*chamberId).wheel())][(*chamberId).station()-1]->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound);
}
} // End of loop over chambers
}
| void DTSegment4DQuality::endJob | ( | void | ) | [virtual] |
Reimplemented from edm::EDAnalyzer.
Definition at line 117 of file DTSegment4DQuality.cc.
{
// Write the histos to file
//theFile->cd();
//h4DHit->Write();
//h4DHit_W0->Write();
//h4DHit_W1->Write();
//h4DHit_W2->Write();
if(doall){
hEff_All->ComputeEfficiency();
hEff_W0->ComputeEfficiency();
hEff_W1->ComputeEfficiency();
hEff_W2->ComputeEfficiency();
}
//hEff_All->Write();
//hEff_W0->Write();
//hEff_W1->Write();
//hEff_W2->Write();
//if ( rootFileName.size() != 0 && dbe_ ) dbe_->save(rootFileName);
//theFile->Close();
}
| void DTSegment4DQuality::endLuminosityBlock | ( | edm::LuminosityBlock const & | lumiSeg, |
| edm::EventSetup const & | c | ||
| ) | [virtual] |
DQMStore* DTSegment4DQuality::dbe_ [private] |
Definition at line 95 of file DTSegment4DQuality.h.
bool DTSegment4DQuality::debug [private] |
Definition at line 67 of file DTSegment4DQuality.h.
bool DTSegment4DQuality::doall [private] |
Definition at line 96 of file DTSegment4DQuality.h.
HRes4DHit* DTSegment4DQuality::h4DHit [private] |
Definition at line 80 of file DTSegment4DQuality.h.
HRes4DHit* DTSegment4DQuality::h4DHit_W0 [private] |
Definition at line 81 of file DTSegment4DQuality.h.
HRes4DHit* DTSegment4DQuality::h4DHit_W1 [private] |
Definition at line 82 of file DTSegment4DQuality.h.
HRes4DHit* DTSegment4DQuality::h4DHit_W2 [private] |
Definition at line 83 of file DTSegment4DQuality.h.
HRes4DHit* DTSegment4DQuality::h4DHitWS[3][4] [private] |
Definition at line 84 of file DTSegment4DQuality.h.
HEff4DHit* DTSegment4DQuality::hEff_All [private] |
Definition at line 86 of file DTSegment4DQuality.h.
HEff4DHit* DTSegment4DQuality::hEff_W0 [private] |
Definition at line 87 of file DTSegment4DQuality.h.
HEff4DHit* DTSegment4DQuality::hEff_W1 [private] |
Definition at line 88 of file DTSegment4DQuality.h.
HEff4DHit* DTSegment4DQuality::hEff_W2 [private] |
Definition at line 89 of file DTSegment4DQuality.h.
HEff4DHit* DTSegment4DQuality::hEffWS[3][4] [private] |
Definition at line 90 of file DTSegment4DQuality.h.
MonitorElement* DTSegment4DQuality::hHitMult[3][4] [private] |
Definition at line 92 of file DTSegment4DQuality.h.
MonitorElement* DTSegment4DQuality::ht0[3][4] [private] |
Definition at line 93 of file DTSegment4DQuality.h.
bool DTSegment4DQuality::local [private] |
Definition at line 97 of file DTSegment4DQuality.h.
std::string DTSegment4DQuality::rootFileName [private] |
Definition at line 69 of file DTSegment4DQuality.h.
Definition at line 72 of file DTSegment4DQuality.h.
double DTSegment4DQuality::sigmaResAlpha [private] |
Definition at line 77 of file DTSegment4DQuality.h.
double DTSegment4DQuality::sigmaResBeta [private] |
Definition at line 78 of file DTSegment4DQuality.h.
double DTSegment4DQuality::sigmaResX [private] |
Definition at line 74 of file DTSegment4DQuality.h.
double DTSegment4DQuality::sigmaResY [private] |
Definition at line 75 of file DTSegment4DQuality.h.
edm::InputTag DTSegment4DQuality::simHitLabel [private] |
Definition at line 71 of file DTSegment4DQuality.h.
1.7.3