CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/Calibration/TkAlCaRecoProducers/src/AlcaBeamSpotManager.cc

Go to the documentation of this file.
00001 
00009 #include "Calibration/TkAlCaRecoProducers/interface/AlcaBeamSpotManager.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 #include "DataFormats/Common/interface/Handle.h"
00012 #include "FWCore/Framework/interface/ESHandle.h"
00013 #include <vector>
00014 #include <math.h>
00015 #include <limits.h>
00016 
00017 using namespace edm;
00018 using namespace reco;
00019 using namespace std;
00020 
00021 //--------------------------------------------------------------------------------------------------
00022 AlcaBeamSpotManager::AlcaBeamSpotManager(void){
00023 }
00024 
00025 //--------------------------------------------------------------------------------------------------
00026 AlcaBeamSpotManager::AlcaBeamSpotManager(const ParameterSet& iConfig) :
00027   beamSpotOutputBase_(iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters").getUntrackedParameter<std::string>("BeamSpotOutputBase")),
00028   beamSpotModuleName_(iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters").getUntrackedParameter<std::string>("BeamSpotModuleName")),
00029   beamSpotLabel_     (iConfig.getParameter<ParameterSet>("AlcaBeamSpotHarvesterParameters").getUntrackedParameter<std::string>("BeamSpotLabel"))
00030 {
00031   LogInfo("AlcaBeamSpotManager") 
00032     << "Output base: " << beamSpotOutputBase_ 
00033     << std::endl;
00034   reset();
00035 }
00036 
00037 //--------------------------------------------------------------------------------------------------
00038 AlcaBeamSpotManager::~AlcaBeamSpotManager(void){
00039 }
00040 
00041 //--------------------------------------------------------------------------------------------------
00042 void AlcaBeamSpotManager::reset(void){
00043   beamSpotMap_.clear();
00044 }
00045 //--------------------------------------------------------------------------------------------------
00046 void AlcaBeamSpotManager::readLumi(const LuminosityBlock& iLumi){
00047 
00048   Handle<BeamSpot> beamSpotHandle;
00049   iLumi.getByLabel(beamSpotModuleName_,beamSpotLabel_, beamSpotHandle);
00050   //iLumi.getByLabel("beamspot","alcaBeamSpot", beamSpotHandle); 
00051 
00052   if(beamSpotHandle.isValid()) { // check the product
00053     beamSpotMap_[iLumi.luminosityBlock()] = *beamSpotHandle;
00054     const BeamSpot* aBeamSpot =  &beamSpotMap_[iLumi.luminosityBlock()];
00055     aBeamSpot = beamSpotHandle.product();
00056     LogInfo("AlcaBeamSpotManager")
00057       << "Lumi: " << iLumi.luminosityBlock() << std::endl;
00058     LogInfo("AlcaBeamSpotManager")
00059       << *aBeamSpot << std::endl;
00060   }
00061   else {
00062     LogInfo("AlcaBeamSpotManager")
00063         << "Lumi: " << iLumi.luminosityBlock() << std::endl;
00064     LogInfo("AlcaBeamSpotManager")
00065         << "   BS is not valid!" << std::endl;
00066   }
00067 
00068 }
00069 
00070 //--------------------------------------------------------------------------------------------------
00071 void AlcaBeamSpotManager::createWeightedPayloads(void){
00072   vector<bsMap_iterator> listToErase;
00073   for(bsMap_iterator it=beamSpotMap_.begin(); it!=beamSpotMap_.end();it++){
00074     if(it->second.type() != BeamSpot::Tracker){
00075       listToErase.push_back(it);
00076     }
00077   }
00078   for(vector<bsMap_iterator>::iterator it=listToErase.begin(); it !=listToErase.end(); it++){
00079     beamSpotMap_.erase(*it);
00080   }
00081   if(beamSpotMap_.size() <= 1){
00082     return;
00083   }
00084   //Return only if lumibased since the collapsing alghorithm requires the next and next to next lumi sections
00085   else if(beamSpotMap_.size() == 2 && beamSpotOutputBase_ == "lumibased"){
00086     return;
00087   }
00088   if(beamSpotOutputBase_ == "lumibased"){
00089 //    bsMap_iterator referenceBS = beamSpotMap_.begin();
00090     bsMap_iterator firstBS     = beamSpotMap_.begin();
00091 //    bsMap_iterator lastBS      = beamSpotMap_.begin();
00092     bsMap_iterator currentBS   = beamSpotMap_.begin();
00093     bsMap_iterator nextBS      = ++beamSpotMap_.begin();
00094     bsMap_iterator nextNextBS  = ++(++(beamSpotMap_.begin()));
00095 
00096     map<LuminosityBlockNumber_t,BeamSpot> tmpBeamSpotMap_;
00097     bool docreate = true;
00098     bool endOfRun = false;//Added
00099     bool docheck = true;//Added
00100     bool foundShift = false;
00101     long countlumi = 0;//Added
00102     string tmprun = "";//Added
00103     long maxNlumis = 60;//Added
00104 //    if weighted:
00105 //        maxNlumis = 999999999
00106 
00107 
00108     unsigned int iteration = 0;
00109     //while(nextNextBS!=beamSpotMap_.end()){
00110     while(nextBS!=beamSpotMap_.end()){
00111       LogInfo("AlcaBeamSpotManager")
00112         << "Iteration: " << iteration << " size: " << beamSpotMap_.size() << "\n"
00113         << "Lumi: " << currentBS->first  << "\n" 
00114         << currentBS->second
00115         << "\n" << nextBS->first     << "\n" << nextBS->second  
00116         << endl;  
00117       if (nextNextBS!=beamSpotMap_.end())
00118         LogInfo("AlcaBeamSpotManager")
00119           << nextNextBS->first << "\n" << nextNextBS->second
00120           << endl;
00121 
00122 
00123       if(docreate){                                                                                                                                                                                                                                                                 
00124         firstBS = currentBS;                                                                                                                                                                                                                                
00125         docreate = false;//Added                                                                                                                                                                                                                                                            
00126       }
00127       //if(iteration >= beamSpotMap_.size()-3){
00128       if(iteration >= beamSpotMap_.size()-2){
00129         LogInfo("AlcaBeamSpotManager")
00130           << "Reached lumi " << currentBS->first
00131           << " now close payload because end of data has been reached.";                                                                                                                                                                                                                    
00132           docreate = true;
00133           endOfRun = true;                                                                                                                                                                                                                                                          
00134       }    
00135       // check we run over the same run                                                                                                                                                                                                                                             
00136 //      if (ibeam->first.Run() != inextbeam->first.Run()){                                                                                                                                                                                                                          
00137 //        LogInfo("AlcaBeamSpotManager")                                                                                                                                                                                                                                            
00138 //          << "close payload because end of run.";                                                                                                                                                                                                                                 
00139 //        docreate = true;                                                                                                                                                                                                                                                          
00140 //      }                                                                                                                                                                                                                                                                           
00141       // check maximum lumi counts
00142       if (countlumi == maxNlumis -1){                                                                                                                                                                                                                                               
00143         LogInfo("AlcaBeamSpotManager")                                                                                                                                                                                                                                              
00144           << "close payload because maximum lumi sections accumulated within run ";                                                                                                                                                                                                 
00145         docreate = true;                                                                                                                                                                                                                                                            
00146         countlumi = 0;                                                                                                                                                                                                                                                              
00147       }
00148 //      if weighted:                                                                                                                                                                                                                                                                
00149 //          docheck = False                                                                                                                                                                                                                                                         
00150       // check offsets                                                                                                                                                                                                                                                              
00151       if(docheck){                                                                                                                                                                                                                                                                  
00152         foundShift = false;
00153         LogInfo("AlcaBeamSpotManager")                                                                                                                                                                                                                                      
00154           << "Checking checking!" << endl;
00155         float limit = 0;                                                                                                                                                                                                                                                            
00156         pair<float,float> adelta1;                                                                                                                                                                                                                                                  
00157         pair<float,float> adelta2;                                                                                                                                                                                                                                                  
00158         pair<float,float> adelta;                                                                                                                                                                                                                                 
00159         pair<float,float> adelta1dxdz;
00160         pair<float,float> adelta2dxdz;
00161         pair<float,float> adelta1dydz;
00162         pair<float,float> adelta2dydz;
00163         pair<float,float> adelta1widthX;
00164         pair<float,float> adelta2widthX;
00165         pair<float,float> adelta1widthY;
00166         pair<float,float> adelta2widthY;
00167         pair<float,float> adelta1z0;
00168         pair<float,float> adelta1sigmaZ;
00169 
00170         // define minimum limit
00171         float min_limit = 0.0025;
00172         
00173         // limit for x and y
00174         limit = currentBS->second.BeamWidthX()/2.;                                                                                                                                                                                                                                  
00175         if(limit < min_limit){
00176           limit = min_limit;
00177         }
00178         
00179         //check movements in X                                                                                                                                                                                                                                                              
00180         adelta1 = delta(currentBS->second.x0(), currentBS->second.x0Error(), nextBS->second.x0(), nextBS->second.x0Error());                                                                                                                                                        
00181         adelta2 = pair<float,float>(0.,1.e9);                                                                                                                                                                                                                                       
00182         if (nextNextBS->second.type() != -1){                                                                                                                                                                                                                                       
00183             adelta2 = delta(nextBS->second.x0(), nextBS->second.x0Error(), nextNextBS->second.x0(), nextNextBS->second.x0Error());                                                                                                                                                  
00184         }                                                                                                                                                                                                                                                                           
00185         bool deltaX = (deltaSig(adelta1.first,adelta1.second) > 3.5 && adelta1.first >= limit)?true:false;                                                                                                                                                                          
00186         if(iteration < beamSpotMap_.size()-2){                                                                                                                                                                                                                                      
00187             if( !deltaX && adelta1.first*adelta2.first > 0. &&  fabs(adelta1.first+adelta2.first) >= limit){                                                                                                                                                                        
00188               LogInfo("AlcaBeamSpotManager")                                                                                                                                                                                                                                        
00189                 << " positive, " << (adelta1.first+adelta2.first) << " limit=" << limit << endl;                                                                                                                                                                                    
00190                 deltaX = true;                                                                                                                                                                                                                                                      
00191             }
00192             else if( deltaX && adelta1.first*adelta2.first < 0 && adelta2.first != 0 && fabs(adelta1.first/adelta2.first) > 0.33 && fabs(adelta1.first/adelta2.first) < 3){                                                                                                                         
00193               LogInfo("AlcaBeamSpotManager")                                                                                                                                                                                                                                        
00194                 << " negative, " << adelta1.first/adelta2.first << endl;                                                                                                                                                                                                            
00195                 deltaX = false;                                                                                                                                                                                                                                                     
00196             }
00197         }                                                                                                                                                                                                                                                         
00198 
00199         //calculating all deltas                                                                                                                                                                                                                                                                    
00200         adelta1dxdz = delta(currentBS->second.dxdz(), currentBS->second.dxdzError(), nextBS->second.dxdz(), nextBS->second.dxdzError());
00201         adelta2dxdz = pair<float,float>(0.,1.e9);
00202         adelta1dydz = delta(currentBS->second.dydz(), currentBS->second.dydzError(), nextBS->second.dydz(), nextBS->second.dydzError());
00203         adelta2dydz = pair<float,float>(0.,1.e9);
00204         adelta1widthX = delta(currentBS->second.BeamWidthX(), currentBS->second.BeamWidthXError(), nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError());
00205         adelta2widthX = pair<float,float>(0.,1.e9);
00206         adelta1widthY = delta(currentBS->second.BeamWidthY(), currentBS->second.BeamWidthYError(), nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError());
00207         adelta2widthY = pair<float,float>(0.,1.e9);
00208         adelta1z0 = delta(currentBS->second.z0(), currentBS->second.z0Error(), nextBS->second.z0(), nextBS->second.z0Error());
00209         adelta1sigmaZ = delta(currentBS->second.sigmaZ(), currentBS->second.sigmaZ0Error(), nextBS->second.sigmaZ(), nextBS->second.sigmaZ0Error());
00210         
00211         //check movements in Y                                                                                                                                                                                                                                                              
00212         adelta1 = delta(currentBS->second.y0(), currentBS->second.y0Error(), nextBS->second.y0(), nextBS->second.y0Error());                                                                                                                                                        
00213         adelta2 = pair<float,float>(0.,1.e9);                                                                                                                                                                                                                                       
00214         if( nextNextBS->second.type() != BeamSpot::Unknown){                                                                                                                                                                                                                        
00215           adelta2 = delta(nextBS->second.y0(), nextBS->second.y0Error(), nextNextBS->second.y0(), nextNextBS->second.y0Error());                                                                                                                                                    
00216           adelta2dxdz = delta(nextBS->second.dxdz(), nextBS->second.dxdzError(), nextNextBS->second.dxdz(), nextNextBS->second.dxdzError());
00217           adelta2dydz = delta(nextBS->second.dydz(), nextBS->second.dydzError(), nextNextBS->second.dydz(), nextNextBS->second.dydzError());
00218           adelta2widthX = delta(nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError(), nextNextBS->second.BeamWidthX(), nextNextBS->second.BeamWidthXError());
00219           adelta2widthY = delta(nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError(), nextNextBS->second.BeamWidthY(), nextNextBS->second.BeamWidthYError());
00220 
00221         }                                                                                                                                                                                                                                                                           
00222         bool deltaY = (deltaSig(adelta1.first,adelta1.second) > 3.5 && adelta1.first >= limit)?true:false;                                                                                                                                                                          
00223         if(iteration < beamSpotMap_.size()-2){                                                                                                                                                                                                                                      
00224           if( !deltaY && adelta1.first*adelta2.first > 0. &&  fabs(adelta1.first+adelta2.first) >= limit){                                                                                                                                                                        
00225             LogInfo("AlcaBeamSpotManager")                                                                                                                                                                                                                                        
00226               << " positive, " << (adelta1.first+adelta2.first) << " limit=" << limit << endl;                                                                                                                                                                                    
00227               deltaY = true;                                                                                                                                                                                                                                                      
00228           }
00229           else if( deltaY && adelta1.first*adelta2.first < 0 && adelta2.first != 0 && fabs(adelta1.first/adelta2.first) > 0.33 && fabs(adelta1.first/adelta2.first) < 3){                                                                                                                                 
00230             LogInfo("AlcaBeamSpotManager")                                                                                                                                                                                                                                        
00231               << " negative, " << adelta1.first/adelta2.first << endl;                                                                                                                                                                                                            
00232               deltaY = false;                                                                                                                                                                                                                                                     
00233           }
00234         }
00235 
00236         limit = currentBS->second.sigmaZ()/2.;                                                                                                                                                                                                                          
00237         bool deltaZ = (deltaSig(adelta1z0.first,adelta1z0.second) > 3.5 && fabs(adelta1z0.first) >= limit)?true:false;                                                
00238         adelta = delta(currentBS->second.sigmaZ(), currentBS->second.sigmaZ0Error(), nextBS->second.sigmaZ(), nextBS->second.sigmaZ0Error());                                                                                                                           
00239         bool deltasigmaZ = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;                                                                                                                                                                                             
00240         bool deltadxdz = false;
00241         bool deltadydz = false;
00242         bool deltawidthX = false;
00243         bool deltawidthY = false;
00244 
00245         if(iteration < beamSpotMap_.size()-2){                                                                                                                                                              
00246 
00247           adelta = delta(currentBS->second.dxdz(), currentBS->second.dxdzError(), nextBS->second.dxdz(), nextBS->second.dxdzError());                                                                                                                                                 
00248           deltadxdz   = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;                                                                                                                                                                                                
00249           if(deltadxdz && (adelta1dxdz.first*adelta2dxdz.first) < 0 && adelta2dxdz.first != 0 && fabs(adelta1dxdz.first/adelta2dxdz.first) > 0.33 && fabs(adelta1dxdz.first/adelta2dxdz.first) < 3){
00250             deltadxdz = false;
00251           }
00252 
00253           adelta = delta(currentBS->second.dydz(), currentBS->second.dydzError(), nextBS->second.dydz(), nextBS->second.dydzError());                                                                                                                                                 
00254           deltadydz   = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;                                                                                                                                                                                                
00255           if(deltadydz && (adelta1dydz.first*adelta2dydz.first) < 0 && adelta2dydz.first != 0 && fabs(adelta1dydz.first/adelta2dydz.first) > 0.33 && fabs(adelta1dydz.first/adelta2dydz.first) < 3){
00256             deltadydz = false;
00257           }
00258 
00259           adelta = delta(currentBS->second.BeamWidthX(), currentBS->second.BeamWidthXError(), nextBS->second.BeamWidthX(), nextBS->second.BeamWidthXError());                                                                                                                   
00260           deltawidthX = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;                                                                                                                                                                                                
00261           if(deltawidthX && (adelta1widthX.first*adelta2widthX.first) < 0 && adelta2widthX.first != 0 && fabs(adelta1widthX.first/adelta2widthX.first) > 0.33 && fabs(adelta1widthX.first/adelta2widthX.first) < 3){
00262             deltawidthX = false;
00263           }
00264 
00265           adelta = delta(currentBS->second.BeamWidthY(), currentBS->second.BeamWidthYError(), nextBS->second.BeamWidthY(), nextBS->second.BeamWidthYError());                                                                                                                   
00266           deltawidthY = (deltaSig(adelta.first,adelta.second) > 5.0)?true:false;                                                                                                                                                                                                
00267           if(deltawidthY && (adelta1widthY.first*adelta2widthY.first) < 0 && adelta2widthY.first != 0 && fabs(adelta1widthY.first/adelta2widthY.first) > 0.33 && fabs(adelta1widthY.first/adelta2widthY.first) < 3){
00268             deltawidthY = false;
00269           }
00270 
00271         }
00272         if (deltaX || deltaY || deltaZ || deltasigmaZ || deltadxdz || deltadydz || deltawidthX || deltawidthY){                                                                                                                                                         
00273           docreate = true;
00274           foundShift = true;
00275           LogInfo("AlcaBeamSpotManager")                                                                                                                                                                                                                                        
00276             << "close payload because of movement in"                                                                                                                                                                                                                   
00277             <<  " X=" << deltaX 
00278             << ", Y=" << deltaY 
00279             << ", Z=" << deltaZ 
00280             << ", sigmaZ=" << deltasigmaZ 
00281             << ", dxdz=" << deltadxdz 
00282             << ", dydz=" << deltadydz
00283             << ", widthX=" << deltawidthX
00284             << ", widthY=" << deltawidthY
00285             << endl;  
00286         }
00287 
00288       }
00289       if(docreate){
00290         if(foundShift){
00291           tmpBeamSpotMap_[firstBS->first] = weight(firstBS,nextBS);
00292           if (endOfRun){ 
00293             //if we're here, then we need to found a shift in the last LS
00294             //We already created a new IOV, now create one just for the last LS
00295             tmpBeamSpotMap_[nextBS->first] = nextBS->second;
00296           }
00297         }
00298         else if(!foundShift && !endOfRun){ //maxLS reached
00299           tmpBeamSpotMap_[firstBS->first] = weight(firstBS,nextBS);
00300         }
00301         else { // end of run with no shift detectred in last LS
00302           tmpBeamSpotMap_[firstBS->first] = weight(firstBS,beamSpotMap_.end());
00303         }
00304         firstBS = nextBS;
00305           countlumi = 0;
00306                                                                                                                                                                                                                                                                     
00307       }
00308       //tmprun = currentBS->second.Run                                                                                                                                                                                                                                              
00309       ++countlumi;                                                                                                                                                                                                                                                                  
00310       
00311       currentBS = nextBS;
00312       nextBS    = nextNextBS;
00313       nextNextBS++;
00314       ++iteration;
00315     }
00316     beamSpotMap_.clear();
00317     beamSpotMap_ = tmpBeamSpotMap_;
00318   }
00319   else if(beamSpotOutputBase_ == "runbased"){
00320     BeamSpot aBeamSpot = weight(beamSpotMap_.begin(),beamSpotMap_.end());
00321     LuminosityBlockNumber_t firstLumi = beamSpotMap_.begin()->first;
00322     beamSpotMap_.clear();
00323     beamSpotMap_[firstLumi] = aBeamSpot;
00324   }
00325   else{
00326     LogInfo("AlcaBeamSpotManager")
00327       << "Unrecognized BeamSpotOutputBase parameter: " << beamSpotOutputBase_
00328       << endl;    
00329   }
00330 }
00331 
00332 //--------------------------------------------------------------------------------------------------
00333 BeamSpot AlcaBeamSpotManager::weight(const bsMap_iterator& begin,
00334                                      const bsMap_iterator& end){
00335   double x,xError = 0;
00336   double y,yError = 0;
00337   double z,zError = 0;
00338   double sigmaZ,sigmaZError = 0;
00339   double dxdz,dxdzError = 0;
00340   double dydz,dydzError = 0;
00341   double widthX,widthXError = 0;
00342   double widthY,widthYError = 0;
00343   LogInfo("AlcaBeamSpotManager")
00344     << "Weighted BeamSpot will span lumi " 
00345     << begin->first << " to " << end->first
00346     << endl;
00347 
00348   BeamSpot::BeamType type = BeamSpot::Unknown;
00349   for(bsMap_iterator it=begin; it!=end; it++){
00350     weight(x     , xError     , it->second.x0()        , it->second.x0Error());
00351     weight(y     , yError     , it->second.y0()        , it->second.y0Error());
00352     weight(z     , zError     , it->second.z0()        , it->second.z0Error());
00353     weight(sigmaZ, sigmaZError, it->second.sigmaZ()    , it->second.sigmaZ0Error());
00354     weight(dxdz  , dxdzError  , it->second.dxdz()      , it->second.dxdzError());
00355     weight(dydz  , dydzError  , it->second.dydz()      , it->second.dydzError());
00356     weight(widthX, widthXError, it->second.BeamWidthX(), it->second.BeamWidthXError());
00357     weight(widthY, widthYError, it->second.BeamWidthY(), it->second.BeamWidthYError());
00358     if(it->second.type() == BeamSpot::Tracker){
00359       type = BeamSpot::Tracker;
00360     }
00361   }
00362   BeamSpot::Point bsPosition(x,y,z);
00363   BeamSpot::CovarianceMatrix error;
00364   error(0,0) = xError*xError;
00365   error(1,1) = yError*yError;
00366   error(2,2) = zError*zError;
00367   error(3,3) = sigmaZError*sigmaZError;
00368   error(4,4) = dxdzError*dxdzError;
00369   error(5,5) = dydzError*dydzError;
00370   error(6,6) = widthXError*widthXError;
00371   BeamSpot weightedBeamSpot(bsPosition,sigmaZ,dxdz,dydz,widthX,error,type);
00372   weightedBeamSpot.setBeamWidthY(widthY);
00373   LogInfo("AlcaBeamSpotManager")
00374     << "Weighted BeamSpot will be:" <<'\n'
00375     << weightedBeamSpot
00376     << endl;
00377   return weightedBeamSpot;
00378 }
00379 
00380 //--------------------------------------------------------------------------------------------------
00381 void AlcaBeamSpotManager::weight(double& mean,double& meanError,const double& val,const double& valError){
00382     double tmpError = 0;
00383     if (meanError < 1e-8){
00384         tmpError = 1/(valError*valError);
00385         mean = val*tmpError;
00386     }
00387     else{
00388         tmpError = 1/(meanError*meanError) + 1/(valError*valError);
00389         mean = mean/(meanError*meanError) + val/(valError*valError);
00390     }
00391     mean = mean/tmpError;
00392     meanError = sqrt(1/tmpError);
00393 }
00394 
00395 //--------------------------------------------------------------------------------------------------
00396 pair<float,float> AlcaBeamSpotManager::delta(const float& x, const float& xError, const float& nextX, const float& nextXError){
00397   return pair<float,float>(x - nextX, sqrt(pow(xError,2) + pow(nextXError,2)) );
00398 }
00399 
00400 //--------------------------------------------------------------------------------------------------
00401 float AlcaBeamSpotManager::deltaSig(const float& num, const float& den){
00402   if(den != 0){
00403     return fabs(num/den);
00404   }
00405   else{
00406     return float(LONG_MAX);
00407   }
00408 }
00409