FMD/scripts/VA1Response.C

   1 //____________________________________________________________________
   2 //
   3 // $Id$
   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=2, Float_t B=5, Float_t dc=.01, Bool_t errors=kFALSE)
  17 {
  18
  19   TF1* response = new TF1("response", "[0] * (1 - exp(-[1] * x))", 0, 1.4);
  20   response->SetParameters(1, B);
  21   response->SetParNames("A", "B");
  22   response->SetLineColor(2);
  23
  24   TF1* fit = new TF1("fit",  "[0] * (1 - exp(-[1] * x))",  0, 1);
  25   fit->SetParameters(.5, B/2);
  26   fit->SetParNames("A", "B");
  27   fit->SetLineColor(3);
  28
  29   TGraph* graph = 0;
  30   if (errors) graph = new TGraphErrors(n);
  31   else        graph = new TGraph(n);
  32   for (Int_t i = 0; i < n; i++) {
  33     Float_t t = Float_t(i + 1) / n;
  34     Float_t c = gRandom->Gaus(response->Eval(t), dc);
  35     graph->SetPoint(i, t, c);
  36     if (errors) ((TGraphErrors*)graph)->SetPointError(i, 0, dc);
  37   }
  38
  39   response->Draw();
  40   response->GetHistogram()->GetYaxis()->SetRangeUser(0, 1.4);
  41   response->GetHistogram()->GetXaxis()->SetRangeUser(0, 1.4);
  42   graph->Draw("P*");
  43   TString fitOpt("E");
  44   if (!errors) fitOpt.Append("W");
  45   graph->Fit("fit", fitOpt.Data());
  46   graph->Fit("fit", fitOpt.Data());
  47
  48   std::cout << "Chi^2/NDF = " << fit->GetChisquare() << "/" << fit->GetNDF()
  49             << " = " << std::flush;
  50   if (fit->GetNDF() == 0)
  51     std::cout << " undefined!" << std::endl;
  52   else
  53     std::cout << (fit->GetChisquare() / fit->GetNDF()) << std::endl;
  54   std::cout << "f(t) = "
  55             << fit->GetParameter(0) << "+/-" << fit->GetParError(0)
  56             << " * (1 - exp("
  57             << fit->GetParameter(1) << "+/-" << fit->GetParError(1)
  58             << " * t))" << std::endl;
  59 }
  60 //____________________________________________________________________
  61 //
  62 // EOF
  63 //