Index: branches/MarsGapdTimeJitter/ceres.rc =================================================================== --- branches/MarsGapdTimeJitter/ceres.rc (revision 18306) +++ branches/MarsGapdTimeJitter/ceres.rc (revision 18331) @@ -134,4 +134,11 @@ ConesTransmission.FileName: resmc/fact/Transmittance_1439Cones_FACT_bearbeitet.txt +# A backward motivated overall PDE adjustment. Introduced by Fabian T. to produce the ICRC2015 Crab spectrum. Fabian, please comment a bit more on this. In case this is only about the additional shielding of our camera support structure, then we schould actually name it like this. + +AdditionalPhotonAcceptance.Function.Name: 0.8 +AdditionalPhotonAcceptance.Function.Npx: 100L +AdditionalPhotonAcceptance.Function.Xmin: 290L +AdditionalPhotonAcceptance.Function.Xmax: 900L + # ------------------------------------------------------------------------- # Setup what in MMCS would be called "additional spot size" @@ -238,4 +245,3 @@ ResidualTimeSpread.Val: 0.0 - # last line comment Index: branches/MarsGapdTimeJitter/fact/analysis/callisto.C =================================================================== --- branches/MarsGapdTimeJitter/fact/analysis/callisto.C (revision 18306) +++ branches/MarsGapdTimeJitter/fact/analysis/callisto.C (revision 18331) @@ -21,4 +21,47 @@ MStatusDisplay *d = new MStatusDisplay; + + // ------------------------------------------------------ + + // ------------------------------------------------------ + + // Calib: 51 / 90 / 197 (20% TH) + // Data: 52 / 64 / 104 (20% TH) + + // Extraction range in slices. It will always(!) contain the full range + // of integration + const int first_slice = 25; // 10ns + const int last_slice = 225; // 125ns + + const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 } + + Long_t max = 0; // All + Long_t max0 = max; // Time marker + Long_t max1 = max; // Light pulser + //Long_t max2 = 3000; // Calibration ratio + Long_t max3 = max; // Pedestal Rndm + Long_t max4 = max; // Pedestal Ext + Long_t max5 = max; // Data + + // ========================= Result ================================== + + //double scale = 0.1; + double scale = 0.1024; + + // ====================================================== + + if (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 1; + } + + // The sequence file which defines the files for the analysis + MSequence seq(seqfile); + if (!seq.IsValid()) + { + gLog << err << "ERROR - Sequence '" << seqfile << "' invalid!" << endl; + return 2; + } MBadPixelsCam badpixels; @@ -38,45 +81,26 @@ // 1393 - // ------------------------------------------------------ - - // ------------------------------------------------------ - - // Calib: 51 / 90 / 197 (20% TH) - // Data: 52 / 64 / 104 (20% TH) - - // Extraction range in slices. It will always(!) contain the full range - // of integration - const int first_slice = 25; // 10ns - const int last_slice = 225; // 125ns - - const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 } - - Long_t max = 0; // All - Long_t max0 = max; // Time marker - Long_t max1 = max; // Light pulser - //Long_t max2 = 3000; // Calibration ratio - Long_t max3 = max; // Pedestal Rndm - Long_t max4 = max; // Pedestal Ext - Long_t max5 = max; // Data - - // ========================= Result ================================== - - //double scale = 0.1; - double scale = 0.1024; - - // ====================================================== - - if (map && gSystem->AccessPathName(map, kFileExists)) - { - gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl; - return 1; - } - - // The sequence file which defines the files for the analysis - MSequence seq(seqfile); - if (!seq.IsValid()) - { - gLog << err << "ERROR - Sequence '" << seqfile << "' invalid!" << endl; - return 2; + // broken drs-board: + // https://www.fact-project.org/logbook/showthread.php?tid=3521 + // in the analysis effect visible 10.1. - 25.5.2015 + // according to logbook: + // first report of problem: + // https://www.fact-project.org/logbook/showthread.php?tid=2772&pid=15756#pid15756 + // repair action: + // https://www.fact-project.org/logbook/showthread.php?tid=3298 + // affected pixels: + // SW: 1193 1194 1195 1391 1392 1393 1304 1305 1306 + if (seq.GetNight().GetNightAsInt()>20141114 && seq.GetNight().GetNightAsInt()<20150526) + { + gLog << "Exclude broken drs-board..." << endl; + badpixels[1193].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1194].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1195].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1391].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1392].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1393].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1304].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1305].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1306].SetUnsuitable(MBadPixelsPix::kUnsuitable); } Index: branches/MarsGapdTimeJitter/fact/processing/buildseqentries.C =================================================================== --- branches/MarsGapdTimeJitter/fact/processing/buildseqentries.C (revision 18306) +++ branches/MarsGapdTimeJitter/fact/processing/buildseqentries.C (revision 18331) @@ -136,4 +136,5 @@ if (arr.GetEntries()>2) fMax = GetId(arr[2]->GetName()); + cout << "----> " << fMin << " " << fMax << endl; } @@ -217,4 +218,5 @@ if (run>0) { + cout << " - - > " << run << endl; if (runfMax) return kFALSE; @@ -915,5 +917,5 @@ } - TString PrepareString(TSQLResult &res, TArrayL &runs) + TString PrepareString(TSQLResult &res, TArrayL &runs, TArrayL &nights) { // Number of result rows @@ -921,4 +923,5 @@ runs.Set(rows); // initialize size of array for run numbers + nights.Set(rows); // initialize size of array for nights TArrayD start(rows); // instantiate array for start times @@ -932,4 +935,5 @@ { runs[idx] = atol((*row)[0]); // run number + nights[idx] = atol((*row)[4]); // night const TString tstart = ((*row)[2]); // start time @@ -1021,5 +1025,5 @@ /*const*/ TString query = - "SELECT fRunID, fRunTypeKEY, fRunStart, fRunStop" + "SELECT fRunID, fRunTypeKEY, fRunStart, fRunStop, fNight*1000 " " FROM RunInfo WHERE fExcludedFDAKEY=1 AND "; query += Form("fRunID BETWEEN %ld AND %ld AND " @@ -1038,5 +1042,6 @@ // and an array with the corresponding run-numbers TArrayL runs; - const TString str = PrepareString(*res, runs); + TArrayL nights; + const TString str = PrepareString(*res, runs, nights); delete res; @@ -1071,5 +1076,5 @@ { // Check if regular expressions matchs - const Ssiz_t len = obj->Match(str, pos, runs[pos]); + const Ssiz_t len = obj->Match(str, pos, nights[pos]); if (len>0) { Index: branches/MarsGapdTimeJitter/fact/processing/camhum.C =================================================================== --- branches/MarsGapdTimeJitter/fact/processing/camhum.C (revision 18331) +++ branches/MarsGapdTimeJitter/fact/processing/camhum.C (revision 18331) @@ -0,0 +1,65 @@ +void camhum(const char *fname, double beg=0, double end=100000) +{ + fits file(fname); + + //file.PrintColumns(); + //file.PrintKeys(); + + Double_t time; + Float_t hum[4]; + file.SetPtrAddress("Time", &time); + file.SetPtrAddress("H", hum); + + UInt_t offset = file.GetUInt("MJDREF"); + if (beg < 30000) + beg+=offset; + if (end < 30000) + end+=offset; + + double avg = 0; + int cnt = 0; + + double last = -1; + double diff = -1; + + while (file.GetNextRow()) + { + time += offset; + + if (time>end) + break; + + int n = 0; + double h = 0; + for (int i=0; i<4; i++) + if (hum[i]>0) + { + h += hum[i]; + n++; + } + + if (n==0) + continue; + + if (timeend) + break; + + if (timeavg*avg ? 0 : sqrt(rms-avg*avg); + + cout << "result " << avg << " " << rms << endl; +} Index: branches/MarsGapdTimeJitter/fact/processing/lastth.C =================================================================== --- branches/MarsGapdTimeJitter/fact/processing/lastth.C (revision 18331) +++ branches/MarsGapdTimeJitter/fact/processing/lastth.C (revision 18331) @@ -0,0 +1,33 @@ +void lastth(const char *fname, double beg=0) +{ + fits file(fname); + + //file.PrintColumns(); + //file.PrintKeys(); + + Double_t time; + UShort_t th; + file.SetPtrAddress("Time", &time); + file.SetPtrAddress("threshold", &th); + + UInt_t offset = file.GetUInt("MJDREF"); + if (beg<30000) + beg += offset; + + UShort_t th_last = 0; + Double_t time_last = -1; + + while (file.GetNextRow()) + { + time += offset; + + if (time>beg && fabs(beg-time)>fabs(time_last)) + break; + + th_last = th; + time_last = beg-time; + } + + if (th_last>0) + cout << "result " << th_last << " " << TMath::Nint(time_last*24*3600) << endl; +} Index: branches/MarsGapdTimeJitter/fact/processing/threshold.C =================================================================== --- branches/MarsGapdTimeJitter/fact/processing/threshold.C (revision 18331) +++ branches/MarsGapdTimeJitter/fact/processing/threshold.C (revision 18331) @@ -0,0 +1,83 @@ +#include + +void threshold(const char *fname, double beg=0, double end=100000) +{ + fits file(fname); + + //file.PrintColumns(); + //file.PrintKeys(); + + Double_t time; + UShort_t th[160]; + file.SetPtrAddress("Time", &time); + file.SetPtrAddress("PatchThresh", th); + + UInt_t offset = file.GetUInt("MJDREF"); + if (beg < 30000) + beg+=offset; + if (end < 30000) + end+=offset; + + double average = 0; + double med_avg = 0; + double med_rms = 0; + UShort_t max_tot = 0; + int cnt = 0; + + double avg_last = -1; + double med_last = -1; + UShort_t max_last = -1; + double diff = -1; + + while (file.GetNextRow()) + { + time += offset; + + if (time>end) + break; + + sort(th, th+160); + + Float_t med = (th[80] + th[81]) / 2; + UShort_t max = th[159]; + + Double_t avg = 0; + for (int i=0; i<160; i++) + avg += th[i]; + avg /= 160; + + if (timemax_tot) + max_tot = max; + + cnt ++; + } + + if (cnt==0) + { + if (diff<5./24/3600) + return; + + cout << "result " << med_last << " 0 " << max_last << " " << avg_last << endl; + return; + } + + average /= cnt; + med_avg /= cnt; + med_rms /= cnt; + med_rms = sqrt(med_rms-med_avg*med_avg); + + cout << "result " << med_avg << " " << med_rms << " " << max_tot << " " << average << endl; +} Index: branches/MarsGapdTimeJitter/fact/resources/sequences.rc =================================================================== --- branches/MarsGapdTimeJitter/fact/resources/sequences.rc (revision 18306) +++ branches/MarsGapdTimeJitter/fact/resources/sequences.rc (revision 18331) @@ -225,4 +225,6 @@ BGBB?TU?PP? 0 PCD+ 0 +BGPP? 20150716 +P?D+ 20150716 #[bB][gG][bB][bB][tT][pP][pP] 0 20120116 Index: branches/MarsGapdTimeJitter/fact/tools/DataFromQLA.C =================================================================== --- branches/MarsGapdTimeJitter/fact/tools/DataFromQLA.C (revision 18331) +++ branches/MarsGapdTimeJitter/fact/tools/DataFromQLA.C (revision 18331) @@ -0,0 +1,222 @@ +#include +#include +#include + +#include +#include +#include +#include + +#include "MSQLMagic.h" +#include "MTime.h" +#include "MMath.h" +#include "MDirIter.h" + +using namespace std; + +void PrintRemarks(ofstream &fout, TString sourcename, TString firstnight, TString lastnight) +{ + MTime now(-1); + + fout << "# " << endl; + fout << "# Please cite the FACT design paper and the QLA webpage when using these data." << endl; + fout << "# FACT design paper: H. Anderhub et al. JINST 8 (2013) P6008 " << endl; + fout << "# http://iopscience.iop.org/1748-0221/8/06/P06008 " << endl; + fout << "# QLA webpage: http://www.fact-project.org/monitoring" << endl; + fout << "# If you intent to use the data, please let us know for reference. " << endl; + fout << "# " << endl; + fout << "# Remarks:" << endl; + fout << "# * These are the results of a fast quick look analyis " << endl; + fout << "# on site, i.e. they are preliminary. " << endl; + fout << "# * The quick look analysis includes all data, " << endl; + fout << "# i.e. no data selection done." << endl; + fout << "# * The given values are not fluxes but excess rates " << endl; + fout << "# (number of excess events per effective ontime), " << endl; + fout << "# i.e. there is a dependence on trigger threshold and " << endl; + fout << "# zenith distance of the observation (with the current " << endl; + fout << "# analysis for zenith distance > 40 degree and trigger " << endl; + fout << "# threshold > 500 DAC counts)." << endl; +// fout << "# * The data are provided with 20 min binning and nightly binning." << endl; + fout << "# * Nights with less than 20 minutes of data are neglected. " << endl; + fout << "# * The QLA results are not reprocessed when a new software " << endl; + fout << "# version is introduced. " << endl; + fout << "# * In case, you need further details about the data or a" << endl; + fout << "# different binning, please do not hesitate to contact us." << endl; + fout << "# * The QLA contains all data since 12.12.2012. " << endl; + fout << "# For older data, please contact us. " << endl; + fout << "# " << endl; + fout << "# Contact: Daniela Dorner dorner@astro.uni-wuerzburg.de " << endl; + fout << "# " << endl; + fout << "# This file was created at " << now.GetString() << endl; + fout << "# Source: " << sourcename << endl; + fout << "# Timerange: " << firstnight << " - " << lastnight << endl; + fout << "# " << endl; + fout << "# start(mjd) stop(mjd) excess-rate(evts/h) error-excess-rate " << endl; + + return; +} + +int DataFromQLA(Int_t sourcekey=1, Int_t nightmin=2011115, Int_t nightmax=20161231) +{ + MSQLServer serv("sql.rc"); + if (!serv.IsConnected()) + { + cout << "ERROR - Connection to database failed." << endl; + return 0; + } + Bool_t dch=kFALSE; + + TString query=Form("SELECT fSourceName FROM Source WHERE fSourceKey=%d", sourcekey); + TSQLResult *res1 = serv.Query(query); + if (!res1) + return 1; + TSQLRow *row1=res1->Next(); + TString sourcename=(*row1)[0]; + delete res1; + sourcename.ReplaceAll(" ", "_"); + + // datacheck + TString datacheck=" "; + //remove data with wrong settings + datacheck+=" AND fNight>20120420 AND NOT fNight IN (20120406,20120410,20120503) AND";//data with different bias voltage + datacheck+=" NOT fNight BETWEEN 20121206 AND 20130110"; // broken bias channel + //datacheck+=" AND NOT (fNight=20120608 AND fRunID=65) "; // something went wrong with tracking? + // 24.6. new coefficients + TString zdparam=" pow(0.753833 * cos(Radians(fZenithDistanceMean)), 7.647435) * exp(-5.753686*pow(Radians(fZenithDistanceMean),2.089609))"; + TString thparam=" pow((if(isnull(fThresholdMinSet),fThresholdMedian,fThresholdMinSet)-329.4203),2) * (-0.0000002044803) "; + TString param=" (fNumEvtsAfterBgCuts/5-fNumSigEvts)/fOnTimeAfterCuts - "+zdparam+" - "+thparam+" "; + datacheck+=" AND -0.085 < ("+param+") "; + datacheck+=" AND 0.25 > ("+param+") "; + + TString select =" SELECT Min(fNight), Max(fNight), Sum(fOnTimeAfterCuts) as bla "; + TString fromjoinwhere=" FROM AnalysisResultsRunLP "; + fromjoinwhere+=" LEFT JOIN RunInfo USING(fNight, fRunID) "; + fromjoinwhere+=Form(" WHERE fSourceKey=%d", sourcekey); + fromjoinwhere+=Form(" AND fNight BETWEEN %d AND %d ", nightmin, nightmax); + fromjoinwhere+=" AND fOnTimeAfterCuts < 1000 "; //exclude runs with wrong/too high ontime + fromjoinwhere+=" AND fOnTimeAfterCuts > 10 "; //exclude runs with wrong/too small ontime + fromjoinwhere+=" AND NOT ISNULL(fNumExcEvts) ";// only where excess was extracted + if (dch) + fromjoinwhere+=datacheck; + query=select+fromjoinwhere; + query+=" GROUP BY fSourceKey "; + query+=" HAVING bla>20*60 "; + cout << "Q: " << query << endl; + TSQLResult *res2 = serv.Query(query); + if (!res2) + return 1; + TSQLRow *row2=res2->Next(); + TString firstnight=(*row2)[0]; + TString lastnight=(*row2)[1]; + delete res2; + + TString filename_nightly="FACT_QuickLookAnalysisResults_NightlyBinning_"+sourcename+".txt"; + TString filename_min="FACT_QuickLookAnalysisResults_20MinuteBinning_"+sourcename+".txt"; + + ofstream nightly(filename_nightly); + ofstream min(filename_min); + if (!nightly) + { + cout << "ERROR - cannot write " << filename_nightly << endl; + return 2; + } + if (!min) + { + cout << "ERROR - cannot write " << filename_min << endl; + return 2; + } + PrintRemarks(nightly, sourcename, firstnight, lastnight); + PrintRemarks(min, sourcename, firstnight, lastnight); + + // query data from AnalysisResultsRunLP + select =" SELECT Sum(fOnTimeAfterCuts) as bla, Sum(fNumExcEvts), Sum(fNumBgEvts), Sum(fNumSigEvts), "; + select+=" Min(fRunStart), Max(fRunStop), fNight, fRunID "; + query=select+fromjoinwhere; + query+=" GROUP BY fNight "; + query+=" HAVING bla>20*60 "; + query+=" ORDER BY fNight "; + //cout << "Q: " << query << endl; + + //variables for calculations and graphs + TString excevts, bgevts, sigevts, ontime, night; + Float_t excevtssum=0; + Float_t bgevtssum=0; + Float_t sigevtssum=0; + Float_t ontimesum=0; + Float_t excerr=0; + MTime start; + MTime stop; + Double_t mjdstart=0; + Double_t mjdstop=0; + //Double_t mjd=0; + //Double_t mjddiff=0; + Int_t counter=0; + + TSQLResult *res3 = serv.Query(query); + if (!res3) + return 1; + TSQLRow *row3=0; + while ((row3=res3->Next())) + { + ontime=(*row3)[0]; + excevts=(*row3)[1]; + bgevts=(*row3)[2]; + sigevts=(*row3)[3]; + start.SetSqlDateTime((*row3)[4]); + stop.SetSqlDateTime((*row3)[5]); + night=(*row3)[6]; + + mjdstart = start.GetMjd(); + mjdstop = stop.GetMjd(); + //mjd = mjdstart+(mjdstop-mjdstart)/2.; + + //significance = MMath::SignificanceLiMaSigned(sigevts.Atof(), bgevts.Atof()*5, 0.2); + excerr = MMath::ErrorExc(sigevts.Atof(), bgevts.Atof()*5, 0.2)/ontime.Atof()*3600.; + nightly.precision(18); + //fout << (*row3)[6] << ": " << start.GetMjd() << " " << stop.GetMjd() << " " << excevts.Atof()/ontime.Atof() << " " << excerr << endl; + nightly << start.GetMjd() << " " << stop.GetMjd() << " " << excevts.Atof()/ontime.Atof()*3600 << " " << excerr << endl; + query="SELECT fOnTimeAfterCuts, fNumExcEvts, fNumBgEvts, fNumSigEvts, "; + query+=" fRunStart, fRunStop "+fromjoinwhere; + query+=" AND fNight="+night; + //cout << "Q: " << query << endl; + TSQLResult *res4 = serv.Query(query); + if (!res4) + return 1; + TSQLRow *row4=0; + counter=0; + while ((row4=res4->Next())) + { + ontime=(*row4)[0]; + excevts=(*row4)[1]; + bgevts=(*row4)[2]; + sigevts=(*row4)[3]; + if(counter==0) + start.SetSqlDateTime((*row4)[4]); + + if (ontimesum+ontime.Atoi()>20*60) + { + excerr = MMath::ErrorExc(sigevtssum, bgevtssum*5, 0.2)/ontimesum*3600.; + min.precision(18); + min << start.GetMjd() << " " << stop.GetMjd() << " " << excevtssum/ontimesum*3600 << " " << excerr << endl; + start.SetSqlDateTime((*row4)[4]); + ontimesum=0; + excevtssum=0; + bgevtssum=0; + sigevtssum=0; + } + counter++; + stop.SetSqlDateTime((*row4)[5]); + ontimesum+=ontime.Atoi(); + excevtssum+=excevts.Atof(); + bgevtssum+=bgevts.Atof(); + sigevtssum+=sigevts.Atof(); + } + } + + delete res3; + return 0; + +} + + + Index: branches/MarsGapdTimeJitter/fact/tools/GetObsTimes.C =================================================================== --- branches/MarsGapdTimeJitter/fact/tools/GetObsTimes.C (revision 18331) +++ branches/MarsGapdTimeJitter/fact/tools/GetObsTimes.C (revision 18331) @@ -0,0 +1,141 @@ +// Be aware that this macro just inserts and not updates the rows + +void GetObsTimes(Bool_t dummy=kTRUE, Int_t start=55866 /*1.11.2011*/, Int_t stop=58848 /*31.12.2019*/) +{ + MSQLMagic serv("sql.rc"); + Bool_t con = serv.IsConnected(); + if (con) + { + cout << "Connected to " << serv.GetName() << endl; + cout << endl; + } + else + { + cout << "Could not connect to " << serv.GetName() << endl; + cout << endl; + return; + } + + + MTime time(-1); + MTime time2, time3, time4, time5; + MTime time6, time7, time8, time9; + Double_t darknight=0; + Int_t night=0; + TString query; + for (int i=start; i<=stop; i++) + { + //GetSolarRst returns the information for the next day + // as JD seems to be rounded to the next value in all + // cases in the libnova function ln_get_solar_rst + //to get the correct values the following would be needed + //Nova::RstTime t1= Nova::GetSolarRst(time.GetJD()-1, -12); + //Nova::RstTime t2= Nova::GetSolarRst(time.GetJD(), -12); + //time2.SetJD(t1.set); + //time3.SetJD(t2.rise); + //however, as time.GetNightAsInt() returns the date of the + // sunrise (MAGIC naming convention), we use: + time.SetMjd(i); + Nova::RstTime t1= Nova::GetSolarRst(time.GetJD(), -12); + Nova::RstTime t2= Nova::GetSolarRst(time.GetJD()+1, -12); + Nova::RstTime t3= Nova::GetSolarRst(time.GetJD(), -18); + Nova::RstTime t4= Nova::GetSolarRst(time.GetJD()+1, -18); + Nova::RstTime t5= Nova::GetLunarRst(time.GetJD()); + Nova::RstTime t6= Nova::GetLunarRst(time.GetJD()+1); + + /* + cout << time.GetNightAsInt() +// << " " << MTime(t1.set).GetSqlDateTime() +// << " " << MTime(t2.rise).GetSqlDateTime() +// << " s" << MTime(t3.set).GetSqlDateTime() +// << " r" << MTime(t4.rise).GetSqlDateTime() + << " s" << MTime(t5.set).GetSqlDateTime() + << " r" << MTime(t5.rise).GetSqlDateTime() + << " s" << MTime(t6.set).GetSqlDateTime() + << " r" << MTime(t6.rise).GetSqlDateTime() + << endl; + */ + + // astronomical twilight (-12 deg) + time2.SetJD(t1.set); + time3.SetJD(t2.rise); + // dark night (-18 deg) + time4.SetJD(t3.set); + time5.SetJD(t4.rise); + // setting rise (time6) and set (time7) of moon + if (t5.setrise + { + if (t5.set>t3.set) + time7.SetJD(t5.set); + else + time7.SetJD(t6.set); + + time6.SetJD(t6.rise); + } + // setting times for dark night + // sunset after moonrise + if (time4>time6) + time8=time7; + else + { + if (time4time7) + time8=time7; + else + time8=time4; + } + // moonset before sunrise + if (time7 r" << time6.GetSqlDateTime() + << " s" << time7.GetSqlDateTime() + << " ==> start" << time8.GetSqlDateTime() + << " stop" << time9.GetSqlDateTime() + << " dark: " << darknight << "h " + << endl; + */ + + query=Form("INSERT ObservationTimes SET fNight=%d, " + "fStartObservation='%s', fStopObservation='%s', " + "fStartDarkTime='%s', fStopDarkTime='%s', " + "fStartDarkNight='%s', fStopDarkNight='%s', " + "fMoonRise='%s', fMoonSet='%s', " + "fNumDarkHours=%f, fMjd=%d", + time.GetNightAsInt(), + time2.GetSqlDateTime().Data(), time3.GetSqlDateTime().Data(),// astronomical twilight + time4.GetSqlDateTime().Data(), time5.GetSqlDateTime().Data(),// darktime + time8.GetSqlDateTime().Data(), time9.GetSqlDateTime().Data(),// darknight + time6.GetSqlDateTime().Data(), time7.GetSqlDateTime().Data(),// moon + darknight, i); //dark hours + + cout << "QUERY: " << query << endl; + + if (dummy) + continue; + TSQLResult *res = serv.Query(query); + if (!res) + { + cout << "problem to insert. " << endl; + return; + } + } + return; + +} Index: branches/MarsGapdTimeJitter/fact/tools/produce_drs_time_fits_file.C =================================================================== --- branches/MarsGapdTimeJitter/fact/tools/produce_drs_time_fits_file.C (revision 18331) +++ branches/MarsGapdTimeJitter/fact/tools/produce_drs_time_fits_file.C (revision 18331) @@ -0,0 +1,165 @@ +// The FACT collaboration +// Dominik Neise and Sebastian Mueller October 2014 +// A callisto like CERN-Root script to process the FACT DRS +// (Domino-Ring-Sampler) time calibration files. +// +// The basic script is taken from the October 2014 analysis framework running on +// La Palma, maintained by Daniela Dorner and Thomas Bretz. There it was called +// "callisto_drs_time.C" + +// Dominik and Sebastian: +// Our first approach was to use the MWriteFitsFile class but we have not been +// able to use it. + +// MWriteFitsFile fitsfile(outfile, +// MWriteFitsFile::kSingleFile, +// "RECREATE", +// "fitsfile", +// "atitle"); +// fitsfile.SetBytesPerSample("Data", 2); +// fitsfile.AddContainer("MDrsCalibrationTime","DANDELION"); + +#include "MLogManip.h" + +int produce_drs_time_fits_file( + const char *drs_time_file, //* path to raw drs-time-calib-file + const char *drs_file, //* path to drs.fits(.gz) file + const char *outfile //* output path like yyyymmdd_rrr.drs.time.fits +){ + gLog.Separator("produce_drs_time_fits_file"); + gLog << all; + gLog << "DRS Timing " << drs_time_file << '\n'; + gLog << "DRS 1024 " << drs_file << '\n'; + gLog << endl; + + MParList plist0; + + MTaskList tlist0; + plist0.AddToList(&tlist0); + + MDrsCalibration drscalib1024; + if (!drscalib1024.ReadFits(drs_file)){ + gLog << "Error while opening drs amplitude calibration file:" + << drs_file << "Aborting" << endl; + return 52; // DN: Why 52 I don't know. + } + plist0.AddToList(&drscalib1024); + + MDrsCalibrationTime timecam; + plist0.AddToList(&timecam); + + MEvtLoop loop0("DetermineTimeCal"); + loop0.SetParList(&plist0); + + factfits drstimeRawFile(drs_time_file); + const int NumberOfChips = drstimeRawFile.GetInt("NPIX")/9; // --> 160 + const int NumberOfCells = drstimeRawFile.GetInt("NROI"); // --> 1024 + drstimeRawFile.close(); + + MRawFitsRead read0(drs_time_file); + tlist0.AddToList(&read0); + MGeomApply apply0; + tlist0.AddToList(&apply0); + MDrsCalibApply drsapply0; + tlist0.AddToList(&drsapply0); + MContinue cont0("MRawEvtHeader.GetTriggerID!=33792", "SelectTim"); + tlist0.AddToList(&cont0); + MFillH fill0("MHDrsCalibrationTime"); + fill0.SetNameTab("DeltaT"); + tlist0.AddToList(&fill0); + + if (!loop0.Eventloop(0)){ + gLog << "Error performing the loop over the drs-time-run:" + << drs_time_file << "Aborting" << endl; + return 8; // DN: why 8 I don't know. + } + + ofits drstimeFile(outfile); + drstimeFile.SetDefaultKeys(); + drstimeFile.AddColumnDouble( + NumberOfChips*NumberOfCells, + "SamplingTimeDeviation", "nominal slices", + "see following COMMENT." + ); + drstimeFile.SetInt("NCHIPS", NumberOfChips, "number of drs chips"); + drstimeFile.SetInt("NCELLS", NumberOfCells, "number of cells of each chip"); + drstimeFile.AddComment( +"The SamplingTimeDeviation specifies the deviation of the actual "); + drstimeFile.AddComment( +"sampling time of a certain cell, from its nominal sampling time."); + drstimeFile.AddComment( +"It is measured in nominal slices."); + drstimeFile.AddComment( +"In order to convert the sample index into the actual sampling time"); + drstimeFile.AddComment( +"just do: "); + drstimeFile.AddComment( +" cell_index + td(cell_index) - td(stop_cell_index)"); + drstimeFile.AddComment( +"where td(i) refers to the SamplingTimeDeviation of cell i of the "); + drstimeFile.AddComment( +"current DRS4 chip."); + drstimeFile.AddComment( +"This gives you the actual sampling time of the cell measured in "); + drstimeFile.AddComment( +"nominal slices [nsl] with respect to the stop cells sampling time."); + drstimeFile.AddComment( +"Convert it to ns, simply by multiplication "); + drstimeFile.AddComment( +"with 0.5 ns/nsl."); + drstimeFile.AddComment( +""); + + drstimeFile.AddComment( +"In FACT it became common practice to measure un-time-calibrated "); + drstimeFile.AddComment( +" DRS4 data in 'time slices' or simply 'samples', while "); + drstimeFile.AddComment( +" time-calibrated DRS4 data usually measured in ns. Since the method,"); + drstimeFile.AddComment( +" that is emplyed to measure the DRS4 time calibration constants has "); + drstimeFile.AddComment( +" no possibility to assert the actual length of a slice is really half"); + drstimeFile.AddComment( +" a nanosecond, this practice is not advisable."); + drstimeFile.AddComment( +" I propose instead to call un-time-calibrated data 'samples',"); + drstimeFile.AddComment( +" since this is what they are. If one wants to stress the fact,"); + drstimeFile.AddComment( +" that no drs time calibration has been applied one should refer to "); + drstimeFile.AddComment( +"'uncalibrated slices'. Since it *should* be common practice to apply"); + drstimeFile.AddComment( +" the drs4 time calibration in all cases, I also propose to use the"); + drstimeFile.AddComment( +" short term of 'slices' or 'sl'. If one wants to stress, that the"); + drstimeFile.AddComment( +" drs4 time calibration has actually been applied to the data,"); + drstimeFile.AddComment( +" the term 'calibrated slices' or 'nominal slices' or short 'nsl'."); + + drstimeFile.WriteTableHeader("DRS_CELL_TIMES"); + + // By drs_sampling_time_deviations we refer to the deviations (measured in nominal slices) of + // the actual sampling time compared to the nominal sampling time of every cell of every chip. + double *drs_sampling_time_deviations = new double[ NumberOfChips * NumberOfCells ]; + for (int chip = 0; chip < NumberOfChips; chip++){ + for (int cell = 0; cell < NumberOfCells; cell++){ + // Dominik and Sebastian: + // We ended with using DrsCalibrateTime.Sum() in order to retrieve + // the contents DrsCalibrateTime.fStat. + // First we wanted to access the member fStat of class DrsCalibrateTime + // which is declared public but it did not work out of the box. + // drs_sampling_time_deviations[chip*NumberOfCells+cell] = + // timecam.fStat[chip*NumberOfCells+cell].first; + drs_sampling_time_deviations[ chip * NumberOfCells + cell ] = + timecam.Sum( chip * NumberOfCells + cell ); + } + } + + drstimeFile.WriteRow(drs_sampling_time_deviations, sizeof(double) * NumberOfChips * NumberOfCells); + drstimeFile.close(); + + return 0; +} Index: branches/MarsGapdTimeJitter/mimage/MImgCleanTime.cc =================================================================== --- branches/MarsGapdTimeJitter/mimage/MImgCleanTime.cc (revision 18306) +++ branches/MarsGapdTimeJitter/mimage/MImgCleanTime.cc (revision 18331) @@ -99,4 +99,7 @@ const Float_t time1 = pix1.GetArrivalTime(); + if (pix1.IsPixelUnmapped()) + return Island(0,time1); + Island island(pix1.GetNumPhotons(), time1); @@ -106,4 +109,6 @@ MSignalPix &pix2 = (*fEvt)[idx2]; + if (pix2.IsPixelUnmapped()) + continue; const Int_t island2 = pix2.GetIdxIsland(); @@ -137,10 +142,11 @@ Int_t MImgCleanTime::Process() { - // The assumption is that all pixels contain valid data + // The assumption is that all mapped pixels contain valid data const UInt_t npix = fEvt->GetNumPixels(); for (UInt_t i=0; i