Reference Manual

Go to the documentation of this file.
00001 
00002 /*
00003  *  See header file for a description of this class.
00004  *
00005  *  $Date: 2012/02/02 13:30:02 $
00006  *  $Revision: 1.4 $
00007  */
00008 
00009 #include "DTResidualCalibration.h"
00010 
00011 // Framework
00012 #include "FWCore/Framework/interface/Event.h"
00013 #include "FWCore/Framework/interface/EventSetup.h"
00014 #include "FWCore/Framework/interface/ESHandle.h"
00015 #include "FWCore/Utilities/interface/InputTag.h"
00016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00017 
00018 //Geometry
00019 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00020 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00021 
00022 //RecHit
00023 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00024 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
00025 
00026 #include "CommonTools/Utils/interface/TH1AddDirectorySentry.h"
00027 #include "CalibMuon/DTCalibration/interface/DTSegmentSelector.h"
00028 #include "CalibMuon/DTCalibration/interface/DTRecHitSegmentResidual.h"
00029 
00030 #include "TFile.h"
00031 #include "TH1F.h"
00032 #include "TH2F.h"
00033 
00034 #include <algorithm>
00035 
00036 DTResidualCalibration::DTResidualCalibration(const edm::ParameterSet& pset):
00037   select_(pset),
00038   segment4DLabel_(pset.getParameter<edm::InputTag>("segment4DLabel")),
00039   rootBaseDir_(pset.getUntrackedParameter<std::string>("rootBaseDir","DT/Residuals")),
00040   detailedAnalysis_(pset.getUntrackedParameter<bool>("detailedAnalysis",false)) {
00041 
00042   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Constructor called.";
00043 
00044   std::string rootFileName = pset.getUntrackedParameter<std::string>("rootFileName","residuals.root");
00045   rootFile_ = new TFile(rootFileName.c_str(), "RECREATE");
00046   rootFile_->cd();
00047 }
00048 
00049 DTResidualCalibration::~DTResidualCalibration() {
00050   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Destructor called.";
00051 }
00052 
00053 void DTResidualCalibration::beginJob() {
00054   TH1::SetDefaultSumw2(true);
00055 }
00056 
00057 void DTResidualCalibration::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
00058   
00059   // get the geometry
00060   edm::ESHandle<DTGeometry> dtGeomH; 
00061   setup.get<MuonGeometryRecord>().get(dtGeomH);
00062   dtGeom_ = dtGeomH.product();
00063 
00064   // Loop over all the chambers
00065   if(histoMapTH1F_.size() == 0) {        
00066      std::vector<DTChamber*>::const_iterator ch_it = dtGeom_->chambers().begin();        
00067      std::vector<DTChamber*>::const_iterator ch_end = dtGeom_->chambers().end();         
00068      for (; ch_it != ch_end; ++ch_it) {          
00069         std::vector<const DTSuperLayer*>::const_iterator sl_it = (*ch_it)->superLayers().begin();        
00070         std::vector<const DTSuperLayer*>::const_iterator sl_end = (*ch_it)->superLayers().end();         
00071         // Loop over the SLs     
00072         for(; sl_it != sl_end; ++sl_it) { 
00073            DTSuperLayerId slId = (*sl_it)->id();
00074            bookHistos(slId);
00075            if(detailedAnalysis_) {
00076               std::vector<const DTLayer*>::const_iterator layer_it = (*sl_it)->layers().begin();         
00077               std::vector<const DTLayer*>::const_iterator layer_end = (*sl_it)->layers().end();
00078               for(; layer_it != layer_end; ++layer_it) { 
00079                  DTLayerId layerId = (*layer_it)->id();
00080                  bookHistos(layerId);
00081               }
00082            }
00083         }
00084      }
00085   }
00086 }
00087 
00088 void DTResidualCalibration::analyze(const edm::Event& event, const edm::EventSetup& setup) {
00089   rootFile_->cd();
00090 
00091   // Get the 4D rechits from the event
00092   edm::Handle<DTRecSegment4DCollection> segment4Ds;
00093   event.getByLabel(segment4DLabel_, segment4Ds);
00094  
00095   // Loop over segments by chamber
00096   DTRecSegment4DCollection::id_iterator chamberIdIt;
00097   for(chamberIdIt = segment4Ds->id_begin(); chamberIdIt != segment4Ds->id_end(); ++chamberIdIt){
00098 
00099      const DTChamber* chamber = dtGeom_->chamber(*chamberIdIt);
00100 
00101      // Get the range for the corresponding ChamberId
00102      DTRecSegment4DCollection::range range = segment4Ds->get((*chamberIdIt));
00103      // Loop over the rechits of this DetUnit
00104      for(DTRecSegment4DCollection::const_iterator segment  = range.first;
00105                                                   segment != range.second; ++segment){
00106 
00107         LogTrace("Calibration") << "Segment local pos (in chamber RF): " << (*segment).localPosition()
00108                                 << "\nSegment global pos: " << chamber->toGlobal((*segment).localPosition());
00109 
00110         if( !select_(*segment, event, setup) ) continue;
00111 
00112         // Get all 1D RecHits at step 3 within the 4D segment
00113         std::vector<DTRecHit1D> recHits1D_S3;
00114   
00115         if( (*segment).hasPhi() ){
00116            const DTChamberRecSegment2D* phiSeg = (*segment).phiSegment();
00117            const std::vector<DTRecHit1D>& phiRecHits = phiSeg->specificRecHits();
00118            std::copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
00119         }
00120 
00121         if( (*segment).hasZed() ){
00122            const DTSLRecSegment2D* zSeg = (*segment).zSegment();
00123            const std::vector<DTRecHit1D>& zRecHits = zSeg->specificRecHits();
00124            std::copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
00125         }
00126 
00127         // Loop over 1D RecHit inside 4D segment
00128         for(std::vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin();
00129                                                     recHit1D != recHits1D_S3.end(); ++recHit1D) {
00130            const DTWireId wireId = recHit1D->wireId();
00131 
00132            float segmDistance = segmentToWireDistance(*recHit1D,*segment);
00133            if(segmDistance > 2.1) LogTrace("Calibration") << "WARNING: segment-wire distance: " << segmDistance;
00134            else                   LogTrace("Calibration") << "segment-wire distance: " << segmDistance;
00135 
00136            float residualOnDistance = DTRecHitSegmentResidual().compute(dtGeom_,*recHit1D,*segment);
00137            LogTrace("Calibration") << "Wire Id " << wireId << " residual on distance: " << residualOnDistance;
00138 
00139            fillHistos(wireId.superlayerId(), segmDistance, residualOnDistance);
00140            if(detailedAnalysis_) fillHistos(wireId.layerId(), segmDistance, residualOnDistance);
00141         }
00142      }
00143   }
00144 
00145 }
00146 
00147 float DTResidualCalibration::segmentToWireDistance(const DTRecHit1D& recHit1D, const DTRecSegment4D& segment){
00148 
00149   // Get the layer and the wire position
00150   const DTWireId wireId = recHit1D.wireId();
00151   const DTLayer* layer = dtGeom_->layer(wireId);
00152   float wireX = layer->specificTopology().wirePosition(wireId.wire());
00153       
00154   // Extrapolate the segment to the z of the wire
00155   // Get wire position in chamber RF
00156   // (y and z must be those of the hit to be coherent in the transf. of RF in case of rotations of the layer alignment)
00157   LocalPoint wirePosInLay(wireX,recHit1D.localPosition().y(),recHit1D.localPosition().z());
00158   GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
00159   const DTChamber* chamber = dtGeom_->chamber(wireId.layerId().chamberId());
00160   LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
00161       
00162   // Segment position at Wire z in chamber local frame
00163   LocalPoint segPosAtZWire = segment.localPosition()    + segment.localDirection()*wirePosInChamber.z()/cos(segment.localDirection().theta());
00164       
00165   // Compute the distance of the segment from the wire
00166   int sl = wireId.superlayer();
00167   float segmDistance = -1;
00168   if(sl == 1 || sl == 3) segmDistance = fabs(wirePosInChamber.x() - segPosAtZWire.x());
00169   else if(sl == 2)       segmDistance =  fabs(segPosAtZWire.y() - wirePosInChamber.y());
00170      
00171   return segmDistance;
00172 }
00173 
00174 void DTResidualCalibration::endJob(){
00175   
00176   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Writing histos to file.";
00177   rootFile_->cd();
00178   rootFile_->Write();
00179   rootFile_->Close();
00180 
00181   /*std::map<DTSuperLayerId, std::vector<TH1F*> >::const_iterator itSlHistos = histoMapTH1F_.begin();
00182   std::map<DTSuperLayerId, std::vector<TH1F*> >::const_iterator itSlHistos_end = histoMapTH1F_.end(); 
00183   for(; itSlHistos != itSlHistos_end; ++itSlHistos){
00184      std::vector<TH1F*>::const_iterator itHistTH1F = (*itSlHistos).second.begin();
00185      std::vector<TH1F*>::const_iterator itHistTH1F_end = (*itSlHistos).second.end();
00186      for(; itHistTH1F != itHistTH1F_end; ++itHistTH1F) (*itHistTH1F)->Write();
00187 
00188      std::vector<TH2F*>::const_iterator itHistTH2F = histoMapTH2F_[(*itSlHistos).first].begin();
00189      std::vector<TH2F*>::const_iterator itHistTH2F_end = histoMapTH2F_[(*itSlHistos).first].end();
00190      for(; itHistTH2F != itHistTH2F_end; ++itHistTH2F) (*itHistTH2F)->Write();
00191   }*/
00192 
00193 }
00194 
00195 void DTResidualCalibration::bookHistos(DTSuperLayerId slId) {
00196   TH1AddDirectorySentry addDir;
00197   rootFile_->cd();
00198 
00199   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Booking histos for SL: " << slId;
00200 
00201   // Compose the chamber name
00202   std::stringstream wheelStr; wheelStr << slId.wheel(); 
00203   std::stringstream stationStr; stationStr << slId.station();   
00204   std::stringstream sectorStr; sectorStr << slId.sector();      
00205   std::stringstream superLayerStr; superLayerStr << slId.superlayer();
00206   // Define the step
00207   int step = 3;
00208   std::stringstream stepStr; stepStr << step;
00209 
00210   std::string slHistoName =
00211     "_STEP" + stepStr.str() +
00212     "_W" + wheelStr.str() +
00213     "_St" + stationStr.str() +
00214     "_Sec" + sectorStr.str() +
00215     "_SL" + superLayerStr.str();
00216   
00217   edm::LogVerbatim("Calibration") << "Accessing " << rootBaseDir_;
00218   TDirectory* baseDir = rootFile_->GetDirectory(rootBaseDir_.c_str());
00219   if(!baseDir) baseDir = rootFile_->mkdir(rootBaseDir_.c_str());
00220   edm::LogVerbatim("Calibration") << "Accessing " << ("Wheel" + wheelStr.str());
00221   TDirectory* wheelDir = baseDir->GetDirectory(("Wheel" + wheelStr.str()).c_str());
00222   if(!wheelDir) wheelDir = baseDir->mkdir(("Wheel" + wheelStr.str()).c_str());
00223   edm::LogVerbatim("Calibration") << "Accessing " << ("Station" + stationStr.str());
00224   TDirectory* stationDir = wheelDir->GetDirectory(("Station" + stationStr.str()).c_str());
00225   if(!stationDir) stationDir = wheelDir->mkdir(("Station" + stationStr.str()).c_str());
00226   edm::LogVerbatim("Calibration") << "Accessing " << ("Sector" + sectorStr.str());
00227   TDirectory* sectorDir = stationDir->GetDirectory(("Sector" + sectorStr.str()).c_str());
00228   if(!sectorDir) sectorDir = stationDir->mkdir(("Sector" + sectorStr.str()).c_str()); 
00229 
00230   /*std::string dirName = rootBaseDir_ + "/Wheel" + wheelStr.str() +
00231                                        "/Station" + stationStr.str() +
00232                                        "/Sector" + sectorStr.str();
00233 
00234   TDirectory* dir = rootFile_->GetDirectory(dirName.c_str());
00235   if(!dir) dir = rootFile_->mkdir(dirName.c_str());
00236   dir->cd();*/
00237   sectorDir->cd();
00238   // Create the monitor elements
00239   std::vector<TH1F*> histosTH1F;
00240   histosTH1F.push_back(new TH1F(("hResDist"+slHistoName).c_str(),
00241                                  "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
00242                                  200, -0.4, 0.4));
00243   std::vector<TH2F*> histosTH2F;
00244   histosTH2F.push_back(new TH2F(("hResDistVsDist"+slHistoName).c_str(),
00245                                  "Residuals on the dist. (cm) from wire (rec_hit - segm_extr) vs dist. (cm)",
00246                                  100, 0, 2.5, 200, -0.4, 0.4));
00247   histoMapTH1F_[slId] = histosTH1F;
00248   histoMapTH2F_[slId] = histosTH2F;
00249 }
00250 
00251 void DTResidualCalibration::bookHistos(DTLayerId layerId) {
00252   TH1AddDirectorySentry addDir;
00253   rootFile_->cd();
00254 
00255   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Booking histos for layer: " << layerId;
00256 
00257   // Compose the chamber name
00258   std::stringstream wheelStr; wheelStr << layerId.wheel();
00259   std::stringstream stationStr; stationStr << layerId.station();
00260   std::stringstream sectorStr; sectorStr << layerId.sector();
00261   std::stringstream superLayerStr; superLayerStr << layerId.superlayer();
00262   std::stringstream layerStr; layerStr << layerId.layer();
00263   // Define the step
00264   int step = 3;
00265   std::stringstream stepStr; stepStr << step;
00266 
00267   std::string layerHistoName =
00268     "_STEP" + stepStr.str() +
00269     "_W" + wheelStr.str() +
00270     "_St" + stationStr.str() +
00271     "_Sec" + sectorStr.str() +
00272     "_SL" + superLayerStr.str() + 
00273     "_Layer" + layerStr.str();
00274   
00275   edm::LogVerbatim("Calibration") << "Accessing " << rootBaseDir_;
00276   TDirectory* baseDir = rootFile_->GetDirectory(rootBaseDir_.c_str());
00277   if(!baseDir) baseDir = rootFile_->mkdir(rootBaseDir_.c_str());
00278   edm::LogVerbatim("Calibration") << "Accessing " << ("Wheel" + wheelStr.str());
00279   TDirectory* wheelDir = baseDir->GetDirectory(("Wheel" + wheelStr.str()).c_str());
00280   if(!wheelDir) wheelDir = baseDir->mkdir(("Wheel" + wheelStr.str()).c_str());
00281   edm::LogVerbatim("Calibration") << "Accessing " << ("Station" + stationStr.str());
00282   TDirectory* stationDir = wheelDir->GetDirectory(("Station" + stationStr.str()).c_str());
00283   if(!stationDir) stationDir = wheelDir->mkdir(("Station" + stationStr.str()).c_str());
00284   edm::LogVerbatim("Calibration") << "Accessing " << ("Sector" + sectorStr.str());
00285   TDirectory* sectorDir = stationDir->GetDirectory(("Sector" + sectorStr.str()).c_str());
00286   if(!sectorDir) sectorDir = stationDir->mkdir(("Sector" + sectorStr.str()).c_str()); 
00287   edm::LogVerbatim("Calibration") << "Accessing " << ("SL" + superLayerStr.str());
00288   TDirectory* superLayerDir = sectorDir->GetDirectory(("SL" + superLayerStr.str()).c_str());
00289   if(!superLayerDir) superLayerDir = sectorDir->mkdir(("SL" + superLayerStr.str()).c_str()); 
00290 
00291   superLayerDir->cd();
00292   // Create histograms
00293   std::vector<TH1F*> histosTH1F;
00294   histosTH1F.push_back(new TH1F(("hResDist"+layerHistoName).c_str(),
00295                                  "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
00296                                  200, -0.4, 0.4));
00297   std::vector<TH2F*> histosTH2F;
00298   histosTH2F.push_back(new TH2F(("hResDistVsDist"+layerHistoName).c_str(),
00299                                  "Residuals on the dist. (cm) from wire (rec_hit - segm_extr) vs dist. (cm)",
00300                                  100, 0, 2.5, 200, -0.4, 0.4));
00301   histoMapPerLayerTH1F_[layerId] = histosTH1F;
00302   histoMapPerLayerTH2F_[layerId] = histosTH2F;
00303 }
00304 
00305 // Fill a set of histograms for a given SL 
00306 void DTResidualCalibration::fillHistos(DTSuperLayerId slId,
00307                                        float distance,
00308                                        float residualOnDistance) {
00309   std::vector<TH1F*> const& histosTH1F = histoMapTH1F_[slId];
00310   std::vector<TH2F*> const& histosTH2F = histoMapTH2F_[slId];                          
00311   histosTH1F[0]->Fill(residualOnDistance);
00312   histosTH2F[0]->Fill(distance, residualOnDistance);
00313 }
00314 
00315 // Fill a set of histograms for a given layer 
00316 void DTResidualCalibration::fillHistos(DTLayerId layerId,
00317                                        float distance,
00318                                        float residualOnDistance) {
00319   std::vector<TH1F*> const& histosTH1F = histoMapPerLayerTH1F_[layerId];
00320   std::vector<TH2F*> const& histosTH2F = histoMapPerLayerTH2F_[layerId];                          
00321   histosTH1F[0]->Fill(residualOnDistance);
00322   histosTH2F[0]->Fill(distance, residualOnDistance);
00323 }
00324
/data/refman/pasoursint/CMSSW_5_3_3/src/CalibMuon/DTCalibration/plugins/DTResidualCalibration.cc
