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
Log Message:	make latexMode a passable parameter
#	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	joshmt	1.16	TString format_nevents(double n,double e, const bool moreDigits=false,const bool latexMode_=true) {
45	joshmt	1.6	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	joshmt	1.7	Double_t weightedMean(Int_t n, Double_t a, Double_t e, const bool returnError=false) {
68	joshmt	1.9	using namespace std;
69	joshmt	1.3
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	joshmt	1.7	if (esq == 0) {cout<<"Error is zero!"<<endl; return -1e9;}
76	joshmt	1.3	numerator += a[i] / esq;
77			denominator += 1.0 / esq;
78			}
79
80			cout<<"mean = "<<numerator / denominator<<" +/- "<<sqrt(1/denominator)<<endl;
81	joshmt	1.7	return returnError ? sqrt(1/denominator) : numerator / denominator;
82	joshmt	1.3	}
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	joshmt	1.9	//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	joshmt	1.15	int signOf(double a) {
123			return (a>= 0.0) ? 1 : -1;
124			}
125
126	joshmt	1.3	// == 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	joshmt	1.9	using namespace std;
145	joshmt	1.10	if (b==0) {
146			cout<<"Warning in errAoverB -- b is zero!"<<endl;
147	joshmt	1.3	return -1;
148			}
149	joshmt	1.10	return (1/b)sqrt( aerraerr + aaberrberr/(bb) );
150	joshmt	1.3	}
151
152			//return error on a*b
153			double errAtimesB(double a, double aerr, double b, double berr) {
154	joshmt	1.10	return sqrt( bbaerraerr + aaberrberr);
155	joshmt	1.3	}
156
157	joshmt	1.11	//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	joshmt	1.2	//======== 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	kreis	1.12
223			void printHist(const TH1D* h, int binlow=1, int binhigh=-1, bool withErrors=true) {
224			if (binhigh == -1) binhigh = h->GetNbinsX();
225	kreis	1.13
226			std::cout << h->GetName() << ": " ;
227	kreis	1.12	for (int i = binlow; i<= binhigh; i++) {
228	kreis	1.13	if(withErrors) std::cout << format_nevents(h->GetBinContent(i), h->GetBinError(i)) << ", ";
229			else std::cout << " " << h->GetBinContent(i);
230	kreis	1.12	}
231			std::cout << std::endl;
232			}
233	kreis	1.13
234	kreis	1.12
235	joshmt	1.2	} //end of namespace
236
237