Changeset 18075 for trunk/Mars/fact
- Timestamp:
- 01/04/15 05:11:16 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/fact/processing/numevents.C
r17977 r18075 8 8 #include <TSQLResult.h> 9 9 #include <TSQLRow.h> 10 #include <TPRegexp.h> 10 11 11 12 #include "MSQLMagic.h" … … 17 18 using namespace std; 18 19 19 //missing: 20 // remove sourcekey in paths 21 // create tables in DB 22 23 int numevents(TString night="20130418", TString inpath="/daq/analysis/1/", TString table= "AnalysisResultsRunLP", Bool_t dummy=kTRUE, Bool_t pernight=kFALSE, Int_t source=0) 20 int 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 198 int numevents(Int_t night, TString inpath="/loc_data/analysis/", TString table= "AnalysisResultsRunLP", Bool_t dummy=kTRUE, Bool_t pernight=kFALSE, Int_t source=0) 24 199 { 25 200 … … 41 216 //and only data runs 42 217 query+=" WHERE fRunTypeKEY=1 "; 43 query+= " AND fNight="+night;218 query+=Form(" AND fNight=%d",night); 44 219 //query+=" AND fRunID=128 "; 45 220 if (pernight) … … 65 240 TString ganymed_fname; 66 241 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; 242 Int_t rc=0; 79 243 80 244 //loop over the data files 81 245 TSQLRow *row=0; 82 246 83 /*84 247 //cout the data names 85 248 cout<<"#evts_after_cleaning" << " " … … 89 252 <<"evts_from_on_region" << " " 90 253 <<"calculated_excess_evts_in_file" << " " << endl; 91 */92 254 93 255 while ((row=res->Next())) … … 115 277 star_fname=inpath+"/star/"+year+"/"+month+"/"+day+"/"+run+"_I.root"; 116 278 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 } 279 } 280 rc=process(serv, ganymed_fname, star_fname, night2, runid, table, pernight, source); 281 if (rc>1) 282 return rc; 283 284 } 285 278 286 279 287 return 1; 280 288 } 289 290 int 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 TracChangeset
for help on using the changeset viewer.