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/06/28 07:59:01 $
00006  *  $Revision: 1.34 $
00007  *  \author G. Cerminara - INFN Torino
00008  *  revised by G. Mila - INFN Torino
00009  */
00010 
00011 #include "DTSegmentAnalysisTask.h"
00012 
00013 // Framework
00014 #include "FWCore/Framework/interface/Event.h"
00015 #include "FWCore/Framework/interface/EventSetup.h"
00016 #include "FWCore/Framework/interface/LuminosityBlock.h"
00017 #include "FWCore/Framework/interface/ESHandle.h"
00018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00019 #include "FWCore/ServiceRegistry/interface/Service.h"
00020 
00021 #include "DQMServices/Core/interface/DQMStore.h"
00022 #include "DQMServices/Core/interface/MonitorElement.h"
00023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00024 
00025 //Geometry
00026 #include "DataFormats/GeometryVector/interface/Pi.h"
00027 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00028 #include "Geometry/DTGeometry/interface/DTTopology.h"
00029 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00030 //RecHit
00031 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00032 
00033 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
00034 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
00035 
00036 #include "DQM/DTMonitorModule/interface/DTTimeEvolutionHisto.h"
00037 
00038 #include <iterator>
00039 #include <TMath.h>
00040 
00041 using namespace edm;
00042 using namespace std;
00043 
00044 DTSegmentAnalysisTask::DTSegmentAnalysisTask(const edm::ParameterSet& pset) : nevents(0) , nEventsInLS(0), hNevtPerLS(0) {
00045 
00046   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "[DTSegmentAnalysisTask] Constructor called!";
00047 
00048   // switch for detailed analysis
00049   detailedAnalysis = pset.getUntrackedParameter<bool>("detailedAnalysis",false);
00050   // the name of the 4D rec hits collection
00051   theRecHits4DLabel = pset.getParameter<string>("recHits4DLabel");
00052   // Get the map of noisy channels
00053   checkNoisyChannels = pset.getUntrackedParameter<bool>("checkNoisyChannels",false);
00054   // # of bins in the time histos
00055   nTimeBins = pset.getUntrackedParameter<int>("nTimeBins",100);
00056   // # of LS per bin in the time histos
00057   nLSTimeBin = pset.getUntrackedParameter<int>("nLSTimeBin",2);
00058   // switch on/off sliding bins in time histos
00059   slideTimeBins = pset.getUntrackedParameter<bool>("slideTimeBins",true);
00060   phiSegmCut = pset.getUntrackedParameter<double>("phiSegmCut",30.);
00061   nhitsCut = pset.getUntrackedParameter<int>("nhitsCut",12);
00062 
00063   // Get the DQM needed services
00064   theDbe = edm::Service<DQMStore>().operator->();
00065 
00066   // top folder for the histograms in DQMStore
00067   topHistoFolder = pset.getUntrackedParameter<string>("topHistoFolder","DT/02-Segments");
00068   // hlt DQM mode
00069   hltDQMMode = pset.getUntrackedParameter<bool>("hltDQMMode",false);
00070 
00071 }
00072 
00073 
00074 DTSegmentAnalysisTask::~DTSegmentAnalysisTask(){
00075   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "[DTSegmentAnalysisTask] Destructor called!";
00076 }
00077 
00078 
00079 void DTSegmentAnalysisTask::beginRun(const Run& run, const edm::EventSetup& context){ 
00080 
00081   // Get the DT Geometry
00082   context.get<MuonGeometryRecord>().get(dtGeom);
00083 
00084   if (!hltDQMMode) {
00085     theDbe->setCurrentFolder("DT/EventInfo/Counters");
00086     nEventMonitor = theDbe->bookFloat("nProcessedEventsSegment");
00087   }
00088 
00089   // loop over all the DT chambers & book the histos
00090   vector<DTChamber*> chambers = dtGeom->chambers();
00091   vector<DTChamber*>::const_iterator ch_it = chambers.begin();
00092   vector<DTChamber*>::const_iterator ch_end = chambers.end();
00093   for (; ch_it != ch_end; ++ch_it) {
00094     bookHistos((*ch_it)->id());
00095   }
00096 
00097   // book sector time-evolution histos
00098   int modeTimeHisto = 0;
00099   if(!slideTimeBins) modeTimeHisto = 1;
00100   for(int wheel = -2; wheel != 3; ++wheel) { // loop over wheels
00101     for(int sector = 1; sector <= 12; ++sector) { // loop over sectors
00102       stringstream wheelstr; wheelstr << wheel; 
00103       stringstream sectorstr; sectorstr << sector;
00104       string sectorHistoName = "NSegmPerEvent_W" + wheelstr.str()
00105         + "_Sec" + sectorstr.str();
00106       string sectorHistoTitle = "# segm. W" + wheelstr.str() + " Sect." + sectorstr.str();
00107 
00108       theDbe->setCurrentFolder(topHistoFolder + "/Wheel" + wheelstr.str() +
00109           "/Sector" + sectorstr.str());
00110       histoTimeEvol[wheel][sector] = new DTTimeEvolutionHisto(&(*theDbe),sectorHistoName,sectorHistoTitle,
00111           nTimeBins,nLSTimeBin,slideTimeBins,modeTimeHisto);
00112     }
00113   }
00114 
00115   if(hltDQMMode) theDbe->setCurrentFolder(topHistoFolder);
00116   else theDbe->setCurrentFolder("DT/EventInfo/");
00117 
00118   hNevtPerLS = new DTTimeEvolutionHisto(&(*theDbe),"NevtPerLS","# evt.",nTimeBins,nLSTimeBin,slideTimeBins,2);
00119 
00120 }
00121 
00122 
00123 void DTSegmentAnalysisTask::endJob() {
00124 
00125   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") <<"[DTSegmentAnalysisTask] endjob called!";
00126 
00127   delete hNevtPerLS;
00128   //theDbe->save("testMonitoring.root");
00129 
00130 }
00131 
00132 
00133 
00134 void DTSegmentAnalysisTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
00135 
00136   nevents++;
00137   nEventMonitor->Fill(nevents);
00138 
00139   nEventsInLS++;
00140   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "[DTSegmentAnalysisTask] Analyze #Run: " << event.id().run()
00141     << " #Event: " << event.id().event();
00142   if(!(event.id().event()%1000))
00143     edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "[DTSegmentAnalysisTask] Analyze #Run: " << event.id().run()
00144       << " #Event: " << event.id().event();
00145 
00146   ESHandle<DTStatusFlag> statusMap;
00147   if(checkNoisyChannels) {
00148     setup.get<DTStatusFlagRcd>().get(statusMap);
00149   } 
00150 
00151 
00152   // -- 4D segment analysis  -----------------------------------------------------
00153 
00154   // Get the 4D segment collection from the event
00155   edm::Handle<DTRecSegment4DCollection> all4DSegments;
00156   event.getByLabel(theRecHits4DLabel, all4DSegments);
00157 
00158   if(!all4DSegments.isValid()) return;
00159 
00160   // Loop over all chambers containing a segment
00161   DTRecSegment4DCollection::id_iterator chamberId;
00162   for (chamberId = all4DSegments->id_begin();
00163       chamberId != all4DSegments->id_end();
00164       ++chamberId){
00165     // Get the range for the corresponding ChamerId
00166     DTRecSegment4DCollection::range  range = all4DSegments->get(*chamberId);
00167 
00168     edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "   Chamber: " << *chamberId << " has " << distance(range.first, range.second)
00169       << " 4D segments";
00170 
00171     // Loop over the rechits of this ChamerId
00172     for (DTRecSegment4DCollection::const_iterator segment4D = range.first;
00173         segment4D!=range.second;
00174         ++segment4D){
00175 
00176       //FOR NOISY CHANNELS////////////////////////////////
00177       bool segmNoisy = false;
00178       if(checkNoisyChannels) {
00179 
00180         if((*segment4D).hasPhi()){
00181           const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment();
00182           vector<DTRecHit1D> phiHits = phiSeg->specificRecHits();
00183           map<DTSuperLayerId,vector<DTRecHit1D> > hitsBySLMap; 
00184           for(vector<DTRecHit1D>::const_iterator hit = phiHits.begin();
00185               hit != phiHits.end(); ++hit) {
00186             DTWireId wireId = (*hit).wireId();
00187 
00188             // Check for noisy channels to skip them
00189             bool isNoisy = false;
00190             bool isFEMasked = false;
00191             bool isTDCMasked = false;
00192             bool isTrigMask = false;
00193             bool isDead = false;
00194             bool isNohv = false;
00195             statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
00196             if(isNoisy) {
00197               edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "Wire: " << wireId << " is noisy, skipping!";
00198               segmNoisy = true;
00199             }      
00200           }
00201         }
00202 
00203         if((*segment4D).hasZed()) {
00204           const DTSLRecSegment2D* zSeg = (*segment4D).zSegment();  // zSeg lives in the SL RF
00205           // Check for noisy channels to skip them
00206           vector<DTRecHit1D> zHits = zSeg->specificRecHits();
00207           for(vector<DTRecHit1D>::const_iterator hit = zHits.begin();
00208               hit != zHits.end(); ++hit) {
00209             DTWireId wireId = (*hit).wireId();
00210             bool isNoisy = false;
00211             bool isFEMasked = false;
00212             bool isTDCMasked = false;
00213             bool isTrigMask = false;
00214             bool isDead = false;
00215             bool isNohv = false;
00216             statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
00217             if(isNoisy) {
00218               edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "Wire: " << wireId << " is noisy, skipping!";
00219               segmNoisy = true;
00220             }     
00221           }
00222         } 
00223 
00224       } // end of switch on noisy channels
00225       if (segmNoisy) {
00226         edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask")<<"skipping the segment: it contains noisy cells";
00227         continue;
00228       }
00229       //END FOR NOISY CHANNELS////////////////////////////////
00230 
00231       int nHits=0;
00232       if((*segment4D).hasPhi())
00233         nHits = (((*segment4D).phiSegment())->specificRecHits()).size();
00234       if((*segment4D).hasZed()) 
00235         nHits = nHits + ((((*segment4D).zSegment())->specificRecHits()).size());
00236 
00237       double anglePhiSegm(0.);
00238       if( (*segment4D).hasPhi() ) {
00239         double xdir = (*segment4D).phiSegment()->localDirection().x();      
00240         double zdir = (*segment4D).phiSegment()->localDirection().z();      
00241 
00242         anglePhiSegm = atan(xdir/zdir)*180./TMath::Pi();
00243       }
00244       if( fabs(anglePhiSegm) > phiSegmCut ) continue;
00245       // If the segment is in Wh+-2/SecX/MB1, get the DT chambers just above and check if there is a segment
00246       // to validate the segment present in MB1
00247       if( fabs((*chamberId).wheel()) == 2 && (*chamberId).station() == 1 ) {
00248 
00249         bool segmOk=false;
00250         int mb(2);
00251         while( mb < 4 ) { 
00252           DTChamberId checkMB((*chamberId).wheel(),mb,(*chamberId).sector());
00253           DTRecSegment4DCollection::range  ckrange = all4DSegments->get(checkMB);
00254 
00255           for (DTRecSegment4DCollection::const_iterator cksegment4D = ckrange.first;
00256               cksegment4D!=ckrange.second;
00257               ++cksegment4D){
00258 
00259             int nHits=0;
00260             if((*cksegment4D).hasPhi())
00261               nHits = (((*cksegment4D).phiSegment())->specificRecHits()).size();
00262             if((*cksegment4D).hasZed()) 
00263               nHits = nHits + ((((*cksegment4D).zSegment())->specificRecHits()).size());
00264 
00265             if( nHits >= nhitsCut ) segmOk=true;
00266 
00267           }
00268           mb++;
00269         }
00270 
00271         if( !segmOk ) continue;
00272 
00273       }
00274       fillHistos(*chamberId,
00275           nHits,
00276           (*segment4D).chi2()/(*segment4D).degreesOfFreedom());
00277     }
00278   }
00279 
00280   // -----------------------------------------------------------------------------
00281 }
00282 
00283 
00284 // Book a set of histograms for a give chamber
00285 void DTSegmentAnalysisTask::bookHistos(DTChamberId chamberId) {
00286 
00287   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "   Booking histos for chamber: " << chamberId;
00288 
00289 
00290   // Compose the chamber name
00291   stringstream wheel; wheel << chamberId.wheel();       
00292   stringstream station; station << chamberId.station(); 
00293   stringstream sector; sector << chamberId.sector();
00294 
00295   string chamberHistoName =
00296     "_W" + wheel.str() +
00297     "_St" + station.str() +
00298     "_Sec" + sector.str();
00299 
00300 
00301   for(int wh=-2; wh<=2; wh++){
00302     stringstream wheel; wheel << wh;
00303     theDbe->setCurrentFolder(topHistoFolder + "/Wheel" + wheel.str());
00304     string histoName =  "numberOfSegments_W" + wheel.str();
00305     if (theDbe->get(topHistoFolder + "/Wheel" + wheel.str() + "/" + histoName))
00306       continue;
00307     summaryHistos[wh] = theDbe->book2D(histoName.c_str(),histoName.c_str(),12,1,13,4,1,5);
00308     summaryHistos[wh]->setAxisTitle("Sector",1);
00309     summaryHistos[wh]->setBinLabel(1,"1",1);
00310     summaryHistos[wh]->setBinLabel(2,"2",1);
00311     summaryHistos[wh]->setBinLabel(3,"3",1);
00312     summaryHistos[wh]->setBinLabel(4,"4",1);
00313     summaryHistos[wh]->setBinLabel(5,"5",1);
00314     summaryHistos[wh]->setBinLabel(6,"6",1);
00315     summaryHistos[wh]->setBinLabel(7,"7",1);
00316     summaryHistos[wh]->setBinLabel(8,"8",1);
00317     summaryHistos[wh]->setBinLabel(9,"9",1);
00318     summaryHistos[wh]->setBinLabel(10,"10",1);
00319     summaryHistos[wh]->setBinLabel(11,"11",1);
00320     summaryHistos[wh]->setBinLabel(12,"12",1);
00321     summaryHistos[wh]->setBinLabel(1,"MB1",2);
00322     summaryHistos[wh]->setBinLabel(2,"MB2",2);
00323     summaryHistos[wh]->setBinLabel(3,"MB3",2);
00324     summaryHistos[wh]->setBinLabel(4,"MB4",2);  
00325   }
00326 
00327 
00328   theDbe->setCurrentFolder(topHistoFolder + "/Wheel" + wheel.str() +
00329       "/Sector" + sector.str() +
00330       "/Station" + station.str());
00331 
00332   // Create the monitor elements
00333   vector<MonitorElement *> histos;
00334   histos.push_back(theDbe->book1D("h4DSegmNHits"+chamberHistoName,
00335         "# of hits per segment",
00336         16, 0.5, 16.5));
00337   if(detailedAnalysis){
00338     histos.push_back(theDbe->book1D("h4DChi2"+chamberHistoName,
00339           "4D Segment reduced Chi2",
00340           20, 0, 20));
00341   }
00342   histosPerCh[chamberId] = histos;
00343 }
00344 
00345 
00346 // Fill a set of histograms for a give chamber 
00347 void DTSegmentAnalysisTask::fillHistos(DTChamberId chamberId,
00348     int nHits,
00349     float chi2) {
00350   int sector = chamberId.sector();
00351   if(chamberId.sector()==13) {
00352     sector = 4;
00353   } else if(chamberId.sector()==14) {
00354     sector = 10;
00355   }
00356 
00357   summaryHistos[chamberId.wheel()]->Fill(sector,chamberId.station());
00358   histoTimeEvol[chamberId.wheel()][sector]->accumulateValueTimeSlot(1);
00359 
00360   vector<MonitorElement *> histos =  histosPerCh[chamberId];                          
00361   histos[0]->Fill(nHits);
00362   if(detailedAnalysis){
00363     histos[1]->Fill(chi2);
00364   }
00365 
00366 }
00367 
00368 
00369 void DTSegmentAnalysisTask::endLuminosityBlock(LuminosityBlock const& lumiSeg, EventSetup const& eSetup) {
00370 
00371   hNevtPerLS->updateTimeSlot(lumiSeg.luminosityBlock(), nEventsInLS);
00372   // book sector time-evolution histos
00373   for(int wheel = -2; wheel != 3; ++wheel) {
00374     for(int sector = 1; sector <= 12; ++sector) {
00375       histoTimeEvol[wheel][sector]->updateTimeSlot(lumiSeg.luminosityBlock(), nEventsInLS);
00376     }
00377   }
00378 }
00379 
00380 
00381 void DTSegmentAnalysisTask::beginLuminosityBlock(LuminosityBlock const& lumiSeg, EventSetup const& eSetup) {
00382   nEventsInLS = 0;
00383 }
00384 
00385
/data/refman/pasoursint/CMSSW_5_2_9/src/DQM/DTMonitorModule/src/DTSegmentAnalysisTask.cc
