#include <DTNoiseComputation.h>

Inheritance diagram for DTNoiseComputation:

Public Member Functions
void	analyze (const edm::Event &event, const edm::EventSetup &setup)
void	beginJob ()
	BeginJob.
void	beginRun (const edm::Run &, const edm::EventSetup &setup)
	DTNoiseComputation (const edm::ParameterSet &ps)
	Constructor.
void	endJob ()
	Endjob.
virtual	~DTNoiseComputation ()
	Destructor.
Private Member Functions
std::string	getChamberName (const DTLayerId &lId) const
	Get the name of the chamber.
std::string	getLayerName (const DTLayerId &lId) const
	Get the name of the layer.
int	getMaxNumBins (const DTChamberId &chId) const
std::string	getSuperLayerName (const DTSuperLayerId &slId) const
	Get the name of the superLayer.
double	getYMaximum (const DTSuperLayerId &slId) const
Private Attributes
std::map< DTChamberId, TH1F * >	AvNoiseIntegratedPerChamber
std::map< DTSuperLayerId, TH1F * >	AvNoiseIntegratedPerSuperLayer
std::map< DTChamberId, TH1F * >	AvNoisePerChamber
std::map< DTSuperLayerId, TH1F * >	AvNoisePerSuperLayer
int	counter
bool	debug
edm::ESHandle< DTGeometry >	dtGeom
bool	fastAnalysis
int	MaxEvents
std::map< std::pair< int, int > , TH1F * >	noisyC
std::map< std::pair< int, int > , TH1F * >	someHowNoisyC
std::map< DTWireId, double >	theAverageNoise
std::map< DTLayerId, std::vector< TH2F * > >	theEvtMap
TFile *	theFile
std::map< DTWireId, TH1F * >	theHistoEvtDistancePerWire
TFile *	theNewFile
std::map< DTWireId, double >	theTimeConstant
std::map< DTLayerId, bool >	toComputeNoiseAverage
std::map< DTWireId, bool >	toDel

Detailed Description

Definition at line 37 of file DTNoiseComputation.h.

Constructor & Destructor Documentation

DTNoiseComputation::DTNoiseComputation ( const edm::ParameterSet & ps )

Constructor.

Definition at line 45 of file DTNoiseComputation.cc.

References gather_cfg::cout, debug, edm::ParameterSet::getUntrackedParameter(), dtT0WireCalibration_cfg::rootFileName, and interactiveExample::theFile.

                                                               {

  cout << "[DTNoiseComputation]: Constructor" <<endl;

  // Get the debug parameter for verbose output
  debug = ps.getUntrackedParameter<bool>("debug");

  // The analysis type
  fastAnalysis = ps.getUntrackedParameter<bool>("fastAnalysis", true);

  // The root file which contain the histos
  string rootFileName = ps.getUntrackedParameter<string>("rootFileName");
  theFile = new TFile(rootFileName.c_str(), "READ");

  // The new root file which contain the histos
  string newRootFileName = ps.getUntrackedParameter<string>("newRootFileName");
  theNewFile = new TFile(newRootFileName.c_str(), "RECREATE");

  // The maximum number of events to analyze
  MaxEvents = ps.getUntrackedParameter<int>("MaxEvents");

}

DTNoiseComputation::~DTNoiseComputation ( ) [virtual]

Destructor.

Definition at line 394 of file DTNoiseComputation.cc.

References interactiveExample::theFile.

                                       {

  theFile->Close();
  theNewFile->Close();

}

Member Function Documentation

void DTNoiseComputation::analyze	(	const edm::Event &	event,
		const edm::EventSetup &	setup
	)		`[inline, virtual]`

Implements edm::EDAnalyzer.

Definition at line 52 of file DTNoiseComputation.h.

{}

void DTNoiseComputation::beginJob ( void ) [inline, virtual]

BeginJob.

Reimplemented from edm::EDAnalyzer.

Definition at line 48 of file DTNoiseComputation.h.

{}

void DTNoiseComputation::beginRun	(	const edm::Run &	,
		const edm::EventSetup &	setup
	)		`[virtual]`

Reimplemented from edm::EDAnalyzer.

Definition at line 68 of file DTNoiseComputation.cc.

