fdpots.C
Go to the documentation of this file.
1 double clamp(const double low, const double x, const double high)
2 {
3  if(x < low) return low;
4  if(x > high) return high;
5  return x;
6 }
7 
8 double rounduptenths(const double x)
9 {
10  return int(x*10 + 1)/10.;
11 }
12 
13 double rounddowntenths(const double x)
14 {
15  return int(x*10)/10.;
16 }
17 
18 void pots(const char * const in, const char * out)
19 {
20  const double textsize = 0.05;
21  const unsigned int week = 7, year = 365, month = 30;
22 
23  TCanvas *c1 = new TCanvas("c1", "c1",0,0,800,500);
24  c1->SetCanvasSize(1000,500);
25  gStyle->SetOptStat(0);
26  gStyle->SetOptTitle(0);
27  c1->SetGrid(0,1);
28  gStyle->SetGridStyle(kSolid);
29  gStyle->SetGridColor(kGray);
30  c1->SetLeftMargin(0.08);
31  c1->SetBottomMargin(0.08);
32  c1->SetTopMargin(0.025);
33  c1->SetRightMargin(0.23);
34  c1->SetTicky(1);
35  c1->SetTickx(1);
36  c1->SetFrameLineWidth(2);
37 
38  ifstream infile(in);
39 
40  vector<string> date;
41  vector<double> uptimefrac, potdelivered, potrecorded, bins;
42 
43  struct datatmp{
44  string s;
45  double a, b, c;
46  };
47  datatmp tmp;
48 
49  int n = 0;
50  while(infile >> tmp.s >> tmp.a >> tmp.b >> tmp.c){
51  date .push_back(tmp.s);
52  uptimefrac .push_back(tmp.a);
53  potrecorded .push_back(tmp.b);
54  potdelivered.push_back(tmp.c);
55  bins .push_back(1 + n++);
56  }
57 
58  double maxpots = 1;
59  for(int i = 0; i < potdelivered.size(); i++)
60  if(potdelivered[i] > maxpots)
61  maxpots = potdelivered[i];
62 
63 
64  TGraph * delivered = new TGraph(date.size(), &bins[0], &potdelivered[0]);
65  delivered->SetName("delivered");
66  TGraph * recorded = new TGraph(date.size(), &bins[0], &potrecorded[0]);
67  recorded->SetName("recorded");
68 
69  enum mode{weekmode, monthmode, yearmode } themode;
70  if(bins.size() < month+1) themode = weekmode;
71  else if(bins.size() < month*4) themode = monthmode;
72  else themode = yearmode;
73  const double highy = rounduptenths(maxpots);
74 
75  const int nbin = themode == weekmode? (week*((bins.size() + week - 1)/ week)+1):
76  themode ==monthmode? (month*((bins.size() + month - 1)/month)+1):
77  (year*((bins.size() + year - 1)/ year)+1);
78 
79  TH2D * dum = new TH2D("dum", "", nbin, 0, nbin, 1, 0, highy);
80 
81  int div = themode == weekmode? 700 + (bins.size() + week - 1)/week:
82  themode ==monthmode? 300 + (bins.size() + month- 1)/month:
83  1200 + (bins.size() + year - 1)/year;
84  TGaxis * xa = new TGaxis(1, 0, nbin, 0, 0, nbin, div, "N");
85  TGaxis * xat= new TGaxis(1, highy, nbin, highy, 0, nbin, div, "N-");
86  xa->SetLabelSize(0);
87  xat->SetLabelSize(0);
88 
89  for(int i = 0; i < bins.size(); i += (themode == weekmode? week :
90  themode == monthmode?month: year))
91  // shift about half a bin right ---v
92  dum->GetXaxis()->SetBinLabel(i+1, Form(" %s", date[i].c_str()));
93 
94  dum->GetXaxis()->SetTitle("");
95  dum->GetYaxis()->SetTitle("10^{18} POT per day");
96  dum->GetYaxis()->SetTitleSize(textsize);
97  dum->GetXaxis()->SetTitleSize(textsize);
98  dum->GetYaxis()->SetLabelSize(textsize);
99  dum->GetXaxis()->SetLabelSize(textsize*1.5);
100  dum->GetXaxis()->SetLabelOffset(0.01);
101  dum->GetYaxis()->SetTitleOffset(0.70);
102  dum->GetXaxis()->SetBit(TAxis::kLabelsHori);
103  dum->GetYaxis()->CenterTitle(1);
104  dum->GetYaxis()->SetDecimals();
105  dum->GetYaxis()->SetTickSize(0.01);
106  dum->GetXaxis()->CenterTitle(1);
107  dum->GetXaxis()->SetTickSize(0);
108  dum->Draw();
109  xa->Draw();
110  xat->Draw();
111 
112  delivered->SetMarkerStyle(kOpenDiamond);
113  delivered->SetMarkerSize(themode == yearmode?0.3: 1.3);
114  delivered->SetMarkerColor(kBlue+1);
115  delivered->Draw("p");
116 
117  recorded->SetMarkerStyle(kMultiply);
118  recorded->SetMarkerSize(themode == yearmode?0.3: 1.3);
119  recorded->SetMarkerColor(kRed);
120  recorded->Draw("p");
121 
122  TGraph * deliveredbyweek = dynamic_cast<TGraph*>(delivered->Clone("deliveredbyweek"));
123 
124  for(int i = 0; i < delivered->GetN(); i+=week){
125  double sum = 0;
126  int n = 0;
127  for(int j = i; j < i+week && j < delivered->GetN(); j++){
128  sum += delivered->GetY()[j];
129  n++;
130  }
131  const double mean = sum/n;
132  for(int j = i; j < i+week && j < delivered->GetN(); j++){
133 
134  // Make transitions sharper
135  if(j == i) deliveredbyweek->GetX()[j] -= 0.25;
136  if(j == i+week-1 || j == delivered->GetN()-1) deliveredbyweek->GetX()[j] += 0.25;
137 
138  deliveredbyweek->GetY()[j] = mean;
139  }
140  }
141 
142  deliveredbyweek->SetLineColor(kBlue+1);
143  deliveredbyweek->SetLineWidth(themode == weekmode?3:2);
144  deliveredbyweek->Draw("l");
145 
146  TGraph * recordedbyweek = dynamic_cast<TGraph*>(recorded->Clone("recordedbyweek"));
147 
148  double lastweeksdelivered = 0, lastweeksrecorded = 0;
149  for(int i = 0; i < recorded->GetN(); i+=week){
150  double sumdem = 0, sumnum = 0, sum = 0;
151  double sumx = 0;
152  int n = 0;
153  for(int j = i; j < i+week && j < recorded->GetN(); j++){
154  sumnum += recorded->GetY()[j];
155  sumdem += delivered->GetY()[j];
156  sumx += recorded->GetX()[j];
157  n++;
158  }
159  const double meaneff = sumnum/sumdem;
160  const double mean = sumnum/n;
161  const double meanx = sumx/n;
162  for(int j = i; j < i+week && j < recorded->GetN(); j++){
163  // Make transitions sharper
164  if(j == i) recordedbyweek->GetX()[j] -= 0.25;
165  if(j == i+week-1 || j == recorded->GetN()-1) recordedbyweek->GetX()[j] += 0.25;
166 
167  recordedbyweek->GetY()[j] = mean;
168  }
169 
170  if(n == 7){ // save last full week
171  lastweeksdelivered = sumdem;
172  lastweeksrecorded = sumnum;
173  }
174 
175  if(recorded->GetN() <= week*9){
176  TText * t = meaneff!=meaneff?new TText(meanx, i%2?highy/8.:highy/6., ""):
177  new TText(meanx, i%2?highy/8.:highy/6., Form("%.1f%%", meaneff*100));
178  t->SetTextSize(textsize);
179  t->SetTextFont(42);
180  t->SetTextAlign(22);
181  t->SetTextColor(kRed);
182  t->Draw();
183  }
184  }
185 
186  recordedbyweek->SetLineColor(kRed);
187  recordedbyweek->SetLineWidth(themode == weekmode?3:2);
188  recordedbyweek->Draw("l");
189 
190  c1->SaveAs(Form("%s-noleg.pdf", out));
191  c1->SaveAs(Form("%s-noleg.eps", out));
192  c1->SaveAs(Form("%s-noleg.png", out));
193 
194  const double legbound = 0.45;
195 
196  TLegend leg(0.77, legbound, 0.995, 0.81);
197  leg.SetTextSize(textsize);
198  leg.SetBorderSize(0);
199  leg.SetMargin(0.19);
200  leg.SetFillStyle(0);
201  leg.AddEntry((TH1D*)NULL, out[0] == 'n'?"Near Detector":"Far Detector", "");
202  leg.AddEntry((TH1D*)NULL, "UTC days", "");
203  leg.AddEntry((TH1D*)NULL, "Sun-Sat weeks", "");
204  leg.AddEntry(delivered, "Delivered daily", "p");
205  leg.AddEntry(recorded, "Recorded daily", "p");
206  leg.AddEntry(deliveredbyweek, "Delivered weekly", "l");
207  leg.AddEntry(recordedbyweek, "Recorded weekly", "l");
208  leg.Draw();
209 
210  TLegend leg2(0.77, 0.24, 0.995, legbound-0.03);
211  leg2.SetTextSize(textsize);
212  leg2.SetBorderSize(0);
213  leg2.SetMargin(0.04);
214  leg2.SetFillStyle(0);
215  leg2.AddEntry((TH1D*)NULL, "Last full week", "");
216  leg2.AddEntry((TH1D*)NULL, Form("%.2f#times10^{18} delivered",lastweeksdelivered), "");
217  leg2.AddEntry((TH1D*)NULL, Form("%.2f#times10^{18} recorded", lastweeksrecorded), "");
218  leg2.Draw();
219 
220  c1->SaveAs(Form("%s.pdf", out));
221  c1->SaveAs(Form("%s.eps", out));
222  c1->SaveAs(Form("%s.png", out));
223 }
224 
225 void ndpots(const char * const more = "")
226 {
227  pots(Form("ndsane%s.dat", more), Form("ndpots%s", more));
228 }
229 
230 void fdpots(const char * const more = "")
231 {
232  pots(Form("fdsane%s.dat", more), Form("fdpots%s", more));
233 }
::xsd::cxx::tree::date< char, simple_type > date
Definition: Database.h:186
void ndpots(const char *const more="")
Definition: fdpots.C:225
double clamp(const double low, const double x, const double high)
Definition: fdpots.C:1
void pots(const char *const in, const char *out)
Definition: fdpots.C:18
double rounduptenths(const double x)
Definition: fdpots.C:8
Float_t tmp
Definition: plot.C:36
const XML_Char * s
Definition: expat.h:262
string infile
const double a
double rounddowntenths(const double x)
Definition: fdpots.C:13
const double j
Definition: BetheBloch.cxx:29
const Binning bins
ifstream in
Definition: comparison.C:7
const hit & b
Definition: hits.cxx:21
void fdpots(const char *const more="")
Definition: fdpots.C:230
c1
Definition: demo5.py:24
Double_t sum
Definition: plot.C:31
static constexpr Double_t year
Definition: Munits.h:185