Changeset 17370 for trunk/Mars/fact
- Timestamp:
- 11/25/13 18:09:42 (11 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/fact/analysis/gain/extract_aux.C
r17326 r17370 9 9 #include "Interpolator2D.h" 10 10 11 #include "MMath.h" 11 12 #include "MHCamera.h" 12 13 #include "MGeomCamFACT.h" 13 14 14 int extract_temp(UInt_t day=20130902, UInt_t run=146, const char *output=0) 15 PixelMap pmap; 16 17 bool extract_current(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g) 15 18 { 16 TString nraw = Form("/fact/raw/%4d/%02d/%02d/%8d_%03d.fits.fz", 17 day/10000, (day/100)%100, day%100, day, run); 18 19 TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FSC_CONTROL_TEMPERATURE.fits", 19 TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FEEDBACK_CALIBRATED_CURRENTS.fits", 20 20 day/10000, (day/100)%100, day%100, day); 21 22 PixelMap pmap;23 if (!pmap.Read("FACTmap111030.txt"))24 {25 cout << "FACTmap111030.txt not found." << endl;26 return 1;27 }28 29 Interpolator2D interpol;30 31 const auto sens = Interpolator2D::ReadGrid("operation/sensor-pos.txt");32 if (sens.size()!=31)33 {34 cout << "Reading sensor-pos.txt failed: " << interpol.getInputGrid().size() << endl;35 return 2;36 }37 38 // =====================================================================39 40 zfits fraw(nraw.Data());41 if (!fraw)42 {43 cout << "Open raw file failed: " << nraw << endl;44 return 1;45 }46 47 double tstart = fraw.GetUInt("TSTARTI") + fraw.GetFloat("TSTARTF");48 double tstop = fraw.GetUInt("TSTOPI") + fraw.GetFloat("TSTOPF");49 21 50 22 // ===================================================================== … … 56 28 { 57 29 cout << "Could not open " << naux << endl; 58 return 1; 30 return false; 31 } 32 33 Double_t time; 34 Float_t I[416]; 35 36 if (!faux.SetPtrAddress("Time", &time)) 37 return false; 38 if (!faux.SetPtrAddress("I", I)) 39 return false; 40 41 TArrayD avg(320); 42 TArrayD rms(320); 43 TArrayI cnt(320); 44 45 while (faux.GetNextRow()) 46 { 47 if (time<tstart || time>tstop) 48 continue; 49 50 TArrayD med(320); 51 52 double mean = 0; 53 54 int j=0; 55 for (int i=0; i<320; i++) 56 { 57 if (i!=66 && i!=191 && i!=193 && i!=17 && i!=206) 58 med[j++] = I[i]; 59 60 avg[i] += I[i]; 61 rms[i] += I[i]*I[i]; 62 cnt[i]++; 63 } 64 65 Double_t median; 66 MMath::MedianDev(315, med.GetArray(), median); 67 g.SetPoint(g.GetN(), time, median); 68 } 69 70 if (cnt[0]==0) 71 return 2; 72 73 for (int i=0; i<320; i++) 74 { 75 avg[i] /= cnt[i]; 76 rms[i] /= cnt[i]; 77 rms[i] -= avg[i]*avg[i]; 78 rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]); 79 } 80 81 double mn = 0; 82 for (int i=0; i<1440; i++) 83 { 84 mn += avg[pmap.index(i).hv()]; 85 havg.SetBinContent(i+1, avg[pmap.index(i).hv()]); 86 hrms.SetBinContent(i+1, rms[pmap.index(i).hv()]); 87 } 88 89 return true; 90 } 91 92 bool extract_dac(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g) 93 { 94 TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.BIAS_CONTROL_DAC.fits", 95 day/10000, (day/100)%100, day%100, day); 96 97 // ===================================================================== 98 99 // Now we read the temperatures of the sensors during the 100 // single pe run and average them 101 fits faux(naux.Data()); 102 if (!faux) 103 { 104 cout << "Could not open " << naux << endl; 105 return false; 106 } 107 108 Double_t time; 109 UShort_t dac[416]; 110 111 if (!faux.SetPtrAddress("Time", &time)) 112 return false; 113 if (!faux.SetPtrAddress("U", dac)) 114 return false; 115 116 TArrayD avg(320); 117 TArrayD rms(320); 118 TArrayI cnt(320); 119 120 while (faux.GetNextRow()) 121 { 122 if (time<tstart || time>tstop) 123 continue; 124 125 TArrayD med(320); 126 127 for (int i=0; i<320; i++) 128 { 129 med[i] = dac[i]; 130 avg[i] += dac[i]; 131 rms[i] += dac[i]*dac[i]; 132 cnt[i]++; 133 } 134 135 Double_t median; 136 MMath::MedianDev(320, med.GetArray(), median); 137 g.SetPoint(g.GetN(), time, median); 138 } 139 140 if (cnt[0]==0) 141 return 2; 142 143 for (int i=0; i<320; i++) 144 { 145 avg[i] /= cnt[i]; 146 rms[i] /= cnt[i]; 147 rms[i] -= avg[i]*avg[i]; 148 rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]); 149 } 150 151 for (int i=0; i<1440; i++) 152 { 153 havg.SetBinContent(i+1, avg[pmap.index(i).hv()]); 154 hrms.SetBinContent(i+1, rms[pmap.index(i).hv()]); 155 } 156 157 return true; 158 } 159 160 bool extract_temp(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g) 161 { 162 TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FSC_CONTROL_TEMPERATURE.fits", 163 day/10000, (day/100)%100, day%100, day); 164 165 Interpolator2D interpol; 166 167 const auto sens = Interpolator2D::ReadGrid("operation/sensor-pos.txt"); 168 if (sens.size()!=31) 169 { 170 cout << "Reading sensor-pos.txt failed: " << interpol.getInputGrid().size() << endl; 171 return false; 172 } 173 174 // ===================================================================== 175 176 // Now we read the temperatures of the sensors during the 177 // single pe run and average them 178 fits faux(naux.Data()); 179 if (!faux) 180 { 181 cout << "Could not open " << naux << endl; 182 return false; 59 183 } 60 184 … … 63 187 64 188 if (!faux.SetPtrAddress("Time", &time)) 65 return -1;189 return false; 66 190 if (!faux.SetPtrAddress("T_sens", temp)) 67 return -1;191 return false; 68 192 69 193 TArrayD avg(31); … … 75 199 if (time<tstart || time>tstop) 76 200 continue; 201 202 double mean = 0; 203 int count = 0; 77 204 78 205 for (int i=0; i<31; i++) … … 83 210 double T = temp[i]; 84 211 212 mean += T; 213 count++; 214 85 215 avg[i] += T; 86 216 rms[i] += T*T; 87 217 cnt[i]++; 88 218 } 89 } 90 91 vector<double> z; 219 220 g.SetPoint(g.GetN(), time, mean/count); 221 222 } 223 224 vector<double> zavg; 225 vector<double> zrms; 92 226 93 227 vector<Interpolator2D::vec> pos; … … 98 232 avg[i] /= cnt[i]; 99 233 rms[i] /= cnt[i]; 100 rms[i] = sqrt(rms[i]-avg[i]*avg[i]); 234 rms[i] -= avg[i]*avg[i]; 235 rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]); 101 236 102 237 pos.push_back(sens[i]); 103 z.push_back(avg[i]); 238 zavg.push_back(avg[i]); 239 zrms.push_back(rms[i]); 104 240 } 105 241 } … … 112 248 { 113 249 cout << "Error setting output grid." << endl; 114 return 3;250 return false; 115 251 } 116 252 if (interpol.getOutputGrid().size()!=320) 117 253 { 118 254 cout << "Reading bias-positions.txt failed: " << interpol.getOutputGrid().size() << endl; 119 return 3;255 return false; 120 256 } 121 257 … … 129 265 sensors.SetPoint(i, input[i].x, input[i].y); 130 266 131 const vector<double> rc = interpol.Interpolate(z); 267 const vector<double> rc_avg = interpol.Interpolate(zavg); 268 const vector<double> rc_rms = interpol.Interpolate(zrms); 132 269 133 270 // ================================================================ 134 271 272 for (int i=0; i<1440; i++) 273 { 274 havg.SetBinContent(i+1, rc_avg[pmap.index(i).hv()]); 275 hrms.SetBinContent(i+1, rc_rms[pmap.index(i).hv()]); 276 } 277 278 return true; 279 280 } 281 282 void extract_aux(UInt_t day=20131114, UInt_t run=136, const char *output=0) 283 { 284 if (!pmap.Read("operation/FACTmap111030.txt")) 285 { 286 cout << "FACTmap111030.txt not found." << endl; 287 return; 288 } 289 290 // ===================================================================== 291 292 TString nraw = Form("/fact/raw/%4d/%02d/%02d/%8d_%03d.fits.fz", 293 day/10000, (day/100)%100, day%100, day, run); 294 295 296 zfits fraw(nraw.Data()); 297 if (!fraw) 298 { 299 cout << "Open raw file failed: " << nraw << endl; 300 return; 301 } 302 303 double tstart = fraw.GetUInt("TSTARTI") + fraw.GetFloat("TSTARTF"); 304 double tstop = fraw.GetUInt("TSTOPI") + fraw.GetFloat("TSTOPF"); 305 306 // ===================================================================== 307 308 MGeomCamFACT geom; 309 310 TGraph gcur; 311 MHCamera hcur_m(geom); 312 MHCamera hcur_r(geom); 313 if (!extract_current(day, tstart, tstop, hcur_m, hcur_r, gcur)) 314 return; 315 316 TGraph gtmp; 317 MHCamera htmp_m(geom); 318 MHCamera htmp_r(geom); 319 if (!extract_temp(day, tstart, tstop, htmp_m, htmp_r, gtmp)) 320 return; 321 322 TGraph gdac; 323 MHCamera hdac_m(geom); 324 MHCamera hdac_r(geom); 325 if (!extract_dac(day, tstart, tstop, hdac_m, hdac_r, gdac)) 326 return; 327 328 // ================================================================ 329 135 330 TCanvas *c = new TCanvas; 136 c->SetName("Temp"); 137 138 MGeomCamFACT geom; 139 140 MHCamera h1(geom); 141 h1.SetName("Temp"); 142 for (int i=0; i<1440; i++) 143 h1.SetBinContent(i+1, rc[pmap.index(i).hv()]); 144 145 h1.SetAllUsed(); 146 h1.DrawCopy(); 147 148 sensors.SetName("Sensors"); 149 sensors.SetMarkerColor(kWhite); 150 sensors.DrawClone("*"); 331 c->Divide(3, 3); 332 333 c->cd(1); 334 hcur_m.SetName("CurrentAvg"); 335 hcur_m.SetAllUsed(); 336 hcur_m.DrawCopy(); 337 c->cd(4); 338 hcur_r.SetName("CurrentRms"); 339 hcur_r.SetAllUsed(); 340 hcur_r.DrawCopy(); 341 c->cd(7); 342 gcur.SetName("CurrentTime"); 343 gcur.SetMarkerStyle(kFullDotMedium); 344 gcur.SetMinimum(0); 345 gcur.DrawClone("AP"); 346 347 c->cd(2); 348 htmp_m.SetName("TempAvg"); 349 htmp_m.SetAllUsed(); 350 htmp_m.DrawCopy(); 351 c->cd(5); 352 htmp_r.SetName("TempRms"); 353 htmp_r.SetAllUsed(); 354 htmp_r.DrawCopy(); 355 c->cd(8); 356 gtmp.SetName("TempTime"); 357 gtmp.SetMarkerStyle(kFullDotMedium); 358 gtmp.SetMinimum(0); 359 gtmp.DrawClone("AP"); 360 361 c->cd(3); 362 hdac_m.SetName("DacAvg"); 363 hdac_m.SetAllUsed(); 364 hdac_m.DrawCopy(); 365 c->cd(6); 366 hdac_r.SetName("DacRms"); 367 hdac_r.SetAllUsed(); 368 hdac_r.DrawCopy(); 369 c->cd(9); 370 gdac.SetName("DacTime"); 371 gdac.SetMarkerStyle(kFullDotMedium); 372 gdac.SetMinimum(0); 373 gdac.DrawClone("AP"); 151 374 152 375 if (output) 153 376 c->SaveAs(output); 154 155 return 0;156 377 }
Note:
See TracChangeset
for help on using the changeset viewer.