FMD/doc/VA1Response.C

   1 //____________________________________________________________________
   2 //
   3 // $Id: VA1Response.C 13249 2006-03-24 16:09:36Z cholm $
   4 //
   5 // Script to try to fit the reponse function of the VA1 signals, based
   6 // on a finite number of ALTRO samples.
   7 //
   8 /** Make Va1 response
   9     @param n
  10     @param B
  11     @param dc
  12     @param errors
  13     @ingroup simple_script
  14 */
  15 void
  16 VA1Response(Int_t n=4, Float_t B=6, Float_t dc=.01, Bool_t errors=kFALSE,
  17             Bool_t doFit=kFALSE)
  18 {
  19   gStyle->SetOptTitle(0);
  20   gStyle->SetOptStat(0);
  21   gStyle->SetOptFit(0);
  22   gStyle->SetLabelFont(132, "xyz");
  23   gStyle->SetTitleFont(132, "xyz");
  24   gStyle->SetTitleSize(0.08, "y");
  25   gStyle->SetTitleOffset(0.5, "y");
  26   gStyle->SetTitleSize(0.06, "x");
  27   gStyle->SetTitleOffset(0.7, "x");
  28
  29   TCanvas* c = new TCanvas("c", "c", 800, 500);
  30   c->SetFillColor(0);
  31   c->SetBorderMode(0);
  32   c->SetBorderSize(0);
  33   c->SetTopMargin(0.05);
  34   c->SetRightMargin(0.05);
  35   c->SetGridx();
  36   c->SetGridy();
  37
  38   TF1* response = new TF1("response", "[0] * (1 - exp(-[1] * x))", 0, 1.4);
  39   response->SetParameters(1, B);
  40   response->SetParNames("A", "B");
  41   response->SetLineColor(2);
  42
  43   TGraph* graph = 0;
  44   if (n >= 2) {
  45     if (errors) graph = new TGraphErrors(n);
  46     else        graph = new TGraph(n);
  47     for (Int_t i = 0; i < n; i++) {
  48       Float_t t = Float_t(i + 1) / n;
  49       Float_t q = gRandom->Gaus(response->Eval(t), dc);
  50       graph->SetPoint(i, t, q);
  51       if (errors) ((TGraphErrors*)graph)->SetPointError(i, 0, dc);
  52     }
  53   }
  54
  55   response->Draw();
  56   response->GetHistogram()->GetYaxis()->SetRangeUser(0, 1.05);
  57   response->GetHistogram()->GetXaxis()->SetRangeUser(0, 1.1);
  58   response->GetHistogram()->GetXaxis()->SetNdivisions(6, kTRUE);
  59   response->GetHistogram()->GetYaxis()->SetNdivisions(10, kTRUE);
  60   response->GetHistogram()->SetXTitle("t");
  61   response->GetHistogram()->SetYTitle(Form("1-e^{-%3.1f t}", B));
  62
  63   if (graph) {
  64
  65
  66     graph->Draw("P*");
  67     TString fitOpt("E");
  68     if (!errors) fitOpt.Append("W");
  69
  70     if (doFit) {
  71       TF1* fit = new TF1("fit",  "[0] * (1 - exp(-[1] * x))",  0, 1);
  72       fit->SetParameters(.5, B/2);
  73       fit->SetParNames("A", "B");
  74       fit->SetLineColor(3);
  75       graph->Fit("fit", fitOpt.Data());
  76       graph->Fit("fit", fitOpt.Data());
  77
  78       std::cout << "Chi^2/NDF = " << fit->GetChisquare() << "/" << fit->GetNDF()
  79                 << " = " << std::flush;
  80       if (fit->GetNDF() == 0)
  81         std::cout << " undefined!" << std::endl;
  82       else
  83         std::cout << (fit->GetChisquare() / fit->GetNDF()) << std::endl;
  84       std::cout << "f(t) = "
  85                 << fit->GetParameter(0) << "+/-" << fit->GetParError(0)
  86                 << " * (1 - exp("
  87                 << fit->GetParameter(1) << "+/-" << fit->GetParError(1)
  88                 << " * t))" << std::endl;
  89     }
  90   }
  91
  92   c->Modified();
  93   c->Update();
  94   c->cd();
  95   c->SaveAs("va1_response.png");
  96 }
  97 //____________________________________________________________________
  98 //
  99 // EOF
 100 //