CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/CalibMuon/DTCalibration/plugins/DTTPAnalyzer.cc

Go to the documentation of this file.
00001 
00008 #include "FWCore/Framework/interface/EDAnalyzer.h"
00009 #include "FWCore/Utilities/interface/InputTag.h"
00010 #include "FWCore/Framework/interface/ESHandle.h"
00011 
00012 #include <string>
00013 #include <map>
00014 
00015 //class DTT0;
00016 class DTLayerId;
00017 class DTWireId;
00018 class DTGeometry;
00019 class DTTTrigBaseSync;
00020 class TFile;
00021 
00022 class DTTPAnalyzer : public edm::EDAnalyzer {
00023 public:
00024   DTTPAnalyzer( const edm::ParameterSet& );
00025   virtual ~DTTPAnalyzer();
00026 
00027   //void beginJob();
00028   void beginRun( const edm::Run& , const edm::EventSetup& );
00029   void analyze( const edm::Event& , const edm::EventSetup& );
00030   void endJob();
00031   
00032 private:
00033   std::string getHistoName( const DTLayerId& ); 
00034 
00035   bool subtractT0_;
00036   edm::InputTag digiLabel_;
00037 
00038   TFile* rootFile_;
00039   //const DTT0* tZeroMap_;
00040   edm::ESHandle<DTGeometry> dtGeom_;
00041   DTTTrigBaseSync* tTrigSync_;
00042 
00043   // Map of the t0 and sigma histos by layer
00044   std::map<DTWireId, int> nDigisPerWire_;
00045   std::map<DTWireId, double> sumWPerWire_;
00046   std::map<DTWireId, double> sumW2PerWire_;
00047   //std::map<DTLayerId, TH1F*> meanHistoMap_;
00048   //std::map<DTLayerId, TH1F*> sigmaHistoMap_;
00049 };
00050 
00051 #include "FWCore/Framework/interface/Event.h"
00052 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00053 #include "FWCore/Framework/interface/MakerMacros.h"
00054 
00055 #include "CalibMuon/DTDigiSync/interface/DTTTrigSyncFactory.h"
00056 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
00057 
00058 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00059 #include "Geometry/Records/interface/MuonNumberingRecord.h"
00060 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00061 
00062 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
00063 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00064 
00065 #include "CondFormats/DTObjects/interface/DTT0.h"
00066 
00067 #include "TH1F.h"
00068 #include "TFile.h"
00069 
00070 DTTPAnalyzer::DTTPAnalyzer(const edm::ParameterSet& pset):
00071   subtractT0_(pset.getParameter<bool>("subtractT0")),
00072   digiLabel_(pset.getParameter<edm::InputTag>("digiLabel")),
00073   tTrigSync_(0) {
00074 
00075   std::string rootFileName = pset.getUntrackedParameter<std::string>("rootFileName");
00076   rootFile_ = new TFile(rootFileName.c_str(), "RECREATE");
00077   rootFile_->cd();
00078 
00079   if(subtractT0_) 
00080      tTrigSync_ = DTTTrigSyncFactory::get()->create(pset.getParameter<std::string>("tTrigMode"),
00081                                                     pset.getParameter<edm::ParameterSet>("tTrigModeConfig"));
00082 
00083 }
00084  
00085 DTTPAnalyzer::~DTTPAnalyzer(){  
00086   rootFile_->Close();
00087 }
00088 
00089 void DTTPAnalyzer::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
00090   // Get the t0 map from the DB
00091   if(subtractT0_){ 
00092      /*ESHandle<DTT0> t0;
00093      setup.get<DTT0Rcd>().get(t0);
00094      tZeroMap_ = &*t0;*/
00095      tTrigSync_->setES(setup);
00096   }
00097   // Get the DT Geometry  
00098   setup.get<MuonGeometryRecord>().get(dtGeom_);
00099 }
00100 
00101 void DTTPAnalyzer::analyze(const edm::Event & event, const edm::EventSetup& setup) {
00102 
00103   // Get the digis from the event
00104   edm::Handle<DTDigiCollection> digis; 
00105   event.getByLabel(digiLabel_, digis);
00106 
00107   // Iterate through all digi collections ordered by LayerId   
00108   DTDigiCollection::DigiRangeIterator dtLayerIt;
00109   for (dtLayerIt = digis->begin();
00110        dtLayerIt != digis->end();
00111        ++dtLayerIt){
00112     // Get the iterators over the digis associated with this LayerId
00113     const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
00114   
00115     // Get the layerId
00116     const DTLayerId layerId = (*dtLayerIt).first; //FIXME: check to be in the right sector
00117 
00118     // Loop over all digis in the given layer
00119     for (DTDigiCollection::const_iterator digi = digiRange.first;
00120          digi != digiRange.second;
00121          digi++) {
00122        const DTWireId wireId( layerId, (*digi).wire() );
00123 
00124        double t0 = (*digi).countsTDC();
00125 
00126        //FIXME: Reject digis not coming from TP
00127 
00128        if(subtractT0_) {
00129           const DTLayer* layer = 0; //fake
00130           const GlobalPoint glPt; //fake
00131           double offset = tTrigSync_->offset(layer, wireId, glPt);
00132           t0 -= offset;
00133        }
00134 
00135        if(nDigisPerWire_.find(wireId) == nDigisPerWire_.end()){
00136           nDigisPerWire_[wireId] = 0;
00137           sumWPerWire_[wireId] = 0.;
00138           sumW2PerWire_[wireId] = 0.;  
00139        }
00140 
00141        ++nDigisPerWire_[wireId]; 
00142        sumWPerWire_[wireId] += t0;
00143        sumW2PerWire_[wireId] += t0*t0;
00144     }
00145 
00146   }
00147 }
00148 
00149 void DTTPAnalyzer::endJob() {
00150   rootFile_->cd();
00151   std::map<DTLayerId, TH1F*> meanHistoMap;
00152   std::map<DTLayerId, TH1F*> sigmaHistoMap; 
00153   for(std::map<DTWireId, int>::const_iterator wireIdIt = nDigisPerWire_.begin();
00154                                               wireIdIt != nDigisPerWire_.end(); ++wireIdIt){
00155      DTWireId wireId((*wireIdIt).first);
00156 
00157      int nDigis = nDigisPerWire_[wireId];
00158      double sumW = sumWPerWire_[wireId];
00159      double sumW2 = sumW2PerWire_[wireId]; 
00160 
00161      double mean = sumW/nDigis;
00162      double rms = sumW2/nDigis - mean*mean;
00163      rms = sqrt(rms);
00164 
00165      DTLayerId layerId = wireId.layerId();
00166      if(meanHistoMap.find(layerId) == meanHistoMap.end()) {
00167         std::string histoName = getHistoName(layerId);
00168         const int firstChannel = dtGeom_->layer(layerId)->specificTopology().firstChannel();
00169         const int nWires = dtGeom_->layer(layerId)->specificTopology().channels();
00170         TH1F* meanHistoTP = new TH1F((histoName + "_tpMean").c_str(),"mean from test pulses by channel", 
00171                                       nWires,firstChannel,(firstChannel + nWires));
00172         TH1F* sigmaHistoTP = new TH1F((histoName + "_tpSigma").c_str(),"sigma from test pulses by channel",
00173                                       nWires,firstChannel,(firstChannel + nWires));
00174         meanHistoMap[layerId] = meanHistoTP;
00175         sigmaHistoMap[layerId] = sigmaHistoTP;
00176      }
00177      // Fill the histograms
00178      int nBin = meanHistoMap[layerId]->GetXaxis()->FindFixBin(wireId.wire());
00179      meanHistoMap[layerId]->SetBinContent(nBin,mean);
00180      sigmaHistoMap[layerId]->SetBinContent(nBin,rms);
00181   }
00182 
00183   for(std::map<DTLayerId, TH1F*>::const_iterator key = meanHistoMap.begin();
00184                                                  key != meanHistoMap.end(); ++key){
00185      meanHistoMap[(*key).first]->Write();
00186      sigmaHistoMap[(*key).first]->Write(); 
00187   } 
00188 
00189 }
00190 
00191 std::string DTTPAnalyzer::getHistoName(const DTLayerId& lId) {
00192   std::string histoName;
00193   std::stringstream theStream;
00194   theStream << "Ch_" << lId.wheel() << "_" << lId.station() << "_" << lId.sector()
00195             << "_SL" << lId.superlayer() << "_L" << lId.layer();
00196   theStream >> histoName;
00197   return histoName;
00198 }
00199 
00200 DEFINE_FWK_MODULE(DTTPAnalyzer);