References newFWLiteAna::bin, DTSuperLayerId::chamberId(), prof2calltree::count, gather_cfg::cout, debug, HcalObjRepresent::Fill(), edm::EventSetup::get(), relativeConstraints::station, DTChamberId::station(), DTLayerId::superlayerId(), interactiveExample::theFile, and DTChamberId::wheel().

{
  // Get the DT Geometry
  setup.get<MuonGeometryRecord>().get(dtGeom);

  static int count = 0;

  if(count == 0){
    string CheckHistoName;
  
    TH1F *hOccHisto;
    TH1F *hAverageNoiseHisto;
    TH1F *hAverageNoiseIntegratedHisto;
    TH1F *hAverageNoiseHistoPerCh;
    TH1F *hAverageNoiseIntegratedHistoPerCh;
    TH2F *hEvtHisto;
    string HistoName;
    string Histo2Name;
    string AverageNoiseName;
    string AverageNoiseIntegratedName;
    string AverageNoiseNamePerCh;
    string AverageNoiseIntegratedNamePerCh;
    TH1F *hnoisyC;
    TH1F *hsomeHowNoisyC;
  
    // Loop over all the chambers        
    vector<DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();       
    vector<DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();        
    // Loop over the SLs         
    for (; ch_it != ch_end; ++ch_it) {   
      DTChamberId ch = (*ch_it)->id();   
      vector<const DTSuperLayer*>::const_iterator sl_it = (*ch_it)->superLayers().begin();       
      vector<const DTSuperLayer*>::const_iterator sl_end = (*ch_it)->superLayers().end();        
      // Loop over the SLs       
      for(; sl_it != sl_end; ++sl_it) {          
        //      DTSuperLayerId sl = (*sl_it)->id();      
        vector<const DTLayer*>::const_iterator l_it = (*sl_it)->layers().begin();        
        vector<const DTLayer*>::const_iterator l_end = (*sl_it)->layers().end();         
        // Loop over the Ls      
        for(; l_it != l_end; ++l_it) {   
          DTLayerId dtLId = (*l_it)->id();
        
          //check if the layer has digi
          theFile->cd();
          CheckHistoName =  "DigiOccupancy_" + getLayerName(dtLId);
          TH1F *hCheckHisto = (TH1F *) theFile->Get(CheckHistoName.c_str());
          if(hCheckHisto){  
            delete hCheckHisto;
            stringstream wheel; wheel << ch.wheel();    
            stringstream station; station << ch.station();
          
            if(someHowNoisyC.find(make_pair(ch.wheel(),ch.station())) == someHowNoisyC.end()) {
              TString histoName_someHowNoisy = "somehowNoisyCell_W"+wheel.str()+"_St"+station.str();
              hsomeHowNoisyC = new TH1F(histoName_someHowNoisy,histoName_someHowNoisy,getMaxNumBins(ch),1,getMaxNumBins(ch)+1);
              someHowNoisyC[make_pair(ch.wheel(),ch.station())]=hsomeHowNoisyC;
            }
          
            if(noisyC.find(make_pair(ch.wheel(),ch.station())) == noisyC.end()) {
              TString histoName_noisy = "noisyCell_W"+wheel.str()+"_St"+station.str();
              hnoisyC = new TH1F(histoName_noisy,histoName_noisy,getMaxNumBins(ch),1,getMaxNumBins(ch)+1);
              noisyC[make_pair(ch.wheel(),ch.station())]=hnoisyC;
            }
          
            //to fill a map with the average noise per wire and fill new noise histo      
            if(AvNoisePerSuperLayer.find(dtLId.superlayerId()) == AvNoisePerSuperLayer.end()) {
              AverageNoiseName = "AverageNoise_" + getSuperLayerName(dtLId.superlayerId());
              hAverageNoiseHisto = new TH1F(AverageNoiseName.c_str(), AverageNoiseName.c_str(), 200, 0, 10000);
              AverageNoiseIntegratedName = "AverageNoiseIntegrated_" + getSuperLayerName(dtLId.superlayerId());
              hAverageNoiseIntegratedHisto = new TH1F(AverageNoiseIntegratedName.c_str(), AverageNoiseIntegratedName.c_str(), 200, 0, 10000);
              AvNoisePerSuperLayer[dtLId.superlayerId()] = hAverageNoiseHisto;
              AvNoiseIntegratedPerSuperLayer[dtLId.superlayerId()] = hAverageNoiseIntegratedHisto;
              if(debug){
                cout << "  New Average Noise Histo per SuperLayer : " << hAverageNoiseHisto->GetName() << endl;
                cout << "  New Average Noise Integrated Histo per SuperLayer : " << hAverageNoiseHisto->GetName() << endl;
              }
            }
            if(AvNoisePerChamber.find(dtLId.superlayerId().chamberId()) == AvNoisePerChamber.end()) {
              AverageNoiseNamePerCh = "AverageNoise_" + getChamberName(dtLId);
              hAverageNoiseHistoPerCh = new TH1F(AverageNoiseNamePerCh.c_str(), AverageNoiseNamePerCh.c_str(), 200, 0, 10000);
              AverageNoiseIntegratedNamePerCh = "AverageNoiseIntegrated_" + getChamberName(dtLId);
              hAverageNoiseIntegratedHistoPerCh = new TH1F(AverageNoiseIntegratedNamePerCh.c_str(), AverageNoiseIntegratedNamePerCh.c_str(), 200, 0, 10000);
              AvNoisePerChamber[dtLId.superlayerId().chamberId()] = hAverageNoiseHistoPerCh;
              AvNoiseIntegratedPerChamber[dtLId.superlayerId().chamberId()] = hAverageNoiseIntegratedHistoPerCh;
              if(debug)
                cout << "  New Average Noise Histo per chamber : " << hAverageNoiseHistoPerCh->GetName() << endl;
            }
          
            HistoName = "DigiOccupancy_" + getLayerName(dtLId);
            theFile->cd();
            hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
            int numBin = hOccHisto->GetXaxis()->GetNbins(); 
            for (int bin=1; bin<=numBin; bin++) {
              DTWireId wireID(dtLId, bin);
              theAverageNoise[wireID]= hOccHisto->GetBinContent(bin);
              if(theAverageNoise[wireID] != 0) {
                AvNoisePerSuperLayer[dtLId.superlayerId()]->Fill(theAverageNoise[wireID]);
                AvNoisePerChamber[dtLId.superlayerId().chamberId()]->Fill(theAverageNoise[wireID]);
              }
            }         
          
            //to compute the average noise per layer (excluding the noisy cells)
            double numCell=0;
            double AvNoise=0;        
            HistoName = "DigiOccupancy_" + getLayerName(dtLId);
            theFile->cd();
            hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
            numBin = hOccHisto->GetXaxis()->GetNbins(); 
            for (int bin=1; bin<=numBin; bin++) {
              DTWireId wireID(dtLId, bin);
              theAverageNoise[wireID]= hOccHisto->GetBinContent(bin);
              if(hOccHisto->GetBinContent(bin)<100){
                numCell++;
                AvNoise += hOccHisto->GetBinContent(bin);
              }
              if(hOccHisto->GetBinContent(bin)>100 && hOccHisto->GetBinContent(bin)<500){
                someHowNoisyC[make_pair(ch.wheel(),ch.station())]->Fill(bin);
                cout<<"filling somehow noisy cell"<<endl;
              }
              if(hOccHisto->GetBinContent(bin)>500){
                noisyC[make_pair(ch.wheel(),ch.station())]->Fill(bin);
                cout<<"filling noisy cell"<<endl;
              }
            }
            AvNoise = AvNoise/numCell;
            cout<<"theAverageNoise for layer "<<getLayerName(dtLId)<<" is : "<<AvNoise << endl;

          
            // book the digi event plots every 1000 events
            int updates = MaxEvents/1000; 
            for(int evt=0; evt<updates; evt++){
              stringstream toAppend; toAppend << evt;
              Histo2Name = "DigiPerWirePerEvent_" + getLayerName(dtLId) + "_" + toAppend.str();
              theFile->cd();
              hEvtHisto = (TH2F *) theFile->Get(Histo2Name.c_str());
              if(hEvtHisto){
                if(debug)
                  cout << "  New Histo with the number of events per evt per wire: " << hEvtHisto->GetName() << endl;
                theEvtMap[dtLId].push_back(hEvtHisto);
              }
            }
          
          }// done if the layer has digi
        }// loop over layers
      }// loop over superlayers
    }// loop over chambers

    count++;
  }

}

