Changeset 17521
- Timestamp:
- 01/20/14 17:42:35 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/fact/plots/quality.C
r15244 r17521 47 47 48 48 TGraph g; 49 g.SetName("Threshold"); 49 50 50 51 while (file.GetNextRow()) … … 54 55 } 55 56 57 g.SetMinimum(281); 56 58 g.SetMarkerStyle(kFullDotMedium); 57 59 g.GetXaxis()->SetTimeDisplay(true); … … 70 72 #include <pair> 71 73 72 vector<pair<double, ln_equ_posn>> vec; 73 74 ln_equ_posn FindPointing(Double_t time) 75 { 76 for (int i=0; i<vec.size(); i++) 77 if (time<vec[i].first) 78 return i==0 ? ln_equ_posn() : vec[i-1].second; 79 80 return vec[vec.size()-1].second; 81 } 82 74 vector<pair<double, Nova::EquPosn>> vecp; 75 76 Nova::EquPosn FindPointing(Double_t time) 77 { 78 for (int i=0; i<vecp.size(); i++) 79 if (time<vecp[i].first) 80 { 81 if (i==0) 82 return Nova::EquPosn(); 83 else 84 return vecp[i-1].second; 85 } 86 87 return vecp[vecp.size()-1].second; 88 } 89 /* 83 90 Float_t Prediction(Double_t time) 84 91 { … … 93 100 94 101 // Distance between source and moon 95 //const double angle = 96 // MAstro::AngularDistance(TMath::DegToRad()*(90-moon.dec), 97 // TMath::DegToRad()*moon.ra, 98 // TMath::DegToRad()*(90-pos.dec), 99 // TMath::DegToRad()*pos.ra); 102 //const double angle = Nova::GetAngularSeparation(moon, hrzm)*TMath::DegToRad(); 100 103 101 104 // Current prediction … … 105 108 double disk = Nova::GetLunarDisk(jd); 106 109 107 double lc = /*cang==0 ? 0 :*/ calt*pow(disk, 2.5);///sqrt(cang); 108 109 return 5+103*lc>4.5 ? 5+103*lc : 4.5; 110 // semi-major axis 111 double lc = calt*pow(disk, 2.2)*pow(Nova::GetLunarEarthDist(jd)/384400, -2);///sqrt(cang); 112 113 return lc>0 ? 4+103*lc : 4; 114 }*/ 115 116 Float_t Prediction(Double_t time) 117 { 118 Double_t jd = time + 40587 + 2400000.5; 119 /* 120 Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01); 121 122 Nova::EquPosn pos = FindPointing(time); 123 124 // get local position of moon 125 Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd); 126 127 // Distance between source and moon 128 const double angle = Nova::GetAngularSeparation(moon, pos); 129 130 // Distance between earth and moon relative to major semi-axis 131 const double dist = Nova::GetLunarEarthDist(jd)/384400; 132 133 // Current prediction 134 double cang = 1-sin(angle*TMath::DegToRad()); 135 double calt = sin(hrzm.alt*TMath::DegToRad()); 136 137 double disk = Nova::GetLunarDisk(jd); 138 139 // light condition 140 double lc = sqrt(calt)*pow(disk, 2.3)*pow(cang+1, 1)*pow(dist, -2); 141 142 return lc>0 ? 7.2 + 69*lc : 7.2; 143 */ 144 145 // Sun properties 146 Nova::EquPosn sun = Nova::GetSolarEquCoords(jd); 147 Nova::ZdAzPosn hrzs = Nova::GetHrzFromEqu(sun, jd); 148 149 // Get source position 150 Nova::EquPosn pos = FindPointing(time); 151 152 // Moon properties 153 Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01); 154 Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd); 155 double disk = Nova::GetLunarDisk(jd); 156 157 // Derived moon properties 158 double angle = Nova::GetAngularSeparation(moon, pos); 159 double edist = Nova::GetLunarEarthDist(jd)/384400; 160 161 // Current prediction 162 double sin_malt = hrzm.alt<0 ? 0 : sin(hrzm.alt*TMath::DegToRad()); 163 double sin_mdist = sin(angle*TMath::DegToRad()); 164 double cos_salt = cos(hrzs.zd*TMath::DegToRad()); 165 166 double c0 = pow(disk, 2.52); 167 double c1 = pow(sin_malt, 0.72); 168 double c2 = pow(edist, -2.00); 169 double c3 = exp(1.46*(1-sin_mdist)*(1-sin_mdist)); 170 double c4 = exp(33.0)*exp(-20.1*(1-cos_salt)*(1-cos_salt)); 171 172 double cur = 6.4 + 96.9*c0*c1*c2*c3 + c4; 173 174 // cout << cur << " " << hrzm.alt << " " << c1 << endl; 175 176 return cur; 177 110 178 } 111 179 … … 131 199 while (file.GetNextRow()) 132 200 { 133 ln_equ_posn p;201 Nova::EquPosn p; 134 202 p.ra = ra*15; 135 203 p.dec = dec; 136 204 137 vec .push_back(make_pair(time, p));205 vecp.push_back(make_pair(time, p)); 138 206 } 139 207 … … 169 237 TGraph g1; 170 238 TGraph g2; 239 g1.SetName("Currents"); 240 g2.SetName("Prediction"); 171 241 172 242 while (file.GetNextRow()) … … 234 304 235 305 TGraph g1, g2; 306 g1.SetName("TriggerRate"); 307 g2.SetName("RelOnTime"); 236 308 237 309 while (file.GetNextRow()) … … 246 318 247 319 g1.SetMinimum(0); 248 g1.SetMaximum(2 00);320 g1.SetMaximum(269); 249 321 g1.SetMarkerStyle(kFullDotMedium); 250 322 g1.GetXaxis()->SetTimeDisplay(true); … … 260 332 gPad->GetCanvas()->cd(4); 261 333 262 gPad->SetGrid x();334 gPad->SetGrid(); 263 335 gPad->SetTopMargin(0); 264 336 gPad->SetBottomMargin(0); … … 278 350 g2.DrawClone("AP"); 279 351 280 281 return 0; 282 } 352 return 0; 353 } 354 355 void PlotRateQC(UInt_t night, MSQLServer &serv) 356 { 357 TString query = 358 "LEFT JOIN AnalysisResultsRunLP USING(fNight, fRunID) " 359 "WHERE fRunTypeKey=1 AND NOT ISNULL (AnalysisResultsRunLP.fNumEvtsAfterCleaning) AND fNight="; 360 query += night; 361 362 TTree *t = serv.GetTree("RunInfo", query); 363 if (!t) 364 return; 365 366 int save = gErrorIgnoreLevel; 367 gErrorIgnoreLevel = kFatal; 368 369 gROOT->SetSelectedPad(0); 370 gPad->GetCanvas()->cd(3); 371 372 t->Draw("AnalysisResultsRunLP.fNumEvtsAfterCleaning/AnalysisResultsRunLP.fOnTimeAfterCuts:(RunInfo.fRunStart+RunInfo.fRunStop)/2+9131*24*3600", "", "same"); 373 TGraph *g = (TGraph*)gPad->GetPrimitive("Graph"); 374 if (g) 375 { 376 g->SetName("CleaningRate"); 377 g->SetMarkerColor(kRed); 378 g->SetMarkerStyle(29);//kFullDotMedium); 379 } 380 381 t->Draw("AnalysisResultsRunLP.fNumEvtsAfterQualCuts/AnalysisResultsRunLP.fOnTimeAfterCuts:(RunInfo.fRunStart+RunInfo.fRunStop)/2+9131*24*3600", "", "same"); 382 g = (TGraph*)gPad->GetPrimitive("Graph"); 383 if (g) 384 { 385 g->SetName("RateAfterQC"); 386 g->SetMarkerColor(kBlue); 387 g->SetMarkerStyle(29);//kFullDotMedium); 388 } 389 390 gErrorIgnoreLevel = save; 391 } 392 283 393 284 394 Int_t PlotPointing(TArrayD **vec, TString fname) … … 304 414 305 415 TGraph g; 416 g.SetName("Zd"); 306 417 307 418 while (file.GetNextRow()) … … 309 420 g.SetPoint(g.GetN(), time*24*3600, 90-zd); 310 421 311 g.SetMinimum( 0);422 g.SetMinimum(1); 312 423 g.SetMaximum(90); 313 424 g.SetMarkerStyle(kFullDotMedium); 314 425 g.GetXaxis()->SetTimeDisplay(true); 315 426 g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT"); 316 g.GetXaxis()->SetLabelSize(0.1 );427 g.GetXaxis()->SetLabelSize(0.12); 317 428 g.GetYaxis()->SetLabelSize(0.1); 318 429 g.GetYaxis()->SetTitle("ELEVATION"); … … 323 434 return 0; 324 435 } 325 /* 326 Int_t Plot Mem(TArrayD **vec, TString fname)327 { 328 fname += ". FAD_CONTROL_STATISTICS1.fits";436 437 Int_t PlotTemperature1(TArrayD **vec, TString fname) 438 { 439 fname += ".TEMPERATURE_DATA.fits"; 329 440 330 441 fits file(fname.Data()); … … 338 449 339 450 Double_t time; 340 UInt_t buf[5]; 341 ULong64_t mem[4]; 342 UInt_t rate[2]; 451 Float_t temp; 343 452 344 453 if (!file.SetPtrAddress("Time", &time)) 345 454 return -1; 346 if (!file.SetPtrAddress("bufferInfo", buf) || 347 !file.SetPtrAddress("memInfo", mem) || 348 !file.SetPtrAddress("rateNew", rate)) 455 if (!file.SetPtrAddress("T", &temp)) 349 456 return -1; 350 457 351 458 TGraph g; 459 g.SetName("ContainerTemp"); 352 460 353 461 while (file.GetNextRow()) 354 462 if (Contains(vec, time)) 355 g.SetPoint(g.GetN(), time*24*3600, mem[1]); 356 463 g.SetPoint(g.GetN(), time*24*3600, temp); 464 465 g.SetMinimum(-5); 466 g.SetMaximum(49); 357 467 g.SetMarkerStyle(kFullDotMedium); 468 g.SetMarkerColor(kRed); 358 469 g.GetXaxis()->SetTimeDisplay(true); 359 g.GetXaxis()->SetTimeFormat("%H:%M"); 360 g.GetXaxis()->SetTimeOffset(0, "gmt"); 361 g.GetXaxis()->SetLabelSize(0.12); 470 g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT"); 471 g.GetXaxis()->SetLabelSize(0.1); 362 472 g.GetYaxis()->SetLabelSize(0.1); 363 //g.GetYaxis()->SetTitle("ELEVATION");473 g.GetYaxis()->SetTitle("TEMP"); 364 474 g.GetYaxis()->SetTitleOffset(0.2); 365 475 g.GetYaxis()->SetTitleSize(0.1); … … 368 478 return 0; 369 479 } 370 */ 371 void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0) 480 481 Int_t PlotTemperature2(TArrayD **vec, TString fname) 482 { 483 fname += ".MAGIC_WEATHER_DATA.fits"; 484 485 fits file(fname.Data()); 486 if (!file) 487 { 488 cerr << fname << ": " << gSystem->GetError() << endl; 489 return -2; 490 } 491 492 //cout << fname << endl; 493 494 Double_t time; 495 Float_t temp; 496 497 if (!file.SetPtrAddress("Time", &time)) 498 return -1; 499 if (!file.SetPtrAddress("T", &temp)) 500 return -1; 501 502 TGraph g; 503 g.SetName("OutsideTemp"); 504 505 while (file.GetNextRow()) 506 if (Contains(vec, time)) 507 g.SetPoint(g.GetN(), time*24*3600, temp); 508 509 g.SetMarkerStyle(kFullDotMedium); 510 g.SetMarkerColor(kBlue); 511 g.DrawClone("P"); 512 513 return 0; 514 } 515 516 Int_t PlotTemperature3(TArrayD **vec, TString fname) 517 { 518 fname += ".FSC_CONTROL_TEMPERATURE.fits"; 519 520 fits file(fname.Data()); 521 if (!file) 522 { 523 cerr << fname << ": " << gSystem->GetError() << endl; 524 return -2; 525 } 526 527 //cout << fname << endl; 528 529 Double_t time; 530 Float_t temp[31]; 531 532 if (!file.SetPtrAddress("Time", &time)) 533 return -1; 534 if (!file.SetPtrAddress("T_sens", temp)) 535 return -1; 536 537 TGraph g, g1, g2; 538 g.SetName("SensorTempAvg"); 539 g1.SetName("SensorTempMin"); 540 g2.SetName("SensorTempMax"); 541 542 while (file.GetNextRow()) 543 if (Contains(vec, time)) 544 { 545 float min = 100; 546 float max = -100; 547 double avg = 0; 548 int num = 0; 549 for (int i=0; i<31; i++) 550 if (temp[i]!=0) 551 { 552 avg += temp[i]; 553 num++; 554 555 min = TMath::Min(min, temp[i]); 556 max = TMath::Max(max, temp[i]); 557 } 558 559 g.SetPoint(g.GetN(), time*24*3600, avg/num); 560 g1.SetPoint(g1.GetN(), time*24*3600, min); 561 g2.SetPoint(g2.GetN(), time*24*3600, max); 562 } 563 564 g.SetMarkerStyle(kFullDotMedium); 565 g.DrawClone("P"); 566 567 /* 568 g1.SetLineWidth(1); 569 g1.DrawClone("L"); 570 571 g2.SetLineWidth(1); 572 g2.DrawClone("L"); 573 */ 574 return 0; 575 } 576 577 Int_t PlotTemperature4(TArrayD **vec, TString fname) 578 { 579 fname += ".FAD_CONTROL_TEMPERATURE.fits"; 580 581 fits file(fname.Data()); 582 if (!file) 583 { 584 cerr << fname << ": " << gSystem->GetError() << endl; 585 return -2; 586 } 587 588 //cout << fname << endl; 589 590 Double_t time; 591 Float_t temp[160]; 592 593 if (!file.SetPtrAddress("Time", &time)) 594 return -1; 595 if (!file.SetPtrAddress("Data1", temp) && 596 !file.SetPtrAddress("temp", temp)) 597 return -1; 598 599 Int_t num = file.GetN("temp")==0 ? file.GetN("Data1") : file.GetN("temp"); 600 Int_t beg = num==82 ? 2 : 0; 601 602 TGraphErrors g1; 603 TGraph g2,g3; 604 605 g1.SetName("FadTempAvg"); 606 g2.SetName("FadTempMin"); 607 g3.SetName("FadTempMax"); 608 609 while (file.GetNextRow()) 610 if (Contains(vec, time)) 611 { 612 double avg = 0; 613 double rms = 0; 614 float min = 100; 615 float max = -100; 616 for (int i=beg; i<num; i++) 617 { 618 avg += temp[i]; 619 rms += temp[i]*temp[i]; 620 621 min = TMath::Min(min, temp[i]); 622 max = TMath::Max(max, temp[i]); 623 } 624 625 avg /= num-beg; 626 rms /= num-beg; 627 628 g1.SetPoint(g1.GetN(), time*24*3600, avg); 629 g1.SetPointError(g1.GetN()-1, 0, sqrt(rms-avg*avg)); 630 631 g2.SetPoint(g2.GetN(), time*24*3600, min); 632 g3.SetPoint(g3.GetN(), time*24*3600, max); 633 } 634 635 g1.SetLineColor(kGreen); 636 g1.DrawClone("[]"); 637 638 g2.SetLineColor(kGreen); 639 g2.SetLineWidth(1); 640 g2.DrawClone("L"); 641 642 g3.SetLineColor(kGreen); 643 g3.SetLineWidth(1); 644 g3.DrawClone("L"); 645 646 return 0; 647 } 648 649 void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0, TString outpath="./") 372 650 { 373 651 // To get correct dates in the histogram you have to add … … 377 655 { 378 656 UInt_t nt = MTime(MTime(-1).GetMjd()-1.5).GetNightAsInt(); 379 y = nt/10000;657 y = nt/10000; 380 658 m = (nt/100)%100; 381 d = nt%100;659 d = nt%100; 382 660 383 661 cout << y << "/" << m << "/" << d << endl; … … 409 687 { 410 688 TString query; 411 query += "SELECT fRunID, fRunStart, fRunStop FROM runinfo";689 query += "SELECT fRunID, fRunStart, fRunStop FROM RunInfo"; 412 690 query += " WHERE fNight="; 413 691 query += night; 414 query += " AND fRunTypeK EY=1 ORDER BY fRunStart";692 query += " AND fRunTypeKey=1 ORDER BY fRunStart"; 415 693 416 694 TSQLResult *res = serv.Query(query); … … 438 716 if (n==0) 439 717 cout << "WARNING - No data runs in db, displaying all data." << endl; 718 else 719 cout << "Num: " << n << "\n" << endl; 440 720 } 441 721 442 722 TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 960); 443 c->Divide(1, 5, 1e-5, 1e-5);723 c->Divide(1, 6, 1e-5, 1e-5); 444 724 445 725 gROOT->SetSelectedPad(0); 726 446 727 c->cd(1); 447 gPad->SetGrid x();728 gPad->SetGrid(); 448 729 gPad->SetTopMargin(0); 449 730 gPad->SetRightMargin(0.001); … … 454 735 gROOT->SetSelectedPad(0); 455 736 c->cd(2); 456 gPad->SetGrid x();737 gPad->SetGrid(); 457 738 gPad->SetTopMargin(0); 458 739 gPad->SetRightMargin(0.001); … … 463 744 gROOT->SetSelectedPad(0); 464 745 c->cd(3); 465 gPad->SetGrid x();746 gPad->SetGrid(); 466 747 gPad->SetTopMargin(0); 467 748 gPad->SetBottomMargin(0); … … 469 750 gPad->SetLeftMargin(0.04); 470 751 cout << PlotRate(runs, fname) << endl; 752 cout << PlotRateQC(night, serv) << endl; 471 753 472 754 gROOT->SetSelectedPad(0); 473 755 c->cd(5); 474 gPad->SetGridx(); 756 gPad->SetGrid(); 757 gPad->SetTopMargin(0); 758 gPad->SetBottomMargin(0); 759 gPad->SetRightMargin(0.001); 760 gPad->SetLeftMargin(0.04); 761 cout << PlotPointing(runs, fname) << endl; 762 763 gROOT->SetSelectedPad(0); 764 c->cd(6); 765 gPad->SetGrid(); 475 766 gPad->SetTopMargin(0); 476 767 gPad->SetRightMargin(0.001); 477 768 gPad->SetLeftMargin(0.04); 478 cout << PlotPointing(runs, fname) << endl; 479 /* 480 gROOT->SetSelectedPad(0); 481 c->cd(6); 482 gPad->SetGridx(); 483 gPad->SetTopMargin(0); 484 gPad->SetRightMargin(0.001); 485 gPad->SetLeftMargin(0.04); 486 cout << PlotMem(runs, fname) << endl; 487 */ 488 c->SaveAs(Form("%04d%02d%02d-quality.pdf", y, m, d)); 769 cout << PlotTemperature1(runs, fname) << endl; 770 cout << PlotTemperature2(runs, fname) << endl; 771 cout << PlotTemperature3(runs, fname) << endl; 772 cout << PlotTemperature4(runs, fname) << endl; 773 774 c->SaveAs(Form("%s/%04d%02d%02d-quality.png", outpath.Data(), y, m, d)); 489 775 } 490 776
Note:
See TracChangeset
for help on using the changeset viewer.