ChanPlots.C
Go to the documentation of this file.
1 #include <TSQLServer.h>
2 #include <TSQLResult.h>
3 #include <iostream.h>
4 #include <fstream>
5 #include <iomanip.h>
6 #include <sstream>
7 #include <string>
8 #include <cstring>
9 #include <algorithm>
10 
11 void ChanPlots(int chan) {
12 
13  gROOT->Reset();
14  gStyle->SetOptStat(0);
15  gStyle->SetOptFit(1);
16  gStyle->SetOptTitle(1);
17  gStyle->SetPalette(1);
18  gStyle->SetFillColor(0);
19  gStyle->SetPadColor(0);
20  gStyle->SetCanvasColor(0);
21  gStyle->SetStatColor(0);
22  gStyle->SetTitleColor(1);
23  gStyle->SetTitleSize(0.06);
24  gStyle->SetPadBorderMode(0);
25  gStyle->SetPadBottomMargin(0.15);
26  gStyle->SetPadLeftMargin(0.12);
27  gStyle->SetPadRightMargin(0.09);
28  gStyle->SetFrameBorderMode(0);
29  gStyle->SetCanvasBorderMode(0);
30 
31  TDatime da(2010,01,01,1,00,00);
32  gStyle->SetTimeOffset(da.Convert());
33  int offtime = da.Convert();
34 
35  TDatime *TNow = new TDatime;//.TDatime();
36  int XNow=TNow->Convert();
37  unsigned int X2DaysAgo =XNow - 2*60*60*24;
38  unsigned int X4WeeksAgo =XNow - 60*60*24*7*4;
39  TDatime *T2d = new TDatime;
40  T2d->Set(X2DaysAgo,false);
41  TDatime *T4w = new TDatime;
42  T4w->Set(X4WeeksAgo,false);
43 
44  TNow->Print();
45  cout<<XNow<<"\t"<<X2DaysAgo<<"\t"<<60*60*24*2<<"\t"<<XNow-X2DaysAgo<<endl;
46  cout<<"DATES: " <<TNow->AsSQLString()<< " " << T2d->AsSQLString()<< " " << T4w->AsSQLString() << endl;
47 
48  //2d = 2 days, hour bins = 48 bins, + underflow and overflow = 50 bins
49 
50  double pe_2d[50] = {0.}; //pe per hit
51  double pecorr_2d[50] = {0.}; // pecorr per hit
52  double nn_2d[50] = {0.}; //nearest neighbor
53  double rate_hi_2d[50] = {0.}; //rate > 100 pe
54  double rate_low_2d[50] = {0.}; // rate < 100 pe
55  int nhit_2d[50] = {0.}; // number of hits
56  int count_2d[50] = {0.}; // number of temperature 'hits'
57  double tave_2d[50] = {0.}; // average temperature
58  double pesq_2d[50] = {0.}; // pe squared per hit (for error)
59  double pecorrsq_2d[50] = {0.}; // pecorr squared per hit (for error)
60  int telapse_2d[50] = {0.};
61 
62  //4w = 4 weeks, day bins = 28 days + under and overflow = 30 bins
63 
64  double pe_4w[30] = {0.};
65  double pecorr_4w[30] = {0.};
66  double nn_4w[30] = {0.};
67  double rate_hi_4w[30] = {0.};
68  double rate_low_4w[30] = {0.};
69  int nhit_4w[30] = {0.};
70  int count_4w[30] = {0.};
71  double tave_4w[30] = {0.};
72  double pesq_4w[30] = {0.};
73  double pecorrsq_4w[30] = {0.};
74  int telapse_4w[30] = {0.};
75 
76  int tnh, ttn, r1, r2, r4, telapse;
77  double tpe, tpec, tnn, trh, trl, ttave, tpesq, tpecorrsq; // temporary variables for DB
78 
79  r1 = XNow-offtime;
80  r2 = X2DaysAgo-offtime;
81  r4 = X4WeeksAgo-offtime;
82 
83  cout << " Time ranges: " << r1 << " " << r2 << " " << r4 << endl;
84 
85  TH1F *pe_time_2d = new TH1F("pe_time_2d", "pe_hist;time;# p.e. per hit", 48, r2, r1);
86  TH1F *pec_time_2d = new TH1F("pec_time_2d", "pecorr_hist;time;corrected # p.e. per hit (if calibrated)", 48, r2, r1);
87  TH1F *nn_time_2d = new TH1F("nn_time_2d", "nn_hist;time;nearest neighbor variable", 48, r2, r1);
88  TH1F *rh_time_2d = new TH1F("rh_time_2d", "ratehi_hist;time;# hits > 100 p.e. per second", 48, r2, r1);
89  TH1F *rl_time_2d = new TH1F("rl_time_2d", "ratelo_hist;time;# hits < 100 p.e. per second", 48, r2, r1);
90  TH1F *apd_temp_2d = new TH1F("apd_temp_2d", "apd_temp;time;degrees C of APD", 48, r2, r1);
91 
92  TH1F *pe_time_4w = new TH1F("pe_time_4w", "pe_hist;time;# p.e. per hit", 28, r4, r1);
93  TH1F *pec_time_4w = new TH1F("pec_time_4w", "pecorr_hist;time;corrected # p.e. per hit (if calibrated)", 28, r4, r1);
94  TH1F *nn_time_4w = new TH1F("nn_time_4w", "nn_hist;time;nearest neighbor variable", 28, r4, r1);
95  TH1F *rh_time_4w = new TH1F("rh_time_4w", "ratehi_hist;time;number of hits > 100 p.e. per second", 28, r4, r1);
96  TH1F *rl_time_4w = new TH1F("rl_time_4w", "ratelo_hist;time;number of hits < 100 p.e. per second", 28, r4, r1);
97  TH1F *apd_temp_4w = new TH1F("apd_temp_4w", "apd_temp;time;degrees C of APD", 28, r4, r1);
98 
99  std::stringstream dbstream;
100 
101  dbstream << "select channel_info.pe,channel_info.pecorr,channel_info.rate_low,channel_info.rate_hi,"
102  << " channel_info.nhits,channel_info.nn_rate,"
103  << " apd_header.time_last,temp_info.temp_ave,apd_header.num_temp,"
104  << " channel_info.pesq,channel_info.pecorrsq,apd_header.time_elapsed from"
105  << " channel_info,daq_info,apd_header,temp_info where "
106  << " daq_info.dcm=temp_info.dcm"
107  << " and daq_info.diblock=temp_info.diblock"
108  << " and daq_info.port=temp_info.port"
109  << " and channel_info.channel=daq_info.channel"
110  << " and channel_info.channel="<<chan
111  << " and apd_header.run=channel_info.run"
112  << " and apd_header.subrun=channel_info.subrun"
113  << " and apd_header.run=temp_info.run"
114  << " and apd_header.subrun=temp_info.subrun"
115  << " and apd_header.time_first> '"<<T2d->AsSQLString()<<"' "
116  << " and apd_header.time_last<'"<<TNow->AsSQLString()<<" ' "
117  << " order by apd_header.run";
118 
119  std::string dbstring = dbstream.str();
120 
121  char * dbcstr;
122  dbcstr = new char [dbstring.size()+1];
123  strcpy(dbcstr, dbstring.c_str());
124 
125  //password removed
126  TSQLServer *db = TSQLServer::Connect("pgsql://131.215.116.171", "", "");
127 
128  std::cout << "Connected? " << db->IsConnected() << std::endl;
129 
130  TSQLStatement *stmt = db->Statement(dbcstr,1000);
131 
132  std::cout << " about to send statement: " << dbcstr << std::endl;
133 
134  if(stmt->Process()){
135  stmt->StoreResult();
136 
137  while (stmt->NextResultRow())
138  {
139  tpe = stmt->GetDouble(0);
140  tpec = stmt->GetDouble(1);
141  trl = stmt->GetDouble(2);
142  trh = stmt->GetDouble(3);
143  tnh = stmt->GetInt(4);
144  tnn = stmt->GetDouble(5);
145  const char* time_last = stmt->GetString(6);
146  ttave = stmt->GetDouble(7);
147  ttn = stmt->GetInt(8);
148  tpesq = stmt->GetDouble(9);
149  tpecorrsq = stmt->GetDouble(10);
150  telapse = stmt->GetInt(11);
151 
152  TDatime *ltime = new TDatime;
153  ltime->Set(time_last);
154  unsigned int tlint = ltime->Convert();
155 
156  int tbin = (tlint-X2DaysAgo) / 3600;
157  if(tbin == 48) tbin = 47;
158  tbin = tbin + 1;
159 
160  // cout << "2d data: " << time_last << " " << tpe << " " << tbin << endl;
161 
162  pe_2d[tbin] += tpe*tnh;
163  pecorr_2d[tbin] += tpec*tnh;
164  nn_2d[tbin] += tnn*telapse;
165  rate_hi_2d[tbin] += trh*telapse;
166  rate_low_2d[tbin] += trl*telapse;
167  nhit_2d[tbin] += tnh;
168  count_2d[tbin] += ttn;
169  tave_2d[tbin] += ttn*ttave;
170  pesq_2d[tbin] += tpesq*tnh;
171  pecorrsq_2d[tbin] += tpecorrsq*tnh;
172  telapse_2d[tbin] += telapse;
173  }
174 
175  double err1[50] = {0.};
176  double err2[50] = {0.};
177  double err3[50] = {0.};
178  double err4[50] = {0.};
179  double err5[50] = {0.};
180  double err6[50] = {0.};
181 
182  for(unsigned int n = 1; n < 49; ++n){
183 
184  if(nhit_2d[n]!=0) {
185 
186  err4[n] = sqrt(rate_hi_2d[n]) / telapse_2d[n];
187  err5[n] = sqrt(rate_low_2d[n]) / telapse_2d[n];
188 
189  pe_2d[n] = pe_2d[n] / nhit_2d[n];
190  pecorr_2d[n] = pecorr_2d[n] / nhit_2d[n];
191  pesq_2d[n] = pesq_2d[n] / nhit_2d[n];
192  pecorrsq_2d[n] = pecorrsq_2d[n] / nhit_2d[n];
193  nn_2d[n] = nn_2d[n]/nhit_2d[n];
194  rate_hi_2d[n] = rate_hi_2d[n]/telapse_2d[n];
195  rate_low_2d[n] = rate_low_2d[n]/telapse_2d[n];
196 
197  err3[n] = ( 1 / sqrt(nhit_2d[n]) ) * nn_2d[n];
198 
199  if(pesq_2d[n]==0) {
200  err1[n] = ( 1 / sqrt(nhit_2d[n]) ) * pe_2d[n];
201  err2[n] = ( 1 / sqrt(nhit_2d[n]) ) * pecorr_2d[n];
202  }
203  else if(pesq_2d[n]>0) {
204  err1[n] = sqrt( pesq_2d[n] - (pe_2d[n]*pe_2d[n]) ) * 1.5 / sqrt(nhit_2d[n]) ;
205  err2[n] = sqrt( pecorrsq_2d[n] - (pecorr_2d[n]*pecorr_2d[n]) ) * 1.5 / sqrt(nhit_2d[n]);
206  }
207  else {
208  std::cout << "Negative squares??? Something is wrong!" << std::endl;
209  }
210 
211  std::cout << "Values: " << pe_2d[n] << " " << pesq_2d[n] << " " << nhit_2d[n] << " " << err1[n] << std::endl;
212 
213  }
214  if(count_2d[n]!=0) {
215  tave_2d[n] = tave_2d[n]/count_2d[n];
216  err6[n] = ( sqrt(count_2d[n]) / count_2d[n] ) * tave_2d[n];
217  }
218 
219  }
220 
221  // cout << "about to try setting content. made it this far????" << endl;
222 
223 
224  pe_time_2d->SetContent(pe_2d);
225  pec_time_2d->SetContent(pecorr_2d);
226  nn_time_2d->SetContent(nn_2d);
227  rh_time_2d->SetContent(rate_hi_2d);
228  rl_time_2d->SetContent(rate_low_2d);
229  apd_temp_2d->SetContent(tave_2d);
230 
231  pe_time_2d->SetError(err1);
232  pec_time_2d->SetError(err2);
233  nn_time_2d->SetError(err3);
234  rh_time_2d->SetError(err4);
235  rl_time_2d->SetError(err5);
236  apd_temp_2d->SetError(err6);
237 
238  delete stmt;
239 
240  }
241 
242  delete db;
243 
244 
245  std::stringstream dbstream2;
246 
247  dbstream2 << "select channel_info.pe,channel_info.pecorr,channel_info.rate_low,channel_info.rate_hi,"
248  << " channel_info.nhits,channel_info.nn_rate,"
249  << " apd_header.time_last,temp_info.temp_ave,apd_header.num_temp,"
250  << " channel_info.pesq,channel_info.pecorrsq,apd_header.time_elapsed from"
251  << " channel_info,daq_info,apd_header,temp_info where "
252  << " daq_info.dcm=temp_info.dcm"
253  << " and daq_info.diblock=temp_info.diblock"
254  << " and daq_info.port=temp_info.port"
255  << " and channel_info.channel=daq_info.channel"
256  << " and channel_info.channel="<<chan
257  << " and apd_header.run=channel_info.run"
258  << " and apd_header.subrun=channel_info.subrun"
259  << " and apd_header.run=temp_info.run"
260  << " and apd_header.subrun=temp_info.subrun"
261  << " and apd_header.time_first> '"<<T4w->AsSQLString()<<"' "
262  << " and apd_header.time_last<'"<<TNow->AsSQLString()<<" ' "
263  << " order by apd_header.run";
264 
265  std::string dbstring2 = dbstream2.str();
266 
267  char * dbcstr2;
268  dbcstr2 = new char [dbstring2.size()+1];
269  strcpy(dbcstr2, dbstring2.c_str());
270 
271  //password removed
272  TSQLServer *db = TSQLServer::Connect("pgsql://131.215.116.171", "", "");
273 
274  std::cout << "Connected? " << db->IsConnected() << std::endl;
275 
276  TSQLStatement *stmt = db->Statement(dbcstr2,1000);
277 
278  std::cout << " about to send statement: " << dbcstr2 << std::endl;
279 
280  if(stmt->Process()){
281  stmt->StoreResult();
282 
283  while (stmt->NextResultRow())
284  {
285  tpe = stmt->GetDouble(0);
286  tpec = stmt->GetDouble(1);
287  trl = stmt->GetDouble(2);
288  trh = stmt->GetDouble(3);
289  tnh = stmt->GetInt(4);
290  tnn = stmt->GetDouble(5);
291  const char* time_last = stmt->GetString(6);
292  ttave = stmt->GetDouble(7);
293  ttn = stmt->GetInt(8);
294  tpesq = stmt->GetDouble(9);
295  tpecorrsq = stmt->GetDouble(10);
296  telapse = stmt->GetInt(11);
297 
298  TDatime *ltime = new TDatime;
299  ltime->Set(time_last);
300  unsigned int tlint = ltime->Convert();
301 
302  int tbin = (tlint-X4WeeksAgo) / (3600*24);
303  if(tbin == 28) tbin = 27;
304  tbin = tbin + 1;
305 
306  // cout << "4w data: " << time_last << " " << tpe << " " << tbin << endl;
307 
308  pe_4w[tbin] += tpe*tnh;
309  pecorr_4w[tbin] += tpec*tnh;
310  nn_4w[tbin] += tnn*telapse;
311  rate_hi_4w[tbin] += trh*telapse;
312  rate_low_4w[tbin] += trl*telapse;
313  nhit_4w[tbin] += tnh;
314  count_4w[tbin] += ttn;
315  tave_4w[tbin] += ttn*ttave;
316  pesq_4w[tbin] += tpesq*tnh;
317  pecorrsq_4w[tbin] += tpecorrsq*tnh;
318  telapse_4w[tbin] += telapse;
319  }
320 
321  double error1[30] = {0.};
322  double error2[30] = {0.};
323  double error3[30] = {0.};
324  double error4[30] = {0.};
325  double error5[30] = {0.};
326  double error6[30] = {0.};
327 
328  for(unsigned int n = 1; n < 29; ++n){
329 
330  if(nhit_4w[n]!=0) {
331 
332  error4[n] = sqrt(rate_hi_4w[n]) / telapse_4w[n];
333  error5[n] = sqrt(rate_low_4w[n]) / telapse_4w[n];
334 
335  pe_4w[n] = pe_4w[n]/nhit_4w[n];
336  pecorr_4w[n] = pecorr_4w[n]/nhit_4w[n];
337  pesq_4w[n] = pesq_4w[n]/nhit_4w[n];
338  pecorrsq_4w[n] = pecorrsq_4w[n]/nhit_4w[n];
339  nn_4w[n] = nn_4w[n]/nhit_4w[n];
340  rate_hi_4w[n] = rate_hi_4w[n]/telapse_4w[n];
341  rate_low_4w[n] = rate_low_4w[n]/telapse_4w[n];
342 
343  error3[n] = ( 1 / sqrt(nhit_4w[n]) ) * nn_4w[n];
344 
345  if(pesq_4w[n]==0) {
346  error1[n] = ( 1 / sqrt(nhit_4w[n]) ) * pe_4w[n];
347  error2[n] = ( 1 / sqrt(nhit_4w[n]) ) * pecorr_4w[n];
348  }
349  else if(pesq_4w[n]>0) {
350  error1[n] = sqrt( fabs(pesq_4w[n] - (pe_4w[n]*pe_4w[n])) ) * 1.5 / sqrt(nhit_4w[n]);
351  error2[n] = sqrt( fabs(pecorrsq_4w[n] - (pecorr_4w[n]*pecorr_4w[n])) * 1.5 / sqrt(nhit_4w[n]) );
352  }
353  else {
354  std::cout << "Negative squares??? Something is wrong!" << std::endl;
355  }
356 
357  std::cout << "Values: " << pe_4w[n] << " " << pesq_4w[n] << " " << nhit_4w[n] << " " << error1[n] << std::endl;
358 
359  }
360 
361  if(count_4w[n]!=0) {
362  tave_4w[n] = tave_4w[n]/count_4w[n];
363  error6[n] = ( sqrt(count_4w[n]) / count_4w[n] ) * tave_4w[n];
364  }
365  }
366 
367  pe_time_4w->SetContent(pe_4w);
368  pec_time_4w->SetContent(pecorr_4w);
369  nn_time_4w->SetContent(nn_4w);
370  rh_time_4w->SetContent(rate_hi_4w);
371  rl_time_4w->SetContent(rate_low_4w);
372  apd_temp_4w->SetContent(tave_4w);
373 
374  pe_time_4w->SetError(error1);
375  pec_time_4w->SetError(error2);
376  nn_time_4w->SetError(error3);
377  rh_time_4w->SetError(error4);
378  rl_time_4w->SetError(error5);
379  apd_temp_4w->SetError(error6);
380 
381  delete stmt;
382 
383  }
384 
385  delete db;
386 
387  pe_time_2d->GetXaxis()->SetTimeDisplay(1);
388  pec_time_2d->GetXaxis()->SetTimeDisplay(1);
389  nn_time_2d->GetXaxis()->SetTimeDisplay(1);
390  rh_time_2d->GetXaxis()->SetTimeDisplay(1);
391  rl_time_2d->GetXaxis()->SetTimeDisplay(1);
392  apd_temp_2d->GetXaxis()->SetTimeDisplay(1);
393 
394  pe_time_2d->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
395  pec_time_2d->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
396  nn_time_2d->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
397  rh_time_2d->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
398  rl_time_2d->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
399  apd_temp_2d->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
400 
401  pe_time_4w->GetXaxis()->SetTimeDisplay(1);
402  pec_time_4w->GetXaxis()->SetTimeDisplay(1);
403  nn_time_4w->GetXaxis()->SetTimeDisplay(1);
404  rh_time_4w->GetXaxis()->SetTimeDisplay(1);
405  rl_time_4w->GetXaxis()->SetTimeDisplay(1);
406  apd_temp_4w->GetXaxis()->SetTimeDisplay(1);
407 
408  pe_time_4w->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
409  pec_time_4w->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
410  nn_time_4w->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
411  rh_time_4w->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
412  rl_time_4w->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
413  apd_temp_4w->GetXaxis()->SetTimeFormat("#splitline{%m\/%d}{%H\:%M}");
414 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TDatime * TNow
Definition: AnaPlotMaker.h:41
T sqrt(T number)
Definition: d0nt_math.hpp:156
OStream cout
Definition: OStream.cxx:6
void ChanPlots(int chan)
Definition: ChanPlots.C:11
enum BeamMode string