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 bool latexMode_=true) {
    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));
  }

  int signOf(double a) {
    return (a>= 0.0) ? 1 : -1;
  }

  // == 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);
  }

  //simple routines for addition in quadrature
  double addInQuad(double a, double b) {return sqrt(a*a + b*b);  }
  double addInQuad(double a, double b, double c) {return sqrt(a*a + b*b + c*c);  }
  double addInQuad(double a, double b, double c, double d) {return sqrt(a*a + b*b + c*c +d*d);  }

  //======== 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;
  }
  
  
  void printHist(const TH1D* h, int binlow=1, int binhigh=-1, bool withErrors=true) {
    if (binhigh == -1) binhigh = h->GetNbinsX();

    std::cout << h->GetName() << ": " ;
    for (int i = binlow; i<= binhigh; i++) {
      if(withErrors) std::cout << format_nevents(h->GetBinContent(i), h->GetBinError(i)) << ", ";
      else std::cout << " " << h->GetBinContent(i);
    }
    std::cout << std::endl;
  }

  
} //end of namespace


Revision:	1.16
Committed:	Fri May 31 09:54:57 2013 UTC (11 years, 11 months ago) by joshmt
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.15:	+1 -2 lines
Error occurred while calculating annotation data.
Log Message:	make latexMode a passable parameter
#	Content
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	#include "TH1D.h"
8	#include "TString.h"
9	#include <iostream>
10	#include "TMath.h"
11
12	#include <cassert>
13
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	cut.ReplaceAll(" ",""); //remove all spaces
21
22	cut.ReplaceAll("==","eq");
23	cut.ReplaceAll(">=","gte");
24	cut.ReplaceAll("<=","lte");
25
26	cut.ReplaceAll(">","gt");
27	cut.ReplaceAll("<","lt");
28
29	cut.ReplaceAll("&&","and");
30	cut.ReplaceAll("\|\|","or");
31
32	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	return cut;
41	}
42
43	//utility function for making output more readable
44	TString format_nevents(double n,double e, const bool moreDigits=false,const bool latexMode_=true) {
45	const TString pm = latexMode_ ? "\\pm ": " +/- ";
46
47	const int eCutoff = moreDigits ? 10 : 1;
48	const int extraDigits = moreDigits ? 1:0;
49
50	TString mathmode = latexMode_ ? "$" : "";
51
52	char out[100];
53	if (e >= eCutoff \|\| e < 0.00001) { //show whole numbers only
54	sprintf(out,"%s%.0f%s%.0f%s",mathmode.Data(),n,pm.Data(),e,mathmode.Data());
55	}
56	else {
57	int nfig = ceil(fabs(log10(e))) + extraDigits;
58	TString form="%s%.";
59	form+=nfig; form+="f%s%.";
60	form+=nfig; form+="f%s";
61	sprintf(out,form.Data(),mathmode.Data(),n,pm.Data(),e,mathmode.Data());
62	}
63	return TString(out);
64	}
65
66
67	Double_t weightedMean(Int_t n, Double_t a, Double_t e, const bool returnError=false) {
68	using namespace std;
69
70	Double_t numerator=0;
71	Double_t denominator=0;
72
73	for (Int_t i=0; i<n; ++i) {
74	Double_t esq = e[i]*e[i];
75	if (esq == 0) {cout<<"Error is zero!"<<endl; return -1e9;}
76	numerator += a[i] / esq;
77	denominator += 1.0 / esq;
78	}
79
80	cout<<"mean = "<<numerator / denominator<<" +/- "<<sqrt(1/denominator)<<endl;
81	return returnError ? sqrt(1/denominator) : numerator / denominator;
82	}
83
84	//implemented in TMath in later versions of root, so i'm using the TMath name
85	bool AreEqualAbs( const double & a, const double & b, const double epsilon=0.001) {
86
87	return ( fabs(a-b) < epsilon ) ;
88
89	}
90
91	//deal with the annoyance of logical booleans that are stored in an ntuple as double
92	bool doubleToBool( double a) {
93	using namespace std;
94
95	int i = TMath::Nint(a);
96	if (i==0) return false;
97	else if (i==1) return true;
98
99	cout<<"[doubleToBool] Something weird going on! "<<a<<"\t"<<i<<endl;
100	if (i>1) return true;
101	return false;
102	}
103
104	//delta phi calculations are pretty commonplace, so put it here
105	double deltaPhi(double phi1, double phi2) {
106	double result = phi1 - phi2;
107	while (result > TMath::Pi()) result -= 2*TMath::Pi();
108	while (result <= -TMath::Pi()) result += 2*TMath::Pi();
109	return fabs(result);
110	}
111
112	double deltaR2(double eta1, double phi1, double eta2, double phi2) {
113	double deta = eta1 - eta2;
114	double dphi = deltaPhi(phi1, phi2);
115	return detadeta + dphidphi;
116	}
117
118	double deltaR(double eta1, double phi1, double eta2, double phi2) {
119	return TMath::Sqrt(deltaR2 (eta1, phi1, eta2, phi2));
120	}
121
122	int signOf(double a) {
123	return (a>= 0.0) ? 1 : -1;
124	}
125
126	// == error propagation ==
127
128	//because not all versions of ROOT have it built in
129	double errOnIntegral(const TH1D* h, int binlow=1, int binhigh=-1) {
130
131	if (binhigh == -1) binhigh = h->GetNbinsX();
132
133	double err=0;
134	for (int i = binlow; i<= binhigh; i++) {
135	err += h->GetBinError(i) * h->GetBinError(i);
136	}
137
138	if (err<0) return err;
139	return sqrt(err);
140	}
141
142	//return error on a/b
143	double errAoverB(double a, double aerr, double b, double berr) {
144	using namespace std;
145	if (b==0) {
146	cout<<"Warning in errAoverB -- b is zero!"<<endl;
147	return -1;
148	}
149	return (1/b)sqrt( aerraerr + aaberrberr/(bb) );
150	}
151
152	//return error on a*b
153	double errAtimesB(double a, double aerr, double b, double berr) {
154	return sqrt( bbaerraerr + aaberrberr);
155	}
156
157	//simple routines for addition in quadrature
158	double addInQuad(double a, double b) {return sqrt(aa + bb); }
159	double addInQuad(double a, double b, double c) {return sqrt(aa + bb + c*c); }
160	double addInQuad(double a, double b, double c, double d) {return sqrt(aa + bb + cc +dd); }
161
162	//======== container for run, lumisection, event number =========
163	class eventID {
164	public:
165	eventID();
166	eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber);
167	~eventID();
168
169	bool operator< (const eventID & id) const;
170	bool operator== (const eventID & id) const;
171	bool operator!= (const eventID & id) const;
172
173	ULong64_t run;
174	ULong64_t ls;
175	ULong64_t ev;
176
177	};
178
179	eventID::eventID() : run(0),ls(0),ev(0) {}
180	eventID::eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber) :
181	run(RunNumber),ls(LumiSection),ev(EventNumber) {}
182	eventID::~eventID() {}
183
184	bool eventID::operator== (const eventID & id) const {
185
186	if (ev == id.ev &&
187	ls == id.ls &&
188	run == id.run) return true;
189
190	return false;
191	}
192
193	bool eventID::operator!= (const eventID & id) const {
194	return !(*this == id);
195	}
196
197	bool eventID::operator< (const eventID & id) const {
198
199	if (run < id.run) return true;
200	else if (run > id.run) return false;
201	else { //if run is equal
202
203	//now compare lumi section
204	if ( ls < id.ls ) return true;
205	else if (ls > id.ls) return false;
206	else { //if ls is equal
207
208	if ( ev < id.ev ) return true;
209	else if (ev > id.ev) return false;
210
211	else { //these are the same event!
212	return false;
213	}
214	}
215	}
216
217	//this shouldn't happen
218	assert(0);
219	return false;
220	}
221
222
223	void printHist(const TH1D* h, int binlow=1, int binhigh=-1, bool withErrors=true) {
224	if (binhigh == -1) binhigh = h->GetNbinsX();
225
226	std::cout << h->GetName() << ": " ;
227	for (int i = binlow; i<= binhigh; i++) {
228	if(withErrors) std::cout << format_nevents(h->GetBinContent(i), h->GetBinError(i)) << ", ";
229	else std::cout << " " << h->GetBinContent(i);
230	}
231	std::cout << std::endl;
232	}
233
234
235	} //end of namespace
236
237