UserCode/joshmt/MiscUtil.cxx

/*
a place for very simple functions that don't depend on anything else

can be easily #included into other code
*/

#include "TH1D.h"
#include "TString.h"
#include <iostream>
#include "TMath.h"

#include <cassert>

namespace jmt {

  //======misc utilities======
  //gets rid of = > < from cuts in order to be better included in file names
  TString fortranize(TString cut) {

    cut.ReplaceAll(" ",""); //remove all spaces

    cut.ReplaceAll("==","eq");
    cut.ReplaceAll(">=","gte");
    cut.ReplaceAll("<=","lte");
    
    cut.ReplaceAll(">","gt");
    cut.ReplaceAll("<","lt");

    cut.ReplaceAll("&&","and");
    cut.ReplaceAll("||","or");

    cut.ReplaceAll("/","Over");
    cut.ReplaceAll("+","Plus");
    cut.ReplaceAll("*","Times");

    //this is pretty ugly
    cut.ReplaceAll("(","L");
    cut.ReplaceAll(")","R");

    return cut;
  }

  //utility function for making output more readable
  TString format_nevents(double n,double e, const bool moreDigits=false) {
    const TString latexMode_=true; //hard code for now
    const TString pm = latexMode_ ?  "\\pm ": " +/- ";

    const int eCutoff = moreDigits ? 10 : 1;
    const int extraDigits = moreDigits ? 1:0;

    TString mathmode = latexMode_ ? "$" : "";
  
    char out[100];
    if (e >= eCutoff || e < 0.00001) { //show whole numbers only
      sprintf(out,"%s%.0f%s%.0f%s",mathmode.Data(),n,pm.Data(),e,mathmode.Data());
    }
    else {
      int nfig = ceil(fabs(log10(e))) + extraDigits;
      TString form="%s%.";
      form+=nfig; form+="f%s%.";
      form+=nfig; form+="f%s";
      sprintf(out,form.Data(),mathmode.Data(),n,pm.Data(),e,mathmode.Data());
    }
    return TString(out);
  }


  Double_t weightedMean(Int_t n, Double_t *a, Double_t *e, const bool returnError=false) {
    using namespace std;

    Double_t numerator=0;
    Double_t denominator=0;

    for (Int_t i=0; i<n; ++i) {
      Double_t esq = e[i]*e[i];
      if (esq == 0) {cout<<"Error is zero!"<<endl; return -1e9;}
      numerator   += a[i] / esq;
      denominator += 1.0  / esq;
    }

    cout<<"mean = "<<numerator / denominator<<" +/- "<<sqrt(1/denominator)<<endl;
    return returnError ? sqrt(1/denominator) : numerator / denominator;
  }

  //implemented in TMath in later versions of root, so i'm using the TMath name
  bool AreEqualAbs( const double & a, const double & b, const double epsilon=0.001) {

    return ( fabs(a-b) < epsilon ) ;

  }

  //deal with the annoyance of logical booleans that are stored in an ntuple as double
  bool doubleToBool( double a) {
    using namespace std;

    int i = TMath::Nint(a);
    if (i==0) return false;
    else if (i==1) return true;
    
    cout<<"[doubleToBool] Something weird going on! "<<a<<"\t"<<i<<endl;
    if (i>1) return true;
    return false;
  }

  //delta phi calculations are pretty commonplace, so put it here
  double deltaPhi(double phi1, double phi2) { 
    double result = phi1 - phi2;
    while (result > TMath::Pi()) result -= 2*TMath::Pi();
    while (result <= -TMath::Pi()) result += 2*TMath::Pi();
    return fabs(result);
  }

  double deltaR2(double eta1, double phi1, double eta2, double phi2) {
    double deta = eta1 - eta2;
    double dphi = deltaPhi(phi1, phi2);
    return deta*deta + dphi*dphi;
  }

  double deltaR(double eta1, double phi1, double eta2, double phi2) {
    return TMath::Sqrt(deltaR2 (eta1, phi1, eta2, phi2));
  }

  // == error propagation ==

  //because not all versions of ROOT have it built in
  double errOnIntegral(const TH1D* h, int binlow=1, int binhigh=-1) {
    
    if (binhigh == -1) binhigh = h->GetNbinsX();
    
    double err=0;
    for (int i = binlow; i<= binhigh; i++) {
      err += h->GetBinError(i) * h->GetBinError(i);
    }
    
    if (err<0) return err;
    return sqrt(err);
  }

  //return error on a/b
  double errAoverB(double a, double aerr, double b, double berr) {
    using namespace std;
    if (b==0) {
      cout<<"Warning in errAoverB -- b is zero!"<<endl;
      return -1;
    }
    return (1/b)*sqrt( aerr*aerr + a*a*berr*berr/(b*b) );
  }

