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>

namespace jmt {

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

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

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

    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) {

    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 ) ;

  }

  // == 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) {
    if (b==0 || a==0) {
      cout<<"Warning in errAoverB -- a or b is zero!"<<endl;
      return -1;
    }
    return (a/b)*sqrt( pow( aerr/a, 2) + pow( berr/b, 2) );
  }

  //return error on a*b
  double errAtimesB(double a, double aerr, double b, double berr) {
    if (b==0 || a==0) {
      cout<<"Warning in errAtimesB -- a or b is zero!"<<endl;
      return -1;
    }
    
    return a*b*sqrt( pow( aerr/a, 2) + pow( berr/b, 2) );
  }

  //======== 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.7
Committed:	Wed Jul 27 13:50:47 2011 UTC (13 years, 9 months ago) by joshmt
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+3 -3 lines
Log Message:	another update
#	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.1
11			namespace jmt {
12
13			//======misc utilities======
14			//gets rid of = > < from cuts in order to be better included in file names
15			TString fortranize(TString cut) {
16
17			cut.ReplaceAll("==","eq");
18			cut.ReplaceAll(">=","gte");
19			cut.ReplaceAll("<=","lte");
20
21			cut.ReplaceAll(">","gt");
22			cut.ReplaceAll("<","lt");
23	joshmt	1.4
24			cut.ReplaceAll("/","Over");
25			cut.ReplaceAll("+","Plus");
26			cut.ReplaceAll("*","Times");
27
28			cut.ReplaceAll("*","Times");
29
30			//this is pretty ugly
31			cut.ReplaceAll("(","L");
32			cut.ReplaceAll(")","R");
33
34	joshmt	1.1	return cut;
35			}
36
37	joshmt	1.6	//utility function for making output more readable
38			TString format_nevents(double n,double e, const bool moreDigits=false) {
39			const TString latexMode_=true; //hard code for now
40			const TString pm = latexMode_ ? "\\pm ": " +/- ";
41
42			const int eCutoff = moreDigits ? 10 : 1;
43			const int extraDigits = moreDigits ? 1:0;
44
45			TString mathmode = latexMode_ ? "$" : "";
46
47			char out[100];
48			if (e >= eCutoff \|\| e < 0.00001) { //show whole numbers only
49			sprintf(out,"%s%.0f%s%.0f%s",mathmode.Data(),n,pm.Data(),e,mathmode.Data());
50			}
51			else {
52			int nfig = ceil(fabs(log10(e))) + extraDigits;
53			TString form="%s%.";
54			form+=nfig; form+="f%s%.";
55			form+=nfig; form+="f%s";
56			sprintf(out,form.Data(),mathmode.Data(),n,pm.Data(),e,mathmode.Data());
57			}
58			return TString(out);
59			}
60
61
62	joshmt	1.7	Double_t weightedMean(Int_t n, Double_t a, Double_t e, const bool returnError=false) {
63	joshmt	1.3
64			Double_t numerator=0;
65			Double_t denominator=0;
66
67			for (Int_t i=0; i<n; ++i) {
68			Double_t esq = e[i]*e[i];
69	joshmt	1.7	if (esq == 0) {cout<<"Error is zero!"<<endl; return -1e9;}
70	joshmt	1.3	numerator += a[i] / esq;
71			denominator += 1.0 / esq;
72			}
73
74			cout<<"mean = "<<numerator / denominator<<" +/- "<<sqrt(1/denominator)<<endl;
75	joshmt	1.7	return returnError ? sqrt(1/denominator) : numerator / denominator;
76	joshmt	1.3	}
77
78			//implemented in TMath in later versions of root, so i'm using the TMath name
79			bool AreEqualAbs( const double & a, const double & b, const double epsilon=0.001) {
80
81			return ( fabs(a-b) < epsilon ) ;
82
83			}
84
85			// == error propagation ==
86
87			//because not all versions of ROOT have it built in
88			double errOnIntegral(const TH1D* h, int binlow=1, int binhigh=-1) {
89
90			if (binhigh == -1) binhigh = h->GetNbinsX();
91
92			double err=0;
93			for (int i = binlow; i<= binhigh; i++) {
94			err += h->GetBinError(i) * h->GetBinError(i);
95			}
96
97			if (err<0) return err;
98			return sqrt(err);
99			}
100
101			//return error on a/b
102			double errAoverB(double a, double aerr, double b, double berr) {
103			if (b==0 \|\| a==0) {
104			cout<<"Warning in errAoverB -- a or b is zero!"<<endl;
105			return -1;
106			}
107			return (a/b)*sqrt( pow( aerr/a, 2) + pow( berr/b, 2) );
108			}
109
110			//return error on a*b
111			double errAtimesB(double a, double aerr, double b, double berr) {
112			if (b==0 \|\| a==0) {
113			cout<<"Warning in errAtimesB -- a or b is zero!"<<endl;
114			return -1;
115			}
116
117			return absqrt( pow( aerr/a, 2) + pow( berr/b, 2) );
118			}
119
120	joshmt	1.2	//======== container for run, lumisection, event number =========
121			class eventID {
122			public:
123			eventID();
124			eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber);
125			~eventID();
126
127			bool operator< (const eventID & id) const;
128			bool operator== (const eventID & id) const;
129			bool operator!= (const eventID & id) const;
130
131			ULong64_t run;
132			ULong64_t ls;
133			ULong64_t ev;
134
135			};
136
137			eventID::eventID() : run(0),ls(0),ev(0) {}
138			eventID::eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber) :
139			run(RunNumber),ls(LumiSection),ev(EventNumber) {}
140			eventID::~eventID() {}
141
142			bool eventID::operator== (const eventID & id) const {
143
144			if (ev == id.ev &&
145			ls == id.ls &&
146			run == id.run) return true;
147
148			return false;
149			}
150
151			bool eventID::operator!= (const eventID & id) const {
152			return !(*this == id);
153			}
154
155			bool eventID::operator< (const eventID & id) const {
156
157			if (run < id.run) return true;
158			else if (run > id.run) return false;
159			else { //if run is equal
160
161			//now compare lumi section
162			if ( ls < id.ls ) return true;
163			else if (ls > id.ls) return false;
164			else { //if ls is equal
165
166			if ( ev < id.ev ) return true;
167			else if (ev > id.ev) return false;
168
169			else { //these are the same event!
170			return false;
171			}
172			}
173			}
174
175			//this shouldn't happen
176			assert(0);
177			return false;
178			}
179
180			} //end of namespace
181
182