Thanks.

Code attached.

Welcome to Geeks to Go - Register now for FREE

**Geeks To Go** is a helpful hub, where thousands of volunteer geeks quickly serve friendly answers and support. Check out the forums and get free advice from the experts. Register now to gain access to all of our features, it's FREE and only takes one minute. Once registered and logged in, you will be able to create topics, post replies to existing threads, give reputation to your fellow members, get your own private messenger, post status updates, manage your profile and so much more.

Started by
quintessenceanx
, Oct 30 2006 10:08 AM

Posted 30 October 2006 - 10:08 AM

Thanks.

Code attached.

Posted 30 October 2006 - 10:11 AM

Since the code doesn't seem to attach, included below:

**EffecFit.C:**

{

gROOT->Reset();

gSystem->Load("libPhysics");

gROOT->LoadMacro("FitEffec2.C");

c1 = new TCanvas("c1","Fitted Effeciency",200,10,700,500);

c1->SetFillColor(42);

c1->SetGrid();

// draw a frame to define the range

TH1F *hr = c1->DrawFrame(1600,0,1850,1);

hr->SetYTitle("Effeciency");

hr->SetXTitle("Voltage");

c1->GetFrame()->SetFillColor(21);

c1->GetFrame()->SetBorderSize(12);

//histogram= new TH1F("histogram","",100,1600,1850,1);

// create first graph

Int_t n1 = 11;

Double_t x1[] = {1600, 1620,1640,1660,1680,1700,1720,1740,1760,1780,1800};

Double_t y1[] = {.313,.380,.488,.592,.682,.780,.825,.865,.900,.884,.892};

Double_t ex1[] = {0,0,0,0,0,0,0,0,0,0,0};

Double_t ey1[] = {.013,.025,.04,.056,0.07,0.08,0.017,0.086,0.078,0.071,0.079};

TF1 *fiteffec1 = new TF1("fiteffec",fiteffec,0.,1850,4);

gr1 = new TGraphErrors(n1,x1,y1,ex1,ey1);

gr1->SetMarkerColor(kBlue);

gr1->SetMarkerStyle(21);

gr1->Draw("P");

gr1->Fit(fiteffec1,"RE");

}

**FitEffec2.C**

Double_t fiteffec(Double_t *x, Double_t *par)

{

if (x[0]<1600) return 0.;

Double_t my_erf = par[0]+par[1]*TMath::Erf(x[0]*par[2]-par[3]);

return my_erf;

// Double_t fit = par[0]*TMath::Exp(-0.5*arg*arg)+par[3]+par[4]*x[0];

void TGraphAsymmErrors::BayesDivide(const TH1 *pass, const TH1 *total, Option_t *option)

{

TString opt = option; opt.ToLower();

Int_t nbins = pass->GetNbinsX();

if (nbins != total->GetNbinsX()){

Error("BayesDivide","Histograms must have the same number of X bins");

return;

}

if (opt.Contains("w")) {

//compare sum of weights with sum of squares of weights

Double_t stats[10];

pass->GetStats(stats);

if (TMath::Abs(stats[0] -stats[1]) > 1e-6) {

Error("BayesDivide","Pass histogram has not been filled with weights = 1");

return;

}

total->GetStats(stats);

if (TMath::Abs(stats[0] -stats[1]) > 1e-6) {

Error("BayesDivide","Total histogram has not been filled with weights = 1");

return;

}

}

//Set the graph to have a number of points equal to the number of histogram bins

Set(nbins);

// Ok, now set the points for each bin

// (Note: the TH1 bin content is shifted to the right by one:

// bin=0 is underflow, bin=nbins+1 is overflow.)

double mode, low, high; //these will hold the result of the Bayes calculation

int npoint=0;//this keeps track of the number of points added to the graph

for (int b=1; b<=nbins; ++b) { // loop through the bins

int t = (int)total->GetBinContent(b);

if (!t) continue; //don't add points for bins with no information

int p = (int)pass->GetBinContent(b);

if (p>t) {

Warning("BayesDivide","Histogram bin %d in pass has more entries than corresponding bin in total! (%d>%d)",b,p,t);

continue; //we may as well go on...

}

//This is the Bayes calculation...

Efficiency(p,t,0.683,mode,low,high);

//These are the low and high error bars

low = mode-low;

high = high-mode;

//If either of the errors are 0, set them to 1/10 of the other error

//so that the fitters don't get confused.

if (low==0.0) low=high/10.;

if (high==0.0) high=low/10.;

//Set the point center and its errors

SetPoint(npoint,pass->GetBinCenter(b),mode);

SetPointError(npoint,

pass->GetBinCenter(b)-pass->GetBinLowEdge(b),

pass->GetBinLowEdge(b)-pass->GetBinCenter(b)+pass->GetBinWidth(b),

low,high);

npoint++;//we have added a point to the graph

}

Set(npoint);//tell the graph how many points we've really added

if (opt.Contains("debug")) {

printf("BayesDivide: made a graph with %d points from %d bins\n",npoint,nbins);

Print();//The debug prints out what we get for each point

}

}

}

