Reference Manual

Go to the documentation of this file.
00001 
00002 /*
00003  *  See header file for a description of this class.
00004  *
00005  *  $Date: 2010/07/20 02:58:24 $
00006  *  $Revision: 1.22 $
00007  *  \author G. Cerminara - INFN Torino
00008  */
00009 
00010 #include "DTResolutionAnalysisTask.h"
00011 
00012 // Framework
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/Framework/interface/EventSetup.h"
00015 #include "FWCore/ServiceRegistry/interface/Service.h"
00016 #include "FWCore/Framework/interface/LuminosityBlock.h"
00017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00018 
00019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00020 
00021 #include "DQMServices/Core/interface/DQMStore.h"
00022 #include "DQMServices/Core/interface/MonitorElement.h"
00023 
00024 //Geometry
00025 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00026 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00027 
00028 //RecHit
00029 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00030 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
00031 
00032 
00033 #include <iterator>
00034 
00035 using namespace edm;
00036 using namespace std;
00037 
00038 DTResolutionAnalysisTask::DTResolutionAnalysisTask(const ParameterSet& pset) {
00039 
00040   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask] Constructor called!" << endl;
00041 
00042   // the name of the 4D rec hits collection
00043   theRecHits4DLabel = pset.getParameter<string>("recHits4DLabel");
00044   
00045   prescaleFactor = pset.getUntrackedParameter<int>("diagnosticPrescale", 1);
00046   resetCycle = pset.getUntrackedParameter<int>("ResetCycle", -1);
00047   // top folder for the histograms in DQMStore
00048   topHistoFolder = pset.getUntrackedParameter<string>("topHistoFolder","DT/02-Segments");
00049 
00050 }
00051 
00052 
00053 DTResolutionAnalysisTask::~DTResolutionAnalysisTask(){
00054 
00055   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask] Destructor called!" << endl;
00056 
00057 }
00058 
00059 void DTResolutionAnalysisTask::beginRun(const Run& run, const EventSetup& setup) {
00060 
00061   // Get the DQM needed services
00062   theDbe = edm::Service<DQMStore>().operator->();
00063   //   theDbe->setCurrentFolder("DT/02-Segments");
00064 
00065   // Get the DT Geometry
00066   setup.get<MuonGeometryRecord>().get(dtGeom);
00067 
00068   // Book the histograms
00069   vector<DTChamber*> chambers = dtGeom->chambers();
00070   for(vector<DTChamber*>::const_iterator chamber = chambers.begin();
00071       chamber != chambers.end(); ++chamber) {  // Loop over all chambers
00072     DTChamberId dtChId = (*chamber)->id();
00073     for(int sl = 1; sl <= 3; ++sl) { // Loop over SLs
00074       if(dtChId.station() == 4 && sl == 2) continue;
00075       const  DTSuperLayerId dtSLId(dtChId,sl);
00076       bookHistos(dtSLId);
00077     }
00078   }
00079 
00080 }
00081 
00082 
00083 
00084 
00085 void DTResolutionAnalysisTask::beginLuminosityBlock(const LuminosityBlock& lumiSeg,
00086                                                     const EventSetup& context) {
00087 
00088   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionTask]: Begin of LS transition"<<endl;
00089   
00090   if(resetCycle != -1 && lumiSeg.id().luminosityBlock() % resetCycle == 0) {
00091     for(map<DTSuperLayerId, vector<MonitorElement*> > ::const_iterator histo = histosPerSL.begin();
00092         histo != histosPerSL.end();
00093         histo++) {
00094       int size = (*histo).second.size();
00095       for(int i=0; i<size; i++){
00096         (*histo).second[i]->Reset();
00097       }
00098     }
00099   }
00100 }
00101 
00102 
00103 void DTResolutionAnalysisTask::endJob(){
00104 
00105  edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask] endjob called!"<<endl;
00106 
00107 }
00108   
00109 
00110 
00111 void DTResolutionAnalysisTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
00112 
00113   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask] Analyze #Run: " << event.id().run()
00114          << " #Event: " << event.id().event() << endl;
00115 
00116   
00117   // Get the 4D segment collection from the event
00118   edm::Handle<DTRecSegment4DCollection> all4DSegments;
00119   event.getByLabel(theRecHits4DLabel, all4DSegments);
00120 
00121   // check the validity of the collection 
00122   if(!all4DSegments.isValid()) return;
00123 
00124   // Loop over all chambers containing a segment
00125   DTRecSegment4DCollection::id_iterator chamberId;
00126   for (chamberId = all4DSegments->id_begin();
00127        chamberId != all4DSegments->id_end();
00128        ++chamberId) {
00129     // Get the range for the corresponding ChamerId
00130     DTRecSegment4DCollection::range  range = all4DSegments->get(*chamberId);
00131     //     int nsegm = distance(range.first, range.second);
00132     //edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "   Chamber: " << *chamberId << " has " << nsegm
00133     //<< " 4D segments" << endl;
00134     // Get the chamber
00135     const DTChamber* chamber = dtGeom->chamber(*chamberId);
00136 
00137     // Loop over the rechits of this ChamerId
00138     for (DTRecSegment4DCollection::const_iterator segment4D = range.first;
00139          segment4D!=range.second;
00140          ++segment4D) {
00141       //edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "   == RecSegment dimension: " << (*segment4D).dimension() << endl;
00142       
00143       // If Statio != 4 skip RecHits with dimension != 4
00144       // For the Station 4 consider 2D RecHits
00145       if((*chamberId).station() != 4 && (*segment4D).dimension() != 4) {
00146         edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask]***Warning: RecSegment dimension is not 4 but "
00147                << (*segment4D).dimension() << "!" << endl;
00148         continue;
00149       } else if((*chamberId).station() == 4 && (*segment4D).dimension() != 2) {
00150         edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask]***Warning: RecSegment dimension is not 2 but "
00151                << (*segment4D).dimension() << "!" << endl;
00152         continue;
00153       }
00154 
00155 
00156       // Get all 1D RecHits at step 3 within the 4D segment
00157       vector<DTRecHit1D> recHits1D_S3;
00158     
00159 
00160       // Get 1D RecHits at Step 3 and select only events with
00161       // 8 hits in phi and 4 hits in theta (if any)
00162 
00163       if((*segment4D).hasPhi()) { // has phi component
00164         const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment();
00165         vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
00166 
00167         if(phiRecHits.size() != 8) {
00168           //edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask] Phi segments has: " << phiRecHits.size()
00169           //<< " hits" << endl; // FIXME: info output
00170           continue;
00171         }
00172         copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
00173       } else {
00174         //edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask] 4D segment has not phi component!" << endl;
00175       }
00176 
00177       if((*segment4D).hasZed()) {
00178         const DTSLRecSegment2D* zSeg = (*segment4D).zSegment();
00179         vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
00180         if(zRecHits.size() != 4) {
00181           //edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "[DTResolutionAnalysisTask] Theta segments has: " << zRecHits.size()
00182           //<< " hits, skipping" << endl; // FIXME: info output
00183           continue;
00184         }
00185         copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
00186       }
00187 
00188       // Loop over 1D RecHit inside 4D segment
00189       for(vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin();
00190           recHit1D != recHits1D_S3.end();
00191           recHit1D++) {
00192         const DTWireId wireId = (*recHit1D).wireId();
00193         
00194         // Get the layer and the wire position
00195         const DTLayer* layer = chamber->superLayer(wireId.superlayerId())->layer(wireId.layerId());
00196         float wireX = layer->specificTopology().wirePosition(wireId.wire());
00197 
00198         // Distance of the 1D rechit from the wire
00199         float distRecHitToWire = fabs(wireX - (*recHit1D).localPosition().x());
00200         
00201         // Extrapolate the segment to the z of the wire
00202         
00203         // Get wire position in chamber RF
00204         LocalPoint wirePosInLay(wireX,(*recHit1D).localPosition().y(),(*recHit1D).localPosition().z());
00205         GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
00206         LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
00207 
00208         // Segment position at Wire z in chamber local frame
00209         LocalPoint segPosAtZWire = (*segment4D).localPosition()
00210           + (*segment4D).localDirection()*wirePosInChamber.z()/cos((*segment4D).localDirection().theta());
00211         
00212         // Compute the distance of the segment from the wire
00213         int sl = wireId.superlayer();
00214   
00215         double distSegmToWire = -1;     
00216         if(sl == 1 || sl == 3) {
00217           // RPhi SL
00218           distSegmToWire = fabs(wirePosInChamber.x() - segPosAtZWire.x());
00219         } else if(sl == 2) {
00220           // RZ SL
00221           distSegmToWire = fabs(wirePosInChamber.y() - segPosAtZWire.y());
00222         }
00223 
00224         if(distSegmToWire > 2.1)
00225           edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "  Warning: dist segment-wire: " << distSegmToWire << endl;
00226 
00227         double residual = distRecHitToWire - distSegmToWire;
00228         // FIXME: Fill the histos
00229         fillHistos(wireId.superlayerId(), distSegmToWire, residual);
00230         
00231         //      edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "     Dist. segment extrapolation - wire (cm): " << distSegmToWire << endl;
00232         //edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "     Dist. RecHit - wire (cm): " << distRecHitToWire << endl;
00233         //edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "     Residual (cm): " << residual << endl;
00234         
00235                           
00236       }// End of loop over 1D RecHit inside 4D segment
00237     }// End of loop over the rechits of this ChamerId
00238   }
00239   // -----------------------------------------------------------------------------
00240 }
00241 
00242 
00243 // Book a set of histograms for a given SL
00244 void DTResolutionAnalysisTask::bookHistos(DTSuperLayerId slId) {
00245   
00246   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "   Booking histos for SL: " << slId << endl;
00247 
00248   // Compose the chamber name
00249   stringstream wheel; wheel << slId.wheel();    
00250   stringstream station; station << slId.station();      
00251   stringstream sector; sector << slId.sector(); 
00252   stringstream superLayer; superLayer << slId.superlayer();     
00253 
00254   
00255   string slHistoName =
00256     "_W" + wheel.str() +
00257     "_St" + station.str() +
00258     "_Sec" + sector.str() +
00259     "_SL" + superLayer.str();
00260   
00261   theDbe->setCurrentFolder(topHistoFolder + "/Wheel" + wheel.str() +
00262                            "/Sector" + sector.str() +
00263                            "/Station" + station.str());
00264   // Create the monitor elements
00265   vector<MonitorElement *> histos;
00266   // Note hte order matters
00267   histos.push_back(theDbe->book1D("hResDist"+slHistoName,
00268                                   "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
00269                                   200, -0.4, 0.4));
00270   //FIXME: 2D plot removed to reduce the # of ME
00271 //   histos.push_back(theDbe->book2D("hResDistVsDist"+slHistoName,
00272 //                                "Residuals on the distance (cm) from wire (rec_hit - segm_extr) vs distance  (cm)",
00273 //                                100, 0, 2.5, 200, -0.4, 0.4));
00274   histosPerSL[slId] = histos;
00275 }
00276 
00277 
00278 // Fill a set of histograms for a given SL 
00279 void DTResolutionAnalysisTask::fillHistos(DTSuperLayerId slId,
00280                                       float distExtr,
00281                                       float residual) {
00282   vector<MonitorElement *> histos =  histosPerSL[slId];                          
00283   histos[0]->Fill(residual);
00284   //FIXME: 2D plot removed to reduce the # of ME
00285   //   histos[1]->Fill(distExtr, residual); 
00286 
00287 }
00288 
00289 
00290
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/DQM/DTMonitorModule/src/DTResolutionAnalysisTask.cc
