CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/CalibMuon/DTCalibration/plugins/DTNoiseComputation.cc

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2012/01/12 15:05:33 $
00005  *  $Revision: 1.7 $
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/MakerMacros.h>
00019 
00020 // Geometry
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 // Digis
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   // Get the debug parameter for verbose output
00050   debug = ps.getUntrackedParameter<bool>("debug");
00051 
00052   // The analysis type
00053   fastAnalysis = ps.getUntrackedParameter<bool>("fastAnalysis", true);
00054 
00055   // The root file which contain the histos
00056   string rootFileName = ps.getUntrackedParameter<string>("rootFileName");
00057   theFile = new TFile(rootFileName.c_str(), "READ");
00058 
00059   // The new root file which contain the histos
00060   string newRootFileName = ps.getUntrackedParameter<string>("newRootFileName");
00061   theNewFile = new TFile(newRootFileName.c_str(), "RECREATE");
00062 
00063   // The maximum number of events to analyze
00064   MaxEvents = ps.getUntrackedParameter<int>("MaxEvents");
00065 
00066 }
00067 
00068 void DTNoiseComputation::beginRun(const edm::Run&, const EventSetup& setup)
00069 {
00070   // Get the DT Geometry
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     // Loop over all the chambers        
00094     vector<DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();       
00095     vector<DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();        
00096     // Loop over the SLs         
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       // Loop over the SLs       
00102       for(; sl_it != sl_end; ++sl_it) {          
00103         //      DTSuperLayerId sl = (*sl_it)->id();      
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         // Loop over the Ls      
00107         for(; l_it != l_end; ++l_it) {   
00108           DTLayerId dtLId = (*l_it)->id();
00109         
00110           //check if the layer has digi
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             //to fill a map with the average noise per wire and fill new noise histo      
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             //to compute the average noise per layer (excluding the noisy cells)
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             // book the digi event plots every 1000 events
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           }// done if the layer has digi
00210         }// loop over layers
00211       }// loop over superlayers
00212     }// loop over chambers
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();//just silence gcc complaining about unused vars
00225   TF1 *funct=0;
00226   TProfile *theNoiseHisto = new TProfile("theNoiseHisto","Time Constant versus Average Noise",100000,0,100000);
00227   
00228 
00229   //compute the time constant
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     //fill the histo with the time constant as a function of the average noise
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   // histos with the integrated noise per layer
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   // histos with the integrated noise per chamber
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   //overimpose the average noise histo
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   // Loop over the chambers
00334   for (; chamber_it != chamber_end; ++chamber_it) {
00335     vector<const DTSuperLayer*>::const_iterator sl_it = (*chamber_it)->superLayers().begin(); 
00336     vector<const DTSuperLayer*>::const_iterator sl_end = (*chamber_it)->superLayers().end();
00337     // Loop over the SLs
00338     for(; sl_it != sl_end; ++sl_it) {
00339       DTSuperLayerId sl = (*sl_it)->id();
00340       vector<const DTLayer*>::const_iterator l_it = (*sl_it)->layers().begin(); 
00341       vector<const DTLayer*>::const_iterator l_end = (*sl_it)->layers().end();
00342 
00343       string canvasName = "c" + getSuperLayerName(sl); 
00344       TCanvas c1(canvasName.c_str(),canvasName.c_str(),600,780);
00345       TLegend *leg=new TLegend(0.5,0.6,0.7,0.8);
00346       for(; l_it != l_end; ++l_it) {
00347         DTLayerId layerId = (*l_it)->id();
00348         string HistoName = "DigiOccupancy_" + getLayerName(layerId);
00349         theFile->cd();
00350         TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
00351         if(hOccHisto){
00352           string TitleHisto = "AverageNoise_" + getSuperLayerName(sl);
00353           cout<<"TitleHisto : "<<TitleHisto<<endl;
00354           hOccHisto->SetTitle(TitleHisto.c_str());
00355           stringstream layer; layer << layerId.layer(); 
00356           string legendHisto = "layer " + layer.str();
00357           leg->AddEntry(hOccHisto,legendHisto.c_str(),"L");
00358           hOccHisto->SetMaximum(getYMaximum(sl));
00359           histo=true;
00360           if(layerId.layer() == 1)
00361             hOccHisto->Draw();
00362           else
00363             hOccHisto->Draw("same");
00364           hOccHisto->SetLineColor(layerId.layer());
00365         }
00366       }
00367       if(histo){
00368         leg->Draw("same");
00369         theNewFile->cd();
00370         c1.Write();
00371       }
00372     }
00373     histo=false;   
00374   }
00375 
00376   //write on file the noisy plots
00377   for(map<pair<int,int>, TH1F*>::const_iterator nCell = noisyC.begin();
00378       nCell != noisyC.end();
00379       nCell++) {
00380     theNewFile->cd();
00381     (*nCell).second->Write();
00382   }
00383   for(map<pair<int,int>, TH1F*>::const_iterator somehownCell = someHowNoisyC.begin();
00384       somehownCell != someHowNoisyC.end();
00385       somehownCell++) {
00386     theNewFile->cd();
00387     (*somehownCell).second->Write();
00388   }
00389     
00390 }
00391 
00392 
00393 DTNoiseComputation::~DTNoiseComputation(){
00394 
00395   theFile->Close();
00396   theNewFile->Close();
00397 
00398 }
00399 
00400 
00401 string DTNoiseComputation::getLayerName(const DTLayerId& lId) const {
00402 
00403   const  DTSuperLayerId dtSLId = lId.superlayerId();
00404   const  DTChamberId dtChId = dtSLId.chamberId(); 
00405   stringstream Layer; Layer << lId.layer();
00406   stringstream superLayer; superLayer << dtSLId.superlayer();
00407   stringstream wheel; wheel << dtChId.wheel();  
00408   stringstream station; station << dtChId.station();    
00409   stringstream sector; sector << dtChId.sector();
00410   
00411   string LayerName = 
00412     "W" + wheel.str()
00413     + "_St" + station.str() 
00414     + "_Sec" + sector.str() 
00415     + "_SL" + superLayer.str()
00416     + "_L" + Layer.str();
00417   
00418   return LayerName;
00419 
00420 }
00421 
00422 
00423 string DTNoiseComputation::getSuperLayerName(const DTSuperLayerId& dtSLId) const {
00424 
00425   const  DTChamberId dtChId = dtSLId.chamberId(); 
00426   stringstream superLayer; superLayer << dtSLId.superlayer();
00427   stringstream wheel; wheel << dtChId.wheel();  
00428   stringstream station; station << dtChId.station();    
00429   stringstream sector; sector << dtChId.sector();
00430   
00431   string SuperLayerName = 
00432     "W" + wheel.str()
00433     + "_St" + station.str() 
00434     + "_Sec" + sector.str() 
00435     + "_SL" + superLayer.str();
00436   
00437   return SuperLayerName;
00438 
00439 }
00440 
00441 
00442 string DTNoiseComputation::getChamberName(const DTLayerId& lId) const {
00443 
00444   const  DTSuperLayerId dtSLId = lId.superlayerId();
00445   const  DTChamberId dtChId = dtSLId.chamberId(); 
00446   stringstream wheel; wheel << dtChId.wheel();  
00447   stringstream station; station << dtChId.station();    
00448   stringstream sector; sector << dtChId.sector();
00449   
00450   string ChamberName = 
00451     "W" + wheel.str()
00452     + "_St" + station.str() 
00453     + "_Sec" + sector.str();
00454   
00455   return ChamberName;
00456 
00457 }
00458 
00459 
00460 int DTNoiseComputation::getMaxNumBins(const DTChamberId& chId) const {
00461   
00462   int maximum=0;
00463   
00464   for(int SL=1; SL<=3; SL++){
00465     if(!(chId.station()==4 && SL==2)){
00466       for (int L=1; L<=4; L++){ 
00467         DTLayerId layerId = DTLayerId(chId, SL, L);
00468         string HistoName = "DigiOccupancy_" + getLayerName(layerId);
00469         theFile->cd();
00470         TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
00471         if(hOccHisto){
00472           if (hOccHisto->GetXaxis()->GetXmax()>maximum) 
00473             maximum = hOccHisto->GetXaxis()->GetNbins();
00474         }
00475       }
00476     }
00477   }
00478   return maximum;
00479 }
00480 
00481 
00482 double DTNoiseComputation::getYMaximum(const DTSuperLayerId& slId) const {
00483   
00484   double maximum=0;
00485   double dummy = pow(10.,10.);
00486 
00487   for (int L=1; L<=4; L++){
00488     DTLayerId layerId = DTLayerId(slId, L);
00489     string HistoName = "DigiOccupancy_" + getLayerName(layerId);
00490     theFile->cd();
00491     TH1F *hOccHisto = (TH1F *) theFile->Get(HistoName.c_str());
00492     if(hOccHisto){
00493       if (hOccHisto->GetMaximum(dummy)>maximum)
00494         maximum = hOccHisto->GetMaximum(dummy);
00495     }
00496   }
00497   return maximum;
00498 }
00499 
00500 
00501