void DTNoiseComputation::endJob ( void ) [virtual]

Endjob.

Reimplemented from edm::EDAnalyzer.

Definition at line 219 of file DTNoiseComputation.cc.

References newFWLiteAna::bin, alignmentValidation::c1, gather_cfg::cout, trackerHits::histo, i, DTLayerId::layer(), fitWZ::par0, and interactiveExample::theFile.

                               {

  cout << "[DTNoiseComputation] endjob called!" <<endl;
  TH1F *hEvtDistance=0;
  TF1 *ExpoFit = new TF1("ExpoFit","expo", 0.5, 1000.5);
  ExpoFit->SetMarkerColor();//just silence gcc complaining about unused vars
  TF1 *funct=0;
  TProfile *theNoiseHisto = new TProfile("theNoiseHisto","Time Constant versus Average Noise",100000,0,100000);
  

  //compute the time constant
  for(map<DTLayerId, vector<TH2F*> >::const_iterator lHisto = theEvtMap.begin();
      lHisto != theEvtMap.end();
      lHisto++) {
    for(int bin=1; bin<(*lHisto).second[0]->GetYaxis()->GetNbins(); bin++){
      int distanceEvt = 1;
      DTWireId wire((*lHisto).first, bin);
      for(int i=0; i<int((*lHisto).second.size()); i++){
        for(int evt=1; evt<=(*lHisto).second[i]->GetXaxis()->GetNbins(); evt++){
          if((*lHisto).second[i]->GetBinContent(evt,bin) == 0) distanceEvt++;
          else { 
            if(toDel.find(wire) == toDel.end()) {
              toDel[wire] = false;
              stringstream toAppend; toAppend << bin;
              string Histo = "EvtDistancePerWire_" + getLayerName((*lHisto).first) + "_" + toAppend.str();
              hEvtDistance = new TH1F(Histo.c_str(),Histo.c_str(), 50000,0.5,50000.5);
            }
            hEvtDistance->Fill(distanceEvt); 
            distanceEvt=1;
          }
        }
      }
      if(toDel.find(wire) != toDel.end()){
        theHistoEvtDistancePerWire[wire] =  hEvtDistance;
        theNewFile->cd();
        theHistoEvtDistancePerWire[wire]->Fit("ExpoFit","R");
        funct = theHistoEvtDistancePerWire[wire]->GetFunction("ExpoFit");
        double par0 = funct->GetParameter(0);
        double par1 = funct->GetParameter(1);
        cout<<"par0: "<<par0<<"  par1: "<<par1<<endl;
        double chi2rid = funct->GetChisquare()/funct->GetNDF();
        if(chi2rid>10)
          theTimeConstant[wire]=1;
        else
          theTimeConstant[wire]=par1;
        toDel[wire] = true;
        theHistoEvtDistancePerWire[wire]->Write();
        delete hEvtDistance;
      }
    }
  }

  if(!fastAnalysis){
    //fill the histo with the time constant as a function of the average noise
    for(map<DTWireId, double>::const_iterator AvNoise = theAverageNoise.begin();
        AvNoise != theAverageNoise.end();
        AvNoise++) {
      DTWireId wire = (*AvNoise).first;
      theNoiseHisto->Fill((*AvNoise).second, theTimeConstant[wire]);
      cout<<"Layer: "<<getLayerName(wire.layerId())<<"  wire: "<<wire.wire()<<endl;
      cout<<"The Average noise: "<<(*AvNoise).second<<endl;
      cout<<"The time constant: "<<theTimeConstant[wire]<<endl;
    }
    theNewFile->cd();
    theNoiseHisto->Write();
  }  


  // histos with the integrated noise per layer
  int numBin;
  double integratedNoise, bin, halfBin, maxBin;
  for(map<DTSuperLayerId, TH1F*>::const_iterator AvNoiseHisto = AvNoisePerSuperLayer.begin();
      AvNoiseHisto != AvNoisePerSuperLayer.end();
      AvNoiseHisto++) {
    integratedNoise=0;
    numBin = (*AvNoiseHisto).second->GetXaxis()->GetNbins();
    maxBin = (*AvNoiseHisto).second->GetXaxis()->GetXmax();
    bin= double(maxBin/numBin);
    halfBin=double(bin/2);
    theNewFile->cd();
    (*AvNoiseHisto).second->Write();
    for(int i=1; i<numBin; i++){
      integratedNoise+=(*AvNoiseHisto).second->GetBinContent(i);
      AvNoiseIntegratedPerSuperLayer[(*AvNoiseHisto).first]->Fill(halfBin,integratedNoise);
      halfBin+=bin;
    }
    theNewFile->cd();
    AvNoiseIntegratedPerSuperLayer[(*AvNoiseHisto).first]->Write(); 
  }
  // histos with the integrated noise per chamber
  for(map<DTChamberId, TH1F*>::const_iterator AvNoiseHisto = AvNoisePerChamber.begin();
      AvNoiseHisto != AvNoisePerChamber.end();
      AvNoiseHisto++) {
    integratedNoise=0;
    numBin = (*AvNoiseHisto).second->GetXaxis()->GetNbins();
    maxBin = (*AvNoiseHisto).second->GetXaxis()->GetXmax();
    bin= maxBin/numBin;
    halfBin=bin/2;
    theNewFile->cd();
    (*AvNoiseHisto).second->Write();
    for(int i=1; i<numBin; i++){
      integratedNoise+=(*AvNoiseHisto).second->GetBinContent(i);
      AvNoiseIntegratedPerChamber[(*AvNoiseHisto).first]->Fill(halfBin,integratedNoise);
      halfBin+=bin;
    } 
    theNewFile->cd();
    AvNoiseIntegratedPerChamber[(*AvNoiseHisto).first]->Write(); 
  }

  
  //overimpose the average noise histo
  bool histo=false;
  vector<DTChamber*>::const_iterator chamber_it = dtGeom->chambers().begin();
  vector<DTChamber*>::const_iterator chamber_end = dtGeom->chambers().end();
  // Loop over the chambers
  for (; chamber_it != chamber_end; ++chamber_it) {
    DTChamberId ch = (*chamber_it)->id();
    vector<const DTSuperLayer*>::const_iterator sl_it = (*chamber_it)->superLayers().begin(); 
    vector<const DTSuperLayer*>::const_iterator sl_end = (*chamber_it)->superLayers().end();
    // Loop over the SLs
    for(; sl_it != sl_end; ++sl_it) {
      DTSuperLayerId sl = (*sl_it)->id();
      vector<const DTLayer*>::const_iterator l_it = (*sl_it)->layers().begin(); 
      vector<const DTLayer*>::const_iterator l_end = (*sl_it)->layers().end();

      string canvasName = "c" + getSuperLayerName(sl); 
      TCanvas c1(canvasName.c_str(),canvasName.c_str(),600,780);
      TLegend *leg=new TLegend(0.5,0.6,0.7,0.8);
      for(; l_it != l_end; ++l_it) {
        DTLayerId layerId = (*l_it)->id();
        string HistoName = "DigiOccupancy_" + getLayerName(layerId);
        theFile->cd();
        TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
        if(hOccHisto){
          string TitleHisto = "AverageNoise_" + getSuperLayerName(sl);
          cout<<"TitleHisto : "<<TitleHisto<<endl;
          hOccHisto->SetTitle(TitleHisto.c_str());
          stringstream layer; layer << layerId.layer(); 
          string legendHisto = "layer " + layer.str();
          leg->AddEntry(hOccHisto,legendHisto.c_str(),"L");
          hOccHisto->SetMaximum(getYMaximum(sl));
          histo=true;
          if(layerId.layer() == 1)
            hOccHisto->Draw();
          else
            hOccHisto->Draw("same");
          hOccHisto->SetLineColor(layerId.layer());
        }
      }
      if(histo){
        leg->Draw("same");
        theNewFile->cd();
        c1.Write();
      }
    }
    histo=false;   
  }

  //write on file the noisy plots
  for(map<pair<int,int>, TH1F*>::const_iterator nCell = noisyC.begin();
      nCell != noisyC.end();
      nCell++) {
    theNewFile->cd();
    (*nCell).second->Write();
  }
  for(map<pair<int,int>, TH1F*>::const_iterator somehownCell = someHowNoisyC.begin();
      somehownCell != someHowNoisyC.end();
      somehownCell++) {
    theNewFile->cd();
    (*somehownCell).second->Write();
  }
    
}

