00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "CalibMuon/DTCalibration/plugins/DTNoiseComputation.h"
00011 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
00012
00013
00014 #include "FWCore/Framework/interface/IOVSyncValue.h"
00015 #include "FWCore/Framework/interface/Event.h"
00016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00017 #include "FWCore/Framework/interface/ESHandle.h"
00018 #include <FWCore/Framework/interface/MakerMacros.h>
00019
00020
00021 #include "Geometry/DTGeometry/interface/DTLayer.h"
00022 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00023 #include "Geometry/DTGeometry/interface/DTTopology.h"
00024 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00025
00026
00027 #include <DataFormats/DTDigi/interface/DTDigi.h>
00028 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
00029 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
00030 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
00031
00032 #include "TH1F.h"
00033 #include "TH2F.h"
00034 #include "TFile.h"
00035 #include "TF1.h"
00036 #include "TProfile.h"
00037 #include "TPostScript.h"
00038 #include "TCanvas.h"
00039 #include "TLegend.h"
00040
00041 using namespace edm;
00042 using namespace std;
00043
00044
00045 DTNoiseComputation::DTNoiseComputation(const edm::ParameterSet& ps){
00046
00047 cout << "[DTNoiseComputation]: Constructor" <<endl;
00048
00049
00050 debug = ps.getUntrackedParameter<bool>("debug");
00051
00052
00053 fastAnalysis = ps.getUntrackedParameter<bool>("fastAnalysis", true);
00054
00055
00056 string rootFileName = ps.getUntrackedParameter<string>("rootFileName");
00057 theFile = new TFile(rootFileName.c_str(), "READ");
00058
00059
00060 string newRootFileName = ps.getUntrackedParameter<string>("newRootFileName");
00061 theNewFile = new TFile(newRootFileName.c_str(), "RECREATE");
00062
00063
00064 MaxEvents = ps.getUntrackedParameter<int>("MaxEvents");
00065
00066 }
00067
00068 void DTNoiseComputation::beginRun(const edm::Run&, const EventSetup& setup)
00069 {
00070
00071 setup.get<MuonGeometryRecord>().get(dtGeom);
00072
00073 static int count = 0;
00074
00075 if(count == 0){
00076 string CheckHistoName;
00077
00078 TH1F *hOccHisto;
00079 TH1F *hAverageNoiseHisto;
00080 TH1F *hAverageNoiseIntegratedHisto;
00081 TH1F *hAverageNoiseHistoPerCh;
00082 TH1F *hAverageNoiseIntegratedHistoPerCh;
00083 TH2F *hEvtHisto;
00084 string HistoName;
00085 string Histo2Name;
00086 string AverageNoiseName;
00087 string AverageNoiseIntegratedName;
00088 string AverageNoiseNamePerCh;
00089 string AverageNoiseIntegratedNamePerCh;
00090 TH1F *hnoisyC;
00091 TH1F *hsomeHowNoisyC;
00092
00093
00094 vector<DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();
00095 vector<DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();
00096
00097 for (; ch_it != ch_end; ++ch_it) {
00098 DTChamberId ch = (*ch_it)->id();
00099 vector<const DTSuperLayer*>::const_iterator sl_it = (*ch_it)->superLayers().begin();
00100 vector<const DTSuperLayer*>::const_iterator sl_end = (*ch_it)->superLayers().end();
00101
00102 for(; sl_it != sl_end; ++sl_it) {
00103
00104 vector<const DTLayer*>::const_iterator l_it = (*sl_it)->layers().begin();
00105 vector<const DTLayer*>::const_iterator l_end = (*sl_it)->layers().end();
00106
00107 for(; l_it != l_end; ++l_it) {
00108 DTLayerId dtLId = (*l_it)->id();
00109
00110
00111 theFile->cd();
00112 CheckHistoName = "DigiOccupancy_" + getLayerName(dtLId);
00113 TH1F *hCheckHisto = (TH1F *) theFile->Get(CheckHistoName.c_str());
00114 if(hCheckHisto){
00115 delete hCheckHisto;
00116 stringstream wheel; wheel << ch.wheel();
00117 stringstream station; station << ch.station();
00118
00119 if(someHowNoisyC.find(make_pair(ch.wheel(),ch.station())) == someHowNoisyC.end()) {
00120 TString histoName_someHowNoisy = "somehowNoisyCell_W"+wheel.str()+"_St"+station.str();
00121 hsomeHowNoisyC = new TH1F(histoName_someHowNoisy,histoName_someHowNoisy,getMaxNumBins(ch),1,getMaxNumBins(ch)+1);
00122 someHowNoisyC[make_pair(ch.wheel(),ch.station())]=hsomeHowNoisyC;
00123 }
00124
00125 if(noisyC.find(make_pair(ch.wheel(),ch.station())) == noisyC.end()) {
00126 TString histoName_noisy = "noisyCell_W"+wheel.str()+"_St"+station.str();
00127 hnoisyC = new TH1F(histoName_noisy,histoName_noisy,getMaxNumBins(ch),1,getMaxNumBins(ch)+1);
00128 noisyC[make_pair(ch.wheel(),ch.station())]=hnoisyC;
00129 }
00130
00131
00132 if(AvNoisePerSuperLayer.find(dtLId.superlayerId()) == AvNoisePerSuperLayer.end()) {
00133 AverageNoiseName = "AverageNoise_" + getSuperLayerName(dtLId.superlayerId());
00134 hAverageNoiseHisto = new TH1F(AverageNoiseName.c_str(), AverageNoiseName.c_str(), 200, 0, 10000);
00135 AverageNoiseIntegratedName = "AverageNoiseIntegrated_" + getSuperLayerName(dtLId.superlayerId());
00136 hAverageNoiseIntegratedHisto = new TH1F(AverageNoiseIntegratedName.c_str(), AverageNoiseIntegratedName.c_str(), 200, 0, 10000);
00137 AvNoisePerSuperLayer[dtLId.superlayerId()] = hAverageNoiseHisto;
00138 AvNoiseIntegratedPerSuperLayer[dtLId.superlayerId()] = hAverageNoiseIntegratedHisto;
00139 if(debug){
00140 cout << " New Average Noise Histo per SuperLayer : " << hAverageNoiseHisto->GetName() << endl;
00141 cout << " New Average Noise Integrated Histo per SuperLayer : " << hAverageNoiseHisto->GetName() << endl;
00142 }
00143 }
00144 if(AvNoisePerChamber.find(dtLId.superlayerId().chamberId()) == AvNoisePerChamber.end()) {
00145 AverageNoiseNamePerCh = "AverageNoise_" + getChamberName(dtLId);
00146 hAverageNoiseHistoPerCh = new TH1F(AverageNoiseNamePerCh.c_str(), AverageNoiseNamePerCh.c_str(), 200, 0, 10000);
00147 AverageNoiseIntegratedNamePerCh = "AverageNoiseIntegrated_" + getChamberName(dtLId);
00148 hAverageNoiseIntegratedHistoPerCh = new TH1F(AverageNoiseIntegratedNamePerCh.c_str(), AverageNoiseIntegratedNamePerCh.c_str(), 200, 0, 10000);
00149 AvNoisePerChamber[dtLId.superlayerId().chamberId()] = hAverageNoiseHistoPerCh;
00150 AvNoiseIntegratedPerChamber[dtLId.superlayerId().chamberId()] = hAverageNoiseIntegratedHistoPerCh;
00151 if(debug)
00152 cout << " New Average Noise Histo per chamber : " << hAverageNoiseHistoPerCh->GetName() << endl;
00153 }
00154
00155 HistoName = "DigiOccupancy_" + getLayerName(dtLId);
00156 theFile->cd();
00157 hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
00158 int numBin = hOccHisto->GetXaxis()->GetNbins();
00159 for (int bin=1; bin<=numBin; bin++) {
00160 DTWireId wireID(dtLId, bin);
00161 theAverageNoise[wireID]= hOccHisto->GetBinContent(bin);
00162 if(theAverageNoise[wireID] != 0) {
00163 AvNoisePerSuperLayer[dtLId.superlayerId()]->Fill(theAverageNoise[wireID]);
00164 AvNoisePerChamber[dtLId.superlayerId().chamberId()]->Fill(theAverageNoise[wireID]);
00165 }
00166 }
00167
00168
00169 double numCell=0;
00170 double AvNoise=0;
00171 HistoName = "DigiOccupancy_" + getLayerName(dtLId);
00172 theFile->cd();
00173 hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
00174 numBin = hOccHisto->GetXaxis()->GetNbins();
00175 for (int bin=1; bin<=numBin; bin++) {
00176 DTWireId wireID(dtLId, bin);
00177 theAverageNoise[wireID]= hOccHisto->GetBinContent(bin);
00178 if(hOccHisto->GetBinContent(bin)<100){
00179 numCell++;
00180 AvNoise += hOccHisto->GetBinContent(bin);
00181 }
00182 if(hOccHisto->GetBinContent(bin)>100 && hOccHisto->GetBinContent(bin)<500){
00183 someHowNoisyC[make_pair(ch.wheel(),ch.station())]->Fill(bin);
00184 cout<<"filling somehow noisy cell"<<endl;
00185 }
00186 if(hOccHisto->GetBinContent(bin)>500){
00187 noisyC[make_pair(ch.wheel(),ch.station())]->Fill(bin);
00188 cout<<"filling noisy cell"<<endl;
00189 }
00190 }
00191 AvNoise = AvNoise/numCell;
00192 cout<<"theAverageNoise for layer "<<getLayerName(dtLId)<<" is : "<<AvNoise << endl;
00193
00194
00195
00196 int updates = MaxEvents/1000;
00197 for(int evt=0; evt<updates; evt++){
00198 stringstream toAppend; toAppend << evt;
00199 Histo2Name = "DigiPerWirePerEvent_" + getLayerName(dtLId) + "_" + toAppend.str();
00200 theFile->cd();
00201 hEvtHisto = (TH2F *) theFile->Get(Histo2Name.c_str());
00202 if(hEvtHisto){
00203 if(debug)
00204 cout << " New Histo with the number of events per evt per wire: " << hEvtHisto->GetName() << endl;
00205 theEvtMap[dtLId].push_back(hEvtHisto);
00206 }
00207 }
00208
00209 }
00210 }
00211 }
00212 }
00213
00214 count++;
00215 }
00216
00217 }
00218
00219 void DTNoiseComputation::endJob(){
00220
00221 cout << "[DTNoiseComputation] endjob called!" <<endl;
00222 TH1F *hEvtDistance=0;
00223 TF1 *ExpoFit = new TF1("ExpoFit","expo", 0.5, 1000.5);
00224 ExpoFit->SetMarkerColor();
00225 TF1 *funct=0;
00226 TProfile *theNoiseHisto = new TProfile("theNoiseHisto","Time Constant versus Average Noise",100000,0,100000);
00227
00228
00229
00230 for(map<DTLayerId, vector<TH2F*> >::const_iterator lHisto = theEvtMap.begin();
00231 lHisto != theEvtMap.end();
00232 lHisto++) {
00233 for(int bin=1; bin<(*lHisto).second[0]->GetYaxis()->GetNbins(); bin++){
00234 int distanceEvt = 1;
00235 DTWireId wire((*lHisto).first, bin);
00236 for(int i=0; i<int((*lHisto).second.size()); i++){
00237 for(int evt=1; evt<=(*lHisto).second[i]->GetXaxis()->GetNbins(); evt++){
00238 if((*lHisto).second[i]->GetBinContent(evt,bin) == 0) distanceEvt++;
00239 else {
00240 if(toDel.find(wire) == toDel.end()) {
00241 toDel[wire] = false;
00242 stringstream toAppend; toAppend << bin;
00243 string Histo = "EvtDistancePerWire_" + getLayerName((*lHisto).first) + "_" + toAppend.str();
00244 hEvtDistance = new TH1F(Histo.c_str(),Histo.c_str(), 50000,0.5,50000.5);
00245 }
00246 hEvtDistance->Fill(distanceEvt);
00247 distanceEvt=1;
00248 }
00249 }
00250 }
00251 if(toDel.find(wire) != toDel.end()){
00252 theHistoEvtDistancePerWire[wire] = hEvtDistance;
00253 theNewFile->cd();
00254 theHistoEvtDistancePerWire[wire]->Fit("ExpoFit","R");
00255 funct = theHistoEvtDistancePerWire[wire]->GetFunction("ExpoFit");
00256 double par0 = funct->GetParameter(0);
00257 double par1 = funct->GetParameter(1);
00258 cout<<"par0: "<<par0<<" par1: "<<par1<<endl;
00259 double chi2rid = funct->GetChisquare()/funct->GetNDF();
00260 if(chi2rid>10)
00261 theTimeConstant[wire]=1;
00262 else
00263 theTimeConstant[wire]=par1;
00264 toDel[wire] = true;
00265 theHistoEvtDistancePerWire[wire]->Write();
00266 delete hEvtDistance;
00267 }
00268 }
00269 }
00270
00271 if(!fastAnalysis){
00272
00273 for(map<DTWireId, double>::const_iterator AvNoise = theAverageNoise.begin();
00274 AvNoise != theAverageNoise.end();
00275 AvNoise++) {
00276 DTWireId wire = (*AvNoise).first;
00277 theNoiseHisto->Fill((*AvNoise).second, theTimeConstant[wire]);
00278 cout<<"Layer: "<<getLayerName(wire.layerId())<<" wire: "<<wire.wire()<<endl;
00279 cout<<"The Average noise: "<<(*AvNoise).second<<endl;
00280 cout<<"The time constant: "<<theTimeConstant[wire]<<endl;
00281 }
00282 theNewFile->cd();
00283 theNoiseHisto->Write();
00284 }
00285
00286
00287
00288 int numBin;
00289 double integratedNoise, bin, halfBin, maxBin;
00290 for(map<DTSuperLayerId, TH1F*>::const_iterator AvNoiseHisto = AvNoisePerSuperLayer.begin();
00291 AvNoiseHisto != AvNoisePerSuperLayer.end();
00292 AvNoiseHisto++) {
00293 integratedNoise=0;
00294 numBin = (*AvNoiseHisto).second->GetXaxis()->GetNbins();
00295 maxBin = (*AvNoiseHisto).second->GetXaxis()->GetXmax();
00296 bin= double(maxBin/numBin);
00297 halfBin=double(bin/2);
00298 theNewFile->cd();
00299 (*AvNoiseHisto).second->Write();
00300 for(int i=1; i<numBin; i++){
00301 integratedNoise+=(*AvNoiseHisto).second->GetBinContent(i);
00302 AvNoiseIntegratedPerSuperLayer[(*AvNoiseHisto).first]->Fill(halfBin,integratedNoise);
00303 halfBin+=bin;
00304 }
00305 theNewFile->cd();
00306 AvNoiseIntegratedPerSuperLayer[(*AvNoiseHisto).first]->Write();
00307 }
00308
00309 for(map<DTChamberId, TH1F*>::const_iterator AvNoiseHisto = AvNoisePerChamber.begin();
00310 AvNoiseHisto != AvNoisePerChamber.end();
00311 AvNoiseHisto++) {
00312 integratedNoise=0;
00313 numBin = (*AvNoiseHisto).second->GetXaxis()->GetNbins();
00314 maxBin = (*AvNoiseHisto).second->GetXaxis()->GetXmax();
00315 bin= maxBin/numBin;
00316 halfBin=bin/2;
00317 theNewFile->cd();
00318 (*AvNoiseHisto).second->Write();
00319 for(int i=1; i<numBin; i++){
00320 integratedNoise+=(*AvNoiseHisto).second->GetBinContent(i);
00321 AvNoiseIntegratedPerChamber[(*AvNoiseHisto).first]->Fill(halfBin,integratedNoise);
00322 halfBin+=bin;
00323 }
00324 theNewFile->cd();
00325 AvNoiseIntegratedPerChamber[(*AvNoiseHisto).first]->Write();
00326 }
00327
00328
00329
00330 bool histo=false;
00331 vector<DTChamber*>::const_iterator chamber_it = dtGeom->chambers().begin();
00332 vector<DTChamber*>::const_iterator chamber_end = dtGeom->chambers().end();
00333
00334 for (; chamber_it != chamber_end; ++chamber_it) {
00335 DTChamberId ch = (*chamber_it)->id();
00336 vector<const DTSuperLayer*>::const_iterator sl_it = (*chamber_it)->superLayers().begin();
00337 vector<const DTSuperLayer*>::const_iterator sl_end = (*chamber_it)->superLayers().end();
00338
00339 for(; sl_it != sl_end; ++sl_it) {
00340 DTSuperLayerId sl = (*sl_it)->id();
00341 vector<const DTLayer*>::const_iterator l_it = (*sl_it)->layers().begin();
00342 vector<const DTLayer*>::const_iterator l_end = (*sl_it)->layers().end();
00343
00344 string canvasName = "c" + getSuperLayerName(sl);
00345 TCanvas c1(canvasName.c_str(),canvasName.c_str(),600,780);
00346 TLegend *leg=new TLegend(0.5,0.6,0.7,0.8);
00347 for(; l_it != l_end; ++l_it) {
00348 DTLayerId layerId = (*l_it)->id();
00349 string HistoName = "DigiOccupancy_" + getLayerName(layerId);
00350 theFile->cd();
00351 TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
00352 if(hOccHisto){
00353 string TitleHisto = "AverageNoise_" + getSuperLayerName(sl);
00354 cout<<"TitleHisto : "<<TitleHisto<<endl;
00355 hOccHisto->SetTitle(TitleHisto.c_str());
00356 stringstream layer; layer << layerId.layer();
00357 string legendHisto = "layer " + layer.str();
00358 leg->AddEntry(hOccHisto,legendHisto.c_str(),"L");
00359 hOccHisto->SetMaximum(getYMaximum(sl));
00360 histo=true;
00361 if(layerId.layer() == 1)
00362 hOccHisto->Draw();
00363 else
00364 hOccHisto->Draw("same");
00365 hOccHisto->SetLineColor(layerId.layer());
00366 }
00367 }
00368 if(histo){
00369 leg->Draw("same");
00370 theNewFile->cd();
00371 c1.Write();
00372 }
00373 }
00374 histo=false;
00375 }
00376
00377
00378 for(map<pair<int,int>, TH1F*>::const_iterator nCell = noisyC.begin();
00379 nCell != noisyC.end();
00380 nCell++) {
00381 theNewFile->cd();
00382 (*nCell).second->Write();
00383 }
00384 for(map<pair<int,int>, TH1F*>::const_iterator somehownCell = someHowNoisyC.begin();
00385 somehownCell != someHowNoisyC.end();
00386 somehownCell++) {
00387 theNewFile->cd();
00388 (*somehownCell).second->Write();
00389 }
00390
00391 }
00392
00393
00394 DTNoiseComputation::~DTNoiseComputation(){
00395
00396 theFile->Close();
00397 theNewFile->Close();
00398
00399 }
00400
00401
00402 string DTNoiseComputation::getLayerName(const DTLayerId& lId) const {
00403
00404 const DTSuperLayerId dtSLId = lId.superlayerId();
00405 const DTChamberId dtChId = dtSLId.chamberId();
00406 stringstream Layer; Layer << lId.layer();
00407 stringstream superLayer; superLayer << dtSLId.superlayer();
00408 stringstream wheel; wheel << dtChId.wheel();
00409 stringstream station; station << dtChId.station();
00410 stringstream sector; sector << dtChId.sector();
00411
00412 string LayerName =
00413 "W" + wheel.str()
00414 + "_St" + station.str()
00415 + "_Sec" + sector.str()
00416 + "_SL" + superLayer.str()
00417 + "_L" + Layer.str();
00418
00419 return LayerName;
00420
00421 }
00422
00423
00424 string DTNoiseComputation::getSuperLayerName(const DTSuperLayerId& dtSLId) const {
00425
00426 const DTChamberId dtChId = dtSLId.chamberId();
00427 stringstream superLayer; superLayer << dtSLId.superlayer();
00428 stringstream wheel; wheel << dtChId.wheel();
00429 stringstream station; station << dtChId.station();
00430 stringstream sector; sector << dtChId.sector();
00431
00432 string SuperLayerName =
00433 "W" + wheel.str()
00434 + "_St" + station.str()
00435 + "_Sec" + sector.str()
00436 + "_SL" + superLayer.str();
00437
00438 return SuperLayerName;
00439
00440 }
00441
00442
00443 string DTNoiseComputation::getChamberName(const DTLayerId& lId) const {
00444
00445 const DTSuperLayerId dtSLId = lId.superlayerId();
00446 const DTChamberId dtChId = dtSLId.chamberId();
00447 stringstream wheel; wheel << dtChId.wheel();
00448 stringstream station; station << dtChId.station();
00449 stringstream sector; sector << dtChId.sector();
00450
00451 string ChamberName =
00452 "W" + wheel.str()
00453 + "_St" + station.str()
00454 + "_Sec" + sector.str();
00455
00456 return ChamberName;
00457
00458 }
00459
00460
00461 int DTNoiseComputation::getMaxNumBins(const DTChamberId& chId) const {
00462
00463 int maximum=0;
00464
00465 for(int SL=1; SL<=3; SL++){
00466 if(!(chId.station()==4 && SL==2)){
00467 for (int L=1; L<=4; L++){
00468 DTLayerId layerId = DTLayerId(chId, SL, L);
00469 string HistoName = "DigiOccupancy_" + getLayerName(layerId);
00470 theFile->cd();
00471 TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
00472 if(hOccHisto){
00473 if (hOccHisto->GetXaxis()->GetXmax()>maximum)
00474 maximum = hOccHisto->GetXaxis()->GetNbins();
00475 }
00476 }
00477 }
00478 }
00479 return maximum;
00480 }
00481
00482
00483 double DTNoiseComputation::getYMaximum(const DTSuperLayerId& slId) const {
00484
00485 double maximum=0;
00486 double dummy = pow(10.,10.);
00487
00488 for (int L=1; L<=4; L++){
00489 DTLayerId layerId = DTLayerId(slId, L);
00490 string HistoName = "DigiOccupancy_" + getLayerName(layerId);
00491 theFile->cd();
00492 TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
00493 if(hOccHisto){
00494 if (hOccHisto->GetMaximum(dummy)>maximum)
00495 maximum = hOccHisto->GetMaximum(dummy);
00496 }
00497 }
00498 return maximum;
00499 }
00500
00501
00502