CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/DQM/EcalCommon/interface/UtilFunctions.h

Go to the documentation of this file.
00001 #ifndef UtilFunctions_H
00002 #define UtilFunctions_H
00003 
00012 #include <cmath>
00013 #include <iostream>
00014 
00015 #include "TH1F.h"
00016 #include "TProfile.h"
00017 #include "TClass.h"
00018 
00019 
00020 namespace ecaldqm {
00021 
00022   // functions implemented here are not universal in the sense that
00023   // the input variables are changed due to efficiency of memory usage.
00024 
00025 
00026   // calculated time intervals and bin locations for time varing profiles
00027 
00028   void calcBins(int binWidth, int divisor, long int start_time, long int last_time, long int current_time,
00029                 long int & binDiff, long int & diff) {
00030 
00031     // changing arguments : binDiff, diff
00032 
00033     // binWidth : time interval
00034     // divisor : time unit - for minute case divisor is 60 and for hour case 3600
00035     // start_time : initial time when the job started
00036     // last_time : the last updated time before calling the current "analyze" function
00037     // current_time : the current time inside "analyze" fucntion
00038     // binDiff : the bin difference for the current time compared to the bin location of the last time
00039     // diff : time difference between the current time and the last time
00040 
00041     long int diff_current_start = current_time - start_time;
00042     long int diff_last_start = last_time - start_time;
00043 
00044     // --------------------------------------------------
00045     // Calculate time interval and bin width
00046     // --------------------------------------------------
00047 
00048     binDiff = diff_current_start/divisor/binWidth - diff_last_start/divisor/binWidth;
00049     diff = (current_time - last_time)/divisor;
00050 
00051     if(diff >= binWidth) {
00052       while(diff >= binWidth) diff -= binWidth;
00053     }
00054 
00055   } // calcBins
00056 
00057 
00058 
00059   // shift bins in TProfile to the right
00060 
00061   void shift2Right(TProfile* p, int bins){
00062 
00063     // bins : how many bins need to be shifted
00064 
00065     if(bins <= 0) return;
00066 
00067     if(!p->GetSumw2()) p->Sumw2();
00068     int nBins = p->GetXaxis()->GetNbins();
00069 
00070     // by shifting n bin to the right, the number of entries are
00071     // reduced by the number in n bins including the overflow bin.
00072     double nentries = p->GetEntries();
00073     for(int i=0; i<bins; i++) nentries -= p->GetBinEntries(i);
00074     p->SetEntries(nentries);
00075   
00076     // the last bin goes to overflow
00077     // each bin moves to the right
00078 
00079     TArrayD* sumw2 = p->GetSumw2();
00080 
00081     for(int i=nBins+1; i>bins; i--) {
00082       // GetBinContent return binContent/binEntries
00083       p->SetBinContent(i, p->GetBinContent(i-bins)*p->GetBinEntries(i-bins));
00084       p->SetBinEntries(i,p->GetBinEntries(i-bins));
00085       sumw2->SetAt(sumw2->GetAt(i-bins),i);
00086     }
00087     
00088   }
00089 
00090 
00091   // shift bins in TProfile to the left
00092 
00093   void shift2Left(TProfile* p, int bins){
00094 
00095     if(bins <= 0) return;
00096 
00097     if(!p->GetSumw2()) p->Sumw2();
00098     int nBins = p->GetXaxis()->GetNbins();
00099 
00100     // by shifting n bin to the left, the number of entries are
00101     // reduced by the number in n bins including the underflow bin.
00102     double nentries = p->GetEntries();
00103     for(int i=0; i<bins; i++) nentries -= p->GetBinEntries(i);
00104     p->SetEntries(nentries);
00105   
00106     // the first bin goes to underflow
00107     // each bin moves to the right
00108 
00109     TArrayD* sumw2 = p->GetSumw2();
00110 
00111     for(int i=0; i<=nBins+1-bins; i++) {
00112       // GetBinContent return binContent/binEntries
00113       p->SetBinContent(i, p->GetBinContent(i+bins)*p->GetBinEntries(i+bins));
00114       p->SetBinEntries(i,p->GetBinEntries(i+bins));
00115       sumw2->SetAt(sumw2->GetAt(i+bins),i);
00116     }
00117 
00118   }
00119 
00120 
00121   // get mean and rms of Y values from TProfile
00122 
00123   void getAverageFromTProfile(TProfile* p, double& mean, double& rms) {
00124 
00125     // changing arguments : mean, rms
00126     mean = rms = 0.0;
00127 
00128     if(!p) return;
00129 
00130     int nbins = p->GetXaxis()->GetNbins();
00131     double y = 0.0;
00132     double y2 = 0.0;
00133     for(int i=1; i<=nbins; i++){
00134       y += p->GetBinContent(i);
00135       y2 += y*y;
00136     }
00137     mean = y/nbins;
00138     rms = std::sqrt(y2/nbins - mean*mean);
00139 
00140   } // getAverageFromTProfile
00141 
00142 
00143   // get mean and rms based on two histograms' difference
00144 
00145   void getMeanRms(TObject* pre, TObject* cur, double& mean, double& rms) {
00146 
00147     // changing arguments : mean, rms
00148 
00149     mean = rms = 0.0;
00150 
00151     if(!cur) return;
00152 
00153     TString name(cur->IsA()->GetName());
00154 
00155     if(name.Contains("TProfile")) {
00156       getAverageFromTProfile((TProfile*)cur,mean,rms);
00157     }
00158     else if(name.Contains("TH2")) {
00159       if(pre) {
00160         mean = ((TH2F*)cur)->GetEntries() - ((TH2F*)pre)->GetEntries();
00161         rms = std::sqrt(mean);
00162       }
00163       else {
00164         mean = ((TH2F*)cur)->GetEntries();
00165         rms = std::sqrt(mean);
00166       }
00167       float nxybins = ((TH2F*)cur)->GetNbinsX()*((TH2F*)cur)->GetNbinsY();
00168       mean /= nxybins;
00169       rms /= nxybins;
00170     }
00171     else if(name.Contains("TH1")) {
00172       if(pre) {
00173         ((TH1F*)pre)->Sumw2();
00174         ((TH1F*)pre)->Add((TH1F*)pre,(TH1F*)cur,-1,1);
00175         mean = ((TH1F*)pre)->GetMean();
00176         rms = ((TH1F*)pre)->GetRMS();
00177       }
00178       else {
00179         mean = ((TH1F*)cur)->GetMean();
00180         rms = ((TH1F*)cur)->GetRMS();
00181       }
00182     }
00183 
00184   } // getMeanRms
00185 
00186   TObject* cloneIt(MonitorElement* me, std::string histo) {
00187 
00188     // The cloned object, ret should be deleted after using it.
00189 
00190     TObject* ret = 0;
00191 
00192     if(!me) return ret;
00193 
00194     std::string title = histo + " Clone";
00195     ret = (TObject*) (me->getRootObject()->Clone(title.c_str()));
00196 
00197     return ret;
00198   }
00199 
00200 
00201 }
00202 
00203 #endif