CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/CalibMuon/DTCalibration/plugins/DTTTrigCalibration.cc

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2011/09/21 08:03:43 $
00005  *  $Revision: 1.6 $
00006  *  \author G. Cerminara - INFN Torino
00007  */
00008 #include "CalibMuon/DTCalibration/plugins/DTTTrigCalibration.h"
00009 #include "CalibMuon/DTCalibration/interface/DTTimeBoxFitter.h"
00010 #include "CalibMuon/DTDigiSync/interface/DTTTrigSyncFactory.h"
00011 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
00012 
00013 #include "FWCore/Framework/interface/EventSetup.h"
00014 #include "FWCore/Framework/interface/Event.h"
00015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00016 #include "FWCore/Framework/interface/ESHandle.h"
00017 
00018 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
00019 #include "DataFormats/MuonDetId/interface/DTWireId.h"
00020 
00021 #include "CondFormats/DTObjects/interface/DTTtrig.h"
00022 
00023 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00024 
00025 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
00026 
00027 
00028 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
00029 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
00030 
00031 
00032 #include "TFile.h"
00033 #include "TH1F.h"
00034 #include "TGraph.h"
00035 
00036 class DTLayer;
00037 
00038 using namespace std;
00039 using namespace edm;
00040 // using namespace cond;
00041 
00042 
00043 
00044 // Constructor
00045 DTTTrigCalibration::DTTTrigCalibration(const edm::ParameterSet& pset) {
00046   // Get the debug parameter for verbose output
00047   debug = pset.getUntrackedParameter<bool>("debug");
00048 
00049   // Get the label to retrieve digis from the event
00050   digiLabel = pset.getUntrackedParameter<string>("digiLabel");
00051 
00052   // Switch on/off the DB writing
00053   findTMeanAndSigma = pset.getUntrackedParameter<bool>("fitAndWrite", false);
00054 
00055   // The TDC time-window (ns)
00056   maxTDCCounts = 5000 * pset.getUntrackedParameter<int>("tdcRescale", 1);
00057   //The maximum number of digis per layer
00058   maxDigiPerLayer = pset.getUntrackedParameter<int>("maxDigiPerLayer", 10);
00059 
00060   // The root file which will contain the histos
00061   string rootFileName = pset.getUntrackedParameter<string>("rootFileName");
00062   theFile = new TFile(rootFileName.c_str(), "RECREATE");
00063   theFile->cd();
00064   theFitter = new DTTimeBoxFitter();
00065   if(debug)
00066     theFitter->setVerbosity(1);
00067   
00068   double sigmaFit = pset.getUntrackedParameter<double>("sigmaTTrigFit",10.);
00069   theFitter->setFitSigma(sigmaFit);
00070 
00071   doSubtractT0 = pset.getUntrackedParameter<bool>("doSubtractT0","false");
00072   // Get the synchronizer
00073   if(doSubtractT0) {
00074     theSync = DTTTrigSyncFactory::get()->create(pset.getUntrackedParameter<string>("tTrigMode"),
00075                                                 pset.getUntrackedParameter<ParameterSet>("tTrigModeConfig"));
00076   } else {
00077     theSync = 0;
00078   }
00079 
00080   checkNoisyChannels = pset.getUntrackedParameter<bool>("checkNoisyChannels","false");
00081 
00082   // the kfactor to be uploaded in the ttrig DB
00083   kFactor =  pset.getUntrackedParameter<double>("kFactor",-0.7);
00084 
00085   if(debug) 
00086     cout << "[DTTTrigCalibration]Constructor called!" << endl;
00087 }
00088 
00089 
00090 
00091 // Destructor
00092 DTTTrigCalibration::~DTTTrigCalibration(){
00093   if(debug) 
00094     cout << "[DTTTrigCalibration]Destructor called!" << endl;
00095 
00096   //   // Delete all histos
00097   //   for(map<DTSuperLayerId, TH1F*>::const_iterator slHisto = theHistoMap.begin();
00098   //       slHisto != theHistoMap.end();
00099   //       slHisto++) {
00100   //     delete (*slHisto).second;
00101   //   }
00102 
00103   theFile->Close();
00104   delete theFitter;
00105 }
00106 
00107 
00108 
00110 void DTTTrigCalibration::analyze(const edm::Event & event, const edm::EventSetup& eventSetup) {
00111 
00112   if(debug) 
00113     cout << "[DTTTrigCalibration] #Event: " << event.id().event() << endl;
00114   
00115   // Get the digis from the event
00116   Handle<DTDigiCollection> digis; 
00117   event.getByLabel(digiLabel, digis);
00118 
00119   ESHandle<DTStatusFlag> statusMap;
00120   if(checkNoisyChannels) {
00121     // Get the map of noisy channels
00122     eventSetup.get<DTStatusFlagRcd>().get(statusMap);
00123   }
00124 
00125   if(doSubtractT0)
00126     theSync->setES(eventSetup);
00127  
00128   //The chambers too noisy in this event
00129   vector<DTChamberId> badChambers;
00130 
00131   // Iterate through all digi collections ordered by LayerId   
00132   DTDigiCollection::DigiRangeIterator dtLayerIt;
00133   for (dtLayerIt = digis->begin();
00134        dtLayerIt != digis->end();
00135        ++dtLayerIt){
00136     // Get the iterators over the digis associated with this LayerId
00137     const DTDigiCollection::Range& digiRange = (*dtLayerIt).second; 
00138     
00139     const DTLayerId layerId = (*dtLayerIt).first;
00140     const DTSuperLayerId slId = layerId.superlayerId();
00141     const DTChamberId chId = slId.chamberId();
00142     bool badChamber=false;
00143 
00144     if(debug)
00145       cout<<"----------- Layer "<<layerId<<" -------------"<<endl;
00146 
00147     //Check if the layer is inside a noisy chamber
00148     for(vector<DTChamberId>::const_iterator chamber = badChambers.begin(); chamber != badChambers.end(); chamber++){
00149       if((*chamber) == chId){
00150         badChamber=true;
00151         break;
00152       }
00153     }
00154     if(badChamber) continue;
00155 
00156     //Check if the layer has too many digis
00157     if((digiRange.second - digiRange.first) > maxDigiPerLayer){
00158       if(debug)
00159         cout<<"Layer "<<layerId<<"has too many digis ("<<(digiRange.second - digiRange.first)<<")"<<endl;
00160       badChambers.push_back(chId);
00161       continue;
00162     }
00163 
00164     // Get the histo from the map
00165     TH1F *hTBox = theHistoMap[slId];
00166     if(hTBox == 0) {
00167       // Book the histogram
00168       theFile->cd();
00169       hTBox = new TH1F(getTBoxName(slId).c_str(), "Time box (ns)", int(0.25*32.0*maxTDCCounts/25.0), 0, maxTDCCounts);
00170       if(debug)
00171         cout << "  New Time Box: " << hTBox->GetName() << endl;
00172       theHistoMap[slId] = hTBox;
00173     }
00174     TH1F *hO = theOccupancyMap[layerId];
00175     if(hO == 0) {
00176       // Book the histogram
00177       theFile->cd();
00178       hO = new TH1F(getOccupancyName(layerId).c_str(), "Occupancy", 100, 0, 100);
00179       if(debug)
00180         cout << "  New Time Box: " << hO->GetName() << endl;
00181       theOccupancyMap[layerId] = hO;
00182     }
00183 
00184     // Loop over all digis in the given range
00185     for (DTDigiCollection::const_iterator digi = digiRange.first;
00186          digi != digiRange.second;
00187          digi++) {
00188       const DTWireId wireId(layerId, (*digi).wire());
00189 
00190       // Check for noisy channels and skip them
00191       if(checkNoisyChannels) {
00192         bool isNoisy = false;
00193         bool isFEMasked = false;
00194         bool isTDCMasked = false;
00195         bool isTrigMask = false;
00196         bool isDead = false;
00197         bool isNohv = false;
00198         statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
00199         if(isNoisy) {
00200           if(debug)
00201             cout << "Wire: " << wireId << " is noisy, skipping!" << endl;
00202           continue;
00203         }      
00204       }
00205       theFile->cd();
00206       double offset = 0;
00207       if(doSubtractT0) {
00208         const DTLayer* layer = 0;//fake
00209         const GlobalPoint glPt;//fake
00210         offset = theSync->offset(layer, wireId, glPt);
00211       }
00212       hTBox->Fill((*digi).time()-offset);
00213       if(debug) {
00214         cout << "   Filling Time Box:   " << hTBox->GetName() << endl;
00215         cout << "           offset (ns): " << offset << endl;
00216         cout << "           time(ns):   " << (*digi).time()-offset<< endl;
00217       }
00218       hO->Fill((*digi).wire());
00219     }
00220   }
00221 }
00222 
00223 
00224 void DTTTrigCalibration::endJob() {
00225   if(debug) 
00226     cout << "[DTTTrigCalibration]Writing histos to file!" << endl;
00227   
00228   // Write all time boxes to file
00229   theFile->cd();
00230   for(map<DTSuperLayerId, TH1F*>::const_iterator slHisto = theHistoMap.begin();
00231       slHisto != theHistoMap.end();
00232       slHisto++) {
00233     (*slHisto).second->Write();
00234   }
00235   for(map<DTLayerId, TH1F*>::const_iterator slHisto = theOccupancyMap.begin();
00236       slHisto != theOccupancyMap.end();
00237       slHisto++) {
00238     (*slHisto).second->Write();
00239   }
00240 
00241   if(findTMeanAndSigma) {
00242       // Create the object to be written to DB
00243       DTTtrig* tTrig = new DTTtrig();
00244 
00245       // Loop over the map, fit the histos and write the resulting values to the DB
00246       for(map<DTSuperLayerId, TH1F*>::const_iterator slHisto = theHistoMap.begin();
00247           slHisto != theHistoMap.end();
00248           slHisto++) {
00249         pair<double, double> meanAndSigma = theFitter->fitTimeBox((*slHisto).second);
00250         tTrig->set((*slHisto).first,
00251                    meanAndSigma.first,
00252                    meanAndSigma.second,
00253                    kFactor,
00254                    DTTimeUnits::ns);
00255     
00256         if(debug) {
00257           cout << " SL: " << (*slHisto).first
00258                << " mean = " << meanAndSigma.first
00259                << " sigma = " << meanAndSigma.second << endl;
00260         }
00261       }
00262 
00263       // Print the ttrig map
00264       dumpTTrigMap(tTrig);
00265   
00266       // Plot the tTrig
00267       plotTTrig(tTrig);
00268 
00269       if(debug) 
00270         cout << "[DTTTrigCalibration]Writing ttrig object to DB!" << endl;
00271 
00272 
00273       // FIXME: to be read from cfg?
00274       string tTrigRecord = "DTTtrigRcd";
00275       
00276       // Write the object to DB
00277       DTCalibDBUtils::writeToDB(tTrigRecord, tTrig);
00278     }  
00279 
00280 }
00281 
00282 
00283 
00284 
00285 string DTTTrigCalibration::getTBoxName(const DTSuperLayerId& slId) const {
00286   string histoName;
00287   stringstream theStream;
00288   theStream << "Ch_" << slId.wheel() << "_" << slId.station() << "_" << slId.sector()
00289             << "_SL" << slId.superlayer() << "_hTimeBox";
00290   theStream >> histoName;
00291   return histoName;
00292 }
00293 
00294 string DTTTrigCalibration::getOccupancyName(const DTLayerId& slId) const {
00295   string histoName;
00296   stringstream theStream;
00297   theStream << "Ch_" << slId.wheel() << "_" << slId.station() << "_" << slId.sector()
00298             << "_SL" << slId.superlayer() << "_L"<< slId.layer() <<"_Occupancy";
00299   theStream >> histoName;
00300   return histoName;
00301 }
00302 
00303 
00304 void DTTTrigCalibration::dumpTTrigMap(const DTTtrig* tTrig) const {
00305   static const double convToNs = 25./32.;
00306   for(DTTtrig::const_iterator ttrig = tTrig->begin();
00307       ttrig != tTrig->end(); ttrig++) {
00308     cout << "Wh: " << (*ttrig).first.wheelId
00309          << " St: " << (*ttrig).first.stationId
00310          << " Sc: " << (*ttrig).first.sectorId
00311          << " Sl: " << (*ttrig).first.slId
00312          << " TTrig mean (ns): " << (*ttrig).second.tTrig * convToNs
00313          << " TTrig sigma (ns): " << (*ttrig).second.tTrms * convToNs<< endl;
00314   }
00315 }
00316 
00317 
00318 void DTTTrigCalibration::plotTTrig(const DTTtrig* tTrig) const {
00319 
00320   TH1F* tTrig_YB1_Se10 = new TH1F("tTrig_YB1_Se10","tTrig YB1_Se10",15,1,16);
00321   TH1F* tTrig_YB2_Se10 = new TH1F("tTrig_YB2_Se10","tTrig YB2_Se10",15,1,16);
00322   TH1F* tTrig_YB2_Se11 = new TH1F("tTrig_YB2_Se11","tTrig YB2_Se11",12,1,13);
00323 
00324   static const double convToNs = 25./32.;
00325   for(DTTtrig::const_iterator ttrig = tTrig->begin();
00326       ttrig != tTrig->end(); ttrig++) {
00327 
00328     // avoid to have wired numbers in the plot
00329     float tTrigValue=0;
00330     float tTrmsValue=0;
00331     if ((*ttrig).second.tTrig * convToNs > 0 &&
00332         (*ttrig).second.tTrig * convToNs < 32000 ) {
00333       tTrigValue = (*ttrig).second.tTrig * convToNs;
00334       tTrmsValue = (*ttrig).second.tTrms * convToNs;
00335     }
00336 
00337     int binx;
00338     string binLabel;
00339     stringstream binLabelStream;
00340     if ((*ttrig).first.sectorId != 14) {
00341       binx = ((*ttrig).first.stationId-1)*3  + (*ttrig).first.slId;
00342       binLabelStream << "MB"<<(*ttrig).first.stationId<<"_SL"<<(*ttrig).first.slId;
00343     }
00344     else {
00345       binx = 12  + (*ttrig).first.slId;
00346       binLabelStream << "MB14_SL"<<(*ttrig).first.slId;
00347     }
00348     binLabelStream >> binLabel;
00349 
00350     if ((*ttrig).first.wheelId == 2) {
00351       if ((*ttrig).first.sectorId == 10 || (*ttrig).first.sectorId == 14) {
00352         tTrig_YB2_Se10->Fill( binx,tTrigValue);
00353         tTrig_YB2_Se10->SetBinError( binx, tTrmsValue);
00354         tTrig_YB2_Se10->GetXaxis()->SetBinLabel(binx,binLabel.c_str());
00355         tTrig_YB2_Se10->GetYaxis()->SetTitle("ns");
00356       }
00357       else {
00358         tTrig_YB2_Se11->Fill( binx,tTrigValue);
00359         tTrig_YB2_Se11->SetBinError( binx,tTrmsValue);
00360         tTrig_YB2_Se11->GetXaxis()->SetBinLabel(binx,binLabel.c_str());
00361         tTrig_YB2_Se11->GetYaxis()->SetTitle("ns");
00362       }
00363     }
00364     else {
00365       tTrig_YB1_Se10->Fill( binx,tTrigValue);
00366       tTrig_YB1_Se10->SetBinError( binx,tTrmsValue);
00367       tTrig_YB1_Se10->GetXaxis()->SetBinLabel(binx,binLabel.c_str());
00368       tTrig_YB1_Se10->GetYaxis()->SetTitle("ns");
00369     }
00370   }
00371 
00372   tTrig_YB1_Se10->Write();
00373   tTrig_YB2_Se10->Write();
00374   tTrig_YB2_Se11->Write();
00375 
00376 }