[BLAST_ANAWARE] GetAsymmetry

From: Jason Seely (seely@mit.edu)
Date: Fri May 21 2004 - 11:52:29 EDT


hey everyone,

for those of you doing vector asymmetries: i am submitting an upgraded version of GetAsymmetry() to ROOT, and i thought that people using it might want the new one, so i've attached it here. 

the only change is that now the function doesn't assume Poisson statistics for the histograms it uses to calculate the asymmetry.  (if things *are* poisson, then no change will be noticed).  this was done mainly just to make it more general, but now it will work correctly if you have, for example, histograms which are some signal+background, where the background is not poisson distributed (?).  anyway, if everything is poisson, then there's no reason to use the new version, but thought i should put it out there.

have fun,

jason

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