source: branches/Mars_MC/fact/processing/numevents.C@ 17054

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