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

Last change on this file since 18688 was 18688, checked in by dorner, 3 years ago
updated to not stop if one file is not found
File size: 9.4 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    Int_t numpb=0;
244
245    //loop over the data files
246    TSQLRow *row=0;
247
248    //cout the data names
249    cout<<"#evts_after_cleaning" << "    "
250        <<"evts_after_qual-cuts" << "    "
251        <<"evts_after_bg-cuts" << "    "
252        <<"mean_evts_from_off_regions"<< "    "
253        <<"evts_from_on_region"  << "    "
254        <<"calculated_excess_evts_in_file" << "    " << endl;
255
256    while ((row=res->Next()))
257    {
258        //use the results from the query
259        night2=(*row)[4];
260        runid=(*row)[5];
261        run=(*row)[0];
262        year=(*row)[1];
263        month=(*row)[2];
264        day=(*row)[3];
265
266        // to be fixed: path for nightly ganymeds
267        if (!pernight)
268            ganymed_fname=inpath+"/ganymed_run/"+year+"/"+month+"/"+day+"/"+run+"-analysis.root";
269        else
270            if (!source)
271                cout << "ERROR - source " << source << " not valid. " << endl;
272            else
273                ganymed_fname=Form("%s/ganymed_night/%d/%s-analysis.root", inpath.Data(), source, night2.Data());
274        cout << "gf => " << ganymed_fname << endl;
275
276        if (!pernight)
277        {
278            star_fname=inpath+"/star/"+year+"/"+month+"/"+day+"/"+run+"_I.root";
279            cout << "sf => " << star_fname << endl;
280        }
281        rc=process(serv, ganymed_fname, star_fname, night2, runid, table, pernight, source);
282        if (rc>1)
283            numpb++;
284            //return rc;
285
286    }
287
288
289    if (numpb>0)
290        return numpb+1;//avoid that 1 (i.e. everything ok) is returned in case of problems
291    else
292        return 1;
293}
294
295int numevents(TString ganymedfile, TString starfile, TString table, Bool_t dummy=kTRUE)
296{
297
298    //connect to mysql server
299    MSQLMagic serv("sql.rc");
300    if (!serv.IsConnected())
301    {
302        cout << "ERROR - Connection to database failed." << endl;
303        return 0;
304    }
305    serv.SetIsDummy(dummy);
306    //get night and run from file name
307    TString basename=gSystem->BaseName(ganymedfile);
308    TPRegexp regexp("20[0-9][0-9][0-2][0-9][0-3][0-9]_[0-9][0-9][0-9]");
309    TString run = basename(regexp);
310    TPRegexp regexp2("20[0-9][0-9][0-2][0-9][0-3][0-9]");
311    TPRegexp regexp3("_[0-9][0-9][0-9]");
312    TPRegexp regexp4("[0-9][0-9][0-9]");
313    TString night = run(regexp2);
314    TString runnum = run(regexp3);
315    TString runid = runnum(regexp4);
316
317    return process(serv, ganymedfile, starfile, night, runid, table, false, 0);
318}
319
Note: See TracBrowser for help on using the repository browser.