Index: /trunk/Mars/fact/processing/calibrate.C =================================================================== --- /trunk/Mars/fact/processing/calibrate.C (revision 15227) +++ /trunk/Mars/fact/processing/calibrate.C (revision 15227) @@ -0,0 +1,288 @@ +/* +#ifndef HAVE_ZLIB +#define HAVE_ZLIB +#endif + +#include + +#include "fits.h" +#include "ofits.h" + +#include "MArrayF.h" +#include "MMath.h" +*/ + +using namespace std; + +bool Read(const char *fname, int *npix) +{ + std::ifstream fin(fname); + + int l = 0; + + while (1) + { + if (l>1439) + break; + + TString line; + line.ReadLine(fin); + if (!fin) + break; + + line = line.Strip(TString::kBoth); + if (line[0]=='#') + continue; + + TObjArray *arr = line.Tokenize(' '); + int cbpx = atoi(((*arr)[1])->GetName()); + int board = atoi(((*arr)[6])->GetName()); + int ch = atoi(((*arr)[7])->GetName()); + delete arr; + + int hv = board*32+ch; + if (hv>=416) + return false; + + int n = cbpx%10>3 ? 5 : 4; + if (npix[hv]>0 && npix[hv]!=n) + return false; + + npix[hv] = n; + l++; + } + + return l==1440; +} + +int calibrate(UInt_t night, TString path="/fact/aux") +{ + int npix[416]; + memset(npix, 0, sizeof(int)*416); + if (!Read("/scratch_nfs/FACTmap111030.txt", npix)) + return 42; + + UInt_t d = night%100; + UInt_t m = (night/100)%100; + UInt_t y = night/10000; + + path += Form("/%04d/%02d/%02d/%d", y, m, d, night); + TString path2 = Form("/scratch_nfs/calibrated_currents/%d", night); + + TString filec = path + ".FEEDBACK_CALIBRATION.fits"; + TString fileu = path + ".BIAS_CONTROL_VOLTAGE.fits"; + TString filei = path + ".BIAS_CONTROL_CURRENT.fits"; + TString fileo = path2 + ".CALIBRATED_CURRENTS.fits"; + + cout << filec << endl; + cout << fileu << endl; + cout << filei << endl; + cout << fileo << endl; + + fits fitsc1(filec.Data()); + fits fitsc2(filec.Data()); + if (!fitsc2) + { + cout << "ERROR - Open " << filec << " failed." << endl; + return 1; + } + + fits fitsu1(fileu.Data()); + fits fitsu2(fileu.Data()); + if (!fitsu2) + { + cout << "ERROR - Open " << fileu << " failed." << endl; + return 2; + } + + fits fitsi1(filei.Data()); + fits fitsi2(filei.Data()); + if (!fitsi2) + { + cout << "ERROR - Open " << filei << " failed." << endl; + return 3; + } + + ofits fout(fileo.Data()); + if (!fout) + { + cout << "ERROR - Open " << fileo << " failed." << endl; + return 4; + } + + //const MTime now; + fout.SetStr( "TELESCOP", "FACT", "Telescope that acquired this data"); + fout.SetStr( "PACKAGE", "MARS", "Package name"); + fout.SetStr( "VERSION", "1.0", "Package description"); + fout.SetFloat("EXTREL", 1.0, "Release Number"); + //fout.SetStr( "COMPILED", __DATE__" "__TIME__, "Compile time"); + fout.SetStr( "ORIGIN", "FACT", "Institution that wrote the file"); + //fout.SetStr( "DATE", now.GetStringFmt("%Y-%m-%dT%H:%M:%S").Data(), "File creation date"); + //fout.SetInt( "NIGHT", now.GetNightAsInt(), "Night as int"); + fout.SetStr( "TIMESYS", "UTC", "Time system"); + fout.SetStr( "TIMEUNIT", "d", "Time given in days w.r.t. to MJDREF"); + fout.SetInt( "MJDREF", 40587, "MJD to UNIX time (seconds since 1970/1/1)"); + + fout.AddColumnDouble("Time", "MJD", "Event time: MJD - MJDREF (= Unix Time)"); + fout.AddColumnDouble("T_cal", "MJD", "Calibration time: MJD - MJDREF (= Unix Time)"); + fout.AddColumnInt( "N", "", "Number of valid values"); + fout.AddColumnFloat( "I_avg", "uA", "Average calibrated current (320 channels)"); + fout.AddColumnFloat( "I_rms", "uA", "RMS of calibrated current (320 channels)"); + fout.AddColumnFloat( "I_med", "uA", "Median of calibrated currents (320 channels)"); + fout.AddColumnFloat( "I_dev", "uA", "Deviation of calibrated currents (320 channels)"); + fout.AddColumnFloat(416, "I", "uA", "My comment"); + + fout.WriteTableHeader("CALIBRATED"); + + Double_t timec, timeu, timei; + Double_t nextc, nextu, nexti; + + if (!fitsc1.SetPtrAddress("Time", &timec)) + { + cout << "ERROR couldn't get pointer to 1st Time for R" << endl; + return 2; + } + if (!fitsu1.SetPtrAddress("Time", &timeu)) + { + cout << "ERROR couldn't get pointer to 1st Time for U" << endl; + return 2; + } + if (!fitsi1.SetPtrAddress("Time", &timei)) + { + cout << "ERROR couldn't get pointer to 1st Time for I" << endl; + return 2; + } + + if (!fitsc2.SetPtrAddress("Time", &nextc)) + { + cout << "ERROR couldn't get pointer to 2nd Time for R" << endl; + return 2; + } + if (!fitsu2.SetPtrAddress("Time", &nextu)) + { + cout << "ERROR couldn't get pointer to 2nd Time for U" << endl; + return 2; + } + if (!fitsi2.SetPtrAddress("Time", &nexti)) + { + cout << "ERROR couldn't get pointer to 2nd Time for I" << endl; + return 2; + } + + Float_t R[416]; + Float_t U[416]; + UShort_t I[416]; + + if (!fitsc1.SetPtrAddress("R", R)) + { + cout << "ERROR couldn't get pointer to R" << endl; + return 2; + } + if (!fitsu1.SetPtrAddress("Uout", U)) + { + cout << "ERROR couldn't get pointer to U" << endl; + return 2; + } + if (!fitsi1.SetPtrAddress("I", I)) + { + cout << "ERROR couldn't get pointer to I" << endl; + return 2; + } + + fitsc2.GetNextRow(); + fitsu2.GetNextRow(); + fitsi2.GetNextRow(); + + if (!fitsc1.GetNextRow()) + return 0; + if (!fitsc2.GetNextRow()) + nextc = 1e50; + + while (1) + { + // First the next two rows (voltage, current) + // with identical time stamp + if (nextu0) + { + cur[i] /= npix[i]; + stat[cnt++] = cur[i]; + } + } + + memcpy(rc.GetArray(), &timeu, sizeof(Double_t)); + memcpy(rc.GetArray()+2, &timec, sizeof(Double_t)); + memcpy(rc.GetArray()+4, &cnt, sizeof(UInt_t)); + + Double_t med; + rc[5] = TMath::Mean(cnt, stat.GetArray()); + rc[6] = TMath::RMS(cnt, stat.GetArray()); + rc[8] = MMath::MedianDev(cnt, stat.GetArray(), med); + rc[7] = med; + + fout.WriteRow(rc.GetArray(), rc.GetSize()*sizeof(Float_t)); + } + + return 0; +}