Changeset 17652 for trunk/Mars/fact
- Timestamp:
- 04/04/14 11:15:26 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/fact/plots/quality.C
r17521 r17652 1 #include <algorithm> 2 #include <functional> 3 1 4 Bool_t Contains(TArrayD **vec, Double_t t0, Double_t range=0) 2 5 { … … 87 90 return vecp[vecp.size()-1].second; 88 91 } 89 /* 92 90 93 Float_t Prediction(Double_t time) 91 94 { 92 95 Double_t jd = time + 40587 + 2400000.5; 93 94 Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);95 96 // Nova::EquPosn pos = FindPointing(time);97 98 // get local position of moon99 Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);100 101 // Distance between source and moon102 //const double angle = Nova::GetAngularSeparation(moon, hrzm)*TMath::DegToRad();103 104 // Current prediction105 //double cang = sin(angle);106 double calt = sin(hrzm.alt*TMath::DegToRad());107 108 double disk = Nova::GetLunarDisk(jd);109 110 // semi-major axis111 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 moon125 Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);126 127 // Distance between source and moon128 const double angle = Nova::GetAngularSeparation(moon, pos);129 130 // Distance between earth and moon relative to major semi-axis131 const double dist = Nova::GetLunarEarthDist(jd)/384400;132 133 // Current prediction134 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 condition140 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 96 145 97 // Sun properties … … 161 113 // Current prediction 162 114 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; 115 double cos_mdist = cos(angle*TMath::DegToRad()); 116 double sin_szd = sin(hrzs.zd*TMath::DegToRad()); 117 118 double c0 = pow(disk, 2.63); 119 double c1 = pow(sin_malt, 0.60); 120 double c2 = pow(edist, -2.00); 121 double c3 = exp(0.67*cos_mdist*cos_mdist*cos_mdist*cos_mdist); 122 double c4 = exp(-97.8+105.8*sin_szd*sin_szd); 123 124 double cur = 6.2 + 95.7*c0*c1*c2*c3 + c4; 175 125 176 126 return cur; 177 178 127 } 179 128 … … 228 177 Double_t time; 229 178 Float_t Imed; 179 Float_t Idev; 180 Float_t I[416]; 230 181 231 182 if (!file.SetPtrAddress("Time", &time)) … … 233 184 234 185 if (!file.SetPtrAddress("I_med", &Imed)) 186 return -1; 187 188 if (!file.SetPtrAddress("I_dev", &Idev)) 189 return -1; 190 191 if (!file.SetPtrAddress("I", I)) 235 192 return -1; 236 193 237 194 TGraph g1; 238 195 TGraph g2; 239 g1.SetName("Currents"); 196 TGraph g3; 197 TGraph g4; 198 TGraph g5; 199 g1.SetName("CurrentsMed"); 240 200 g2.SetName("Prediction"); 201 g3.SetName("CurrentsDev"); 202 g4.SetName("CurrentsMax-4"); 203 g5.SetName("CurrentsMax"); 241 204 242 205 while (file.GetNextRow()) 243 206 if (Contains(vec, time)) 244 207 { 208 // crazy pixels 209 I[66] = 0; 210 I[191] = 0; 211 I[193] = 0; 212 213 sort(I, I+320); 214 245 215 g1.SetPoint(g1.GetN(), time*24*3600, Imed); 246 216 g2.SetPoint(g2.GetN(), time*24*3600, Prediction(time)); 217 g3.SetPoint(g3.GetN(), time*24*3600, Imed+Idev); 218 g4.SetPoint(g4.GetN(), time*24*3600, I[315]); 219 g5.SetPoint(g5.GetN(), time*24*3600, I[319]); 247 220 } 248 221 249 222 g1.SetMinimum(0); 250 g1.SetMaximum( 99);223 g1.SetMaximum(149); 251 224 252 225 g1.SetMarkerStyle(kFullDotMedium); … … 260 233 g1.DrawClone("AP"); 261 234 235 g5.SetMarkerColor(kGray); 236 g5.DrawClone("P"); 237 238 g4.SetMarkerColor(kGray+1); 239 g4.DrawClone("P"); 240 241 g3.SetMarkerColor(kGray+2); 242 g3.DrawClone("P"); 243 262 244 g2.SetMarkerColor(kBlue); 263 245 g2.SetMarkerStyle(kFullDotMedium); 264 g2.GetXaxis()->SetTimeDisplay(true);265 g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");266 g2.GetXaxis()->SetLabelSize(0.12);267 g2.GetYaxis()->SetLabelSize(0.1);268 g2.GetYaxis()->SetTitle("CURRENT");269 g2.GetYaxis()->SetTitleOffset(0.2);270 g2.GetYaxis()->SetTitleSize(0.1);271 246 g2.DrawClone("P"); 272 247 … … 647 622 } 648 623 649 void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0, TString outpath="./")624 void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0, const char *outpath="quality") 650 625 { 651 626 // To get correct dates in the histogram you have to add … … 772 747 cout << PlotTemperature4(runs, fname) << endl; 773 748 774 c->SaveAs(Form("%s/%04d%02d%02d-quality.png", outpath.Data(), y, m, d)); 775 } 776 777 //20130314_141 749 c->SaveAs(Form("%s/%04d%02d%02d.png", outpath, y, m, d)); 750 }
Note:
See TracChangeset
for help on using the changeset viewer.