Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:26 

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Mila - INFN Torino
0005  */
0006 
0007 #include "CalibMuon/DTCalibration/plugins/DTNoiseComputation.h"
0008 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
0009 
0010 // Framework
0011 #include "FWCore/Framework/interface/IOVSyncValue.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 
0017 // Geometry
0018 #include "Geometry/DTGeometry/interface/DTLayer.h"
0019 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0020 #include "Geometry/DTGeometry/interface/DTTopology.h"
0021 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0022 
0023 // Digis
0024 #include "DataFormats/DTDigi/interface/DTDigi.h"
0025 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
0026 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
0027 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
0028 
0029 #include "TH1F.h"
0030 #include "TH2F.h"
0031 #include "TFile.h"
0032 #include "TF1.h"
0033 #include "TProfile.h"
0034 #include "TPostScript.h"
0035 #include "TCanvas.h"
0036 #include "TLegend.h"
0037 
0038 using namespace edm;
0039 using namespace std;
0040 
0041 DTNoiseComputation::DTNoiseComputation(const edm::ParameterSet &ps) : dtGeomToken_(esConsumes()) {
0042   cout << "[DTNoiseComputation]: Constructor" << endl;
0043 
0044   // Get the debug parameter for verbose output
0045   debug = ps.getUntrackedParameter<bool>("debug");
0046 
0047   // The analysis type
0048   fastAnalysis = ps.getUntrackedParameter<bool>("fastAnalysis", true);
0049 
0050   // The root file which contain the histos
0051   string rootFileName = ps.getUntrackedParameter<string>("rootFileName");
0052   theFile = new TFile(rootFileName.c_str(), "READ");
0053 
0054   // The new root file which contain the histos
0055   string newRootFileName = ps.getUntrackedParameter<string>("newRootFileName");
0056   theNewFile = new TFile(newRootFileName.c_str(), "RECREATE");
0057 
0058   // The maximum number of events to analyze
0059   MaxEvents = ps.getUntrackedParameter<int>("MaxEvents");
0060 }
0061 
0062 void DTNoiseComputation::beginRun(const edm::Run &, const EventSetup &setup) {
0063   // Get the DT Geometry
0064   dtGeom = setup.getHandle(dtGeomToken_);
0065 
0066   static int count = 0;
0067 
0068   if (count == 0) {
0069     string CheckHistoName;
0070 
0071     TH1F *hOccHisto;
0072     TH1F *hAverageNoiseHisto;
0073     TH1F *hAverageNoiseIntegratedHisto;
0074     TH1F *hAverageNoiseHistoPerCh;
0075     TH1F *hAverageNoiseIntegratedHistoPerCh;
0076     TH2F *hEvtHisto;
0077     string HistoName;
0078     string Histo2Name;
0079     string AverageNoiseName;
0080     string AverageNoiseIntegratedName;
0081     string AverageNoiseNamePerCh;
0082     string AverageNoiseIntegratedNamePerCh;
0083     TH1F *hnoisyC;
0084     TH1F *hsomeHowNoisyC;
0085 
0086     // Loop over all the chambers
0087     vector<const DTChamber *>::const_iterator ch_it = dtGeom->chambers().begin();
0088     vector<const DTChamber *>::const_iterator ch_end = dtGeom->chambers().end();
0089     // Loop over the SLs
0090     for (; ch_it != ch_end; ++ch_it) {
0091       DTChamberId ch = (*ch_it)->id();
0092       vector<const DTSuperLayer *>::const_iterator sl_it = (*ch_it)->superLayers().begin();
0093       vector<const DTSuperLayer *>::const_iterator sl_end = (*ch_it)->superLayers().end();
0094       // Loop over the SLs
0095       for (; sl_it != sl_end; ++sl_it) {
0096         //      DTSuperLayerId sl = (*sl_it)->id();
0097         vector<const DTLayer *>::const_iterator l_it = (*sl_it)->layers().begin();
0098         vector<const DTLayer *>::const_iterator l_end = (*sl_it)->layers().end();
0099         // Loop over the Ls
0100         for (; l_it != l_end; ++l_it) {
0101           DTLayerId dtLId = (*l_it)->id();
0102 
0103           //check if the layer has digi
0104           theFile->cd();
0105           CheckHistoName = "DigiOccupancy_" + getLayerName(dtLId);
0106           TH1F *hCheckHisto = (TH1F *)theFile->Get(CheckHistoName.c_str());
0107           if (hCheckHisto) {
0108             delete hCheckHisto;
0109             string wheel = std::to_string(ch.wheel());
0110             string station = std::to_string(ch.station());
0111 
0112             if (someHowNoisyC.find(make_pair(ch.wheel(), ch.station())) == someHowNoisyC.end()) {
0113               TString histoName_someHowNoisy = "somehowNoisyCell_W" + wheel + "_St" + station;
0114               hsomeHowNoisyC =
0115                   new TH1F(histoName_someHowNoisy, histoName_someHowNoisy, getMaxNumBins(ch), 1, getMaxNumBins(ch) + 1);
0116               someHowNoisyC[make_pair(ch.wheel(), ch.station())] = hsomeHowNoisyC;
0117             }
0118 
0119             if (noisyC.find(make_pair(ch.wheel(), ch.station())) == noisyC.end()) {
0120               TString histoName_noisy = "noisyCell_W" + wheel + "_St" + station;
0121               hnoisyC = new TH1F(histoName_noisy, histoName_noisy, getMaxNumBins(ch), 1, getMaxNumBins(ch) + 1);
0122               noisyC[make_pair(ch.wheel(), ch.station())] = hnoisyC;
0123             }
0124 
0125             //to fill a map with the average noise per wire and fill new noise histo
0126             if (AvNoisePerSuperLayer.find(dtLId.superlayerId()) == AvNoisePerSuperLayer.end()) {
0127               AverageNoiseName = "AverageNoise_" + getSuperLayerName(dtLId.superlayerId());
0128               hAverageNoiseHisto = new TH1F(AverageNoiseName.c_str(), AverageNoiseName.c_str(), 200, 0, 10000);
0129               AverageNoiseIntegratedName = "AverageNoiseIntegrated_" + getSuperLayerName(dtLId.superlayerId());
0130               hAverageNoiseIntegratedHisto =
0131                   new TH1F(AverageNoiseIntegratedName.c_str(), AverageNoiseIntegratedName.c_str(), 200, 0, 10000);
0132               AvNoisePerSuperLayer[dtLId.superlayerId()] = hAverageNoiseHisto;
0133               AvNoiseIntegratedPerSuperLayer[dtLId.superlayerId()] = hAverageNoiseIntegratedHisto;
0134               if (debug) {
0135                 cout << "  New Average Noise Histo per SuperLayer : " << hAverageNoiseHisto->GetName() << endl;
0136                 cout << "  New Average Noise Integrated Histo per SuperLayer : " << hAverageNoiseHisto->GetName()
0137                      << endl;
0138               }
0139             }
0140             if (AvNoisePerChamber.find(dtLId.superlayerId().chamberId()) == AvNoisePerChamber.end()) {
0141               AverageNoiseNamePerCh = "AverageNoise_" + getChamberName(dtLId);
0142               hAverageNoiseHistoPerCh =
0143                   new TH1F(AverageNoiseNamePerCh.c_str(), AverageNoiseNamePerCh.c_str(), 200, 0, 10000);
0144               AverageNoiseIntegratedNamePerCh = "AverageNoiseIntegrated_" + getChamberName(dtLId);
0145               hAverageNoiseIntegratedHistoPerCh = new TH1F(
0146                   AverageNoiseIntegratedNamePerCh.c_str(), AverageNoiseIntegratedNamePerCh.c_str(), 200, 0, 10000);
0147               AvNoisePerChamber[dtLId.superlayerId().chamberId()] = hAverageNoiseHistoPerCh;
0148               AvNoiseIntegratedPerChamber[dtLId.superlayerId().chamberId()] = hAverageNoiseIntegratedHistoPerCh;
0149               if (debug)
0150                 cout << "  New Average Noise Histo per chamber : " << hAverageNoiseHistoPerCh->GetName() << endl;
0151             }
0152 
0153             HistoName = "DigiOccupancy_" + getLayerName(dtLId);
0154             theFile->cd();
0155             hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0156             int numBin = hOccHisto->GetXaxis()->GetNbins();
0157             for (int bin = 1; bin <= numBin; bin++) {
0158               DTWireId wireID(dtLId, bin);
0159               theAverageNoise[wireID] = hOccHisto->GetBinContent(bin);
0160               if (theAverageNoise[wireID] != 0) {
0161                 AvNoisePerSuperLayer[dtLId.superlayerId()]->Fill(theAverageNoise[wireID]);
0162                 AvNoisePerChamber[dtLId.superlayerId().chamberId()]->Fill(theAverageNoise[wireID]);
0163               }
0164             }
0165 
0166             //to compute the average noise per layer (excluding the noisy cells)
0167             double numCell = 0;
0168             double AvNoise = 0;
0169             HistoName = "DigiOccupancy_" + getLayerName(dtLId);
0170             theFile->cd();
0171             hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0172             numBin = hOccHisto->GetXaxis()->GetNbins();
0173             for (int bin = 1; bin <= numBin; bin++) {
0174               DTWireId wireID(dtLId, bin);
0175               theAverageNoise[wireID] = hOccHisto->GetBinContent(bin);
0176               if (hOccHisto->GetBinContent(bin) < 100) {
0177                 numCell++;
0178                 AvNoise += hOccHisto->GetBinContent(bin);
0179               }
0180               if (hOccHisto->GetBinContent(bin) > 100 && hOccHisto->GetBinContent(bin) < 500) {
0181                 someHowNoisyC[make_pair(ch.wheel(), ch.station())]->Fill(bin);
0182                 cout << "filling somehow noisy cell" << endl;
0183               }
0184               if (hOccHisto->GetBinContent(bin) > 500) {
0185                 noisyC[make_pair(ch.wheel(), ch.station())]->Fill(bin);
0186                 cout << "filling noisy cell" << endl;
0187               }
0188             }
0189             AvNoise = AvNoise / numCell;
0190             cout << "theAverageNoise for layer " << getLayerName(dtLId) << " is : " << AvNoise << endl;
0191 
0192             // book the digi event plots every 1000 events
0193             int updates = MaxEvents / 1000;
0194             for (int evt = 0; evt < updates; evt++) {
0195               Histo2Name = "DigiPerWirePerEvent_" + getLayerName(dtLId) + "_" + std::to_string(evt);
0196               theFile->cd();
0197               hEvtHisto = (TH2F *)theFile->Get(Histo2Name.c_str());
0198               if (hEvtHisto) {
0199                 if (debug)
0200                   cout << "  New Histo with the number of events per evt per wire: " << hEvtHisto->GetName() << endl;
0201                 theEvtMap[dtLId].push_back(hEvtHisto);
0202               }
0203             }
0204 
0205           }  // done if the layer has digi
0206         }    // loop over layers
0207       }      // loop over superlayers
0208     }        // loop over chambers
0209 
0210     count++;
0211   }
0212 }
0213 
0214 void DTNoiseComputation::endJob() {
0215   cout << "[DTNoiseComputation] endjob called!" << endl;
0216   TH1F *hEvtDistance = nullptr;
0217   TF1 *ExpoFit = new TF1("ExpoFit", "expo", 0.5, 1000.5);
0218   ExpoFit->SetMarkerColor();  //just silence gcc complaining about unused vars
0219   TProfile *theNoiseHisto = new TProfile("theNoiseHisto", "Time Constant versus Average Noise", 100000, 0, 100000);
0220 
0221   //compute the time constant
0222   for (map<DTLayerId, vector<TH2F *> >::const_iterator lHisto = theEvtMap.begin(); lHisto != theEvtMap.end();
0223        ++lHisto) {
0224     for (int bin = 1; bin < (*lHisto).second[0]->GetYaxis()->GetNbins(); bin++) {
0225       int distanceEvt = 1;
0226       DTWireId wire((*lHisto).first, bin);
0227       for (int i = 0; i < int((*lHisto).second.size()); i++) {
0228         for (int evt = 1; evt <= (*lHisto).second[i]->GetXaxis()->GetNbins(); evt++) {
0229           if ((*lHisto).second[i]->GetBinContent(evt, bin) == 0)
0230             distanceEvt++;
0231           else {
0232             if (toDel.find(wire) == toDel.end()) {
0233               toDel[wire] = false;
0234               string Histo = "EvtDistancePerWire_" + getLayerName((*lHisto).first) + "_" + std::to_string(bin);
0235               hEvtDistance = new TH1F(Histo.c_str(), Histo.c_str(), 50000, 0.5, 50000.5);
0236             }
0237             hEvtDistance->Fill(distanceEvt);
0238             distanceEvt = 1;
0239           }
0240         }
0241       }
0242       if (toDel.find(wire) != toDel.end()) {
0243         theHistoEvtDistancePerWire[wire] = hEvtDistance;
0244         theNewFile->cd();
0245         theHistoEvtDistancePerWire[wire]->Fit("ExpoFit", "R");
0246         TF1 *funct = theHistoEvtDistancePerWire[wire]->GetFunction("ExpoFit");
0247         double par0 = funct->GetParameter(0);
0248         double par1 = funct->GetParameter(1);
0249         cout << "par0: " << par0 << "  par1: " << par1 << endl;
0250         double chi2rid = funct->GetChisquare() / funct->GetNDF();
0251         if (chi2rid > 10)
0252           theTimeConstant[wire] = 1;
0253         else
0254           theTimeConstant[wire] = par1;
0255         toDel[wire] = true;
0256         theHistoEvtDistancePerWire[wire]->Write();
0257         delete hEvtDistance;
0258       }
0259     }
0260   }
0261 
0262   if (!fastAnalysis) {
0263     //fill the histo with the time constant as a function of the average noise
0264     for (map<DTWireId, double>::const_iterator AvNoise = theAverageNoise.begin(); AvNoise != theAverageNoise.end();
0265          ++AvNoise) {
0266       DTWireId wire = (*AvNoise).first;
0267       theNoiseHisto->Fill((*AvNoise).second, theTimeConstant[wire]);
0268       cout << "Layer: " << getLayerName(wire.layerId()) << "  wire: " << wire.wire() << endl;
0269       cout << "The Average noise: " << (*AvNoise).second << endl;
0270       cout << "The time constant: " << theTimeConstant[wire] << endl;
0271     }
0272     theNewFile->cd();
0273     theNoiseHisto->Write();
0274   }
0275 
0276   // histos with the integrated noise per layer
0277   int numBin;
0278   double integratedNoise, bin, halfBin, maxBin;
0279   for (map<DTSuperLayerId, TH1F *>::const_iterator AvNoiseHisto = AvNoisePerSuperLayer.begin();
0280        AvNoiseHisto != AvNoisePerSuperLayer.end();
0281        ++AvNoiseHisto) {
0282     integratedNoise = 0;
0283     numBin = (*AvNoiseHisto).second->GetXaxis()->GetNbins();
0284     maxBin = (*AvNoiseHisto).second->GetXaxis()->GetXmax();
0285     bin = double(maxBin / numBin);
0286     halfBin = double(bin / 2);
0287     theNewFile->cd();
0288     (*AvNoiseHisto).second->Write();
0289     for (int i = 1; i < numBin; i++) {
0290       integratedNoise += (*AvNoiseHisto).second->GetBinContent(i);
0291       AvNoiseIntegratedPerSuperLayer[(*AvNoiseHisto).first]->Fill(halfBin, integratedNoise);
0292       halfBin += bin;
0293     }
0294     theNewFile->cd();
0295     AvNoiseIntegratedPerSuperLayer[(*AvNoiseHisto).first]->Write();
0296   }
0297   // histos with the integrated noise per chamber
0298   for (map<DTChamberId, TH1F *>::const_iterator AvNoiseHisto = AvNoisePerChamber.begin();
0299        AvNoiseHisto != AvNoisePerChamber.end();
0300        ++AvNoiseHisto) {
0301     integratedNoise = 0;
0302     numBin = (*AvNoiseHisto).second->GetXaxis()->GetNbins();
0303     maxBin = (*AvNoiseHisto).second->GetXaxis()->GetXmax();
0304     bin = maxBin / numBin;
0305     halfBin = bin / 2;
0306     theNewFile->cd();
0307     (*AvNoiseHisto).second->Write();
0308     for (int i = 1; i < numBin; i++) {
0309       integratedNoise += (*AvNoiseHisto).second->GetBinContent(i);
0310       AvNoiseIntegratedPerChamber[(*AvNoiseHisto).first]->Fill(halfBin, integratedNoise);
0311       halfBin += bin;
0312     }
0313     theNewFile->cd();
0314     AvNoiseIntegratedPerChamber[(*AvNoiseHisto).first]->Write();
0315   }
0316 
0317   //overimpose the average noise histo
0318   bool histo = false;
0319   vector<const DTChamber *>::const_iterator chamber_it = dtGeom->chambers().begin();
0320   vector<const DTChamber *>::const_iterator chamber_end = dtGeom->chambers().end();
0321   // Loop over the chambers
0322   for (; chamber_it != chamber_end; ++chamber_it) {
0323     vector<const DTSuperLayer *>::const_iterator sl_it = (*chamber_it)->superLayers().begin();
0324     vector<const DTSuperLayer *>::const_iterator sl_end = (*chamber_it)->superLayers().end();
0325     // Loop over the SLs
0326     for (; sl_it != sl_end; ++sl_it) {
0327       DTSuperLayerId sl = (*sl_it)->id();
0328       vector<const DTLayer *>::const_iterator l_it = (*sl_it)->layers().begin();
0329       vector<const DTLayer *>::const_iterator l_end = (*sl_it)->layers().end();
0330 
0331       string canvasName = "c" + getSuperLayerName(sl);
0332       TCanvas c1(canvasName.c_str(), canvasName.c_str(), 600, 780);
0333       TLegend *leg = new TLegend(0.5, 0.6, 0.7, 0.8);
0334       for (; l_it != l_end; ++l_it) {
0335         DTLayerId layerId = (*l_it)->id();
0336         string HistoName = "DigiOccupancy_" + getLayerName(layerId);
0337         theFile->cd();
0338         TH1F *hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0339         if (hOccHisto) {
0340           string TitleHisto = "AverageNoise_" + getSuperLayerName(sl);
0341           cout << "TitleHisto : " << TitleHisto << endl;
0342           hOccHisto->SetTitle(TitleHisto.c_str());
0343           string legendHisto = "layer " + std::to_string(layerId.layer());
0344           leg->AddEntry(hOccHisto, legendHisto.c_str(), "L");
0345           hOccHisto->SetMaximum(getYMaximum(sl));
0346           histo = true;
0347           if (layerId.layer() == 1)
0348             hOccHisto->Draw();
0349           else
0350             hOccHisto->Draw("same");
0351           hOccHisto->SetLineColor(layerId.layer());
0352         }
0353       }
0354       if (histo) {
0355         leg->Draw("same");
0356         theNewFile->cd();
0357         c1.Write();
0358       }
0359     }
0360     histo = false;
0361   }
0362 
0363   //write on file the noisy plots
0364   for (map<pair<int, int>, TH1F *>::const_iterator nCell = noisyC.begin(); nCell != noisyC.end(); ++nCell) {
0365     theNewFile->cd();
0366     (*nCell).second->Write();
0367   }
0368   for (map<pair<int, int>, TH1F *>::const_iterator somehownCell = someHowNoisyC.begin();
0369        somehownCell != someHowNoisyC.end();
0370        ++somehownCell) {
0371     theNewFile->cd();
0372     (*somehownCell).second->Write();
0373   }
0374 }
0375 
0376 DTNoiseComputation::~DTNoiseComputation() {
0377   theFile->Close();
0378   theNewFile->Close();
0379 }
0380 
0381 string DTNoiseComputation::getLayerName(const DTLayerId &lId) const {
0382   const DTSuperLayerId dtSLId = lId.superlayerId();
0383   const DTChamberId dtChId = dtSLId.chamberId();
0384   string layerName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0385                      std::to_string(dtChId.sector()) + "_SL" + std::to_string(dtSLId.superlayer()) + "_L" +
0386                      std::to_string(lId.layer());
0387 
0388   return layerName;
0389 }
0390 
0391 string DTNoiseComputation::getSuperLayerName(const DTSuperLayerId &dtSLId) const {
0392   const DTChamberId dtChId = dtSLId.chamberId();
0393 
0394   string superLayerName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0395                           std::to_string(dtChId.sector()) + "_SL" + std::to_string(dtSLId.superlayer());
0396 
0397   return superLayerName;
0398 }
0399 
0400 string DTNoiseComputation::getChamberName(const DTLayerId &lId) const {
0401   const DTSuperLayerId dtSLId = lId.superlayerId();
0402   const DTChamberId dtChId = dtSLId.chamberId();
0403   string chamberName = "W" + std::to_string(dtChId.wheel()) + "_St" + std::to_string(dtChId.station()) + "_Sec" +
0404                        std::to_string(dtChId.sector());
0405 
0406   return chamberName;
0407 }
0408 
0409 int DTNoiseComputation::getMaxNumBins(const DTChamberId &chId) const {
0410   int maximum = 0;
0411 
0412   for (int SL = 1; SL <= 3; SL++) {
0413     if (!(chId.station() == 4 && SL == 2)) {
0414       for (int L = 1; L <= 4; L++) {
0415         DTLayerId layerId = DTLayerId(chId, SL, L);
0416         string HistoName = "DigiOccupancy_" + getLayerName(layerId);
0417         theFile->cd();
0418         TH1F *hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0419         if (hOccHisto) {
0420           if (hOccHisto->GetXaxis()->GetXmax() > maximum)
0421             maximum = hOccHisto->GetXaxis()->GetNbins();
0422         }
0423       }
0424     }
0425   }
0426   return maximum;
0427 }
0428 
0429 double DTNoiseComputation::getYMaximum(const DTSuperLayerId &slId) const {
0430   double maximum = 0;
0431   double dummy = pow(10., 10.);
0432 
0433   for (int L = 1; L <= 4; L++) {
0434     DTLayerId layerId = DTLayerId(slId, L);
0435     string HistoName = "DigiOccupancy_" + getLayerName(layerId);
0436     theFile->cd();
0437     TH1F *hOccHisto = (TH1F *)theFile->Get(HistoName.c_str());
0438     if (hOccHisto) {
0439       if (hOccHisto->GetMaximum(dummy) > maximum)
0440         maximum = hOccHisto->GetMaximum(dummy);
0441     }
0442   }
0443   return maximum;
0444 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!