Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PS_timing/scripts/fits.C
Go to the documentation of this file.
1 using namespace RooFit;
2 TH1* GetHistogram(TFile *f, TString hname, TString htype, int counter) {
3  stringstream ss; ss << counter;
4  TH2I *h2 = (TH2I*)f->Get("PS_timing/"+hname);
5  if (htype != "PSC" && htype != "PS" && htype != "TDCADC") {
6  cerr << "Unsupported histogram type: Choose 'PS', 'PSC', or 'TDCADC'" << endl;
7  }
8  TH1 *h = h2->ProjectionY(TString(h2->GetName())+"_"+TString(ss.str()),counter,counter);
9  TString title = (htype == "TDCADC") ? "PSC counter "+TString(ss.str()) : htype+" counter "+TString(ss.str());
10  h->SetTitle(title);
11  h->GetXaxis()->SetTitle(h2->GetYaxis()->GetTitle());
12  return h;
13 }
14 TH1* GetTaggerHistogram(TFile *f) {
15  TH2I *h2 = (TH2I*)f->Get("PS_timing/TAGHRF_tdcTimeDiffVsID");
16  const int N = 127; // Use upstream counters (relative to microscope)
17  TH1 *h = h2->ProjectionY(TString(h2->GetName()),1,N+1);
18  h->GetXaxis()->SetTitle("time(TDC) - time(RF) [ns]");
19  return h;
20 }
21 double GetMode(TH1 *h) {
22  if (h->GetEntries() == 0.0) return 999.0;
23  int max_bin = h->GetMaximumBin();
24  double max = h->GetBinContent(max_bin);
25  if (max < 7.0) return 999.0;
26  return h->GetBinCenter(max_bin);
27 }
28 void WriteFitResults(ofstream &fout, TH1 *h, TString htype, int counter) {
29  TString sep = ",";
30  double max = h->GetBinContent(h->GetMaximumBin());
31  double mode = GetMode(h);
32  if (h->GetEntries() < 100.0 || mode == 999.0) {
33  fout << counter << sep << max << sep << 0.0 << sep << 0.0 << sep << 0.0 << endl;
34  return;
35  }
36  double w = 0.75;
37  double xlow = mode - w; double xhigh = mode + w;
38  if (htype == "TDCADC") {
39  xlow = h->GetBinCenter(h->FindFirstBinAbove(0.5*max)) - w;
40  xhigh = h->GetBinCenter(h->FindLastBinAbove(0.5*max)) + w;
41  }
42  RooRealVar x("TDC time difference","TDC time difference [ns]",xlow,xhigh);
43  RooDataHist data("data","data",RooArgList(x),h);
44  // model: gaussian
45  RooRealVar mean("mean","mean",mode,xlow,xhigh);
46  RooRealVar sigma("sigma","sigma",0.2,0.01,0.5);
47  RooGaussian gauss("gauss","gauss",x,mean,sigma);
48  gauss.fitTo(data,PrintLevel(-1));
49  // Make plot
50  TCanvas *canvas = new TCanvas("c","c",800,500);
51  RooPlot *plot = x.frame();
52  plot->SetTitle(h->GetTitle());
53  plot->SetXTitle(h->GetXaxis()->GetTitle());
54  TString ytitle = (htype == "PS") ? "PS hits" : "PSC hits";
55  plot->SetYTitle(ytitle);
56  plot->SetTitleSize(0.0474,"XY");
57  data.plotOn(plot);
58  gauss.plotOn(plot,LineColor(kRed));
59  plot->Draw();
60  gauss.paramOn(plot,Layout(0.58,0.9,0.9),Format("NEU",AutoPrecision(1)));
61  plot->Draw();
62  stringstream ss; ss << counter;
63  system(TString("mkdir -p fits_"+htype).Data());
64  canvas->Print("fits_"+htype+"/counter_"+TString(ss.str())+".gif");
65  delete canvas; delete plot;
66  fout << counter << sep << max << sep << mean.getVal() << sep << mean.getError() << sep << sigma.getVal() << endl;
67 }
68 void WriteFitResults2(ofstream &fout, TH1 *h, TString htype, int counter) {
69  TString sep = ",";
70  double max = h->GetBinContent(h->GetMaximumBin());
71  double mode = GetMode(h);
72  if (h->GetEntries() < 100.0 || mode == 999.0) {
73  if (counter > 0) fout << counter << sep << max << sep << 0.0 << sep << 0.0 << sep << 0.0 << endl;
74  return;
75  }
76  double w = (htype == "TAGH") ? 1.0 : 0.75;
77  double xlow = mode - w;
78  double xhigh = mode + w;
79  RooRealVar x("TDC time difference","TDC time difference [ns]",xlow,xhigh);
80  RooDataHist data("data","data",RooArgList(x),h);
81  // model: add a narrow and wide gaussian with same mean
82  RooRealVar mean("mean","mean",mode,xlow,xhigh);
83  RooRealVar sigma1("sigma1","sigma1",0.2,0.01,0.4);//0.01,0.6
84  RooGaussian gauss1("gauss1","gauss1",x,mean,sigma1);
85  RooRealVar sigma2("sigma2","sigma2",0.7,0.3,2.5);//0.1,2.5
86  RooGaussian gauss2("gauss2","gauss2",x,mean,sigma2);
87  // f1: fraction of entries in first gaussian
88  RooRealVar f1("f1","f1",0.75,0.01,1.);
89  RooAddPdf doubleGauss("doubleGauss","doubleGauss",RooArgList(gauss1,gauss2),RooArgList(f1));
90  doubleGauss.fitTo(data,PrintLevel(-1));
91  // Make plot
92  TCanvas *canvas = new TCanvas("c","c",800,500);
93  RooPlot *plot = x.frame();
94  plot->SetTitle(h->GetTitle());
95  plot->SetXTitle(h->GetXaxis()->GetTitle());
96  TString ytitle = (htype == "PS") ? "PS hits" : "PSC hits";
97  if (htype == "TAGH") ytitle = "TAGH hits";
98  plot->SetYTitle(ytitle);
99  plot->SetTitleSize(0.0474,"XY");
100  data.plotOn(plot);
101  doubleGauss.plotOn(plot,LineColor(kRed));
102  plot->Draw();
103  doubleGauss.plotOn(plot,Components("gauss1"),LineColor(kBlue),LineStyle(kDashed));
104  doubleGauss.plotOn(plot,Components("gauss2"),LineColor(kGreen),LineStyle(kDashed));
105  doubleGauss.paramOn(plot,Layout(0.15,0.45,0.9),Format("NEU",AutoPrecision(1)));
106  plot->Draw();
107  stringstream ss; ss << counter;
108  system(TString("mkdir -p fits_"+htype).Data());
109  canvas->Print("fits_"+htype+"/counter_"+TString(ss.str())+".gif");
110  delete canvas; delete plot;
111  if (counter > 0) fout << counter << sep << max << sep << mean.getVal() << sep << mean.getError() << sep << sigma1.getVal() << endl;
112 }
113 int fits(TString rootFile, bool doAllFits) {
114  TFile *f = new TFile(rootFile,"read");
115  TString hnames[] = {"PSC_tdcadcTimeDiffVsID","PSCRF_tdcTimeDiffVsID","PSRF_adcTimeDiffVsID"};
116  TString htypes[] = {"TDCADC","PSC","PS"};
117  int N_htypes = 2;
118  if (doAllFits) N_htypes = 3;
119  for (int htype = 0; htype < N_htypes; htype++) {
120  system("mkdir -p fits-csv");
121  ofstream fout; fout.open("fits-csv/results_"+htypes[htype]+".txt");
122  const int N = (htypes[htype] == "PS") ? 290 : 16;
123  for (int i = 1; i <= N; i++) { // counters
124  TH1 *h = GetHistogram(f,hnames[htype],htypes[htype],i);
125  if (htypes[htype] != "PSC") WriteFitResults(fout,h,htypes[htype],i);
126  else WriteFitResults2(fout,h,htypes[htype],i);
127  }
128  fout.close();
129  }
130  ofstream fout; fout.open("fits-csv/results_base.txt");
131  TH1 *h = GetTaggerHistogram(f);
132  WriteFitResults2(fout,h,"TAGH",1);
133  fout.close();
134  return 0;
135 }
int fits(TString rootFile, bool doAllFits)
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
void WriteFitResults2(ofstream &fout, TH1 *h, TString htype, int counter)
TF1 * f
Definition: FitGains.C:21
double counter
Definition: FitGains.C:151
TFile f1(fnam)
Double_t sigma[NCHANNELS]
Definition: st_tw_resols.C:37
TH1 * GetTaggerHistogram(TFile *f)
double GetMode(TH1 *h)
TH1 * GetHistogram(TFile *f, TString hname, TString htype, int counter)
void WriteFitResults(ofstream &fout, TH1 *h, TString htype, int counter)