CMSSW Reference Manual

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