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

Last change on this file since 19267 was 19267, checked in by tbretz, 6 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.