{

gROOT->Reset();

gSystem->Load("libPhysics");

gROOT->LoadMacro("FitEffec2.C");

c1 = new TCanvas("c1","Fitted Effeciency",200,10,700,500);

c1->SetFillColor(42);

c1->SetGrid();

// draw a frame to define the range

TH1F *hr = c1->DrawFrame(1600,0,1850,1);

hr->SetYTitle("Effeciency");

hr->SetXTitle("Voltage");

c1->GetFrame()->SetFillColor(21);

c1->GetFrame()->SetBorderSize(12);

//histogram= new TH1F("histogram","",100,1600,1850,1);

// create first graph

Int_t n1 = 11;

Double_t x1[] = {1600, 1620,1640,1660,1680,1700,1720,1740,1760,1780,1800};

Double_t y1[] = {.313,.380,.488,.592,.682,.780,.825,.865,.900,.884,.892};

Double_t ex1[] = {0,0,0,0,0,0,0,0,0,0,0};

Double_t ey1[] = {.013,.025,.04,.056,0.07,0.08,0.017,0.086,0.078,0.071,0.079};

TF1 *fiteffec1 = new TF1("fiteffec",fiteffec,0.,1850,4);

gr1 = new TGraphErrors(n1,x1,y1,ex1,ey1);

gr1->SetMarkerColor(kBlue);

gr1->SetMarkerStyle(21);

gr1->Draw("P");

gr1->Fit(fiteffec1,"RE");

}

Double_t fiteffec(Double_t *x, Double_t *par)

{

if (x[0]<1600) return 0.;

Double_t my_erf = par[0]+par[1]*TMath::Erf(x[0]*par[2]-par[3]);

return my_erf;

// Double_t fit = par[0]*TMath::Exp(-0.5*arg*arg)+par[3]+par[4]*x[0];

void TGraphAsymmErrors::BayesDivide(const TH1 *pass, const TH1 *total, Option_t *option)

{

TString opt = option; opt.ToLower();

Int_t nbins = pass->GetNbinsX();

if (nbins != total->GetNbinsX()){

Error("BayesDivide","Histograms must have the same number of X bins");

return;

}

if (opt.Contains("w")) {

//compare sum of weights with sum of squares of weights

Double_t stats[10];

pass->GetStats(stats);

if (TMath::Abs(stats[0] -stats[1]) > 1e-6) {

Error("BayesDivide","Pass histogram has not been filled with weights = 1");

return;

}

total->GetStats(stats);

if (TMath::Abs(stats[0] -stats[1]) > 1e-6) {

Error("BayesDivide","Total histogram has not been filled with weights = 1");

return;

}

}

//Set the graph to have a number of points equal to the number of histogram bins

Set(nbins);

// Ok, now set the points for each bin

// (Note: the TH1 bin content is shifted to the right by one:

// bin=0 is underflow, bin=nbins+1 is overflow.)

double mode, low, high; //these will hold the result of the Bayes calculation

int npoint=0;//this keeps track of the number of points added to the graph

for (int b=1; b<=nbins; ++b) { // loop through the bins

int t = (int)total->GetBinContent(b);

if (!t) continue; //don't add points for bins with no information

int p = (int)pass->GetBinContent(b);

if (p>t) {

Warning("BayesDivide","Histogram bin %d in pass has more entries than corresponding bin in total! (%d>%d)",b,p,t);

continue; //we may as well go on...

}

//This is the Bayes calculation...

Efficiency(p,t,0.683,mode,low,high);

//These are the low and high error bars

low = mode-low;

high = high-mode;

//If either of the errors are 0, set them to 1/10 of the other error

//so that the fitters don't get confused.

if (low==0.0) low=high/10.;

if (high==0.0) high=low/10.;

//Set the point center and its errors

SetPoint(npoint,pass->GetBinCenter(b),mode);

SetPointError(npoint,

pass->GetBinCenter(b)-pass->GetBinLowEdge(b),

pass->GetBinLowEdge(b)-pass->GetBinCenter(b)+pass->GetBinWidth(b),

low,high);

npoint++;//we have added a point to the graph

}

Set(npoint);//tell the graph how many points we've really added

if (opt.Contains("debug")) {

printf("BayesDivide: made a graph with %d points from %d bins\n",npoint,nbins);

Print();//The debug prints out what we get for each point

}

}

}

0 members, 0 guests, 0 anonymous users

Community Forum Software by IP.Board

Licensed to: Geeks to Go, Inc.