| 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.3 |
void weightedMean(Int_t n, Double_t *a, Double_t *e) {
|
| 38 |
|
|
|
| 39 |
|
|
Double_t numerator=0;
|
| 40 |
|
|
Double_t denominator=0;
|
| 41 |
|
|
|
| 42 |
|
|
for (Int_t i=0; i<n; ++i) {
|
| 43 |
|
|
Double_t esq = e[i]*e[i];
|
| 44 |
|
|
if (esq == 0) {cout<<"Error is zero!"<<endl; return;}
|
| 45 |
|
|
numerator += a[i] / esq;
|
| 46 |
|
|
denominator += 1.0 / esq;
|
| 47 |
|
|
}
|
| 48 |
|
|
|
| 49 |
|
|
cout<<"mean = "<<numerator / denominator<<" +/- "<<sqrt(1/denominator)<<endl;
|
| 50 |
|
|
|
| 51 |
|
|
}
|
| 52 |
|
|
|
| 53 |
|
|
//implemented in TMath in later versions of root, so i'm using the TMath name
|
| 54 |
|
|
bool AreEqualAbs( const double & a, const double & b, const double epsilon=0.001) {
|
| 55 |
|
|
|
| 56 |
|
|
return ( fabs(a-b) < epsilon ) ;
|
| 57 |
|
|
|
| 58 |
|
|
}
|
| 59 |
|
|
|
| 60 |
|
|
// == error propagation ==
|
| 61 |
|
|
|
| 62 |
|
|
//because not all versions of ROOT have it built in
|
| 63 |
|
|
double errOnIntegral(const TH1D* h, int binlow=1, int binhigh=-1) {
|
| 64 |
|
|
|
| 65 |
|
|
if (binhigh == -1) binhigh = h->GetNbinsX();
|
| 66 |
|
|
|
| 67 |
|
|
double err=0;
|
| 68 |
|
|
for (int i = binlow; i<= binhigh; i++) {
|
| 69 |
|
|
err += h->GetBinError(i) * h->GetBinError(i);
|
| 70 |
|
|
}
|
| 71 |
|
|
|
| 72 |
|
|
if (err<0) return err;
|
| 73 |
|
|
return sqrt(err);
|
| 74 |
|
|
}
|
| 75 |
|
|
|
| 76 |
|
|
//return error on a/b
|
| 77 |
|
|
double errAoverB(double a, double aerr, double b, double berr) {
|
| 78 |
|
|
if (b==0 || a==0) {
|
| 79 |
|
|
cout<<"Warning in errAoverB -- a or b is zero!"<<endl;
|
| 80 |
|
|
return -1;
|
| 81 |
|
|
}
|
| 82 |
|
|
return (a/b)*sqrt( pow( aerr/a, 2) + pow( berr/b, 2) );
|
| 83 |
|
|
}
|
| 84 |
|
|
|
| 85 |
|
|
//return error on a*b
|
| 86 |
|
|
double errAtimesB(double a, double aerr, double b, double berr) {
|
| 87 |
|
|
if (b==0 || a==0) {
|
| 88 |
|
|
cout<<"Warning in errAtimesB -- a or b is zero!"<<endl;
|
| 89 |
|
|
return -1;
|
| 90 |
|
|
}
|
| 91 |
|
|
|
| 92 |
|
|
return a*b*sqrt( pow( aerr/a, 2) + pow( berr/b, 2) );
|
| 93 |
|
|
}
|
| 94 |
|
|
|
| 95 |
joshmt |
1.2 |
//======== container for run, lumisection, event number =========
|
| 96 |
|
|
class eventID {
|
| 97 |
|
|
public:
|
| 98 |
|
|
eventID();
|
| 99 |
|
|
eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber);
|
| 100 |
|
|
~eventID();
|
| 101 |
|
|
|
| 102 |
|
|
bool operator< (const eventID & id) const;
|
| 103 |
|
|
bool operator== (const eventID & id) const;
|
| 104 |
|
|
bool operator!= (const eventID & id) const;
|
| 105 |
|
|
|
| 106 |
|
|
ULong64_t run;
|
| 107 |
|
|
ULong64_t ls;
|
| 108 |
|
|
ULong64_t ev;
|
| 109 |
|
|
|
| 110 |
|
|
};
|
| 111 |
|
|
|
| 112 |
|
|
eventID::eventID() : run(0),ls(0),ev(0) {}
|
| 113 |
|
|
eventID::eventID(ULong64_t RunNumber, ULong64_t LumiSection, ULong64_t EventNumber) :
|
| 114 |
|
|
run(RunNumber),ls(LumiSection),ev(EventNumber) {}
|
| 115 |
|
|
eventID::~eventID() {}
|
| 116 |
|
|
|
| 117 |
|
|
bool eventID::operator== (const eventID & id) const {
|
| 118 |
|
|
|
| 119 |
|
|
if (ev == id.ev &&
|
| 120 |
|
|
ls == id.ls &&
|
| 121 |
|
|
run == id.run) return true;
|
| 122 |
|
|
|
| 123 |
|
|
return false;
|
| 124 |
|
|
}
|
| 125 |
|
|
|
| 126 |
|
|
bool eventID::operator!= (const eventID & id) const {
|
| 127 |
|
|
return !(*this == id);
|
| 128 |
|
|
}
|
| 129 |
|
|
|
| 130 |
|
|
bool eventID::operator< (const eventID & id) const {
|
| 131 |
|
|
|
| 132 |
|
|
if (run < id.run) return true;
|
| 133 |
|
|
else if (run > id.run) return false;
|
| 134 |
|
|
else { //if run is equal
|
| 135 |
|
|
|
| 136 |
|
|
//now compare lumi section
|
| 137 |
|
|
if ( ls < id.ls ) return true;
|
| 138 |
|
|
else if (ls > id.ls) return false;
|
| 139 |
|
|
else { //if ls is equal
|
| 140 |
|
|
|
| 141 |
|
|
if ( ev < id.ev ) return true;
|
| 142 |
|
|
else if (ev > id.ev) return false;
|
| 143 |
|
|
|
| 144 |
|
|
else { //these are the same event!
|
| 145 |
|
|
return false;
|
| 146 |
|
|
}
|
| 147 |
|
|
}
|
| 148 |
|
|
}
|
| 149 |
|
|
|
| 150 |
|
|
//this shouldn't happen
|
| 151 |
|
|
assert(0);
|
| 152 |
|
|
return false;
|
| 153 |
|
|
}
|
| 154 |
|
|
|
| 155 |
|
|
} //end of namespace
|
| 156 |
|
|
|
| 157 |
|
|
|