string DTNoiseComputation::getChamberName ( const DTLayerId & lId ) const [private]

Get the name of the chamber.

Definition at line 443 of file DTNoiseComputation.cc.

References DTSuperLayerId::chamberId(), DTChamberId::sector(), relativeConstraints::station, DTChamberId::station(), DTLayerId::superlayerId(), and DTChamberId::wheel().

                                                                    {

  const  DTSuperLayerId dtSLId = lId.superlayerId();
  const  DTChamberId dtChId = dtSLId.chamberId(); 
  stringstream wheel; wheel << dtChId.wheel();  
  stringstream station; station << dtChId.station();    
  stringstream sector; sector << dtChId.sector();
  
  string ChamberName = 
    "W" + wheel.str()
    + "_St" + station.str() 
    + "_Sec" + sector.str();
  
  return ChamberName;

}

string DTNoiseComputation::getLayerName ( const DTLayerId & lId ) const [private]

Get the name of the layer.

Definition at line 402 of file DTNoiseComputation.cc.

References DTSuperLayerId::chamberId(), DTLayerId::layer(), DTChamberId::sector(), relativeConstraints::station, DTChamberId::station(), DTSuperLayerId::superlayer(), DTLayerId::superlayerId(), and DTChamberId::wheel().

                                                                  {

  const  DTSuperLayerId dtSLId = lId.superlayerId();
  const  DTChamberId dtChId = dtSLId.chamberId(); 
  stringstream Layer; Layer << lId.layer();
  stringstream superLayer; superLayer << dtSLId.superlayer();
  stringstream wheel; wheel << dtChId.wheel();  
  stringstream station; station << dtChId.station();    
  stringstream sector; sector << dtChId.sector();
  
  string LayerName = 
    "W" + wheel.str()
    + "_St" + station.str() 
    + "_Sec" + sector.str() 
    + "_SL" + superLayer.str()
    + "_L" + Layer.str();
  
  return LayerName;

}

