CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/DQM/DTMonitorModule/src/DTChamberEfficiencyTask.cc

Go to the documentation of this file.
00001 
00002 
00003 /*
00004  *  See header file for a description of this class.
00005  *
00006  *  $Date: 2011/06/14 08:53:04 $
00007  *  $Revision: 1.15 $
00008  *  \author G. Mila - INFN Torino
00009  */
00010 
00011 
00012 
00013 #include "DTChamberEfficiencyTask.h"
00014 
00015 //Framework
00016 #include "FWCore/Framework/interface/Event.h"
00017 #include "FWCore/Framework/interface/EventSetup.h"
00018 #include "FWCore/ServiceRegistry/interface/Service.h"
00019 
00020 #include "DQMServices/Core/interface/DQMStore.h"
00021 #include "DQMServices/Core/interface/MonitorElement.h"
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023 
00024 //Geometry
00025 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00026 
00027 
00028 #include <iterator>
00029 #include <iostream>
00030 #include <cmath>
00031 using namespace edm;
00032 using namespace std;
00033 
00034 
00035 
00036 DTChamberEfficiencyTask::DTChamberEfficiencyTask(const ParameterSet& pset) {
00037 
00038   debug = pset.getUntrackedParameter<bool>("debug",false);
00039 
00040   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "[DTChamberEfficiencyTask] Constructor called!";
00041 
00042   // Get the DQM needed services
00043   theDbe = edm::Service<DQMStore>().operator->();
00044 
00045   theDbe->setCurrentFolder("DT/DTChamberEfficiencyTask");
00046 
00047   parameters = pset;
00048 
00049 }
00050 
00051 
00052 DTChamberEfficiencyTask::~DTChamberEfficiencyTask(){
00053   
00054   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "[DTChamberEfficiencyTask] Destructor called!";
00055 }  
00056 
00057 
00058 void DTChamberEfficiencyTask::beginJob(){
00059 
00060   // the name of the 4D rec hits collection
00061   theRecHits4DLabel = parameters.getParameter<string>("recHits4DLabel");
00062 
00063   // parameters to use for the segment quality check
00064   theMinHitsSegment = static_cast<unsigned int>(parameters.getParameter<int>("minHitsSegment"));
00065   theMinChi2NormSegment = parameters.getParameter<double>("minChi2NormSegment");
00066   // parameter to use for the exstrapolated segment check
00067   theMinCloseDist = parameters.getParameter<double>("minCloseDist");
00068 
00069   // the running modality
00070   onlineMonitor = parameters.getUntrackedParameter<bool>("onlineMonitor");
00071 
00072   // the analysis mode
00073   detailedAnalysis = parameters.getUntrackedParameter<bool>("detailedAnalysis");
00074 
00075 }
00076 
00077 
00078 void DTChamberEfficiencyTask::beginLuminosityBlock(LuminosityBlock const& lumiSeg, EventSetup const& context) {
00079 
00080   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTChamberEfficiencyTask")<<"[DTChamberEfficiencyTask]: Begin of LS transition";
00081   
00082   if(lumiSeg.id().luminosityBlock()%parameters.getUntrackedParameter<int>("ResetCycle", 3) == 0 && onlineMonitor) {
00083     for(map<DTChamberId, vector<MonitorElement*> > ::const_iterator histo = histosPerCh.begin();
00084         histo != histosPerCh.end();
00085         histo++) {
00086       int size = (*histo).second.size();
00087       for(int i=0; i<size; i++){
00088         (*histo).second[i]->Reset();
00089       }
00090     }
00091   }
00092   
00093 }
00094 
00095 
00096 void DTChamberEfficiencyTask::beginRun(const edm::Run& run, const edm::EventSetup& setup){
00097   
00098   // Get the DT Geometry
00099   setup.get<MuonGeometryRecord>().get(dtGeom);
00100 
00101   // Loop over all the chambers
00102   vector<DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();
00103   vector<DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();
00104   for (; ch_it != ch_end; ++ch_it) {
00105     // histo booking
00106     bookHistos((*ch_it)->id());
00107   }
00108 
00109 }
00110 
00111 
00112 
00113 void DTChamberEfficiencyTask::endJob(){
00114 
00115     edm::LogVerbatim ("DTDQM|DTMonitorModule|DTChamberEfficiencyTask")<<"[DTChamberEfficiencyTask] endjob called!";
00116 }
00117   
00118 
00119 
00120 // Book a set of histograms for a given Layer
00121 void DTChamberEfficiencyTask::bookHistos(DTChamberId chId) {
00122   
00123   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "   Booking histos for CH : " << chId;
00124   
00125   // Compose the chamber name
00126   stringstream wheel; wheel << chId.wheel();    
00127   stringstream station; station << chId.station();      
00128   stringstream sector; sector << chId.sector(); 
00129   
00130   string HistoName =
00131     "_W" + wheel.str() +
00132     "_St" + station.str() +
00133     "_Sec" + sector.str();
00134   
00135   theDbe->setCurrentFolder("DT/01-DTChamberEfficiency/Task/Wheel" + wheel.str() +
00136                            "/Sector" + sector.str() +
00137                            "/Station" + station.str());
00138                 
00139   // Create the monitor elements
00140   vector<MonitorElement *> histos;
00141 
00142   //efficiency selection cuts
00143   // a- number of segments of the top chamber > 0 && number of segments of the bottom chamber > 0
00144   // b- number of segments of the middle chamber > 0
00145   // c- check of the top and bottom segment quality 
00146   // d- check if interpolation falls inside the middle chamber
00147   // e- check of the middle segment quality 
00148   // f- check if the distance between the reconstructed and the exstrapolated segments is ok
00149 
00150 
00151   // histo for efficiency with cuts a-/c-/d-
00152   histos.push_back(theDbe->book2D("hEffGoodSegVsPosDen"+HistoName,"Eff vs local position (good) ",25,-250.,250., 25,-250.,250.));
00153   // histo for efficiency with cuts a-/b-/c-/d-/e-/f-
00154   histos.push_back(theDbe->book2D("hEffGoodCloseSegVsPosNum"+HistoName, "Eff vs local position (good and close segs) ", 25,-250.,250., 25,-250.,250.));
00155   if(detailedAnalysis){
00156     histos.push_back(theDbe->book1D("hDistSegFromExtrap"+HistoName, "Distance segments from extrap position ",200,0.,200.));
00157     // histo for efficiency from segment counting
00158     histos.push_back(theDbe->book1D("hNaiveEffSeg"+HistoName, "Naive eff ",10,0.,10.));
00159     // histo for efficiency with cuts a-/c-
00160   histos.push_back(theDbe->book2D("hEffSegVsPosDen"+HistoName,"Eff vs local position (all) ",25,-250.,250., 25,-250.,250.));    
00161     // histo for efficiency with cuts a-/b-/c-/d-
00162     histos.push_back(theDbe->book2D("hEffSegVsPosNum"+HistoName, "Eff vs local position ",25,-250.,250., 25,-250.,250.));
00163     // histo for efficiency with cuts a-/b-/c-/d-/e-
00164     histos.push_back(theDbe->book2D("hEffGoodSegVsPosNum"+HistoName, "Eff vs local position (good segs) ", 25,-250.,250., 25,-250.,250.));
00165   }
00166   histosPerCh[chId] = histos;
00167 }
00168 
00169 
00170 void DTChamberEfficiencyTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
00171 
00172   edm::LogVerbatim ("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "[DTChamberEfficiencyTask] Analyze #Run: " << event.id().run()
00173          << " #Event: " << event.id().event();
00174   
00175   // Get the 4D rechit collection from the event
00176   event.getByLabel(theRecHits4DLabel, segs);
00177 
00178   int bottom=0, top=0;
00179 
00180 
00181   // Loop over all the chambers
00182   vector<DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();
00183   vector<DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();
00184   for (; ch_it != ch_end; ++ch_it) {
00185     
00186     DTChamberId ch = (*ch_it)->id();
00187     int wheel =  ch.wheel();
00188     int sector = ch.sector();
00189     int station = ch.station();
00190 
00191 
00192     DTChamberId MidId(wheel, station, sector);
00193     
00194     // get efficiency for MB1 using MB2 and MB3 
00195     if( station == 1 ) {
00196       bottom = 2;
00197       top = 3;
00198     }
00199     
00200     // get efficiency for MB2 using MB1 and MB3
00201     if( station == 2 ) {
00202       bottom = 1;
00203       top = 3;
00204     }
00205 
00206     // get efficiency for MB2 using MB2 and MB4
00207     if( station == 3 ) {
00208       bottom = 2;
00209       top = 4;
00210     }
00211     
00212     // get efficiency for MB4 using MB2 and MB3
00213     if( station == 4 ) {
00214       bottom = 2;
00215       top = 3;
00216     }
00217 
00218     // Select events with (good) segments in Bot and Top
00219     DTChamberId BotId(wheel, bottom, sector);
00220     DTChamberId TopId(wheel, top, sector);      
00221     
00222     // Get segments in the bottom chambers (if any)
00223     DTRecSegment4DCollection::range segsBot= segs->get(BotId);
00224     int nSegsBot=segsBot.second-segsBot.first;
00225     // check if any segments is there
00226     if (nSegsBot==0) continue;
00227 
00228     vector<MonitorElement *> histos =  histosPerCh[MidId];  
00229     
00230     // Get segments in the top chambers (if any)
00231     DTRecSegment4DCollection::range segsTop= segs->get(TopId);
00232     int nSegsTop=segsTop.second-segsTop.first;
00233     
00234     // Select one segment for the bottom chamber
00235     const DTRecSegment4D& bestBotSeg= getBestSegment(segsBot);
00236     
00237     // Select one segment for the top chamber
00238     DTRecSegment4D* pBestTopSeg=0;
00239     if (nSegsTop>0) 
00240       pBestTopSeg = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop));
00241     //if top chamber is MB4 sector 10, consider also sector 14
00242     if (TopId.station() == 4 && TopId.sector() == 10) {
00243       DTChamberId TopId14(wheel, top, 14);
00244       DTRecSegment4DCollection::range segsTop14= segs->get(TopId14);
00245       int nSegsTop14=segsTop14.second-segsTop14.first;
00246       nSegsTop+=nSegsTop14;
00247       if (nSegsTop14) {
00248         DTRecSegment4D* pBestTopSeg14 = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop14));
00249         // get best between sector 10 and 14
00250         pBestTopSeg = const_cast<DTRecSegment4D*>(getBestSegment(pBestTopSeg, pBestTopSeg14));
00251       }
00252     }
00253     if (!pBestTopSeg) continue;
00254     const DTRecSegment4D& bestTopSeg= *pBestTopSeg;
00255     
00256     DTRecSegment4DCollection::range segsMid= segs->get(MidId);
00257     int nSegsMid=segsMid.second-segsMid.first;
00258     
00259     if(detailedAnalysis){
00260       // very trivial efficiency, just count segments
00261       histos[3]->Fill(0);
00262       if (nSegsMid>0) histos[3]->Fill(1);
00263     }
00264     
00265     // get position at Mid by interpolating the position (not direction) of best
00266     // segment in Bot and Top to Mid surface
00267     LocalPoint posAtMid = interpolate(bestBotSeg, bestTopSeg, MidId);
00268     
00269     // is best segment good enough?
00270     if (isGoodSegment(bestBotSeg) && isGoodSegment(bestTopSeg)) {
00271       if(detailedAnalysis)
00272         histos[4]->Fill(posAtMid.x(),posAtMid.y());
00273       //check if interpolation fall inside middle chamber
00274       if ((dtGeom->chamber(MidId))->surface().bounds().inside(posAtMid)) {
00275         histos[0]->Fill(posAtMid.x(),posAtMid.y());         
00276         if (nSegsMid>0) {
00277 
00278           if(detailedAnalysis){
00279             histos[3]->Fill(2);
00280             histos[5]->Fill(posAtMid.x(),posAtMid.y());
00281           }
00282 
00283           const DTRecSegment4D& bestMidSeg= getBestSegment(segsMid);
00284           // check if middle segments is good enough
00285           if (isGoodSegment(bestMidSeg)) {
00286         
00287             if(detailedAnalysis)
00288               histos[6]->Fill(posAtMid.x(),posAtMid.y());
00289             LocalPoint midSegPos=bestMidSeg.localPosition();
00290             
00291             // check if middle segments is also close enough
00292             double dist;
00293             if (bestMidSeg.hasPhi()) { 
00294               if (bestTopSeg.hasZed() && bestBotSeg.hasZed() && bestMidSeg.hasZed()) { 
00295                 dist = (midSegPos-posAtMid).mag();
00296               } else {
00297                     dist = fabs((midSegPos-posAtMid).x());
00298               }
00299             } else {
00300               dist = fabs((midSegPos-posAtMid).y());
00301             }
00302             if (dist < theMinCloseDist ) {
00303               histos[1]->Fill(posAtMid.x(),posAtMid.y());
00304             }
00305             if(detailedAnalysis)
00306               histos[2]->Fill(dist);
00307           }
00308         }
00309       }
00310     }
00311   }// loop over stations
00312 
00313 }
00314 
00315 
00316 
00317         
00318 // requirements : max number of hits and min chi2
00319 const DTRecSegment4D& DTChamberEfficiencyTask::getBestSegment(const DTRecSegment4DCollection::range& segs) const{
00320   DTRecSegment4DCollection::const_iterator bestIter;
00321   unsigned int nHitBest=0;
00322   double chi2Best=99999.;
00323   for (DTRecSegment4DCollection::const_iterator seg=segs.first ;
00324        seg!=segs.second ; ++seg ) {
00325     unsigned int nHits= ((*seg).hasPhi() ? (*seg).phiSegment()->recHits().size() : 0 ) ;
00326     nHits+= ((*seg).hasZed() ?  (*seg).zSegment()->recHits().size() : 0 );
00327 
00328     if (nHits==nHitBest) {
00329       if ((*seg).chi2()/(*seg).degreesOfFreedom() < chi2Best ) {
00330         chi2Best=(*seg).chi2()/(*seg).degreesOfFreedom();
00331         bestIter = seg;
00332       }
00333     }
00334     else if (nHits>nHitBest) {
00335       nHitBest=nHits;
00336       bestIter = seg;
00337     }
00338   }
00339   return *bestIter;
00340 }
00341 
00342 const DTRecSegment4D* DTChamberEfficiencyTask::getBestSegment(const DTRecSegment4D* s1,
00343                                                     const DTRecSegment4D* s2) const{
00344 
00345   if (!s1) return s2;
00346   if (!s2) return s1;
00347   unsigned int nHits1= (s1->hasPhi() ? s1->phiSegment()->recHits().size() : 0 ) ;
00348   nHits1+= (s1->hasZed() ?  s1->zSegment()->recHits().size() : 0 );
00349 
00350   unsigned int nHits2= (s2->hasPhi() ? s2->phiSegment()->recHits().size() : 0 ) ;
00351   nHits2+= (s2->hasZed() ?  s2->zSegment()->recHits().size() : 0 );
00352 
00353   if (nHits1==nHits2) {
00354     if (s1->chi2()/s1->degreesOfFreedom() < s2->chi2()/s2->degreesOfFreedom() ) 
00355       return s1;
00356     else
00357       return s2;
00358   }
00359   else if (nHits1>nHits2) return s1;
00360   return s2;
00361 }
00362 
00363 
00364 LocalPoint DTChamberEfficiencyTask::interpolate(const DTRecSegment4D& seg1,
00365                                             const DTRecSegment4D& seg3,
00366                                             const DTChamberId& id2) const {
00367   // Get GlobalPoition of Seg in MB1
00368   GlobalPoint gpos1=(dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localPosition());
00369 
00370   // Get GlobalPoition of Seg in MB3
00371   GlobalPoint gpos3=(dtGeom->chamber(seg3.chamberId()))->toGlobal(seg3.localPosition());
00372 
00373   // interpolate
00374   // get all in MB2 frame
00375   LocalPoint pos1=(dtGeom->chamber(id2))->toLocal(gpos1);
00376   LocalPoint pos3=(dtGeom->chamber(id2))->toLocal(gpos3);
00377 
00378   // case 1: 1 and 3 has both projection. No problem
00379 
00380   // case 2: one projection is missing for one of the segments. Keep the other's segment position
00381   if (!seg1.hasZed()) pos1=LocalPoint(pos1.x(),pos3.y(),pos1.z());
00382   if (!seg3.hasZed()) pos3=LocalPoint(pos3.x(),pos1.y(),pos3.z());
00383 
00384   if (!seg1.hasPhi()) pos1=LocalPoint(pos3.x(),pos1.y(),pos1.z());
00385   if (!seg3.hasPhi()) pos3=LocalPoint(pos1.x(),pos3.y(),pos3.z());
00386 
00387   // direction
00388   LocalVector dir = (pos3-pos1).unit(); // z points inward!
00389   LocalPoint pos2 = pos1+dir*pos1.z()/(-dir.z());
00390 
00391   return pos2;
00392 }
00393 
00394 
00395 bool DTChamberEfficiencyTask::isGoodSegment(const DTRecSegment4D& seg) const {
00396   if (seg.chamberId().station()!=4 && !seg.hasZed()) return false;
00397   unsigned int nHits= (seg.hasPhi() ? seg.phiSegment()->recHits().size() : 0 ) ;
00398   nHits+= (seg.hasZed() ?  seg.zSegment()->recHits().size() : 0 );
00399   return ( nHits >= theMinHitsSegment &&
00400            seg.chi2()/seg.degreesOfFreedom() < theMinChi2NormSegment );
00401 }