  //return error on a*b
  double errAtimesB(double a, double aerr, double b, double berr) {
    return sqrt( b*b*aerr*aerr + a*a*berr*berr);
  }

  //======== container for run, lumisection, event number =========
  class eventID {
  public:
    eventID();    
    eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber);
    ~eventID();
    
    bool operator< (const eventID & id) const;
    bool operator== (const eventID & id) const;
    bool operator!= (const eventID & id) const;
    
    ULong64_t run;
    ULong64_t ls;
    ULong64_t ev;
    
  };
  
  eventID::eventID() : run(0),ls(0),ev(0) {}
  eventID::eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber) : 
    run(RunNumber),ls(LumiSection),ev(EventNumber) {}
  eventID::~eventID() {}
  
  bool eventID::operator== (const eventID & id) const {

    if (ev == id.ev &&
        ls == id.ls &&
        run == id.run) return true;
    
    return false;
  }

  bool eventID::operator!= (const eventID & id) const {
    return !(*this == id);
  }

  bool eventID::operator< (const eventID & id) const {
    
    if (run < id.run) return true;
    else if (run > id.run) return false;    
    else { //if run is equal

      //now compare lumi section
      if ( ls < id.ls ) return true;
      else if (ls > id.ls) return false;    
      else { //if ls is equal
        
        if ( ev < id.ev ) return true;
        else if (ev > id.ev) return false;    
        
        else { //these are the same event!
          return false;
        }
      }
    }

    //this shouldn't happen
    assert(0);
    return false;
  }
  
} //end of namespace