int DTNoiseComputation::getMaxNumBins ( const DTChamberId & chId ) const [private]

Definition at line 461 of file DTNoiseComputation.cc.

References DTLayerId, dttmaxenums::L, DTChamberId::station(), and interactiveExample::theFile.

                                                                   {
  
  int maximum=0;
  
  for(int SL=1; SL<=3; SL++){
    if(!(chId.station()==4 && SL==2)){
      for (int L=1; L<=4; L++){ 
        DTLayerId layerId = DTLayerId(chId, SL, L);
        string HistoName = "DigiOccupancy_" + getLayerName(layerId);
        theFile->cd();
        TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
        if(hOccHisto){
          if (hOccHisto->GetXaxis()->GetXmax()>maximum) 
            maximum = hOccHisto->GetXaxis()->GetNbins();
        }
      }
    }
  }
  return maximum;
}

string DTNoiseComputation::getSuperLayerName ( const DTSuperLayerId & slId ) const [private]

Get the name of the superLayer.

Definition at line 424 of file DTNoiseComputation.cc.

References DTSuperLayerId::chamberId(), DTChamberId::sector(), relativeConstraints::station, DTChamberId::station(), DTSuperLayerId::superlayer(), and DTChamberId::wheel().

                                                                               {

  const  DTChamberId dtChId = dtSLId.chamberId(); 
  stringstream superLayer; superLayer << dtSLId.superlayer();
  stringstream wheel; wheel << dtChId.wheel();  
  stringstream station; station << dtChId.station();    
  stringstream sector; sector << dtChId.sector();
  
  string SuperLayerName = 
    "W" + wheel.str()
    + "_St" + station.str() 
    + "_Sec" + sector.str() 
    + "_SL" + superLayer.str();
  
  return SuperLayerName;

}

double DTNoiseComputation::getYMaximum ( const DTSuperLayerId & slId ) const [private]

Definition at line 483 of file DTNoiseComputation.cc.

References DTLayerId, dttmaxenums::L, funct::pow(), and interactiveExample::theFile.

                                                                       {
  
  double maximum=0;
  double dummy = pow(10.,10.);

  for (int L=1; L<=4; L++){
    DTLayerId layerId = DTLayerId(slId, L);
    string HistoName = "DigiOccupancy_" + getLayerName(layerId);
    theFile->cd();
    TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
    if(hOccHisto){
      if (hOccHisto->GetMaximum(dummy)>maximum)
        maximum = hOccHisto->GetMaximum(dummy);
    }
  }
  return maximum;
}