|
|
|
#include <DTVDriftCalibration.h>
Classes | |
| class | cellInfo |
Public Member Functions | |
| void | analyze (const edm::Event &event, const edm::EventSetup &eventSetup) |
| DTVDriftCalibration (const edm::ParameterSet &pset) | |
| Constructor. | |
| void | endJob () |
| virtual | ~DTVDriftCalibration () |
| Destructor. | |
Private Types | |
| typedef DTTMax::TMax | TMax |
| enum | TMaxGranularity { byChamber, bySL, byPartition } |
Private Attributes | |
| bool | checkNoisyChannels |
| bool | debug |
| std::string | digiLabel |
| bool | findVDriftAndT0 |
| h2DSegm * | h2DSegmRPhi |
| h2DSegm * | h2DSegmRZ |
| h4DSegm * | h4DSegmAllCh |
| std::string | theCalibChamber |
| edm::ParameterSet | theCalibFilePar |
| TFile * | theFile |
| DTMeanTimerFitter * | theFitter |
| TMaxGranularity | theGranularity |
| double | theMaxChi2 |
| double | theMaxPhiAngle |
| double | theMaxZAngle |
| std::string | theRecHits4DLabel |
| DTTTrigBaseSync * | theSync |
| std::string | theVDriftOutputFile |
| std::map< DTWireId, cellInfo * > | theWireIdAndCellMap |
No description available.
Definition at line 34 of file DTVDriftCalibration.h.
typedef DTTMax::TMax DTVDriftCalibration::TMax [private] |
Definition at line 53 of file DTVDriftCalibration.h.
enum DTVDriftCalibration::TMaxGranularity [private] |
Definition at line 86 of file DTVDriftCalibration.h.
{byChamber, bySL, byPartition};
| DTVDriftCalibration::DTVDriftCalibration | ( | const edm::ParameterSet & | pset | ) |
Constructor.
Definition at line 49 of file DTVDriftCalibration.cc.
References bookHistos(), ExpressReco_HICollisions_FallBack::checkNoisyChannels, gather_cfg::cout, debug, reco::get(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hChi2, interactiveExample::theFile, and ExpressReco_HICollisions_FallBack::theMaxChi2.
{
// the root file which will contain the histos
string rootFileName = pset.getUntrackedParameter<string>("rootFileName");
theFile = new TFile(rootFileName.c_str(), "RECREATE");
theFile->cd();
debug = pset.getUntrackedParameter<bool>("debug", "false");
theFitter = new DTMeanTimerFitter(theFile);
if(debug)
theFitter->setVerbosity(1);
hChi2 = new TH1F("hChi2","Chi squared tracks",100,0,100);
h2DSegmRPhi = new h2DSegm("SLRPhi");
h2DSegmRZ = new h2DSegm("SLRZ");
h4DSegmAllCh = new h4DSegm("AllCh");
bookHistos();
findVDriftAndT0 = pset.getUntrackedParameter<bool>("fitAndWrite", "false");
// Chamber/s to calibrate
theCalibChamber = pset.getUntrackedParameter<string>("calibChamber", "All");
// the name of the 4D rec hits collection
theRecHits4DLabel = pset.getUntrackedParameter<string>("recHits4DLabel");
// the txt file which will contain the calibrated constants
theVDriftOutputFile = pset.getUntrackedParameter<string>("vDriftFileName");
//get the switch to check the noisy channels
checkNoisyChannels = pset.getUntrackedParameter<bool>("checkNoisyChannels","false");
// Get the synchronizer
theSync = DTTTrigSyncFactory::get()->create(pset.getUntrackedParameter<string>("tTrigMode"),
pset.getUntrackedParameter<ParameterSet>("tTrigModeConfig"));
// get parameter set for DTCalibrationMap constructor
theCalibFilePar = pset.getUntrackedParameter<ParameterSet>("calibFileConfig");
// get maximum chi2 value
theMaxChi2 = pset.getParameter<double>("maxChi2");
// Maximum incidence angle for Phi SL
theMaxPhiAngle = pset.getParameter<double>("maxAnglePhi");
// Maximum incidence angle for Theta SL
theMaxZAngle = pset.getParameter<double>("maxAngleZ");
// the granularity to be used for tMax
string tMaxGranularity = pset.getUntrackedParameter<string>("tMaxGranularity","bySL");
// Initialize correctly the enum which specify the granularity for the calibration
if(tMaxGranularity == "bySL") {
theGranularity = bySL;
} else if(tMaxGranularity == "byChamber"){
theGranularity = byChamber;
} else if(tMaxGranularity== "byPartition") {
theGranularity = byPartition;
} else {
cout << "[DTVDriftCalibration]###Warning: Check parameter tMaxGranularity: "
<< tMaxGranularity << " options not available!" << endl;
}
if(debug)
cout << "[DTVDriftCalibration]Constructor called!" << endl;
}
| DTVDriftCalibration::~DTVDriftCalibration | ( | ) | [virtual] |
Destructor.
Definition at line 116 of file DTVDriftCalibration.cc.
References gather_cfg::cout, debug, and interactiveExample::theFile.
| void DTVDriftCalibration::analyze | ( | const edm::Event & | event, |
| const edm::EventSetup & | eventSetup | ||
| ) | [virtual] |
Implements edm::EDAnalyzer.
Definition at line 125 of file DTVDriftCalibration.cc.
References DTVDriftCalibration::cellInfo::add(), ExpressReco_HICollisions_FallBack::checkNoisyChannels, gather_cfg::cout, debug, DTChamberId, edm::EventSetup::get(), DTTMax::getTMax(), hChi2, DTRecSegment2D::localDirection(), DTRecSegment2D::localPosition(), AlCaRecoCosmics_cfg::name, Geom::pi(), pi, DTChamberId::sector(), DTRecSegment2D::specificRecHits(), DTChamberId::station(), crabStatusFromReport::statusMap, DTSuperLayerId::superLayer(), DTChamber::superLayer(), DTSLRecSegment2D::superLayerId(), DTLayerId::superlayerId(), interactiveExample::theFile, ExpressReco_HICollisions_FallBack::theMaxChi2, PV3DBase< T, PVType, FrameType >::theta(), GeomDet::toGlobal(), GeomDet::toLocal(), DTVDriftCalibration::cellInfo::update(), DTChamberId::wheel(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
{
cout << endl<<"--- [DTVDriftCalibration] Event analysed #Run: " << event.id().run()
<< " #Event: " << event.id().event() << endl;
theFile->cd();
DTChamberId chosenChamberId;
if(theCalibChamber != "All") {
stringstream linestr;
int selWheel, selStation, selSector;
linestr << theCalibChamber;
linestr >> selWheel >> selStation >> selSector;
chosenChamberId = DTChamberId(selWheel, selStation, selSector);
cout << "chosen chamber " << chosenChamberId << endl;
}
// Get the DT Geometry
ESHandle<DTGeometry> dtGeom;
eventSetup.get<MuonGeometryRecord>().get(dtGeom);
// Get the rechit collection from the event
Handle<DTRecSegment4DCollection> all4DSegments;
event.getByLabel(theRecHits4DLabel, all4DSegments);
// Get the map of noisy channels
ESHandle<DTStatusFlag> statusMap;
if(checkNoisyChannels) {
eventSetup.get<DTStatusFlagRcd>().get(statusMap);
}
// Set the event setup in the Synchronizer
theSync->setES(eventSetup);
// Loop over segments by chamber
DTRecSegment4DCollection::id_iterator chamberIdIt;
for (chamberIdIt = all4DSegments->id_begin();
chamberIdIt != all4DSegments->id_end();
++chamberIdIt){
// Get the chamber from the setup
const DTChamber* chamber = dtGeom->chamber(*chamberIdIt);
if(debug)
cout << "Chamber Id: " << *chamberIdIt << endl;
// Calibrate just the chosen chamber/s
if((theCalibChamber != "All") && ((*chamberIdIt) != chosenChamberId))
continue;
// Get the range for the corresponding ChamberId
DTRecSegment4DCollection::range range = all4DSegments->get((*chamberIdIt));
// Loop over the rechits of this DetUnit
for (DTRecSegment4DCollection::const_iterator segment = range.first;
segment!=range.second; ++segment){
if(debug) {
cout << "Segment local pos (in chamber RF): " << (*segment).localPosition() << endl;
cout << "Segment global pos: " << chamber->toGlobal((*segment).localPosition()) << endl;;
}
//get the segment chi2
double chiSquare = ((*segment).chi2()/(*segment).degreesOfFreedom());
// cut on the segment chi2
if(chiSquare > theMaxChi2) continue;
// get the Phi 2D segment and plot the angle in the chamber RF
if(!((*segment).phiSegment())){
cout<<"No phi segment"<<endl;
}
LocalPoint phiSeg2DPosInCham;
LocalVector phiSeg2DDirInCham;
bool segmNoisy = false;
map<DTSuperLayerId,vector<DTRecHit1D> > hitsBySLMap;
if((*segment).hasPhi()){
const DTChamberRecSegment2D* phiSeg = (*segment).phiSegment(); // phiSeg lives in the chamber RF
phiSeg2DPosInCham = phiSeg->localPosition();
phiSeg2DDirInCham = phiSeg->localDirection();
vector<DTRecHit1D> phiHits = phiSeg->specificRecHits();
for(vector<DTRecHit1D>::const_iterator hit = phiHits.begin();
hit != phiHits.end(); ++hit) {
DTWireId wireId = (*hit).wireId();
DTSuperLayerId slId = wireId.superlayerId();
hitsBySLMap[slId].push_back(*hit);
// Check for noisy channels to skip them
if(checkNoisyChannels) {
bool isNoisy = false;
bool isFEMasked = false;
bool isTDCMasked = false;
bool isTrigMask = false;
bool isDead = false;
bool isNohv = false;
statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
if(isNoisy) {
if(debug)
cout << "Wire: " << wireId << " is noisy, skipping!" << endl;
segmNoisy = true;
}
}
}
}
// get the Theta 2D segment and plot the angle in the chamber RF
LocalVector zSeg2DDirInCham;
LocalPoint zSeg2DPosInCham;
if((*segment).hasZed()) {
const DTSLRecSegment2D* zSeg = (*segment).zSegment(); // zSeg lives in the SL RF
const DTSuperLayer* sl = chamber->superLayer(zSeg->superLayerId());
zSeg2DPosInCham = chamber->toLocal(sl->toGlobal((*zSeg).localPosition()));
zSeg2DDirInCham = chamber->toLocal(sl->toGlobal((*zSeg).localDirection()));
hitsBySLMap[zSeg->superLayerId()] = zSeg->specificRecHits();
// Check for noisy channels to skip them
vector<DTRecHit1D> zHits = zSeg->specificRecHits();
for(vector<DTRecHit1D>::const_iterator hit = zHits.begin();
hit != zHits.end(); ++hit) {
DTWireId wireId = (*hit).wireId();
if(checkNoisyChannels) {
bool isNoisy = false;
bool isFEMasked = false;
bool isTDCMasked = false;
bool isTrigMask = false;
bool isDead = false;
bool isNohv = false;
statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
if(isNoisy) {
if(debug)
cout << "Wire: " << wireId << " is noisy, skipping!" << endl;
segmNoisy = true;
}
}
}
}
if (segmNoisy) continue;
LocalPoint segment4DLocalPos = (*segment).localPosition();
LocalVector segment4DLocalDir = (*segment).localDirection();
if(fabs(atan(segment4DLocalDir.y()/segment4DLocalDir.z())* 180./Geom::pi()) > theMaxZAngle) continue; // cut on the angle
if(fabs(atan(segment4DLocalDir.x()/segment4DLocalDir.z())* 180./Geom::pi()) > theMaxPhiAngle) continue; // cut on the angle
hChi2->Fill(chiSquare);
if((*segment).hasPhi())
h2DSegmRPhi->Fill(phiSeg2DPosInCham.x(), phiSeg2DDirInCham.x()/phiSeg2DDirInCham.z());
if((*segment).hasZed())
h2DSegmRZ->Fill(zSeg2DPosInCham.y(), zSeg2DDirInCham.y()/zSeg2DDirInCham.z());
if((*segment).hasZed() && (*segment).hasPhi())
h4DSegmAllCh->Fill(segment4DLocalPos.x(),
segment4DLocalPos.y(),
atan(segment4DLocalDir.x()/segment4DLocalDir.z())* 180./Geom::pi(),
atan(segment4DLocalDir.y()/segment4DLocalDir.z())* 180./Geom::pi(),
180 - segment4DLocalDir.theta()* 180./Geom::pi());
else if((*segment).hasPhi())
h4DSegmAllCh->Fill(segment4DLocalPos.x(),
atan(segment4DLocalDir.x()/segment4DLocalDir.z())* 180./Geom::pi());
else if((*segment).hasZed())
cout<<"4d segment with only Z"<<endl;
else{
cout<<"ERROR: 4D segment without Z and Phi. Aborting!"<<endl;
abort();
}
//loop over the segments
for(map<DTSuperLayerId,vector<DTRecHit1D> >::const_iterator slIdAndHits = hitsBySLMap.begin(); slIdAndHits != hitsBySLMap.end(); ++slIdAndHits) {
if (slIdAndHits->second.size() < 3) continue;
DTSuperLayerId slId = slIdAndHits->first;
// Create the DTTMax, that computes the 4 TMax
DTTMax slSeg(slIdAndHits->second, *(chamber->superLayer(slIdAndHits->first)),chamber->toGlobal((*segment).localDirection()), chamber->toGlobal((*segment).localPosition()), theSync);
if(theGranularity == bySL) {
vector<const TMax*> tMaxes = slSeg.getTMax(slId);
DTWireId wireId(slId, 0, 0);
theFile->cd();
cellInfo* cell = theWireIdAndCellMap[wireId];
if (cell==0) {
TString name = (((((TString) "TMax"+(long) slId.wheel()) +(long) slId.station())
+(long) slId.sector())+(long) slId.superLayer());
cell = new cellInfo(name);
theWireIdAndCellMap[wireId] = cell;
}
cell->add(tMaxes);
cell->update(); // FIXME to reset the counter to avoid triple counting, which actually is not used...
}
else {
cout << "[DTVDriftCalibration]###Warning: the chosen granularity is not implemented yet, only bySL available!" << endl;
}
// to be implemented: granularity different from bySL
// else if (theGranularity == byPartition) {
// // Use the custom granularity defined by partition();
// // in this case, add() should be called once for each Tmax of each layer
// // and triple counting should be avoided within add()
// vector<cellInfo*> cells;
// for (int i=1; i<=4; i++) {
// const DTTMax::InfoLayer* iLayer = slSeg.getInfoLayer(i);
// if(iLayer == 0) continue;
// cellInfo * cell = partition(iLayer->idWire);
// cells.push_back(cell);
// vector<const TMax*> tMaxes = slSeg.getTMax(iLayer->idWire);
// cell->add(tMaxes);
// }
// //reset the counter to avoid triple counting
// for (vector<cellInfo*>::const_iterator i = cells.begin();
// i!= cells.end(); i++) {
// (*i)->update();
// }
// }
}
}
}
}
| void DTVDriftCalibration::endJob | ( | void | ) | [virtual] |
Reimplemented from edm::EDAnalyzer.
Definition at line 342 of file DTVDriftCalibration.cc.
References DTCalibrationMap::addCell(), DTVelocityUnits::cm_per_ns, gather_cfg::cout, debug, DTCalibrationMap::getConsts(), DTVDriftCalibration::cellInfo::getHists(), DTCalibrationMap::getKey(), hChi2, DTWireId::layerId(), MultiGaussianStateTransform::N, record, DTChamberId::sector(), DTMtime::set(), DTChamberId::station(), DTSuperLayerId::superLayer(), interactiveExample::theFile, DTChamberId::wheel(), hTMaxCell::Write(), and DTCalibrationMap::writeConsts().
{
theFile->cd();
gROOT->GetList()->Write();
h2DSegmRPhi->Write();
h2DSegmRZ->Write();
h4DSegmAllCh->Write();
hChi2->Write();
// Instantiate a DTCalibrationMap object if you want to calculate the calibration constants
DTCalibrationMap calibValuesFile(theCalibFilePar);
// Create the object to be written to DB
DTMtime* mTime = new DTMtime();
// write the TMax histograms of each SL to the root file
if(theGranularity == bySL) {
for(map<DTWireId, cellInfo*>::const_iterator wireCell = theWireIdAndCellMap.begin();
wireCell != theWireIdAndCellMap.end(); wireCell++) {
cellInfo* cell= theWireIdAndCellMap[(*wireCell).first];
hTMaxCell* cellHists = cell->getHists();
theFile->cd();
cellHists->Write();
if(findVDriftAndT0) { // if TRUE: evaluate calibration constants from TMax hists filled in this job
// evaluate v_drift and sigma from the TMax histograms
DTWireId wireId = (*wireCell).first;
vector<float> newConstants;
TString N=(((((TString) "TMax"+(long) wireId.wheel()) +(long) wireId.station())
+(long) wireId.sector())+(long) wireId.superLayer());
vector<float> vDriftAndReso = theFitter->evaluateVDriftAndReso(N);
// Don't write the constants for the SL if the vdrift was not computed
if(vDriftAndReso.front() == -1)
continue;
const DTCalibrationMap::CalibConsts* oldConstants = calibValuesFile.getConsts(wireId);
if(oldConstants != 0) {
newConstants.push_back((*oldConstants)[0]);
newConstants.push_back((*oldConstants)[1]);
newConstants.push_back((*oldConstants)[2]);
} else {
newConstants.push_back(-1);
newConstants.push_back(-1);
newConstants.push_back(-1);
}
for(int ivd=0; ivd<=5;ivd++) {
// 0=vdrift, 1=reso, 2=(3deltat0-2deltat0), 3=(2deltat0-1deltat0),
// 4=(1deltat0-0deltat0), 5=deltat0 from hists with max entries,
newConstants.push_back(vDriftAndReso[ivd]);
}
calibValuesFile.addCell(calibValuesFile.getKey(wireId), newConstants);
// vdrift is cm/ns , resolution is cm
mTime->set((wireId.layerId()).superlayerId(),
vDriftAndReso[0],
vDriftAndReso[1],
DTVelocityUnits::cm_per_ns);
if(debug) {
cout << " SL: " << (wireId.layerId()).superlayerId()
<< " vDrift = " << vDriftAndReso[0]
<< " reso = " << vDriftAndReso[1] << endl;
}
}
}
}
// to be implemented: granularity different from bySL
// if(theGranularity == "byChamber") {
// const vector<DTChamber*> chambers = dMap.chambers();
// // Loop over all chambers
// for(vector<MuBarChamber*>::const_iterator chamber = chambers.begin();
// chamber != chambers.end(); chamber ++) {
// MuBarChamberId chamber_id = (*chamber)->id();
// MuBarDigiParameters::Key wire_id(chamber_id, 0, 0, 0);
// vector<float> newConstants;
// vector<float> vDriftAndReso = evaluateVDriftAndReso(wire_id, f);
// const CalibConsts* oldConstants = digiParams.getConsts(wire_id);
// if(oldConstants !=0) {
// newConstants = *oldConstants;
// newConstants.push_back(vDriftAndReso[0]);
// newConstants.push_back(vDriftAndReso[1]);
// newConstants.push_back(vDriftAndReso[2]);
// newConstants.push_back(vDriftAndReso[3]);
// } else {
// newConstants.push_back(-1);
// newConstants.push_back(-1);
// newConstants.push_back(vDriftAndReso[0]);
// newConstants.push_back(vDriftAndReso[1]);
// newConstants.push_back(vDriftAndReso[2]);
// newConstants.push_back(vDriftAndReso[3]);
// }
// digiParams.addCell(wire_id, newConstants);
// }
// }
//write values to a table
calibValuesFile.writeConsts(theVDriftOutputFile);
if(debug)
cout << "[DTVDriftCalibration]Writing vdrift object to DB!" << endl;
// Write the vdrift object to DB
string record = "DTMtimeRcd";
DTCalibDBUtils::writeToDB<DTMtime>(record, mTime);
}
bool DTVDriftCalibration::checkNoisyChannels [private] |
Definition at line 115 of file DTVDriftCalibration.h.
bool DTVDriftCalibration::debug [private] |
Definition at line 93 of file DTVDriftCalibration.h.
std::string DTVDriftCalibration::digiLabel [private] |
Definition at line 96 of file DTVDriftCalibration.h.
bool DTVDriftCalibration::findVDriftAndT0 [private] |
Definition at line 106 of file DTVDriftCalibration.h.
h2DSegm* DTVDriftCalibration::h2DSegmRPhi [private] |
Definition at line 79 of file DTVDriftCalibration.h.
h2DSegm* DTVDriftCalibration::h2DSegmRZ [private] |
Definition at line 78 of file DTVDriftCalibration.h.
h4DSegm* DTVDriftCalibration::h4DSegmAllCh [private] |
Definition at line 80 of file DTVDriftCalibration.h.
std::string DTVDriftCalibration::theCalibChamber [private] |
Definition at line 133 of file DTVDriftCalibration.h.
Definition at line 121 of file DTVDriftCalibration.h.
TFile* DTVDriftCalibration::theFile [private] |
Definition at line 99 of file DTVDriftCalibration.h.
DTMeanTimerFitter* DTVDriftCalibration::theFitter [private] |
Definition at line 102 of file DTVDriftCalibration.h.
Definition at line 87 of file DTVDriftCalibration.h.
double DTVDriftCalibration::theMaxChi2 [private] |
Definition at line 124 of file DTVDriftCalibration.h.
double DTVDriftCalibration::theMaxPhiAngle [private] |
Definition at line 127 of file DTVDriftCalibration.h.
double DTVDriftCalibration::theMaxZAngle [private] |
Definition at line 130 of file DTVDriftCalibration.h.
std::string DTVDriftCalibration::theRecHits4DLabel [private] |
Definition at line 90 of file DTVDriftCalibration.h.
DTTTrigBaseSync* DTVDriftCalibration::theSync [private] |
Definition at line 118 of file DTVDriftCalibration.h.
std::string DTVDriftCalibration::theVDriftOutputFile [private] |
Definition at line 109 of file DTVDriftCalibration.h.
std::map<DTWireId, cellInfo*> DTVDriftCalibration::theWireIdAndCellMap [private] |
Definition at line 112 of file DTVDriftCalibration.h.
1.7.3