source: branches/Corsika7500Compatibility/fact/processing/DataFromQLA.C@ 19004

Last change on this file since 19004 was 18177, checked in by Daniela Dorner, 10 years ago
macro to retrieve official data from QLA
File size: 8.6 KB
Line 
1#include <iostream>
2#include <iomanip>
3#include <fstream>
4
5#include <TFile.h>
6#include <TSQLResult.h>
7#include <TSQLRow.h>
8#include <TPRegexp.h>
9
10#include "MSQLMagic.h"
11#include "MTime.h"
12#include "MMath.h"
13#include "MDirIter.h"
14
15using namespace std;
16
17void PrintRemarks(ofstream &fout, TString sourcename, TString firstnight, TString lastnight)
18{
19 MTime now(-1);
20
21 fout << "# " << endl;
22 fout << "# Please cite the FACT design paper and the QLA webpage when using these data." << endl;
23 fout << "# FACT design paper: H. Anderhub et al. JINST 8 (2013) P6008 " << endl;
24 fout << "# http://iopscience.iop.org/1748-0221/8/06/P06008 " << endl;
25 fout << "# QLA webpage: http://www.fact-project.org/monitoring" << endl;
26 fout << "# If you intent to use the data, please let us know for reference. " << endl;
27 fout << "# " << endl;
28 fout << "# Remarks:" << endl;
29 fout << "# * These are the results of a fast quick look analyis " << endl;
30 fout << "# on site, i.e. they are preliminary. " << endl;
31 fout << "# * The quick look analysis includes all data, " << endl;
32 fout << "# i.e. no data selection done." << endl;
33 fout << "# * The given values are not fluxes but excess rates " << endl;
34 fout << "# (number of excess events per effective ontime), " << endl;
35 fout << "# i.e. there is a dependence on trigger threshold and " << endl;
36 fout << "# zenith distance of the observation (with the current " << endl;
37 fout << "# analysis for zenith distance > 40 degree and trigger " << endl;
38 fout << "# threshold > 500 DAC counts)." << endl;
39// fout << "# * The data are provided with 20 min binning and nightly binning." << endl;
40 fout << "# * Nights with less than 20 minutes of data are neglected. " << endl;
41 fout << "# * The QLA results are not reprocessed when a new software " << endl;
42 fout << "# version is introduced. " << endl;
43 fout << "# * In case, you need further details about the data or a" << endl;
44 fout << "# different binning, please do not hesitate to contact us." << endl;
45 fout << "# * The QLA contains all data since 12.12.2012. " << endl;
46 fout << "# For older data, please contact us. " << endl;
47 fout << "# " << endl;
48 fout << "# Contact: Daniela Dorner dorner@astro.uni-wuerzburg.de " << endl;
49 fout << "# " << endl;
50 fout << "# This file was created at " << now.GetString() << endl;
51 fout << "# Source: " << sourcename << endl;
52 fout << "# Timerange: " << firstnight << " - " << lastnight << endl;
53 fout << "# " << endl;
54 fout << "# start(mjd) stop(mjd) excess-rate(evts/h) error-excess-rate " << endl;
55
56 return;
57}
58
59int DataFromQLA(Int_t sourcekey=1, Int_t nightmin=2011115, Int_t nightmax=20161231)
60{
61 MSQLServer serv("sql.rc");
62 if (!serv.IsConnected())
63 {
64 cout << "ERROR - Connection to database failed." << endl;
65 return 0;
66 }
67 Bool_t dch=kFALSE;
68
69 TString query=Form("SELECT fSourceName FROM Source WHERE fSourceKey=%d", sourcekey);
70 TSQLResult *res1 = serv.Query(query);
71 if (!res1)
72 return 1;
73 TSQLRow *row1=res1->Next();
74 TString sourcename=(*row1)[0];
75 delete res1;
76 sourcename.ReplaceAll(" ", "_");
77
78 // datacheck
79 TString datacheck=" ";
80 //remove data with wrong settings
81 datacheck+=" AND fNight>20120420 AND NOT fNight IN (20120406,20120410,20120503) AND";//data with different bias voltage
82 datacheck+=" NOT fNight BETWEEN 20121206 AND 20130110"; // broken bias channel
83 //datacheck+=" AND NOT (fNight=20120608 AND fRunID=65) "; // something went wrong with tracking?
84 // 24.6. new coefficients
85 TString zdparam=" pow(0.753833 * cos(Radians(fZenithDistanceMean)), 7.647435) * exp(-5.753686*pow(Radians(fZenithDistanceMean),2.089609))";
86 TString thparam=" pow((if(isnull(fThresholdMinSet),fThresholdMedian,fThresholdMinSet)-329.4203),2) * (-0.0000002044803) ";
87 TString param=" (fNumEvtsAfterBgCuts/5-fNumSigEvts)/fOnTimeAfterCuts - "+zdparam+" - "+thparam+" ";
88 datacheck+=" AND -0.085 < ("+param+") ";
89 datacheck+=" AND 0.25 > ("+param+") ";
90
91 TString select =" SELECT Min(fNight), Max(fNight), Sum(fOnTimeAfterCuts) as bla ";
92 TString fromjoinwhere=" FROM AnalysisResultsRunLP ";
93 fromjoinwhere+=" LEFT JOIN RunInfo USING(fNight, fRunID) ";
94 fromjoinwhere+=Form(" WHERE fSourceKey=%d", sourcekey);
95 fromjoinwhere+=Form(" AND fNight BETWEEN %d AND %d ", nightmin, nightmax);
96 fromjoinwhere+=" AND fOnTimeAfterCuts < 1000 "; //exclude runs with wrong/too high ontime
97 fromjoinwhere+=" AND fOnTimeAfterCuts > 10 "; //exclude runs with wrong/too small ontime
98 fromjoinwhere+=" AND NOT ISNULL(fNumExcEvts) ";// only where excess was extracted
99 if (dch)
100 fromjoinwhere+=datacheck;
101 query=select+fromjoinwhere;
102 query+=" GROUP BY fSourceKey ";
103 query+=" HAVING bla>20*60 ";
104 cout << "Q: " << query << endl;
105 TSQLResult *res2 = serv.Query(query);
106 if (!res2)
107 return 1;
108 TSQLRow *row2=res2->Next();
109 TString firstnight=(*row2)[0];
110 TString lastnight=(*row2)[1];
111 delete res2;
112
113 TString filename_nightly="FACT_QuickLookAnalysisResults_NightlyBinning_"+sourcename+".txt";
114 TString filename_min="FACT_QuickLookAnalysisResults_20MinuteBinning_"+sourcename+".txt";
115
116 ofstream nightly(filename_nightly);
117 ofstream min(filename_min);
118 if (!nightly)
119 {
120 cout << "ERROR - cannot write " << filename_nightly << endl;
121 return 2;
122 }
123 if (!min)
124 {
125 cout << "ERROR - cannot write " << filename_min << endl;
126 return 2;
127 }
128 PrintRemarks(nightly, sourcename, firstnight, lastnight);
129 PrintRemarks(min, sourcename, firstnight, lastnight);
130
131 // query data from AnalysisResultsRunLP
132 select =" SELECT Sum(fOnTimeAfterCuts) as bla, Sum(fNumExcEvts), Sum(fNumBgEvts), Sum(fNumSigEvts), ";
133 select+=" Min(fRunStart), Max(fRunStop), fNight, fRunID ";
134 query=select+fromjoinwhere;
135 query+=" GROUP BY fNight ";
136 query+=" HAVING bla>20*60 ";
137 query+=" ORDER BY fNight ";
138 //cout << "Q: " << query << endl;
139
140 //variables for calculations and graphs
141 TString excevts, bgevts, sigevts, ontime, night;
142 Float_t excevtssum=0;
143 Float_t bgevtssum=0;
144 Float_t sigevtssum=0;
145 Float_t ontimesum=0;
146 Float_t excerr=0;
147 MTime start;
148 MTime stop;
149 Double_t mjdstart=0;
150 Double_t mjdstop=0;
151 //Double_t mjd=0;
152 //Double_t mjddiff=0;
153 Int_t counter=0;
154
155 TSQLResult *res3 = serv.Query(query);
156 if (!res3)
157 return 1;
158 TSQLRow *row3=0;
159 while ((row3=res3->Next()))
160 {
161 ontime=(*row3)[0];
162 excevts=(*row3)[1];
163 bgevts=(*row3)[2];
164 sigevts=(*row3)[3];
165 start.SetSqlDateTime((*row3)[4]);
166 stop.SetSqlDateTime((*row3)[5]);
167 night=(*row3)[6];
168
169 mjdstart = start.GetMjd();
170 mjdstop = stop.GetMjd();
171 //mjd = mjdstart+(mjdstop-mjdstart)/2.;
172
173 //significance = MMath::SignificanceLiMaSigned(sigevts.Atof(), bgevts.Atof()*5, 0.2);
174 excerr = MMath::ErrorExc(sigevts.Atof(), bgevts.Atof()*5, 0.2)/ontime.Atof()*3600.;
175 nightly.precision(18);
176 //fout << (*row3)[6] << ": " << start.GetMjd() << " " << stop.GetMjd() << " " << excevts.Atof()/ontime.Atof() << " " << excerr << endl;
177 nightly << start.GetMjd() << " " << stop.GetMjd() << " " << excevts.Atof()/ontime.Atof()*3600 << " " << excerr << endl;
178 query="SELECT fOnTimeAfterCuts, fNumExcEvts, fNumBgEvts, fNumSigEvts, ";
179 query+=" fRunStart, fRunStop "+fromjoinwhere;
180 query+=" AND fNight="+night;
181 //cout << "Q: " << query << endl;
182 TSQLResult *res4 = serv.Query(query);
183 if (!res4)
184 return 1;
185 TSQLRow *row4=0;
186 counter=0;
187 while ((row4=res4->Next()))
188 {
189 ontime=(*row4)[0];
190 excevts=(*row4)[1];
191 bgevts=(*row4)[2];
192 sigevts=(*row4)[3];
193 if(counter==0)
194 start.SetSqlDateTime((*row4)[4]);
195
196 if (ontimesum+ontime.Atoi()>20*60)
197 {
198 excerr = MMath::ErrorExc(sigevtssum, bgevtssum*5, 0.2)/ontimesum*3600.;
199 min.precision(18);
200 min << start.GetMjd() << " " << stop.GetMjd() << " " << excevtssum/ontimesum*3600 << " " << excerr << endl;
201 start.SetSqlDateTime((*row4)[4]);
202 ontimesum=0;
203 excevtssum=0;
204 bgevtssum=0;
205 sigevtssum=0;
206 }
207 counter++;
208 stop.SetSqlDateTime((*row4)[5]);
209 ontimesum+=ontime.Atoi();
210 excevtssum+=excevts.Atof();
211 bgevtssum+=bgevts.Atof();
212 sigevtssum+=sigevts.Atof();
213 }
214 }
215
216 delete res3;
217 return 0;
218
219}
220
221
222
Note: See TracBrowser for help on using the repository browser.