source: branches/Corsika7500Compatibility/fact/processing/numevents.C@ 18454

Last change on this file since 18454 was 18075, checked in by Daniela Dorner, 10 years ago
added the possibility to call macro also for 1 ganymed file
File size: 9.2 KB
Line 
1#include <iostream>
2
3#include <TCanvas.h>
4#include <TString.h>
5#include <TTree.h>
6#include <TFile.h>
7#include <TSystem.h>
8#include <TSQLResult.h>
9#include <TSQLRow.h>
10#include <TPRegexp.h>
11
12#include "MSQLMagic.h"
13#include "MAlphaFitter.h"
14#include "MHThetaSq.h"
15#include "MStatusArray.h"
16
17
18using namespace std;
19
20int process(MSQLMagic &serv, TString gfile, TString sfile, TString night2, TString runid, TString table, Bool_t pernight, Int_t source)
21{
22 TString where;
23 TString vars;
24 //the variables to save into the db
25 Int_t fNumEvtsAfterCleaning = 0;
26 Int_t fNumEvtsAfterQualCuts = 0;
27 Int_t fNumEvtsAfterBgCuts = 0;
28 Float_t fNumBgEvts = 0;
29 Float_t fNumSigEvts = 0;
30 Float_t fNumExcEvts = 0;
31 Float_t fNumIslandsMean = 0;
32 Float_t fOnTimeAfterCuts = 0;
33
34 if (!pernight)
35 {
36 //check star file status
37 if (gSystem->AccessPathName(sfile))
38 {
39 cout << "ERROR - file " << sfile << " does not exist." << endl;
40 return 2;
41 }
42
43 //try to open the star file
44 TFile star_file(sfile, "READ");
45 if (!star_file.IsOpen())
46 {
47 cout << "ERROR - Could not open file " << sfile << endl;
48 return 2;
49 }
50
51 //Get the events tree
52 TTree *star_tree = (TTree *)star_file.Get("Events");
53 if (!star_tree)
54 {
55 cout << "ERROR - Could not read tree " << sfile << endl;
56 return 2;
57 }
58
59 //the number of events after cleaning
60 fNumEvtsAfterCleaning = star_tree->GetEntries();
61
62 //get the mean number of islands in the _I.root file on runbasis
63 star_tree->Draw("MImagePar.fNumIslands>>IslandHisto","","");
64 TH1F *HistJo = (TH1F*)gDirectory->Get("IslandHisto");
65 if (!HistJo)
66 {
67 cout << "ERROR - Reading of IslandHisto failed " << HistJo << endl;
68 return 2;
69 }
70 fNumIslandsMean = HistJo->GetMean();
71 }
72
73 //check ganymed file status
74 if (gSystem->AccessPathName(gfile))
75 {
76 cout << "ERROR - file " << gfile << " does not exist." << endl;
77 return 2;
78 }
79
80 //try to open the ganymed file
81 TFile ganymed_file(gfile, "READ");
82 if (!ganymed_file.IsOpen())
83 {
84 cout << "ERROR - Could not open file " << gfile << endl;
85 return 2;
86 }
87
88 //Get the status display contents
89 MStatusArray arr;
90 if (arr.Read()<=0)
91 {
92 cout << "ERROR - could not read MStatusArray, file " << gfile << endl;
93 return 2;
94 }
95
96 TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("OnTime", "TH1D", "OnTime");
97 if (!vstime)
98 {
99 cout << "WARNING - Reading of Theta failed." << endl;
100 return 2;
101 }
102 fOnTimeAfterCuts = vstime->Integral();
103
104 //Get number of events after quality cuts
105 //from the number of entries in a histogram
106 TH1F *precut =
107 (TH1F*)arr.FindCanvas("PreCut")->FindObject("SizeSame");
108 if (!precut)
109 {
110 cout << "WARNING - Reading of PreCut canvas failed." << endl;
111 return 2;
112 }
113 fNumEvtsAfterQualCuts = precut->GetEntries();
114
115 //Get number of events after background rejection cuts
116 //from the number of entries in a histogram
117 TH1F *postcut =
118 (TH1F*)arr.FindCanvas("PostCut")->FindObject("SizeSame");
119 if (!postcut)
120 {
121 cout << "WARNING - Reading of PostCut canvas failed." << endl;
122 return 2;
123 }
124 fNumEvtsAfterBgCuts = postcut->GetEntries();
125
126 //Get signal, bg and excess events
127 MHThetaSq *halphaon = (MHThetaSq*)arr.FindObjectInCanvas("MHThetaSq", "Hist");
128 if (!halphaon)
129 {
130 cout << "WARNING - Reading of MHThetaSq failed, file " << gfile << endl;
131 return 2;
132 }
133 const MAlphaFitter &fit = halphaon->GetAlphaFitter();
134 if (!&fit)
135 {
136 cout << "WARNING - Reading of MAlphaFitter failed, file " << gfile << endl;
137 return 2;
138 }
139
140 //check if the excess events number is small and close to 0
141 //due to the computer precision error when subtracting sig and bg events
142 fNumBgEvts = fit.GetEventsBackground();
143 fNumSigEvts = fit.GetEventsSignal();
144 fNumExcEvts = fit.GetEventsExcess();
145 fNumExcEvts = fabs(fNumExcEvts)<1e-5 ? 0 : fNumExcEvts;
146
147 /*
148 cout<< fNumEvtsAfterCleaning<< " "
149 << fNumEvtsAfterQualCuts<< " "
150 << fNumEvtsAfterBgCuts << " "
151 << fNumBgEvts << " "
152 << fNumSigEvts << " "
153 << fNumExcEvts << endl;
154 */
155
156 //inserting or updating the information in the database
157 if (!pernight)
158 vars = Form("fRunID=%s, fNight=%s,", runid.Data(), night2.Data());
159 else
160 vars = Form("fNight=%s, fSourceKey=%d, ", night2.Data(), source);
161 vars += Form(" fNumEvtsAfterQualCuts=%d, "
162 " fNumEvtsAfterBgCuts=%d, "
163 " fNumBgEvts=%6.1f, "
164 " fNumSigEvts=%6.1f, "
165 " fNumExcEvts=%6.1f, "
166 " fOnTimeAfterCuts=%7.2f ",
167 fNumEvtsAfterQualCuts,
168 fNumEvtsAfterBgCuts,
169 fNumBgEvts,
170 fNumSigEvts,
171 fNumExcEvts,
172 fOnTimeAfterCuts
173 );
174 if (!pernight)
175 vars += Form(", fNumIslandsMean=%6.2f, fNumEvtsAfterCleaning=%d ",
176 fNumIslandsMean,
177 fNumEvtsAfterCleaning
178 );
179
180
181 if (pernight)
182 where = Form("fNight=%s AND fSourceKey=%d", night2.Data(), source);
183 else
184 where = Form("fRunID=%s AND fNight=%s", runid.Data(), night2.Data());
185
186 cout << "vars: " << vars << endl;
187 cout << "where: " << where << endl;
188
189 if (!serv.InsertUpdate(table, vars, where))
190 {
191 cout << "ERROR - insert/update to DB failed (vars:" << vars << " where: " << where << ")." << endl;
192 return 2;
193 }
194 return 1;
195}
196
197
198int numevents(Int_t night, TString inpath="/loc_data/analysis/", TString table= "AnalysisResultsRunLP", Bool_t dummy=kTRUE, Bool_t pernight=kFALSE, Int_t source=0)
199{
200
201 //connect to mysql server
202 //MSQLServer serv("sql.rc");
203 MSQLMagic serv("sql.rc");
204 if (!serv.IsConnected())
205 {
206 cout << "ERROR - Connection to database failed." << endl;
207 return 0;
208 }
209 serv.SetIsDummy(dummy);
210
211 //select runs, where star finished, format the list nicely
212 TString query;
213 query = " SELECT CONCAT(fNight, '_', LPAD(fRunID, 3, 0)), SUBSTRING(fNight, 1,4), SUBSTRING(fNight, 5,2), SUBSTRING(fNight, 7,2), ";
214 query += " fNight, fRunID FROM RunInfo";
215
216 //and only data runs
217 query+=" WHERE fRunTypeKEY=1 ";
218 query+=Form(" AND fNight=%d",night);
219 //query+=" AND fRunID=128 ";
220 if (pernight)
221 query+=" GROUP BY fNight";
222
223 cout << "Q: " << query << endl;
224
225 //post the query
226 TSQLResult *res = serv.Query(query);
227 if (!res)
228 {
229 cout << "ERROR - no result from query " << query << endl;
230 return 2;
231 }
232
233 //allocate variables to use in the loop
234 TString night2;
235 TString runid;
236 TString run;
237 TString year;
238 TString month;
239 TString day;
240 TString ganymed_fname;
241 TString star_fname;
242 Int_t rc=0;
243
244 //loop over the data files
245 TSQLRow *row=0;
246
247 //cout the data names
248 cout<<"#evts_after_cleaning" << " "
249 <<"evts_after_qual-cuts" << " "
250 <<"evts_after_bg-cuts" << " "
251 <<"mean_evts_from_off_regions"<< " "
252 <<"evts_from_on_region" << " "
253 <<"calculated_excess_evts_in_file" << " " << endl;
254
255 while ((row=res->Next()))
256 {
257 //use the results from the query
258 night2=(*row)[4];
259 runid=(*row)[5];
260 run=(*row)[0];
261 year=(*row)[1];
262 month=(*row)[2];
263 day=(*row)[3];
264
265 // to be fixed: path for nightly ganymeds
266 if (!pernight)
267 ganymed_fname=inpath+"/ganymed_run/"+year+"/"+month+"/"+day+"/"+run+"-analysis.root";
268 else
269 if (!source)
270 cout << "ERROR - source " << source << " not valid. " << endl;
271 else
272 ganymed_fname=Form("%s/ganymed_night/%d/%s-analysis.root", inpath.Data(), source, night2.Data());
273 cout << "gf => " << ganymed_fname << endl;
274
275 if (!pernight)
276 {
277 star_fname=inpath+"/star/"+year+"/"+month+"/"+day+"/"+run+"_I.root";
278 cout << "sf => " << star_fname << endl;
279 }
280 rc=process(serv, ganymed_fname, star_fname, night2, runid, table, pernight, source);
281 if (rc>1)
282 return rc;
283
284 }
285
286
287 return 1;
288}
289
290int numevents(TString ganymedfile, TString starfile, TString table, Bool_t dummy=kTRUE)
291{
292
293 //connect to mysql server
294 MSQLMagic serv("sql.rc");
295 if (!serv.IsConnected())
296 {
297 cout << "ERROR - Connection to database failed." << endl;
298 return 0;
299 }
300 serv.SetIsDummy(dummy);
301 //get night and run from file name
302 TString basename=gSystem->BaseName(ganymedfile);
303 TPRegexp regexp("20[0-9][0-9][0-2][0-9][0-3][0-9]_[0-9][0-9][0-9]");
304 TString run = basename(regexp);
305 TPRegexp regexp2("20[0-9][0-9][0-2][0-9][0-3][0-9]");
306 TPRegexp regexp3("_[0-9][0-9][0-9]");
307 TPRegexp regexp4("[0-9][0-9][0-9]");
308 TString night = run(regexp2);
309 TString runnum = run(regexp3);
310 TString runid = runnum(regexp4);
311
312 return process(serv, ganymedfile, starfile, night, runid, table, false, 0);
313}
314
Note: See TracBrowser for help on using the repository browser.