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(" ",""); //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) {

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