| 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 |
|
|
| 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 |
< |
|
| 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 |
< |
void weightedMean(Int_t n, Double_t *a, Double_t *e) { |
| 43 |
> |
//utility function for making output more readable |
| 44 |
> |
TString format_nevents(double n,double e, const bool moreDigits=false) { |
| 45 |
> |
const TString latexMode_=true; //hard code for now |
| 46 |
> |
const TString pm = latexMode_ ? "\\pm ": " +/- "; |
| 47 |
> |
|
| 48 |
> |
const int eCutoff = moreDigits ? 10 : 1; |
| 49 |
> |
const int extraDigits = moreDigits ? 1:0; |
| 50 |
> |
|
| 51 |
> |
TString mathmode = latexMode_ ? "$" : ""; |
| 52 |
> |
|
| 53 |
> |
char out[100]; |
| 54 |
> |
if (e >= eCutoff || e < 0.00001) { //show whole numbers only |
| 55 |
> |
sprintf(out,"%s%.0f%s%.0f%s",mathmode.Data(),n,pm.Data(),e,mathmode.Data()); |
| 56 |
> |
} |
| 57 |
> |
else { |
| 58 |
> |
int nfig = ceil(fabs(log10(e))) + extraDigits; |
| 59 |
> |
TString form="%s%."; |
| 60 |
> |
form+=nfig; form+="f%s%."; |
| 61 |
> |
form+=nfig; form+="f%s"; |
| 62 |
> |
sprintf(out,form.Data(),mathmode.Data(),n,pm.Data(),e,mathmode.Data()); |
| 63 |
> |
} |
| 64 |
> |
return TString(out); |
| 65 |
> |
} |
| 66 |
> |
|
| 67 |
> |
|
| 68 |
> |
Double_t weightedMean(Int_t n, Double_t *a, Double_t *e, const bool returnError=false) { |
| 69 |
> |
using namespace std; |
| 70 |
|
|
| 71 |
|
Double_t numerator=0; |
| 72 |
|
Double_t denominator=0; |
| 73 |
|
|
| 74 |
|
for (Int_t i=0; i<n; ++i) { |
| 75 |
|
Double_t esq = e[i]*e[i]; |
| 76 |
< |
if (esq == 0) {cout<<"Error is zero!"<<endl; return;} |
| 76 |
> |
if (esq == 0) {cout<<"Error is zero!"<<endl; return -1e9;} |
| 77 |
|
numerator += a[i] / esq; |
| 78 |
|
denominator += 1.0 / esq; |
| 79 |
|
} |
| 80 |
|
|
| 81 |
|
cout<<"mean = "<<numerator / denominator<<" +/- "<<sqrt(1/denominator)<<endl; |
| 82 |
< |
|
| 82 |
> |
return returnError ? sqrt(1/denominator) : numerator / denominator; |
| 83 |
|
} |
| 84 |
|
|
| 85 |
|
//implemented in TMath in later versions of root, so i'm using the TMath name |
| 89 |
|
|
| 90 |
|
} |
| 91 |
|
|
| 92 |
+ |
//deal with the annoyance of logical booleans that are stored in an ntuple as double |
| 93 |
+ |
bool doubleToBool( double a) { |
| 94 |
+ |
using namespace std; |
| 95 |
+ |
|
| 96 |
+ |
int i = TMath::Nint(a); |
| 97 |
+ |
if (i==0) return false; |
| 98 |
+ |
else if (i==1) return true; |
| 99 |
+ |
|
| 100 |
+ |
cout<<"[doubleToBool] Something weird going on! "<<a<<"\t"<<i<<endl; |
| 101 |
+ |
if (i>1) return true; |
| 102 |
+ |
return false; |
| 103 |
+ |
} |
| 104 |
+ |
|
| 105 |
+ |
//delta phi calculations are pretty commonplace, so put it here |
| 106 |
+ |
double deltaPhi(double phi1, double phi2) { |
| 107 |
+ |
double result = phi1 - phi2; |
| 108 |
+ |
while (result > TMath::Pi()) result -= 2*TMath::Pi(); |
| 109 |
+ |
while (result <= -TMath::Pi()) result += 2*TMath::Pi(); |
| 110 |
+ |
return fabs(result); |
| 111 |
+ |
} |
| 112 |
+ |
|
| 113 |
+ |
double deltaR2(double eta1, double phi1, double eta2, double phi2) { |
| 114 |
+ |
double deta = eta1 - eta2; |
| 115 |
+ |
double dphi = deltaPhi(phi1, phi2); |
| 116 |
+ |
return deta*deta + dphi*dphi; |
| 117 |
+ |
} |
| 118 |
+ |
|
| 119 |
+ |
double deltaR(double eta1, double phi1, double eta2, double phi2) { |
| 120 |
+ |
return TMath::Sqrt(deltaR2 (eta1, phi1, eta2, phi2)); |
| 121 |
+ |
} |
| 122 |
+ |
|
| 123 |
|
// == error propagation == |
| 124 |
|
|
| 125 |
|
//because not all versions of ROOT have it built in |
| 138 |
|
|
| 139 |
|
//return error on a/b |
| 140 |
|
double errAoverB(double a, double aerr, double b, double berr) { |
| 141 |
+ |
using namespace std; |
| 142 |
|
if (b==0 || a==0) { |
| 143 |
|
cout<<"Warning in errAoverB -- a or b is zero!"<<endl; |
| 144 |
|
return -1; |
| 148 |
|
|
| 149 |
|
//return error on a*b |
| 150 |
|
double errAtimesB(double a, double aerr, double b, double berr) { |
| 151 |
+ |
using namespace std; |
| 152 |
|
if (b==0 || a==0) { |
| 153 |
|
cout<<"Warning in errAtimesB -- a or b is zero!"<<endl; |
| 154 |
|
return -1; |
