StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Plot6.cc
1 /*
2  Author:David Kapukchyan
3  @[August 26, 2022]
4  > First instance
5 
6  A macro to plot the results from testing the StFcsWaveformFitMaker methods using StFcsWaveformFitMaker set to test level 6. The "prefix" is to set the output save file names prefix.
7 
8  How to determine the sum8 vs. Fit scale factor.
9 
10  1. Run StFcsWaveformFitMaker on a microdst file `root4star -b -q testWff6MuDst.C'(nevents,"OutFileName.root","InFile.MuDst.root")'`
11  - Before running check to make sure that the line that sets the energy sum scale reads *wff->setEnergySumScale(1.0,1.0,1.0)*
12  2. Take output root file "OutFileName.root" and run this macro giving it the file name and a prefix for the saved image files
13  - `root -b -q Plot6.C'("OutFileName.root","TestScale1")'`
14  3. Look at *TestScale1_sum8vfit.png* to check what values to scale by. slope of fitted line is the value to use
15  4. Change the line in *testWff6MuDst.C* that sets the energy sum scale to the values you get from the slopes of the fit.
16  - *wff->setEnergySumScale(slopecal,slopehcal,slopepres)*
17  5. Run StFcsWaveformFitMaker on the same microdst file with same number of events `root4star -b -q testWff6MuDst.C'(nevents,"OutFileNameNew.root","InFile.MuDst.root")'`
18  - You can change the output file name if you want
19  6. Go to step 2 and repeat for new output file
20 
21 */
22 
23 //#include "Rtools.h"
24 //#include "HistColl2.h"
25 
26 void Plot6(std::string filename="test.root",std::string prefix="test")
27 {
28  gStyle->SetOptStat(111111);
29  gStyle->SetOptDate(0);
30 
31  TFile* infile = TFile::Open(filename.c_str());
32  if( infile==0 || infile->IsZombie() ){ std::cout << "ERROR:Unable to open: " << filename << std::endl; return 0; }
33  else{ std::cout << "Opened file:"<<infile << " Named:"<<infile->GetName() << std::endl; }
34 
35  TH1F* npeaks[7] = {0};
36  TH1F* npeaksfiltered[7] = {0};
37  TH1F* res0[7] = {0};
38  TH1F* res0zoom[7] = {0};
39  TH1F* sum8res0[7] = {0};
40  TH1F* sum8res0zoom[7] = {0};
41  TH1F* fitres0[7] = {0};
42  TH1F* fitres0zoom[7] = {0};
43  TH2F* sum8vfit[7] = {0};
44  for( UInt_t i=0; i<7; ++i ){
45  char name[100];
46 
47  sprintf(name,"H1_NPeaks_%d",i);
48  npeaks[i] = (TH1F*)infile->Get(name);
49  npeaks[i]->SetTitle(";number of peaks");
50 
51  sprintf(name,"H1_NPeaksFiltered_%d",i);
52  npeaksfiltered[i] = (TH1F*)infile->Get(name);
53  npeaksfiltered[i]->SetTitle(";number of peaks");
54 
55  sprintf(name,"H1_Res0_%d",i);
56  res0[i] = (TH1F*)infile->Get(name);
57  res0[i]->SetTitle(";ADC sum");
58 
59  sprintf(name,"H1_Res0Zoom_%d",i);
60  res0zoom[i] = (TH1F*)infile->Get(name);
61  res0zoom[i]->SetTitle(";ADC sum");
62 
63  sprintf(name,"H1_Sum8Res0_%d",i);
64  sum8res0[i] = (TH1F*)infile->Get(name);
65  sum8res0[i]->SetTitle(";ADC sum");
66 
67  sprintf(name,"H1_Sum8Res0Zoom_%d",i);
68  sum8res0zoom[i] = (TH1F*)infile->Get(name);
69  sum8res0zoom[i]->SetTitle(";ADC sum");
70 
71  sprintf(name,"H1_FitRes0_%d",i);
72  fitres0[i] = (TH1F*)infile->Get(name);
73  fitres0[i]->SetTitle(";ADC sum");
74 
75  sprintf(name,"H1_FitRes0Zoom_%d",i);
76  fitres0zoom[i] = (TH1F*)infile->Get(name);
77  fitres0zoom[i]->SetTitle(";ADC sum");
78 
79  sprintf(name,"H2_Sum8vFit_%d",i);
80  sum8vfit[i] = (TH2F*)infile->Get(name);
81  sum8vfit[i]->SetTitle(";Sum8;FitSum");
82 
83  res0[i] ->SetLineColor(kBlack);
84  sum8res0[i] ->SetLineColor(kBlue);
85  fitres0[i] ->SetLineColor(kGreen);
86  res0zoom[i] ->SetLineColor(kBlack);
87  sum8res0zoom[i]->SetLineColor(kBlue);
88  fitres0zoom[i] ->SetLineColor(kGreen);
89  res0[i] ->SetStats(0);
90  sum8res0[i] ->SetStats(0);
91  fitres0[i] ->SetStats(0);
92  res0zoom[i] ->SetStats(0);
93  sum8res0zoom[i]->SetStats(0);
94  fitres0zoom[i] ->SetStats(0);
95  }
96  res0[6]->SetTitle("");
97  sum8res0[6]->SetTitle("");
98  fitres0[6]->SetTitle("");
99  res0zoom[6]->SetTitle("");
100  sum8res0zoom[6]->SetTitle("");
101  fitres0zoom[6]->SetTitle("");
102 
103  //Rtools::HistColl1F* npeaks = new Rtools::HistColl1F(infile,"NPeaks",";number of peaks");
104  //Rtools::HistColl1F* npeaksfiltered = new Rtools::HistColl1F(infile,"NPeaksFiltered",";number of peaks");
105  //Rtools::HistColl1F* res0 = new Rtools::HistColl1F(infile,"Res0",";ADC sum");
106  //Rtools::HistColl1F* res0zoom = new Rtools::HistColl1F(infile,"Res0Zoom", ";ADC sum");
107  //Rtools::HistColl1F* sum8res0 = new Rtools::HistColl1F(infile,"Sum8Res0", ";ADC sum");
108  //Rtools::HistColl1F* sum8res0zoom = new Rtools::HistColl1F(infile,"Sum8Res0Zoom",";ADC sum");
109  //Rtools::HistColl1F* fitres0 = new Rtools::HistColl1F(infile,"FitRes0",";ADC sum");
110  //Rtools::HistColl1F* fitres0zoom = new Rtools::HistColl1F(infile,"FitRes0Zoom",";ADC sum");
111 
112  //Rtools::HistColl2F* sum8vfit = new Rtools::HistColl2F(infile,"Sum8vFit",";Sum8;FitSum");
113  //sum8vfit->Print();
114 
115  TH1F* timefit = (TH1F*)infile->Get("H1_TimeFitPulse");
116  TH1F* measuretime = (TH1F*)infile->Get("FitTime");
117 
118  int width = 1280;
119  int height = 960;
120 
121  //Rtools::DrawTools* dt0 = new Rtools::DrawTools("dt0","",width,height);
122  TCanvas* canv = new TCanvas("canv","",width,height);
123  npeaks[6]->SetTitle("");
124  npeaks[6]->SetStats(0);
125  npeaks[6]->GetXaxis()->SetTitle("Found number of peaks");
126  npeaks[6]->Draw("hist e");
127  canv->SaveAs( (prefix+"_Npeaks.png").c_str() );
128 
129  canv->Clear();
130  npeaksfiltered[6]->Draw("hist e");
131  canv->SaveAs( (prefix+"_NpeaksFiltered.png").c_str() );
132 
133  canv->Clear();
134  canv->Divide(3,2);
135  for( UInt_t i=0; i<6; ++i ){
136  canv->cd(i+1);
137  npeaks[i]->Draw("hist e");
138  }
139  canv->SaveAs( (prefix+"_detNpeaks.png").c_str() );
140 
141  canv->Clear();
142  canv->Divide(3,2);
143  for( UInt_t i=0; i<6; ++i ){
144  canv->cd(i+1);
145  npeaksfiltered[i]->Draw("hist e");
146  }
147  canv->SaveAs( (prefix+"_detNpeaksFiltered.png").c_str() );
148 
149  //Draw sum8 result, fitted result, and final result on top of each other
150  canv->Clear();
151  //TPad* pad = dt0->DrawCanv(width,height);
152  canv->Divide(2,1);
153  canv->cd(1)->SetLogy(1);
154  res0[6]->Draw("hist e");
155  res0[6]->GetXaxis()->SetTitle("Signal Integral");
156  sum8res0[6]->Draw("hist e sames");
157  fitres0[6]->Draw("hist e sames");
158  canv->cd(2)->SetLogy(1);
159  res0zoom[6]->Draw("hist e");
160  res0zoom[6]->GetXaxis()->SetTitle("Signal Integral");
161  sum8res0zoom[6]->Draw("hist e sames");
162  fitres0zoom[6]->Draw("hist e sames");
163  TLegend* legres = 0;//new TLegend(0.5,0.5,0.99,0.99,"Sum Comparison","nbNDC");
164  //Rtools::AddHistStats(legres,res0->At(6));
165  //Rtools::AddHistStats(legres,sum8res0->At(6));
166  //Rtools::AddHistStats(legres,fitres0->At(6));
167  //legres->Draw();
168  canv->SaveAs( (prefix+"_res0.png").c_str() );
169 
170  canv->Clear();
171  canv->Divide(3,2);
172  for( UInt_t i=0; i<6; ++i ){
173  canv->cd(i+1)->SetLogy(1);
174  res0[i]->Draw("hist e");
175  sum8res0[i]->Draw("hist e same");
176  fitres0[i]->Draw("hist e same");
177  }
178  canv->SaveAs( (prefix+"_detres0.png").c_str() );
179 
180  canv->Clear();
181  canv->Divide(3,2);
182  for( UInt_t i=0; i<6; ++i ){
183  canv->cd(i+1)->SetLogy(1);
184  res0zoom[i]->Draw("hist e");
185  sum8res0zoom[i]->Draw("hist e same");
186  fitres0zoom[i]->Draw("hist e same");
187  }
188  canv->SaveAs( (prefix+"_detres0zoom.png").c_str() );
189 
190  //Draw just sum8 and fitted result on top of each other
191  canv->Clear();
192  canv->Divide(2,1);
193  canv->cd(1)->SetLogy(1);
194  sum8res0[6]->Draw("hist e");
195  sum8res0[6]->GetXaxis()->SetTitle("Signal Integral");
196  fitres0[6]->Draw("hist e same");
197  fitres0[6]->GetXaxis()->SetTitle("Signal Integral");
198  canv->cd(2)->SetLogy(1);
199  sum8res0zoom[6]->Draw("hist e");
200  sum8res0zoom[6]->GetXaxis()->SetTitle("Signal Integral");
201  fitres0zoom[6]->Draw("hist e same");
202  fitres0zoom[6]->GetXaxis()->SetTitle("Signal Integral");
203  TLegend* legres2 = 0;//new TLegend(0.5,0.5,0.99,0.99,"Comparing Sum Methods","nbNDC");
204  //Rtools::AddHistStats(legres2,sum8res0->At(6));
205  //Rtools::AddHistStats(legres2,fitres0->At(6));
206  //legres2->Draw();
207  canv->SaveAs( (prefix+"_s8fit.png").c_str() );
208 
209  canv->Clear();
210  TF1* linfit[7] = {0};
211  linfit[6] = new TF1("linfit","[0]+[1]*x",0,2000);
212  canv->SetLogz(1);
213  sum8vfit[6]->Fit(linfit[6]);
214  sum8vfit[6]->Draw("colz");
215  canv->SaveAs( (prefix+"_sum8vfit.png").c_str() );
216 
217  canv->Clear();
218  canv->Divide(3,2);
219  for( UInt_t i=0; i<6; ++i ){
220  linfit[i] = new TF1("linfit","[0]+[1]*x",0,2000);
221  canv->cd(i+1)->SetLogz(1);
222  sum8vfit[i]->Fit(linfit[i]);
223  sum8vfit[i]->Draw("colz");
224  }
225  canv->SaveAs( (prefix+"_detsum8vfit.png").c_str() );
226 
227 
228  canv->Clear();
229  timefit->SetLineColor(kBlue);
230  timefit->GetXaxis()->SetRangeUser(0,100);
231  timefit->GetXaxis()->SetTitle("time (ms)");
232  timefit->Draw("hist e");
233  double timeintegral = timefit->Integral(1,timefit->GetNbinsX());
234  gPad->Update();//h1 should be on the current pad otherwise won't work
235  TPaveStats* TPS_timefit = timefit->FindObject("stats");
236  TPS_timefit->SetName("timefitstats");
237  //TPaveStats* TPS_timefit = Rtools::GetStats(timefit,111111);
238  //std::cout<< "|TPS:"<<TPS_timefit<< "|TPS_Name:"<<TPS_timefit->GetName() << "|S:"<< TPS_timefit->GetListOfLines()->GetEntries() << std::endl;
239  timefit->SetStats(0);
240  std::stringstream ss_tinttext;
241  ss_tinttext << "Integral = " << timeintegral;
242  TLatex* ttimeint = new TLatex(0,0, ss_tinttext.str().c_str() );
243  //ttimeint->SetTextFont(42);
244  ttimeint->SetTextSize(0.03);
245  TPS_timefit->GetListOfLines()->Add(ttimeint);
246  TPS_timefit->Draw();
247  canv->SaveAs( (prefix+"_timefit.png").c_str() );
248 
249  if( measuretime!=0 ){
250  canv->Clear();
251  canv->SetLogy();
252  measuretime->SetLineColor(kBlue);
253  //measuretime->GetXaxis()->SetRangeUser(0,100);
254  measuretime->GetXaxis()->SetTitle("time (ms)");
255  measuretime->Draw("hist e");
256  double totaltime = measuretime->Integral(1,measuretime->GetNbinsX());
257  gPad->Update();//h1 should be on the current pad otherwise won't work
258  TPaveStats* TPS_measuretime = measuretime->FindObject("stats");
259  TPS_measuretime->SetName("measuretimestats");
260  measuretime->SetStats(0);
261  std::stringstream ss_minttext;
262  ss_minttext << "Integral = " << totaltime;
263  TLatex* mtimeint = new TLatex(0,0, ss_minttext.str().c_str() );
264  //mtimeint->SetTextFont(42);
265  mtimeint->SetTextSize(0.03);
266  TPS_measuretime->GetListOfLines()->Add(mtimeint);
267  TPS_measuretime->Draw();
268  canv->SaveAs( (prefix+"_measuretime.png").c_str() );
269  }
270 
271  /*
272  delete npeaks;
273  delete npeaksfiltered;
274  delete res0;
275  delete res0zoom;
276  delete sum8res0;
277  delete sum8res0zoom;
278  delete fitres0;
279  delete fitres0zoom;
280  delete sum8vfit;
281  delete timefit;
282  delete measuretime;
283  delete legres;
284  delete legres2;
285  */
286  infile->Close();
287  delete infile;
288  delete canv;
289 
290 }
291 
Generated by   doxygen 1.8.5