Revision:	1.10
Committed:	Wed Feb 15 15:13:44 2012 UTC (13 years, 2 months ago) by joshmt
Content type:	text/plain
Branch:	MAIN
Changes since 1.9:	+4 -10 lines
Log Message:	more robust against 0
#	User	Rev	Content
1	joshmt	1.1	/*
2			a place for very simple functions that don't depend on anything else
3
4			can be easily #included into other code
5			*/
6
7	joshmt	1.5	#include "TH1D.h"
8	joshmt	1.1	#include "TString.h"
9	joshmt	1.3	#include <iostream>
10	joshmt	1.9	#include "TMath.h"
11
12			#include <cassert>
13	joshmt	1.1
14			namespace jmt {
15
16			//======misc utilities======
17			//gets rid of = > < from cuts in order to be better included in file names
18			TString fortranize(TString cut) {
19
20	joshmt	1.8	cut.ReplaceAll(" ",""); //remove all spaces
21
22	joshmt	1.1	cut.ReplaceAll("==","eq");
23			cut.ReplaceAll(">=","gte");
24			cut.ReplaceAll("<=","lte");
25
26			cut.ReplaceAll(">","gt");
27			cut.ReplaceAll("<","lt");
28	joshmt	1.4
29	joshmt	1.8	cut.ReplaceAll("&&","and");
30			cut.ReplaceAll("\|\|","or");
31
32	joshmt	1.4	cut.ReplaceAll("/","Over");
33			cut.ReplaceAll("+","Plus");
34			cut.ReplaceAll("*","Times");
35
36			//this is pretty ugly
37			cut.ReplaceAll("(","L");
38			cut.ReplaceAll(")","R");
39
40	joshmt	1.1	return cut;
41			}
42
43	joshmt	1.6	//utility function for making output more readable
44			TString format_nevents(double n,double e, const bool moreDigits=false) {
45			const TString latexMode_=true; //hard code for now
46			const TString pm = latexMode_ ? "\\pm ": " +/- ";
47
48			const int eCutoff = moreDigits ? 10 : 1;
49			const int extraDigits = moreDigits ? 1:0;
50
51			TString mathmode = latexMode_ ? "$" : "";
52
53			char out[100];
54			if (e >= eCutoff \|\| e < 0.00001) { //show whole numbers only
55			sprintf(out,"%s%.0f%s%.0f%s",mathmode.Data(),n,pm.Data(),e,mathmode.Data());
56			}
57			else {
58			int nfig = ceil(fabs(log10(e))) + extraDigits;
59			TString form="%s%.";
60			form+=nfig; form+="f%s%.";
61			form+=nfig; form+="f%s";
62			sprintf(out,form.Data(),mathmode.Data(),n,pm.Data(),e,mathmode.Data());
63			}
64			return TString(out);
65			}
66
67
68	joshmt	1.7	Double_t weightedMean(Int_t n, Double_t a, Double_t e, const bool returnError=false) {
69	joshmt	1.9	using namespace std;
70	joshmt	1.3
71			Double_t numerator=0;
72			Double_t denominator=0;
73
74			for (Int_t i=0; i<n; ++i) {
75			Double_t esq = e[i]*e[i];
76	joshmt	1.7	if (esq == 0) {cout<<"Error is zero!"<<endl; return -1e9;}
77	joshmt	1.3	numerator += a[i] / esq;
78			denominator += 1.0 / esq;
79			}
80
81			cout<<"mean = "<<numerator / denominator<<" +/- "<<sqrt(1/denominator)<<endl;
82	joshmt	1.7	return returnError ? sqrt(1/denominator) : numerator / denominator;
83	joshmt	1.3	}
84
85			//implemented in TMath in later versions of root, so i'm using the TMath name
86			bool AreEqualAbs( const double & a, const double & b, const double epsilon=0.001) {
87
88			return ( fabs(a-b) < epsilon ) ;
89
90			}
91
92	joshmt	1.9	//deal with the annoyance of logical booleans that are stored in an ntuple as double
93			bool doubleToBool( double a) {
94			using namespace std;
95
96			int i = TMath::Nint(a);
97			if (i==0) return false;
98			else if (i==1) return true;
99
100			cout<<"[doubleToBool] Something weird going on! "<<a<<"\t"<<i<<endl;
101			if (i>1) return true;
102			return false;
103			}
104
105			//delta phi calculations are pretty commonplace, so put it here
106			double deltaPhi(double phi1, double phi2) {
107			double result = phi1 - phi2;
108			while (result > TMath::Pi()) result -= 2*TMath::Pi();
109			while (result <= -TMath::Pi()) result += 2*TMath::Pi();
110			return fabs(result);
111			}
112
113			double deltaR2(double eta1, double phi1, double eta2, double phi2) {
114			double deta = eta1 - eta2;
115			double dphi = deltaPhi(phi1, phi2);
116			return detadeta + dphidphi;
117			}
118
119			double deltaR(double eta1, double phi1, double eta2, double phi2) {
120			return TMath::Sqrt(deltaR2 (eta1, phi1, eta2, phi2));
121			}
122
123	joshmt	1.3	// == error propagation ==
124
125			//because not all versions of ROOT have it built in
126			double errOnIntegral(const TH1D* h, int binlow=1, int binhigh=-1) {
127
128			if (binhigh == -1) binhigh = h->GetNbinsX();
129
130			double err=0;
131			for (int i = binlow; i<= binhigh; i++) {
132			err += h->GetBinError(i) * h->GetBinError(i);
133			}
134
135			if (err<0) return err;
136			return sqrt(err);
137			}
138
139			//return error on a/b
140			double errAoverB(double a, double aerr, double b, double berr) {
141	joshmt	1.9	using namespace std;
142	joshmt	1.10	if (b==0) {
143			cout<<"Warning in errAoverB -- b is zero!"<<endl;
144	joshmt	1.3	return -1;
145			}
146	joshmt	1.10	return (1/b)sqrt( aerraerr + aaberrberr/(bb) );
147	joshmt	1.3	}
148
149			//return error on a*b
150			double errAtimesB(double a, double aerr, double b, double berr) {
151	joshmt	1.10	return sqrt( bbaerraerr + aaberrberr);
152	joshmt	1.3	}
153
154	joshmt	1.2	//======== container for run, lumisection, event number =========
155			class eventID {
156			public:
157			eventID();
158			eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber);
159			~eventID();
160
161			bool operator< (const eventID & id) const;
162			bool operator== (const eventID & id) const;
163			bool operator!= (const eventID & id) const;
164
165			ULong64_t run;
166			ULong64_t ls;
167			ULong64_t ev;
168
169			};
170
171			eventID::eventID() : run(0),ls(0),ev(0) {}
172			eventID::eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber) :
173			run(RunNumber),ls(LumiSection),ev(EventNumber) {}
174			eventID::~eventID() {}
175
176			bool eventID::operator== (const eventID & id) const {
177
178			if (ev == id.ev &&
179			ls == id.ls &&
180			run == id.run) return true;
181
182			return false;
183			}
184
185			bool eventID::operator!= (const eventID & id) const {
186			return !(*this == id);
187			}
188
189			bool eventID::operator< (const eventID & id) const {
190
191			if (run < id.run) return true;
192			else if (run > id.run) return false;
193			else { //if run is equal
194
195			//now compare lumi section
196			if ( ls < id.ls ) return true;
197			else if (ls > id.ls) return false;
198			else { //if ls is equal
199
200			if ( ev < id.ev ) return true;
201			else if (ev > id.ev) return false;
202
203			else { //these are the same event!
204			return false;
205			}
206			}
207			}
208
209			//this shouldn't happen
210			assert(0);
211			return false;
212			}
213
214			} //end of namespace
215
216