-- -------------------------------------------------------------------- jason seely 26.650.b massachusetts institute of technology 77 massachusetts avenue cambridge, ma 02139-4307 email: seely@mit.edu phone: 617.253.4798/6734 html: web.mit.edu/seely/www --------------------------------------------------------------------
TH1 *TH1::GetAsymmetry(TH1* h2, Double_t c2, Double_t dc2)
{
// return an histogram containing the asymmetry of this histogram with h2,
// where the asymmetry is defined as:
//
// Asymmetry = (h1 - c2*h2)/(h1 + c2*h2) where h1 = this
//
// works for 1D, 2D, etc. histograms
//
// c2 is an optional argument that gives a relative weight between the two
// histograms (default is 1), and dc2 is the error on this weight (default is 0).
// This is useful, for example, when forming an asymmetry between two histograms
// from 2 different data sets that need to be normalized to each other in some way.
//
// example:
// assuming 'h1' and 'h2' are already filled, and 'h2'
//
// h3 = h1->GetAsymmetry(h2,c2,dc2)
//
// then 'h3' is created and filled with the asymmetry between 'h1' and 'h2';
// h1 and h2 are left intact.
//
// h3 = (h1 - a*h2)/(h1 + a*h2)
//
// Note that this function no longer assumes Poissson errors for 'h1' and 'h2'.
// So, for example, if 'h1' and 'h2' are sums of histograms ('counts + background'
// or something), then the asymmetry will now compute the proper uncertainty on
// the asymmetry. (Thanks to Christos Leonidopoulos for suggesting this fix.)
//
// Note that it is the user's responsibility to manage the created histogram.
//
// code proposed by Jason Seely (seely@mit.edu) and adapted by R.Brun
// clone the histograms so top and bottom will have the
// correct dimensions:
// Sumw2 just makes sure the errors will be computed properly
// when we form sums and ratios below.
Bool_t addStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
TH1 *asym = (TH1*)Clone();
asym->Sumw2();
TH1 *top = (TH1*)asym->Clone();
TH1 *bottom = (TH1*)asym->Clone();
TH1::AddDirectory(addStatus);
// form the top and bottom of the asymmetry, and then divide:
TH1 *h1 = this;
top->Add(h1,h2,1,-c2);
bottom->Add(h1,h2,1,c2);
asym->Divide(top,bottom);
Int_t xmax = asym->GetNbinsX();
Int_t ymax = asym->GetNbinsY();
Int_t zmax = asym->GetNbinsZ();
Float_t as, bot, error, a, b, da, db;
// now loop over bins to calculate the correct errors
// the reason this error calculation looks complex is because of c2
for(Int_t i=1; i<= xmax; i++){
for(Int_t j=1; j<= ymax; j++){
for(Int_t k=1; k<= zmax; k++){
// here some bin contents are written into variables to make the error
// calculation a little more legible:
a = h1->GetBinContent(i,j,k);
b = h2->GetBinContent(i,j,k);
da = h1->GetBinError(i,j,k);
db = h2->GetBinError(i,j,k);
as = asym->GetBinContent(i,j,k);
bot = bottom->GetBinContent(i,j,k);
// make sure there are some events, if not, then the errors are set = 0
// automatically.
if(bot < 1){}
else{
// remember: da = sqrt(a), and form the error:
error = TMath::Sqrt((1-as)*(1-as)*da*da
+(1+as)*(1+as)*(c2*c2*db*db+b*b*dc2*dc2))/bot;
asym->SetBinError(i,j,k,error);
}
}
}
}
delete top;
delete bottom;
return asym;
}
This archive was generated by hypermail 2.1.2 : Mon Feb 24 2014 - 14:07:31 EST