source: trunk/Mars/fact/processing/numevents.C @ 19267

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