Index: branches/Mars_McMismatchStudy/fact/analysis/callisto.C =================================================================== --- branches/Mars_McMismatchStudy/fact/analysis/callisto.C (revision 17961) +++ branches/Mars_McMismatchStudy/fact/analysis/callisto.C (revision 17962) @@ -1,5 +1,5 @@ #include "MLogManip.h" -int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output") +int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output", bool use_delays=false) { // ====================================================== @@ -12,182 +12,8 @@ // map file to use (get that from La Palma!) - const char *map = usemap ? "/home/fact/FACT++/FACTmap111030.txt" : NULL; -/* - Bool_t maximum = kTRUE; - - const char *lp_template = maximum ? - "template-lp-extractor-maximum.root" : - "template-lp-extractor-leading-edge.root"; - - const char *pulse_template = "template-pulse.root"; -*/ +// const char *map = usemap ? "/home/fact/FACT++/FACTmap111030.txt" : NULL; + const char *map = usemap ? "/scratch/fact/FACTmap111030.txt" : NULL; + // ------------------------------------------------------ - - //bool use_delays=false; - - int spike_removal=3; - - // The gain as extracted from our dark count spectra - double gain = 258; - - // ------------------------------------------------------ - - // 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 - - // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel: - // - // IntegralFixed: - // * fRiseTime: Number of slices left from arrival time - // * fFallTime: Number of slices right from arrival time - // IntegralRel: - // * fRiseTime: Number of slices left from maximum time - // * fFallTime: Number of slices right from maximum time - // -/* - const int rise_time_cal = maximum ? 40 : 10; // was 13; 5ns - const int fall_time_cal = maximum ? 120 : 160; // was 23; 80ns - - const int rise_time_dat = maximum ? 10 : 2; // was 13; was 10; 1ns - const int fall_time_dat = maximum ? 40 : 48; // was 23; was 40; 24ns - - // Extraction type: Extract integral and half leading edge - - const int type = maximum ? (MExtralgoSpline::kAmplitudeRel) : (MExtralgoSpline::kIntegralFixed); - //const int type = MExtralgoSpline::kIntegralFixed; - - - 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 - - // ====================================================== - - 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; - } - - // -------------------------------------------------------------------------------- - - gLog.Separator("Callisto"); - gLog << all << "Calibrate data of sequence '" << seq.GetFileName() << "'" << endl; - gLog << endl; - - // ------------------------------------------------------ - - ostringstream drsname; - drsname << gSystem->DirName(seqfile) << "/"; - drsname << seq.GetNight().GetNightAsInt() << "_"; - drsname << Form("%03d", seq.GetDrsSequence()) << ".drs.seq"; - - MSequence drs(drsname.str().c_str()); - if (!drs.IsValid()) - { - gLog << err << "ERROR - DRS sequence invalid!" << endl; - return 3; - } - - gLog << all << "DRS sequence file: " << drsname.str() << '\n' << endl; - - TString drsfile = seq.GetFileName(0, MSequence::kRawDrs); - if (drsfile.IsNull()) - { - gLog << err << "No DRS file available in sequence." << endl; - return 4; - } - - TString timfile = drs.GetFileName(0, MSequence::kFitsDat); - TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs); - TString calfile = seq.GetFileName(0, MSequence::kFitsCal); - //TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); - - gLog << all; - gLog << "DRS calib 300: " << drsfile << '\n'; - gLog << "DRS calib 1024: " << drs1024 << "\n\n"; - - MDrsCalibration drscalib300; - if (!drscalib300.ReadFits(drsfile.Data())) - return 5; - - MDrsCalibration drscalib1024; - if (!drscalib1024.ReadFits(drs1024.Data())) - return 6; - - gLog << all; - gLog << "Time calibration : " << timfile << '\n'; - gLog << "Light Pulser file: " << calfile << '\n' << endl; - //gLog << "Pedestal file: " << pedfile << '\n'; - - // ------------------------------------------------------ - - MDirIter iter; - if (seq.GetRuns(iter, MSequence::kFitsDat)<=0) - { - gLog << err << "ERROR - Sequence valid but without files." << endl; - return 7; - } - iter.Print(); - - // ====================================================== - - MStatusArray arrt, arrp; - - TFile ft(lp_template); - if (arrt.Read()<=0) - { - gLog << err << "ERROR - Reading LP template from " << lp_template << endl; - return 100; - } - - MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam"); - if (!lpref) - { - gLog << err << "ERROR - LP Template not found in " << lp_template << endl; - return 101; - } - lpref->SetDirectory(0); - - MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain"); - if (!gain) - { - gLog << err << "ERROR - Gain not found in " << lp_template << endl; - return 101; - } - gain->SetDirectory(0); - - TFile fp(pulse_template); - if (arrp.Read()<=0) - { - gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl; - return 102; - } - - TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0"); - if (!hpulse) - { - gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl; - return 103; - } - hpulse->SetDirectory(0); - - // ====================================================== MStatusDisplay *d = new MStatusDisplay; @@ -201,5 +27,4 @@ badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable); - // Twin pixel // 113 @@ -210,5 +35,152 @@ // 1393 - MDrsCalibrationTime timecam; + // ------------------------------------------------------ + + // ------------------------------------------------------ + + // 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; + } + + // -------------------------------------------------------------------------------- + + gLog.Separator("Callisto"); + gLog << all << "Calibrate data of sequence '" << seq.GetFileName() << "'" << endl; + gLog << endl; + + // ------------------------------------------------------ + + ostringstream drsname; + drsname << gSystem->DirName(seqfile) << "/"; + drsname << seq.GetNight().GetNightAsInt() << "_"; + drsname << Form("%03d", seq.GetDrsSequence()) << ".drs.seq"; + + MSequence drs(drsname.str().c_str()); + if (!drs.IsValid()) + { + gLog << err << "ERROR - DRS sequence invalid!" << endl; + return 3; + } + + gLog << all << "DRS sequence file: " << drsname.str() << '\n' << endl; + + TString drsfile = seq.GetFileName(0, MSequence::kRawDrs); + if (drsfile.IsNull()) + { + gLog << err << "No DRS file available in sequence." << endl; + return 4; + } + + TString timfile = drs.GetFileName(0, MSequence::kFitsDat); + TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs); + TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); + TString calfile = seq.GetFileName(0, MSequence::kFitsCal); + + gLog << all; + gLog << "DRS calib 300: " << drsfile << '\n'; + gLog << "DRS calib 1024: " << drs1024 << "\n\n"; + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile.Data())) + return 5; + + MDrsCalibration drscalib1024; + if (!drscalib1024.ReadFits(drs1024.Data())) + return 6; + + gLog << all; + gLog << "Time calibration : " << timfile << '\n'; + gLog << "Pedestal file: " << pedfile << '\n'; + gLog << "Light Pulser file: " << calfile << '\n' << endl; + + // ------------------------------------------------------ + + MDirIter iter; + if (seq.GetRuns(iter, MSequence::kFitsDat)<=0) + { + gLog << err << "ERROR - Sequence valid but without files." << endl; + return 7; + } + iter.Print(); + + // ====================================================== + +/* + MStatusArray arrt, arrp; + + TFile ft(lp_template); + if (arrt.Read()<=0) + { + gLog << err << "ERROR - Reading LP template from " << lp_template << endl; + return 100; + } + + MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam"); + if (!lpref) + { + gLog << err << "ERROR - LP Template not found in " << lp_template << endl; + return 101; + } + lpref->SetDirectory(0); + + MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain"); + if (!gain) + { + gLog << err << "ERROR - Gain not found in " << lp_template << endl; + return 101; + } + gain->SetDirectory(0); + + TFile fp(pulse_template); + if (arrp.Read()<=0) + { + gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl; + return 102; + } + + TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0"); + if (!hpulse) + { + gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl; + return 103; + } + hpulse->SetDirectory(0); +*/ + // ====================================================== // Plot the trigger pattern rates vs. run-number @@ -223,9 +195,5 @@ // hrate.DefaultLabelY("ERROR"); - Bool_t isinteg = - (type&MExtralgoSpline::kIntegral) || - (type&MExtralgoSpline::kFixedWidth) || - (type&MExtralgoSpline::kDynWidth) - ? kTRUE : kFALSE; + MDrsCalibrationTime timecam; gStyle->SetOptFit(kTRUE); @@ -241,5 +209,4 @@ plist0.AddToList(&tlist0); plist0.AddToList(&drscalib1024); - plist0.AddToList(&badpixels); plist0.AddToList(&timecam); @@ -257,5 +224,5 @@ MDrsCalibApply drsapply0; - drsapply0.SetRemoveSpikes(spike_removal); + //drsapply0.SetRemoveSpikes(4); MFillH fill0("MHDrsCalibrationTime"); @@ -281,5 +248,5 @@ // ====================================================== -/* + gLog << endl; gLog.Separator("Processing external light pulser run"); @@ -308,35 +275,13 @@ MDrsCalibApply drsapply1; - drsapply1.SetRemoveSpikes(spike_removal); + //drsapply1.SetRemoveSpikes(4); + + MFilterData filterdata1; // --- - MExtractTimeAndChargeSpline extractor1b("ExtractPulse"); + MExtractFACT extractor1b("ExtractPulse"); extractor1b.SetRange(first_slice, last_slice); - extractor1b.SetRiseTimeHiGain(rise_time_cal); - extractor1b.SetFallTimeHiGain(fall_time_cal); - extractor1b.SetHeightTm(heighttm); - extractor1b.SetChargeType(type); - extractor1b.SetSaturationLimit(600000); extractor1b.SetNoiseCalculation(kFALSE); - - MExtractTimeAndChargeSpline extractor1c("ExtractAmplitude"); - extractor1c.SetRange(first_slice, last_slice); - extractor1c.SetChargeType(MExtralgoSpline::kAmplitude); - extractor1c.SetSaturationLimit(600000); - extractor1c.SetNoiseCalculation(kFALSE); - extractor1c.SetNameSignalCam("Amplitude"); - extractor1c.SetNameTimeCam("AmplitudePos"); - - // --- - - MHCamEvent evt1a(5, "CalRatio", "Ratio per slice between integrated signal and amplitude;; r [1/n]"); - evt1a.SetNameSub("Amplitude", kTRUE); - MFillH fill1a(&evt1a, "MExtractedSignalCam", "FillRatio"); - fill1a.SetDrawOption("gaus"); - - MParameterD ratio1a; - ratio1a.SetVal(1./(fall_time_cal+rise_time_cal)); - fill1a.SetWeight(&ratio1a); // --- @@ -391,10 +336,6 @@ tlist1.AddToList(&drsapply1); tlist1.AddToList(&cont1); + tlist1.AddToList(&filterdata1); tlist1.AddToList(&extractor1b); - if (isinteg) - { - tlist1.AddToList(&extractor1c); - tlist1.AddToList(&fill1a); - } tlist1.AddToList(&calctm1a); tlist1.AddToList(&treat1); @@ -416,123 +357,6 @@ timecam.SetDelays(*evt1h.GetHist()); - // ========================= Result ================================== - - Double_t avgS = evt1f.GetHist()->GetMean(); - Double_t medS = evt1f.GetHist()->GetMedian(); - Double_t rmsS = evt1f.GetHist()->GetRMS(); - Double_t maxS = evt1f.GetHist()->GetMaximum(); - - MArrayF der1(hpulse->GetNbinsX()); - MArrayF der2(hpulse->GetNbinsX()); - - MExtralgoSpline spline(hpulse->GetArray()+1, hpulse->GetNbinsX(), - der1.GetArray(), der2.GetArray()); - spline.SetRiseFallTime(rise_time_dat, fall_time_dat); - spline.SetExtractionType(type); - spline.SetHeightTm(heighttm); - - spline.Extract(hpulse->GetMaximumBin()-1); - - // The pulser signal is most probably around 400mV/9.5mV - // IntegraFixed 2/48 corresponds to roughly 215mV*50slices - Double_t scale = 1./spline.GetSignal(); - - MArrayD calib(1440); - for (int i=0; i<1440; i++) - { - Double_t g = gain->GetBinContent(i+1)>0.5 ? gain->GetBinContent(i+1) : 1; - if (evt1f.GetHist()->GetBinContent(i+1)>0 && !badpixels[i].IsUnsuitable()) - calib[i] = lpref->GetBinContent(i+1) / evt1f.GetHist()->GetBinContent(i+1) / g; - } - - gROOT->SetSelectedPad(0); - d->AddTab("PulseTemp"); - gPad->SetGrid(); - hpulse->SetNameTitle("Pulse", "Single p.e. pulse template"); - hpulse->SetDirectory(0); - hpulse->SetLineColor(kBlack); - hpulse->DrawCopy(); - - TAxis *ax = hpulse->GetXaxis(); - - Double_t w = hpulse->GetBinWidth(1); - Double_t T = w*(spline.GetTime()+0.5) +ax->GetXmin(); - Double_t H = w*(hpulse->GetMaximumBin()+0.5)+ax->GetXmin(); - - TLine line; - line.SetLineColor(kRed); - line.DrawLine(T-rise_time_dat*w, spline.GetHeight(), - T+fall_time_dat*w, spline.GetHeight()); - line.DrawLine(T, spline.GetHeight()/4, T, 3*spline.GetHeight()/4); - line.DrawLine(T-rise_time_dat*w, 0, - T-rise_time_dat*w, spline.GetHeight()); - line.DrawLine(T+fall_time_dat*w, 0, - T+fall_time_dat*w, spline.GetHeight()); - - TGraph gg; - for (int ix=1; ix<=hpulse->GetNbinsX(); ix++) - for (int i=0; i<10; i++) - { - Double_t x = hpulse->GetBinLowEdge(ix)+i*hpulse->GetBinWidth(ix)/10.; - gg.SetPoint(gg.GetN(), x+w/2, spline.EvalAt(ix-1+i/10.)); - } - - gg.SetLineColor(kBlue); - gg.SetMarkerColor(kBlue); - gg.SetMarkerStyle(kFullDotMedium); - gg.DrawClone("L"); - - gROOT->SetSelectedPad(0); - d->AddTab("CalConst"); - MGeomCamFACT fact; - MHCamera hcalco(fact); - hcalco.SetName("CalConst"); - hcalco.SetTitle(Form("Relative calibration constant [%.0f/pe]", 1./scale)); - hcalco.SetCamContent(calib); - hcalco.SetAllUsed(); - //hcalco.Scale(scale); - hcalco.DrawCopy(); -*/ - - TF1 f("f", "[0]*(1-1/(1+exp(x/[1])))*exp(-x/[2])", -15, 150); - f.SetParameter(0, gain*0.05143); - f.SetParameter(1, 1.075); - f.SetParameter(2, 19.3); - f.SetNpx(2*165); //2GHz - - // Convert to graph - TGraph g(&f); - - // Convert to float - MArrayF x(g.GetN()); - MArrayF y(g.GetN()); - for (int i=0; iGetPathInfo(outfile, fstat); + int rc = gSystem->GetPathInfo(outpath, fstat); bool isdir = !rc || R_ISDIR(fstat.fMode); - const char *buf = gSystem->ConcatFileName(outfile, "callisto.root"); - outfile = buf; + TString filename = datfile + "_callisto.root"; + filename.Replace(0, filename.Last('/')+1, ""); + const char *buf = gSystem->ConcatFileName(outpath, filename); + TString outfile = buf; delete [] buf; @@ -119,116 +121,7 @@ // map file to use (get that from La Palma!) - const char *pmap = usemap ? "TestForThomas/FACT/FACTmap111030.txt" : NULL; - - Bool_t maximum = kTRUE; - - //const char *lp_template = maximum ? - // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-maximum.root" : - // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-leading-edge.root"; - - const char *pulse_template = "TestForThomas/FACT/template-pulse.root"; + const char *pmap = usemap ? "/home/isdc/toscanos/FACT/Mars_svn/resources/FACTmap111030.txt" : NULL; // ------------------------------------------------------ - - // 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 = 20; // 10ns - const int last_slice = 250; // 125ns - - // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel: - // - // IntegralFixed: - // * fRiseTime: Number of slices left from arrival time - // * fFallTime: Number of slices right from arrival time - // IntegralRel: - // * fRiseTime: Number of slices left from maximum time - // * fFallTime: Number of slices right from maximum time - // - const int rise_time_cal = maximum ? 40 : 10; // was 13; 5ns - const int fall_time_cal = maximum ? 120 : 160; // was 23; 80ns - - const int rise_time_dat = maximum ? 10 : 2; // was 13; was 10; 1ns - const int fall_time_dat = maximum ? 40 : 48; // was 23; was 40; 24ns - - // Extraction type: Extract integral and half leading edge - - const MExtralgoSpline::ExtractionType_t type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed); - //const int type = MExtralgoSpline::kIntegralFixed; - - - const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 } - - Long_t max = 0; // All - Long_t max3 = max; // Pedestal Rndm - Long_t max4 = max; // Pedestal Ext - - // ====================================================== - - if (pmap && gSystem->AccessPathName(pmap, kFileExists)) - { - gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl; - return 1; - } - - gLog.Separator("Callisto"); - gLog << all; - gLog << "Data File: " << datfile << '\n'; - gLog << "DRS calib 300: " << drsfile << endl;; - - MDrsCalibration drscalib300; - if (!drscalib300.ReadFits(drsfile.Data())) { - gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl; - return 5; - } - gLog << all; - gLog << "Pedestal file: " << pedfile << '\n'; - gLog << "Output file: " << outfile << '\n'; - gLog << "Display file: " << displayfile << '\n'; - gLog << "Display title: " << displaytitle << endl; - - // ------------------------------------------------------ - MStatusArray arrt, arrp; - - // TFile ft(lp_template); - // if (arrt.Read()<=0) - // { - // gLog << err << "ERROR - Reading LP template from " << lp_template << endl; - // return 100; - // } - - // MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam"); - // if (!lpref) - // { - // gLog << err << "ERROR - LP Template not found in " << lp_template << endl; - // return 101; - // } - // lpref->SetDirectory(0); - - // MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain"); - // if (!gain) - // { - // gLog << err << "ERROR - Gain not found in " << lp_template << endl; - // return 101; - // } - // gain->SetDirectory(0); - - TFile fp(pulse_template); - if (arrp.Read()<=0) - { - gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl; - return 102; - } - - TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0"); - if (!hpulse) - { - gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl; - return 103; - } - hpulse->SetDirectory(0); - // ====================================================== MStatusDisplay *d = new MStatusDisplay; @@ -236,5 +129,4 @@ MBadPixelsCam badpixels; badpixels.InitSize(1440); - /* badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable); @@ -243,5 +135,4 @@ badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable); - */ // Twin pixel // 113 @@ -251,4 +142,91 @@ // 1195 // 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 + + Long_t max = 0; // All + Long_t max3 = max; // Pedestal Rndm + Long_t max4 = max; // Pedestal Ext + + // ====================================================== + + //double scale = 0.1; + double scale = 0.1024; + + // ====================================================== + + if (pmap && gSystem->AccessPathName(pmap, kFileExists)) + { + gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl; + return 1; + } + + gLog.Separator("Callisto"); + gLog << all; + gLog << "Data File: " << datfile << '\n'; + gLog << "DRS calib 300: " << drsfile << endl;; + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile.Data())) { + gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl; + return 5; + } + gLog << all; + gLog << "Pedestal file: " << pedfile << '\n'; + gLog << "Output file: " << outfile << '\n'; + gLog << "Display file: " << displayfile << '\n'; + gLog << "Display title: " << displaytitle << endl; + + // ------------------------------------------------------ +/* + MStatusArray arrt, arrp; + + TFile ft(lp_template); + if (arrt.Read()<=0) + { + gLog << err << "ERROR - Reading LP template from " << lp_template << endl; + return 100; + } + + MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam"); + if (!lpref) + { + gLog << err << "ERROR - LP Template not found in " << lp_template << endl; + return 101; + } + lpref->SetDirectory(0); + + MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain"); + if (!gain) + { + gLog << err << "ERROR - Gain not found in " << lp_template << endl; + return 101; + } + gain->SetDirectory(0); + + TFile fp(pulse_template); + if (arrp.Read()<=0) + { + gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl; + return 102; + } + + TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0"); + if (!hpulse) + { + gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl; + return 103; + } + hpulse->SetDirectory(0); + */ + // ====================================================== MDrsCalibrationTime timecam; @@ -268,5 +246,5 @@ // ========================= Result ================================== - +/* //~ Double_t avgS = evt1f.GetHist()->GetMean(); //~ Double_t medS = evt1f.GetHist()->GetMedian(); @@ -340,5 +318,5 @@ //hcalco.Scale(scale); hcalco.DrawCopy(); - +*/ // ====================================================== @@ -372,34 +350,15 @@ MDrsCalibApply drsapply3; + MFilterData filterdata3; + //--- - MExtractTimeAndChargeSpline extractor3; + MExtractFACT extractor3; extractor3.SetRange(first_slice, last_slice); - extractor3.SetRiseTimeHiGain(rise_time_dat); - extractor3.SetFallTimeHiGain(fall_time_dat); - extractor3.SetHeightTm(heighttm); - extractor3.SetChargeType(type); - extractor3.SetSaturationLimit(600000); extractor3.SetNoiseCalculation(kTRUE); - -// MHCamEvent evt2a(0, "PedRdm", "Extracted Pedestal Signal;;S"); - -// MFillH fill2a(&evt2a, "MExtractedSignalCam", "FillPedRndm"); - - // Use this for data, but not for calibration events -// evt2a.SetErrorSpread(kFALSE); - - /* - MCalibrateData conv3; - conv3.SetCalibrationMode(MCalibrateData::kNone); - conv3.SetPedestalFlag(MCalibrateData::kNo); - conv3.SetCalibConvMinLimit(0); - conv3.SetCalibConvMaxLimit(10000); - conv3.SetScaleFactor(scale); - */ MCalibrateFact conv3; conv3.SetScale(scale); - conv3.SetCalibConst(calib); + //conv3.SetCalibConst(calib); MBadPixelsTreat treat3; @@ -420,4 +379,5 @@ tlist3.AddToList(&drsapply3); tlist3.AddToList(&cont3); + tlist3.AddToList(&filterdata3); tlist3.AddToList(&extractor3); // tlist3.AddToList(&fill3a); @@ -461,30 +421,13 @@ MDrsCalibApply drsapply4; - MExtractTimeAndChargeSpline extractor4; + MFilterData filterdata4; + + MExtractFACT extractor4; extractor4.SetRange(first_slice, last_slice); - extractor4.SetRiseTimeHiGain(rise_time_dat); - extractor4.SetFallTimeHiGain(fall_time_dat); - extractor4.SetHeightTm(heighttm); - extractor4.SetChargeType(type); - extractor4.SetSaturationLimit(600000); extractor4.SetNoiseCalculation(kFALSE); - // MHCamEvent evt3a(0, "PedExt", "Extracted Pedestal Signal;;S"); - - // MFillH fill3a(&evt3a, "MExtractedSignalCam", "FillPedExt"); - - // Use this for data, but not for calibration events -// evt3a.SetErrorSpread(kFALSE); -/* - MCalibrateData conv4; - conv4.SetCalibrationMode(MCalibrateData::kNone); - conv4.SetPedestalFlag(MCalibrateData::kNo); - conv4.SetCalibConvMinLimit(0); - conv4.SetCalibConvMaxLimit(10000); - conv4.SetScaleFactor(scale); -*/ MCalibrateFact conv4; conv4.SetScale(scale); - conv4.SetCalibConst(calib); + //conv4.SetCalibConst(calib); MBadPixelsTreat treat4; @@ -504,6 +447,6 @@ tlist4.AddToList(&drsapply4); tlist4.AddToList(&cont4); + tlist4.AddToList(&filterdata4); tlist4.AddToList(&extractor4); -// tlist4.AddToList(&fill4a); tlist4.AddToList(&conv4); tlist4.AddToList(&treat4); @@ -550,4 +493,8 @@ MDrsCalibApply drsapply5; + MTreatSaturation treatsat5; + + MFilterData filterdata5; + MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0", "SelectDat"); MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal"); @@ -559,29 +506,14 @@ // --- - MExtractTimeAndChargeSpline extractor5dat; + MExtractFACT extractor5dat; extractor5dat.SetRange(first_slice, last_slice); - extractor5dat.SetRiseTimeHiGain(rise_time_dat); - extractor5dat.SetFallTimeHiGain(fall_time_dat); - extractor5dat.SetHeightTm(heighttm); - extractor5dat.SetChargeType(type); - extractor5dat.SetSaturationLimit(600000); extractor5dat.SetNoiseCalculation(kFALSE); - MExtractTimeAndChargeSpline extractor5cal; + MExtractFACT extractor5cal; extractor5cal.SetRange(first_slice, last_slice); - extractor5cal.SetRiseTimeHiGain(rise_time_cal); - extractor5cal.SetFallTimeHiGain(fall_time_cal); - extractor5cal.SetHeightTm(heighttm); - extractor5cal.SetChargeType(type); - extractor5cal.SetSaturationLimit(600000); extractor5cal.SetNoiseCalculation(kFALSE); - MExtractTimeAndChargeSpline extractor5tm("ExtractTM"); + MExtractFACT extractor5tm("ExtractTM"); extractor5tm.SetRange(last_slice, 294); - extractor5tm.SetRiseTimeHiGain(1); - extractor5tm.SetFallTimeHiGain(1); - extractor5tm.SetHeightTm(0.5); - extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel); - extractor5tm.SetSaturationLimit(600000); extractor5tm.SetNoiseCalculation(kFALSE); extractor5tm.SetNameSignalCam("TimeMarkerAmplitude"); @@ -593,15 +525,7 @@ // --- -/* - MCalibrateData conv5; - conv5.SetCalibrationMode(MCalibrateData::kNone); - conv5.SetPedestalFlag(MCalibrateData::kNo); - conv5.SetCalibConvMinLimit(0); - conv5.SetCalibConvMaxLimit(10000); - conv5.SetScaleFactor(scale); -*/ MCalibrateFact conv5; conv5.SetScale(scale); - conv5.SetCalibConst(calib); + //conv5.SetCalibConst(calib); MCalibrateDrsTimes calctm5; @@ -689,4 +613,6 @@ tlist5.AddToList(&filterped); tlist5.AddToList(&fill5a); + tlist5.AddToList(&treatsat5); + tlist5.AddToList(&filterdata5); tlist5.AddToList(&extractor5dat); tlist5.AddToList(&extractor5cal); Index: branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C.GOOD =================================================================== --- branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C.GOOD (revision 17962) +++ branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C.GOOD (revision 17962) @@ -0,0 +1,721 @@ +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#if !defined(__CINT__) || defined(__MAKECINT__) + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.h" + +#include "../mcore/DrsCalib.h" +#include "MDrsCalibration.h" +#include "MExtralgoSpline.h" +#include "MSequence.h" +#include "MStatusArray.h" +#include "MHCamera.h" +#include "MJob.h" +#include "MWriteRootFile.h" +#include "MHCamera.h" +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" +#include "MDirIter.h" +#include "MTaskList.h" +#include "MFDataPhrase.h" +#include "MArrayF.h" +#include "MBadPixelsTreat.h" +#include "MCalibrateDrsTimes.h" +#include "MHSectorVsTime.h" +#include "MHCamEvent.h" +#include "MExtractTimeAndChargeSpline.h" +#include "MFillH.h" +#include "MDrsCalibApply.h" +#include "MGeomApply.h" +#include "MContinue.h" +#include "MRawFitsRead.h" +#include "MReadMarsFile.h" +#include "MEvtLoop.h" +#include "MParList.h" +#include "MStatusDisplay.h" +#include "MDrsCalibrationTime.h" +#include "MH3.h" +#include "MGeomCamFACT.h" +#include "MCalibrateFact.h" +#include "MParameters.h" +#include "MWriteAsciiFile.h" + +#endif + +using namespace std; + +/* Maybe you wanna use this macro like this: + * + * 0.) ---- call root ---- + * root -b + * + * 1.) ---- compile the stuff ---- + * .L fact/analysis/callisto_buildable_no_sequence_file.C++ + * + * + * 2.) ---- you can call it then ---- + * Therefore you need to specify all the paths ... see below. + * + * When you wanna call the stuff directly from the bash make sure to + * escape the bracets and quotes correctly. + * + * your can do: + * root -b -q callisto_buildable_no_sequence_file.C++'("path1","path2",...)' + * or: + * root -b -q callisto_buildable_no_sequence_file.C++(\"path1\",\"$HOME\",...) + * using bash enviroment variables like $HOME is not possible in the upper variant. + */ + +int callisto(const TString drsfile="test300samples2.drs.fits.gz", + const TString pedfile="00000001.001_P_MonteCarlo000_Events.fits", + const TString datfile="00000003.387_D_MonteCarlo010_Events.fits", + TString outpath = "", + TString displayfile = "", TString displaytitle = "") +{ + + // ====================================================== + + if (displaytitle.IsNull()) + displaytitle = gSystem->BaseName(datfile); + + FileStat_t fstat; + int rc = gSystem->GetPathInfo(outpath, fstat); + bool isdir = !rc || R_ISDIR(fstat.fMode); + + TString filename = datfile + "_callisto.root"; + filename.Replace(0, filename.Last('/')+1, ""); + const char *buf = gSystem->ConcatFileName(outpath, filename); + TString outfile = buf; + delete [] buf; + + if (displayfile.IsNull()) + { + displayfile = outfile; + displayfile.Insert(displayfile.Last('.'), "-display"); + } + else + { + if (isdir && gSystem->DirName(displayfile)==TString(".")) + { + buf = gSystem->ConcatFileName(outfile, displayfile); + displayfile = buf; + delete [] buf; + } + } + + // ====================================================== + + // true: Display correctly mapped pixels in the camera displays + // but the value-vs-index plot is in software/spiral indices + // false: Display pixels in hardware/linear indices, + // but the order is the camera display is distorted + bool usemap = true; + + // map file to use (get that from La Palma!) + const char *pmap = usemap ? "/home/isdc/toscanos/FACT/Mars_svn/resources/FACTmap111030.txt" : NULL; + + Bool_t maximum = kTRUE; + + //const char *lp_template = maximum ? + // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-maximum.root" : + // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-leading-edge.root"; + + // const char *lp_template = "/home/isdc/toscanos/FACT/Mars_svn/resources/template-lp.root"; + + //const char *pulse_template = "/cm/shared/apps/fact/Mars_svn_LP/template-pulse.root"; + + const char *pulse_template = "/home/isdc/toscanos/FACT/Mars_svn/resources/template-pulse.root"; + + // ------------------------------------------------------ + + // 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 = 20; // 10ns + const int last_slice = 250; // 125ns + + // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel: + // + // IntegralFixed: + // * fRiseTime: Number of slices left from arrival time + // * fFallTime: Number of slices right from arrival time + // IntegralRel: + // * fRiseTime: Number of slices left from maximum time + // * fFallTime: Number of slices right from maximum time + // + const int rise_time_cal = maximum ? 40 : 10; // was 13; 5ns + const int fall_time_cal = maximum ? 120 : 160; // was 23; 80ns + + const int rise_time_dat = maximum ? 10 : 2; // was 13; was 10; 1ns + const int fall_time_dat = maximum ? 40 : 48; // was 23; was 40; 24ns + + // Extraction type: Extract integral and half leading edge + + const MExtralgoSpline::ExtractionType_t type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed); + //const int type = MExtralgoSpline::kIntegralFixed; + + + const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 } + + Long_t max = 0; // All + Long_t max3 = max; // Pedestal Rndm + Long_t max4 = max; // Pedestal Ext + + // ====================================================== + + if (pmap && gSystem->AccessPathName(pmap, kFileExists)) + { + gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl; + return 1; + } + + gLog.Separator("Callisto"); + gLog << all; + gLog << "Data File: " << datfile << '\n'; + gLog << "DRS calib 300: " << drsfile << endl;; + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile.Data())) { + gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl; + return 5; + } + gLog << all; + gLog << "Pedestal file: " << pedfile << '\n'; + gLog << "Output file: " << outfile << '\n'; + gLog << "Display file: " << displayfile << '\n'; + gLog << "Display title: " << displaytitle << endl; + + // ------------------------------------------------------ + MStatusArray arrt, arrp; + + // TFile ft(lp_template); + // if (arrt.Read()<=0) + // { + // gLog << err << "ERROR - Reading LP template from " << lp_template << endl; + // return 100; + // } + + // MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam"); + // if (!lpref) + // { + // gLog << err << "ERROR - LP Template not found in " << lp_template << endl; + // return 101; + // } + // lpref->SetDirectory(0); + + // MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain"); + // if (!gain) + // { + // gLog << err << "ERROR - Gain not found in " << lp_template << endl; + // return 101; + // } + // gain->SetDirectory(0); + + TFile fp(pulse_template); + if (arrp.Read()<=0) + { + gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl; + return 102; + } + + TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0"); + if (!hpulse) + { + gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl; + return 103; + } + hpulse->SetDirectory(0); + // ====================================================== + + MStatusDisplay *d = new MStatusDisplay; + + MBadPixelsCam badpixels; + badpixels.InitSize(1440); + /* + badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable); + */ + // Twin pixel + // 113 + // 115 + // 354 + // 423 + // 1195 + // 1393 + + MDrsCalibrationTime timecam; + + // Plot the trigger pattern rates vs. run-number + MH3 hrate("MRawRunHeader.GetFileID", "MRawEvtHeader.GetTriggerID&0xff00"); + hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)"); + hrate.SetName("Rate"); + hrate.SetTitle("Event rate [Hz];File Id;Trigger Type;"); + hrate.InitLabels(MH3::kLabelsXY); + hrate.DefineLabelY( 0, "Data"); // What if TriggerID==0 already??? + hrate.DefineLabelY(0x100, "Cal"); + hrate.DefineLabelY(0x400, "Ped"); + // hrate.DefaultLabelY("ERROR"); + gStyle->SetOptFit(kTRUE); + + + // ========================= Result ================================== + + //~ Double_t avgS = evt1f.GetHist()->GetMean(); + //~ Double_t medS = evt1f.GetHist()->GetMedian(); + //~ Double_t rmsS = evt1f.GetHist()->GetRMS(); + //~ Double_t maxS = evt1f.GetHist()->GetMaximum(); + + MArrayF der1(hpulse->GetNbinsX()); + MArrayF der2(hpulse->GetNbinsX()); + + MExtralgoSpline spline(hpulse->GetArray()+1, hpulse->GetNbinsX(), + der1.GetArray(), der2.GetArray()); + spline.SetRiseFallTime(rise_time_dat, fall_time_dat); + spline.SetExtractionType(type); + spline.SetHeightTm(heighttm); + + spline.Extract(hpulse->GetMaximumBin()-1); + + // The pulser signal is most probably around 400mV/9.5mV + // IntegraFixed 2/48 corresponds to roughly 215mV*50slices + Double_t scale = 1./spline.GetSignal(); + + MArrayD calib(1440); + for (int i=0; i<1440; i++) + calib[i] =1.; + + gROOT->SetSelectedPad(0); + d->AddTab("PulseTemp"); + gPad->SetGrid(); + hpulse->SetNameTitle("Pulse", "Single p.e. pulse template"); + hpulse->SetDirectory(0); + hpulse->SetLineColor(kBlack); + hpulse->DrawCopy(); + + TAxis *ax = hpulse->GetXaxis(); + + Double_t w = hpulse->GetBinWidth(1); + Double_t T = w*(spline.GetTime()+0.5) +ax->GetXmin(); + //~ Double_t H = w*(hpulse->GetMaximumBin()+0.5)+ax->GetXmin(); + + TLine line; + line.SetLineColor(kRed); + line.DrawLine(T-rise_time_dat*w, spline.GetHeight(), + T+fall_time_dat*w, spline.GetHeight()); + line.DrawLine(T, spline.GetHeight()/4, T, 3*spline.GetHeight()/4); + line.DrawLine(T-rise_time_dat*w, 0, + T-rise_time_dat*w, spline.GetHeight()); + line.DrawLine(T+fall_time_dat*w, 0, + T+fall_time_dat*w, spline.GetHeight()); + + TGraph gg; + for (int ix=1; ix<=hpulse->GetNbinsX(); ix++) + for (int i=0; i<10; i++) + { + Double_t x = hpulse->GetBinLowEdge(ix)+i*hpulse->GetBinWidth(ix)/10.; + gg.SetPoint(gg.GetN(), x+w/2, spline.EvalAt(ix-1+i/10.)); + } + + gg.SetLineColor(kBlue); + gg.SetMarkerColor(kBlue); + gg.SetMarkerStyle(kFullDotMedium); + gg.DrawClone("L"); + + gROOT->SetSelectedPad(0); + d->AddTab("CalConst"); + MGeomCamFACT fact; + MHCamera hcalco(fact); + hcalco.SetName("CalConst"); + hcalco.SetTitle(Form("Relative calibration constant [%.0f/pe]", 1./scale)); + hcalco.SetCamContent(calib); + hcalco.SetAllUsed(); + //hcalco.Scale(scale); + hcalco.DrawCopy(); + + // ====================================================== + + gLog << endl; + gLog.Separator("Extracting random pedestal"); + + MTaskList tlist3; + + MParList plist3; + plist3.AddToList(&tlist3); + plist3.AddToList(&drscalib300); + plist3.AddToList(&badpixels); + plist3.AddToList(&timecam); + + MEvtLoop loop3("DetermineRndmPed"); + loop3.SetDisplay(d); + loop3.SetParList(&plist3); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read3; + read3.LoadMap(pmap); + read3.AddFile(pedfile); + + MFillH fill3a(&hrate); + + MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); + + MGeomApply apply3; + + MDrsCalibApply drsapply3; + + //MFilterData filterdata5; + + //--- + + MExtractTimeAndChargeSpline extractor3; + extractor3.SetRange(first_slice, last_slice); + extractor3.SetRiseTimeHiGain(rise_time_dat); + extractor3.SetFallTimeHiGain(fall_time_dat); + extractor3.SetHeightTm(heighttm); + extractor3.SetChargeType(type); + extractor3.SetSaturationLimit(600000); + extractor3.SetNoiseCalculation(kTRUE); + +// MHCamEvent evt2a(0, "PedRdm", "Extracted Pedestal Signal;;S"); + +// MFillH fill2a(&evt2a, "MExtractedSignalCam", "FillPedRndm"); + + // Use this for data, but not for calibration events +// evt2a.SetErrorSpread(kFALSE); + + /* + MCalibrateData conv3; + conv3.SetCalibrationMode(MCalibrateData::kNone); + conv3.SetPedestalFlag(MCalibrateData::kNo); + conv3.SetCalibConvMinLimit(0); + conv3.SetCalibConvMaxLimit(10000); + conv3.SetScaleFactor(scale); + */ + + MCalibrateFact conv3; + conv3.SetScale(scale); + conv3.SetCalibConst(calib); + + MBadPixelsTreat treat3; + treat3.SetProcessPedestalRun(kFALSE); + treat3.SetProcessPedestalEvt(kFALSE); + treat3.SetProcessTimes(kFALSE); + + MHCamEvent evt3b(0, "PedRdm","Interpolated random pedestal;;Signal [~phe]"); + //evt2b.SetErrorSpread(kFALSE); + + MFillH fill3b(&evt3b, "MSignalCam", "FillPedRdm"); + fill3b.SetDrawOption("gaus"); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist3.AddToList(&read3); + tlist3.AddToList(&apply3); + tlist3.AddToList(&drsapply3); + tlist3.AddToList(&cont3); + tlist3.AddToList(&extractor3); +// tlist3.AddToList(&fill3a); + tlist3.AddToList(&conv3); + tlist3.AddToList(&treat3); + tlist3.AddToList(&fill3b); + + if (!loop3.Eventloop(max3)) + return 14; + + if (!loop3.GetDisplay()) + return 15; + + // ====================================================== + + gLog << endl; + gLog.Separator("Extracting pedestal"); + + MTaskList tlist4; + + MParList plist4; + plist4.AddToList(&tlist4); + plist4.AddToList(&drscalib300); + plist4.AddToList(&badpixels); + plist4.AddToList(&timecam); + + MEvtLoop loop4("DetermineExtractedPed"); + loop4.SetDisplay(d); + loop4.SetParList(&plist4); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read4; + read4.LoadMap(pmap); + read4.AddFile(pedfile); + + MContinue cont4("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); + + MGeomApply apply4; + + MDrsCalibApply drsapply4; + + MExtractTimeAndChargeSpline extractor4; + extractor4.SetRange(first_slice, last_slice); + extractor4.SetRiseTimeHiGain(rise_time_dat); + extractor4.SetFallTimeHiGain(fall_time_dat); + extractor4.SetHeightTm(heighttm); + extractor4.SetChargeType(type); + extractor4.SetSaturationLimit(600000); + extractor4.SetNoiseCalculation(kFALSE); + + // MHCamEvent evt3a(0, "PedExt", "Extracted Pedestal Signal;;S"); + + // MFillH fill3a(&evt3a, "MExtractedSignalCam", "FillPedExt"); + + // Use this for data, but not for calibration events +// evt3a.SetErrorSpread(kFALSE); +/* + MCalibrateData conv4; + conv4.SetCalibrationMode(MCalibrateData::kNone); + conv4.SetPedestalFlag(MCalibrateData::kNo); + conv4.SetCalibConvMinLimit(0); + conv4.SetCalibConvMaxLimit(10000); + conv4.SetScaleFactor(scale); +*/ + MCalibrateFact conv4; + conv4.SetScale(scale); + conv4.SetCalibConst(calib); + + MBadPixelsTreat treat4; + treat4.SetProcessPedestalRun(kFALSE); + treat4.SetProcessPedestalEvt(kFALSE); + + MHCamEvent evt4b(0, "PedExt","Interpolated extracted pedestal;;Signal [~phe]"); + //evt4b.SetErrorSpread(kFALSE); + + MFillH fill4b(&evt4b, "MSignalCam", "FillPedExt"); + fill4b.SetDrawOption("gaus"); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist4.AddToList(&read4); + tlist4.AddToList(&apply4); + tlist4.AddToList(&drsapply4); + tlist4.AddToList(&cont4); + tlist4.AddToList(&extractor4); +// tlist4.AddToList(&fill4a); + tlist4.AddToList(&conv4); + tlist4.AddToList(&treat4); + tlist4.AddToList(&fill4b); + + if (!loop4.Eventloop(max4)) + return 15; + + if (!loop4.GetDisplay()) + return 16; + + // =================================================================== + + gLog << endl; + gLog.Separator("Extracting and calibration data"); + + MTaskList tlist5; + + MParList plist5; + plist5.AddToList(&tlist5); + plist5.AddToList(&drscalib300); + plist5.AddToList(&badpixels); + plist5.AddToList(&timecam); + + MEvtLoop loop5("CalibratingData"); + loop5.SetDisplay(d); + loop5.SetParList(&plist5); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read5a; + MReadMarsFile read5b("Events"); + read5a.LoadMap(pmap); + read5a.AddFile(datfile); + read5b.DisableAutoScheme(); + read5b.AddFile(datfile); + + MRead &read5 = datfile.EndsWith(".root") ? static_cast(read5b) : static_cast(read5a); + + MFillH fill5a(&hrate); + + MGeomApply apply5; + + MDrsCalibApply drsapply5; + + MTreatSaturation treatsat5; + + MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0", "SelectDat"); + MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal"); + MFDataPhrase filterped("(MRawEvtHeader.GetTriggerID&0xff00)==0x400", "SelectPed"); + MFDataPhrase filterncl("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectNonCal"); + + //MContinue cont4("MRawEvtHeader.GetTriggerID!=4", "SelectData"); + + // --- + + MExtractTimeAndChargeSpline extractor5dat; + extractor5dat.SetRange(first_slice, last_slice); + extractor5dat.SetRiseTimeHiGain(rise_time_dat); + extractor5dat.SetFallTimeHiGain(fall_time_dat); + extractor5dat.SetHeightTm(heighttm); + extractor5dat.SetChargeType(type); + extractor5dat.SetSaturationLimit(600000); + extractor5dat.SetNoiseCalculation(kFALSE); + + MExtractTimeAndChargeSpline extractor5cal; + extractor5cal.SetRange(first_slice, last_slice); + extractor5cal.SetRiseTimeHiGain(rise_time_cal); + extractor5cal.SetFallTimeHiGain(fall_time_cal); + extractor5cal.SetHeightTm(heighttm); + extractor5cal.SetChargeType(type); + extractor5cal.SetSaturationLimit(600000); + extractor5cal.SetNoiseCalculation(kFALSE); + + MExtractTimeAndChargeSpline extractor5tm("ExtractTM"); + extractor5tm.SetRange(last_slice, 294); + extractor5tm.SetRiseTimeHiGain(1); + extractor5tm.SetFallTimeHiGain(1); + extractor5tm.SetHeightTm(0.5); + extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel); + extractor5tm.SetSaturationLimit(600000); + extractor5tm.SetNoiseCalculation(kFALSE); + extractor5tm.SetNameSignalCam("TimeMarkerAmplitude"); + extractor5tm.SetNameTimeCam("TimeMarkerTime"); + + extractor5dat.SetFilter(&filterncl); + extractor5cal.SetFilter(&filtercal); + //extractor4tm.SetFilter(&filtercal); + + // --- +/* + MCalibrateData conv5; + conv5.SetCalibrationMode(MCalibrateData::kNone); + conv5.SetPedestalFlag(MCalibrateData::kNo); + conv5.SetCalibConvMinLimit(0); + conv5.SetCalibConvMaxLimit(10000); + conv5.SetScaleFactor(scale); +*/ + MCalibrateFact conv5; + conv5.SetScale(scale); + conv5.SetCalibConst(calib); + + MCalibrateDrsTimes calctm5; + calctm5.SetNameUncalibrated("UncalibratedTimes"); + + MCalibrateDrsTimes calctm5tm("CalibrateTimeMarker"); + calctm5tm.SetNameArrivalTime("TimeMarkerTime"); + calctm5tm.SetNameUncalibrated("UncalTimeMarker"); + calctm5tm.SetNameCalibrated("TimeMarker"); + calctm5tm.SetTimeMarker(); + //calctm4tm.SetFilter(&filtercal); + + MBadPixelsTreat treat5; + treat5.SetProcessPedestalRun(kFALSE); + treat5.SetProcessPedestalEvt(kFALSE); + + MHCamEvent evt5b(0, "ExtSig", "Extracted signal;;S [mV·sl]"); + MHCamEvent evt5c(0, "CalSig", "Calibrated and interpolated signal;;S [~phe]"); + MHCamEvent evt5d(4, "ExtSigTm", "Extracted time;;T [sl]"); + MHCamEvent evt5e(6, "CalSigTm", "Calibrated and interpolated time;;T [ns]"); + + MFillH fill5b(&evt5b, "MExtractedSignalCam", "FillExtSig"); + MFillH fill5c(&evt5c, "MSignalCam", "FillCalSig"); + MFillH fill5d(&evt5d, "MArrivalTimeCam", "FillExtTm"); + MFillH fill5e(&evt5e, "MSignalCam", "FillCalTm"); + + fill5c.SetDrawOption("gaus"); + fill5d.SetDrawOption("gaus"); + fill5e.SetDrawOption("gaus"); + + /* + fill4b.SetFilter(&filterdat); + fill4c.SetFilter(&filterdat); + fill4d.SetFilter(&filterdat); + fill4e.SetFilter(&filterdat); + */ + + //MFSoftwareTrigger swtrig; + //MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger"); + //contsw.SetInverted(); + + // The second rule is for the case reading raw-files! + + MWriteRootFile write5(outfile, "RECREATE", "Calibrated Data", 2); + write5.AddContainer("MRawRunHeader", "RunHeaders"); + write5.AddContainer("MGeomCam", "RunHeaders"); + write5.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); + write5.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE); + write5.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); + + // Common events + write5.AddContainer("MCorsikaEvtHeader", "Events", kFALSE); + write5.AddContainer("MMcEvt", "Events", kFALSE); + write5.AddContainer("IncidentAngle", "Events", kFALSE); + write5.AddContainer("MPointingPos", "Events", kFALSE); + write5.AddContainer("MSignalCam", "Events"); + write5.AddContainer("MTime", "Events", kFALSE); + write5.AddContainer("MRawEvtHeader", "Events"); + //write.AddContainer("MTriggerPattern", "Events"); + + // ------------------ Setup histograms and fill tasks ---------------- + + MContinue test; + test.SetFilter(&filterncl); + + MTaskList tlist5tm; + tlist5tm.AddToList(&extractor5tm); + tlist5tm.AddToList(&calctm5tm); + tlist5tm.SetFilter(&filtercal); + + MTaskList tlist5dat; + tlist5dat.AddToList(&fill5b); + tlist5dat.AddToList(&fill5c); + tlist5dat.AddToList(&fill5d); + tlist5dat.AddToList(&fill5e); + tlist5dat.SetFilter(&filterdat); + + tlist5.AddToList(&read5); + tlist5.AddToList(&apply5); + tlist5.AddToList(&drsapply5); + tlist5.AddToList(&treatsat5); + tlist5.AddToList(&filterncl); + //tlist5.AddToList(&test); + tlist5.AddToList(&filterdat); + tlist5.AddToList(&filtercal); + tlist5.AddToList(&filterped); + tlist5.AddToList(&fill5a); + tlist5.AddToList(&extractor5dat); + tlist5.AddToList(&extractor5cal); + tlist5.AddToList(&calctm5); + tlist5.AddToList(&tlist5tm); + tlist5.AddToList(&conv5); + tlist5.AddToList(&treat5); + tlist5.AddToList(&tlist5dat); + tlist5.AddToList(&write5); + + if (!loop5.Eventloop(max4)) + return 18; + + if (!loop5.GetDisplay()) + return 19; + + d->SetTitle(displaytitle, kFALSE); + d->SaveAs(displayfile); + + return 0; +} Index: branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C =================================================================== --- branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C (revision 17961) +++ branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C (revision 17962) @@ -169,4 +169,5 @@ MEnergyEstimate energyest; energyest.SetRule("MHillas.fSize^0.815*pow(10, 0.0199*(MHillasSrc.fDist*1.17e-2 - (-6.95))^2*(MHillasSrc.fDist*1.17e-2>(-6.95)) + 0.5997 *(MHillasSrc.fDist*1.17e-2 - 0.7568)^2 * (MHillasSrc.fDist*1.17e-2>0.7568) )"); + // energyest.SetRule("1737 * (MHillas.fSize / 300.)^0.810"); // Cuts before writing Index: branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C.ORIGINAL =================================================================== --- branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C.ORIGINAL (revision 17962) +++ branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C.ORIGINAL (revision 17962) @@ -0,0 +1,525 @@ +void SetupWriter(MWriteRootFile *write, const char *name, Bool_t ismc) const +{ + if (!write) + return; + + write->SetName(name); + write->AddContainer("MHillas", "Events"); + write->AddContainer("MHillasSrc", "Events"); + write->AddContainer("MHillasExt", "Events"); + write->AddContainer("MPointingPos", "Events"); + write->AddContainer("MHillasSrcAnti", "Events", kFALSE); + write->AddContainer("MImagePar", "Events", kFALSE); + write->AddContainer("MNewImagePar", "Events", kFALSE); + write->AddContainer("MNewImagePar2", "Events", kFALSE); + write->AddContainer("Hadronness", "Events", kFALSE); + write->AddContainer("MSrcPosCam", "Events", kFALSE); + write->AddContainer("MSrcPosAnti", "Events", kFALSE); + write->AddContainer("ThetaSquared", "Events", kFALSE); + write->AddContainer("OpticalAxis", "Events", kFALSE); + write->AddContainer("Disp", "Events", kFALSE); + write->AddContainer("Ghostbuster", "Events", kFALSE); + write->AddContainer("MEnergyEst", "Events", kFALSE); + write->AddContainer("MTime", "Events", kFALSE); + write->AddContainer("MMcEvt", "Events", kFALSE); + write->AddContainer("MEnergyEst", "Events", kFALSE); + write->AddContainer("DataType", "Events"); + write->AddContainer("Weight", "Events"); + write->AddContainer("FileId", "Events"); + write->AddContainer("EvtNumber", "Events"); + + write->AddContainer("MMcRunHeader", "RunHeaders", kFALSE); + write->AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); + write->AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE); + + //if (ismc) + // write->AddCopySource("OriginalMC"); +} + +int ganymed(const char *dataset, const char *outfile, Float_t ra, Float_t dec, Bool_t ismc=false) +{ + MDirIter files; + if (files.ReadList(dataset)<=0) + return 2; + files.Sort(); + + // ============================================================= + + // Crab (05h34'32" 22d00'52") + MPointingPos source("MSourcePos"); + //source.SetSkyPosition(5.5755, 22.0144); // ra[h], dec[deg] + //source.SetSkyPosition(12.356089, 30.176897); // ra[h], dec[deg]//1218 + //source.SetSkyPosition(16.897867,39.760201);// mrk501 + //source.SetSkyPosition(11.0742, 38.2089); // ra[h], dec[deg]//mrk421 + + if (ismc) + source.SetLocalPosition(ra, dec); // RA is ZD, DEC, is AZ!!!! + else + source.SetSkyPosition(ra, dec); // ra[h], dec[deg]//mrk421 + + UInt_t fNumOffSourcePos = 5; + + // ============================================================= + + // FIXME: If fPathIn read cuts and energy estimator from file! + MContinue contq("MImagePar.fNumIslands>3.5" + "|| MNewImagePar.fNumUsedPixels<5.5" + "|| MNewImagePar.fLeakage1>0.1", + "CutQ"); + + MContinue cont0("MHillas.fSize<0" + "|| MHillasExt.fSlopeSpreadWeighted>0.68" + "|| MHillasExt.fSlopeSpreadWeighted<0.18" + "|| log10(MHillas.GetArea*MGeomCam.fConvMm2Deg^2)>(log10(MHillas.fSize)-2)*1.1-1.55" + "|| MNewImagePar.fConcCore < 0.13" + "|| MNewImagePar.fConcCOG < 0.15" + , "Cut0"); + + // Cuts before plots + MContinue cont2("0" + "|| MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg<-0.07" + "|| log10(MNewImagePar.fConc1)<-log10(MHillas.fSize)*0.57+0.13" + //"|| log10(MNewImagePar.fConc1)>-log10(MHillas.fSize)*0.57+0.5" + //"|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>log10(MHillas.fSize)*0.42+0.29" + "|| MHillas.GetArea*MGeomCam.fConvMm2Deg^2<0.01" + "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>log10(MHillas.fSize)*0.5+0.095" + "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg<0.4" + "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>1.2", + "Cut2"); + + TArrayD param(12); + // Parametrization of Disp + param[0] = 1.15136; // constant + param[8] = 0.0681437; // slope + param[9] = 2.62932; // leak + param[10] = 1.51279; // size-offset + param[11] = 0.0507821; // size-slope + // Parametrization for sign of disp (m3long, slope) + param[5] = -0.07; + param[6] = 7.2; + param[7] = 0.5; + // ThetaSq-Cut + param[1] = 0.11; // 0.215 + // Area-Cut + param[2] = 0.215468; + param[3] = 5.63973; + param[4] = 0.0836169; + + // -------------------------------------------------------------------------------- + + // Setup fitter and histograms + MAlphaFitter fit; + fit.SetBackgroundFitMin(0.09); + fit.SetBackgroundFitMax(0.25); + fit.SetPolynomOrder(1); + fit.SetSignalFunction(MAlphaFitter::kThetaSq); + fit.SetScaleUser(1./fNumOffSourcePos); // includes fit.SetScaleMode(MAlphaFitter::kUserScale); + + // -------------------------------------------------------------------------------- + + MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin"); + MBinning binsD( 25, 0, 2.5, "BinningDist"); + MBinning binsC( 25, 0, 2.5, "BinningDistC"); + MBinning binsR( 60, 0, 300, "BinningRate"); + + // -------------------------------------------------------------------------------- + + MParList plist; + plist.AddToList(&source); + plist.AddToList(&fit); + + // Initialize default binnings + plist.AddToList(&bins3); + plist.AddToList(&binsD); + plist.AddToList(&binsC); + plist.AddToList(&binsR); + + MParameterI par("DataType"); + MParameterD weight("Weight"); + plist.AddToList(&par); + plist.AddToList(&weight); + + // Setup Tasklist + MTaskList tlist; + plist.AddToList(&tlist); + + // La Palma Magic1 + MObservatory obs; + plist.AddToList(&obs); + + // -------------------------------------------------------------------------------- + + MParameterD scale; + scale.SetVal(1./fNumOffSourcePos); + + MHThetaSq *halphaoff = plist.FindCreateObj("MHThetaSq", "HistOff"); + halphaoff->ForceUsingSize(); + + MFillH falpha(halphaoff, "", "FillHist"); + MFillH fenergy1("HistEOff [MHThetaSq]", "", "FillHistE"); + + MFMagicCuts cuts; + cuts.SetHadronnessCut(MFMagicCuts::kNoCut); + cuts.SetVariables(param); + + // Cuts before thetasq + MContinue cont1(&cuts, "Cut1"); + + // Energy estimation + MEnergyEstimate energyest; + energyest.SetRule("MHillas.fSize^0.815*pow(10, 0.0199*(MHillasSrc.fDist*1.17e-2 - (-6.95))^2*(MHillasSrc.fDist*1.17e-2>(-6.95)) + 0.5997 *(MHillasSrc.fDist*1.17e-2 - 0.7568)^2 * (MHillasSrc.fDist*1.17e-2>0.7568) )"); + + // Cuts before writing + MContinue cont3("", "Cut3"); + contq.SetAllowEmpty(); + cont0.SetAllowEmpty(); + cont1.SetAllowEmpty(); + cont2.SetAllowEmpty(); + cont3.SetAllowEmpty(); + + cont1.SetInverted(); + + // Filter for VsSize + MFDataPhrase ftheta("ThetaSquared.fVal < ThetaSquaredCut.fVal", "CutT"); + + // ------------- Loop Off Data -------------------- + MRead *readoffdata = NULL; + if (ismc) + { + MReadMarsFile *read = new MReadMarsFile("Events"); + read->DisableAutoScheme(); + readoffdata = read; + } + else + { + MReadReports *read = new MReadReports; + read->AddTree("Events", "MTime.", MReadReports::kMaster); + read->AddTree("Drive",MReadReports::kRequired); + read->AddTree("Rates",MReadReports::kRequired); + readoffdata = read; + } + + readoffdata->AddFiles(files); + + TString fname0 = Form("%s-summary.root", outfile); + TString fname1 = Form("%s-analysis.root", outfile); + + MWriteRootFile dummy0(fname0.Data(), "RECREATE"); + MWriteRootFile dummy1(fname1.Data(), "RECREATE"); + + MWriteRootFile *write0 = &dummy0; + MWriteRootFile *write1 = &dummy1; + SetupWriter(write0, "WriteAfterCut0", ismc); + SetupWriter(write1, "WriteAfterCut3", ismc); + + MWriteRootFile write2(fname1.Data(), "RECRATE"); + write2.AddContainer("MReportRates", "Rates"); + write2.AddContainer("MTimeRates", "Rates"); + + MParameterCalc setevtnum("MRawEvtHeader.fDAQEvtNumber", "SetEvtNumber"); + setevtnum.SetNameParameter("EvtNumber"); + + MParameterCalc setrunnum("MRawRunHeader.GetFileID", "SetFileId"); + setrunnum.SetNameParameter("FileId"); + + MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre"); + MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost"); + MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost"); + MFillH fill3c("MHVsSizeOffTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta"); + MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost"); + MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost"); + MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost"); + MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost"); + //MFillH fill9a("MHEffOffTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime"); + fill1a.SetNameTab("PreCut"); + fill2a.SetNameTab("PostCut"); + fill3a.SetNameTab("VsSize"); + fill3c.SetNameTab("CutT"); + fill4a.SetNameTab("HilExt"); + fill5a.SetNameTab("HilSrc"); + fill6a.SetNameTab("ImgPar"); + fill7a.SetNameTab("NewPar"); + //fill9a.SetNameTab("EffOffT"); + + fill3c.SetFilter(&ftheta); + + //MFDataMember fbin("Bin.fVal", '>', 0); + //fill9a.SetFilter(&fbin); + + MTaskList tlist2; + tlist2.SetNumPasses(fNumOffSourcePos); + fill2a.SetWeight(&scale); + fill3a.SetWeight(&scale); + fill3c.SetWeight(&scale); + fill4a.SetWeight(&scale); + fill5a.SetWeight(&scale); + fill6a.SetWeight(&scale); + fill7a.SetWeight(&scale); + + // How to get source position from off- and on-data? + MSrcPosCalc scalc; + scalc.SetMode(MSrcPosCalc::kWobble); + scalc.SetCallback(&tlist2); + + MSrcPosCorrect scor; + + MHillasCalc hcalc; + MHillasCalc hcalc2("MHillasCalcAnti"); + hcalc.SetFlags(MHillasCalc::kCalcHillasSrc); + hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc); + hcalc2.SetNameHillasSrc("MHillasSrcAnti"); + hcalc2.SetNameSrcPosCam("MSrcPosAnti"); + + MSrcPosRndm srcrndm; + + tlist2.AddToList(&scalc); + //tlist2.AddToList(&scor); + tlist2.AddToList(&srcrndm); + tlist2.AddToList(&hcalc); + tlist2.AddToList(&hcalc2); + tlist2.AddToList(&setrunnum); + tlist2.AddToList(&setevtnum); + tlist2.AddToList(&cuts); + tlist2.AddToList(&energyest); + tlist2.AddToList(write0); + tlist2.AddToList(&cont0); // Post/Pre shown cuts + tlist2.AddToList(&cont1); // Calculate ThetaSq + tlist2.AddToList(&fill2a); + tlist2.AddToList(&ftheta); // Calculate theta Cut + if (!ismc) + { + tlist2.AddToList(&falpha); + tlist2.AddToList(&fenergy1); + } + tlist2.AddToList(&fill3c); // With theta cut + tlist2.AddToList(&cont2); + tlist2.AddToList(&fill3a); + tlist2.AddToList(&fill4a); + tlist2.AddToList(&fill5a); + tlist2.AddToList(&fill6a); + tlist2.AddToList(&fill7a); + tlist2.AddToList(&cont3); + tlist2.AddToList(write1); + + MPointingPosCalc pcalc; + MPointingDevCalc devcalc; + + tlist.AddToList(readoffdata); + if (!ismc) + tlist.AddToList(&pcalc, "Drive"); + tlist.AddToList(&contq, "Events"); + tlist.AddToList(&fill1a, "Events"); + tlist.AddToList(&tlist2, "Events"); + + // by setting it here it is distributed to all consecutive tasks + tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); + + par.SetVal(0); + weight.SetVal(-1./fNumOffSourcePos); + + MStatusDisplay *d = new MStatusDisplay; + + // Create and setup the eventloop + MEvtLoop evtloop("Off"); + evtloop.SetParList(&plist); + evtloop.SetDisplay(d); + + // Execute first analysis + if (!evtloop.Eventloop()) + return 2; + + if (!evtloop.GetDisplay()) + return 2; + + // ------------- Loop On Data -------------------- + MRead *readondata = NULL; + if (ismc) + { + MReadMarsFile *read = new MReadMarsFile("Events"); + read->DisableAutoScheme(); + readondata = read; + } + else + { + MReadReports *read = new MReadReports; + read->AddTree("Events", "MTime.", MReadReports::kMaster); + read->AddTree("Drive", MReadReports::kRequired); + read->AddTree("Rates", MReadReports::kRequired); + read->AddTree("Humidity", MReadReports::kRequired); + read->AddTree("Temperatures", MReadReports::kRequired); + readondata = read; + } + readondata->AddFiles(files); + + scalc.SetMode(MSrcPosCalc::kDefault); + scalc.SetNumRandomOffPositions(fNumOffSourcePos); + + MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre"); + MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost"); + MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost"); + MFillH fill3d("MHVsSizeOnTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta"); + MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost"); + MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost"); + MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost"); + MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost"); + //MFillH fill9b("MHEffOnTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime"); + fill1b.SetNameTab("PreCut"); + fill2b.SetNameTab("PostCut"); + fill3b.SetNameTab("VsSize"); + fill3d.SetNameTab("CutT"); + fill4b.SetNameTab("HilExt"); + fill5b.SetNameTab("HilSrc"); + fill6b.SetNameTab("ImgPar"); + fill7b.SetNameTab("NewPar"); + //fill9b.SetNameTab("EffOnT"); + fill1b.SetDrawOption("same"); + fill2b.SetDrawOption("same"); + fill3b.SetDrawOption("same"); + fill3d.SetDrawOption("same"); + fill4b.SetDrawOption("same"); + fill5b.SetDrawOption("same"); + fill6b.SetDrawOption("same"); + fill7b.SetDrawOption("same"); + //fill9b.SetFilter(&fbin); + + fill3d.SetFilter(&ftheta); + + MHThetaSq *halphaon = plist.FindCreateObj("MHThetaSq", "Hist"); + halphaon->ForceUsingSize(); + + MFillH falpha2(halphaon, "", "FillHist"); + MFillH fenergy2("HistE [MHThetaSq]", "", "FillHistE"); + + tlist2.SetNumPasses(); + + tlist.Replace(readondata); + + MHSrcPosCam hsrcpos(kTRUE); + + MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam"); + fillsrc.SetNameTab("SrcPos"); + + tlist2.AddToListBefore(&fillsrc, &hcalc); + + MH3 hvs1("MPointingPos.fZd"); + MH3 hvs2("MPointingPos.fZd"); + MHVsTime hvstM("MTimeDrive.GetMjd"); + MHVsTime hvst0("MPointingPos.fZd"); + MHVsTime hvst1("MReportRates.fTriggerRate"); + MH3 hvst2("MPointingPos.fZd", "MReportRates.fTriggerRate"); + MHVsTime hvst("MReportTemperatures.GetTempMean"); + MHVsTime hvsh("MReportHumidity.GetMean"); + hvs1.SetName("ObsTime;Theta"); + hvs2.SetName("OnTime;Theta"); + hvstM.SetName("MJD"); + hvst0.SetName("ZdTime"); + hvst1.SetName("Rate"); + hvst2.SetName("RateZd;Theta;Rate"); + hvst.SetName("CamTemp"); + hvsh.SetName("CamHum"); + hvs1.SetTitle("Observation time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]"); + hvs2.SetTitle("On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]"); + hvstM.SetTitle("MJD vs time;Time;MJD"); + hvst0.SetTitle("Zenith distance vs time;Time;Zd [\\circ]"); + hvst1.SetTitle("Camera trigger rate vs time;Time;Rate [Hz]"); + hvst2.SetTitle("Camera trigger rate vs Zd; Zd [\\circ];Rate [Hz]"); + hvst.SetTitle("Average sensor compartment temperature;Time;T [\\circ C]"); + hvsh.SetTitle("Average camera humidity;Time;Hum [%]"); + + MFillH fillvs1(&hvs1, "", "FillObsTime"); + MFillH fillvs2(&hvs2, "", "FillOnTime"); + MFillH fillvstM(&hvstM, "MTime", "FillMjdVsTime"); + MFillH fillvst0(&hvst0, "MTimeDrive", "FillZdVsTime"); + MFillH fillvst1(&hvst1, "MTimeRates", "FillRate"); + MFillH fillvst2(&hvst2, "", "FillRateVsZd"); + MFillH fillvst(&hvst, "MTimeTemperatures", "FillCamTemp"); + MFillH fillvsh(&hvsh, "MTimeHumidity", "FillCamHum"); + fillvs1.SetNameTab("ObsTime"); + fillvs2.SetNameTab("OnTime"); + fillvst1.SetNameTab("Rate"); + fillvst2.SetNameTab("RateZd"); + fillvst.SetNameTab("CamTemp"); + fillvsh.SetNameTab("CamHum"); + + hvs1.SetWeight("MReportRates.fElapsedTime"); + hvs2.SetWeight("MReportRates.fElapsedOnTime"); + + hvst1.SetMinimum(0); + hvst2.SetDrawOption("colz"); + + MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/); + hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)"); + hrate.SetName("RateVsRun"); + hrate.SetTitle("Event rate after quality cuts;File Id;Event Rate [Hz];"); + hrate.InitLabels(MH3::kLabelsX); + hrate.GetHist().SetMinimum(0); + + MFillH fillr(&hrate, "", "FillRateVsRun"); + tlist.AddToListAfter(&fillr, &contq, "Events"); + + tlist.Replace(&fill1b); + + tlist2.Replace(&fill2b); + tlist2.Replace(&fill3b); + tlist2.Replace(&fill3d); + tlist2.Replace(&fill4b); + tlist2.Replace(&fill5b); + tlist2.Replace(&fill6b); + tlist2.Replace(&fill7b); + tlist2.Replace(&falpha2); + tlist2.Replace(&fenergy2); + if (!ismc) + { + tlist.AddToList(&fillvs1, "Rates"); + tlist.AddToList(&fillvs2, "Rates"); + tlist.AddToList(&fillvst1, "Rates"); + tlist.AddToList(&fillvst2, "Rates"); + tlist.AddToList(&write2, "Rates"); + tlist.AddToList(&fillvst, "Temperatures"); + tlist.AddToList(&fillvsh, "Humidity"); + tlist.AddToList(&fillvst0, "Drive"); + tlist.AddToList(&fillvstM, "Drive"); + } + + // by setting it here it is distributed to all consecutive tasks + tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); + + par.SetVal(1); + weight.SetVal(1); + + // Execute first analysis + if (!evtloop.Eventloop()) + return 2; + + if (!evtloop.GetDisplay()) + return 2; + + delete readondata; + delete readoffdata; + + TObjArray cont; + //cont.Add(const_cast(GetEnv())); + //cont.Add(const_cast(&fSequence)); + + //TNamed cmdline("CommandLine", fCommandLine.Data()); + //cont.Add(&cmdline); + + if (d) + { + TString title = "-- Analysis #"; + title += dataset; + title += " --"; + d->SetTitle(title, kFALSE); + + cont.Add(d); + TString path; + path += fname1; + + d->Write(); + } + +// if (!WriteContainer(cont, fname1, "RECREATE")) +// return; + + return 0; +} Index: branches/Mars_McMismatchStudy/fact/analysis/mc/star.C.NEWCLEANING =================================================================== --- branches/Mars_McMismatchStudy/fact/analysis/mc/star.C.NEWCLEANING (revision 17962) +++ branches/Mars_McMismatchStudy/fact/analysis/mc/star.C.NEWCLEANING (revision 17962) @@ -0,0 +1,420 @@ +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#if !defined(__CINT__) || defined(__MAKECINT__) + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.h" + +#include "../mcore/DrsCalib.h" +#include "MDrsCalibration.h" +//#include "MLogManip.h" +#include "MExtralgoSpline.h" +#include "MSequence.h" +#include "MStatusArray.h" +#include "MHCamera.h" +#include "MJob.h" +#include "MWriteRootFile.h" +#include "MHCamera.h" +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" +#include "MDirIter.h" +#include "MTaskList.h" +#include "MFDataPhrase.h" +#include "MArrayF.h" +#include "MBadPixelsTreat.h" +#include "MCalibrateDrsTimes.h" +#include "MHSectorVsTime.h" +#include "MHCamEvent.h" +#include "MExtractTimeAndChargeSpline.h" +#include "MFillH.h" +#include "MDrsCalibApply.h" +#include "MGeomApply.h" +#include "MContinue.h" +#include "MRawFitsRead.h" +#include "MEvtLoop.h" +#include "MParList.h" +#include "MStatusDisplay.h" +#include "MDrsCalibrationTime.h" +#include "MH3.h" +#include "MGeomCamFACT.h" +#include "MCalibrateFact.h" +#include "MParameters.h" +#include "MWriteAsciiFile.h" +#include "MFMagicCuts.h" + +#include "MMuonSetup.h" +#include "MReadMarsFile.h" +#include "MHillasCalc.h" +#include "MHn.h" +#include "MMuonSearchParCalc.h" +#include "MMuonCalibParCalc.h" +#include "MBinning.h" +#include "MImgCleanTime.h" + + +#endif + +using namespace std; + +int star(const TString callistofile, TString outfile, + TString displayfile, TString displaytitle, + const TString txtfile, + Double_t lvl1=4.0, Double_t lvl2=2.5, Double_t deltat = 17.5); + + +int star(const TString infile = "callisto.root", TString outfile = "", + Double_t lvl1=4.0, Double_t lvl2=2.5, Double_t deltat = 17.5) +{ + return star(infile, outfile, "", "", "", lvl1, lvl2, deltat); +} + +int star(const TString callistofile, TString outfile, + TString displayfile, TString displaytitle, + const TString txtfile, + Double_t lvl1, Double_t lvl2, Double_t deltat) +{ + + //Double_t lvl1=7.8, + //Double_t lvl2=3.9, + // lvl1 = 2.5; + // lvl2 = 0.5; + // ------------------------------------------------------ + + if (displaytitle.IsNull()) + displaytitle = gSystem->BaseName(callistofile); + + FileStat_t fstat; + int rc = gSystem->GetPathInfo(outfile, fstat); + bool isdir = !rc || R_ISDIR(fstat.fMode); + + char *buf = gSystem->ConcatFileName(outfile, "star.root"); + outfile = buf; + delete [] buf; + + if (displayfile.IsNull()) + { + displayfile = outfile; + displayfile.Insert(displayfile.Last('.'), "-display"); + } + else + { + if (isdir && gSystem->DirName(displayfile)==TString(".")) + { + buf = gSystem->ConcatFileName(outfile, displayfile); + displayfile = buf; + delete [] buf; + } + } + + // ------------------------------------------------------ + + gLog.Separator("Star"); + gLog << all << "Calculate image parameters of sequence "; + gLog << endl; + gLog << "Callisto file: " << callistofile << endl; + gLog << "Output file: " << outfile << endl; + gLog << "Display file: " << displayfile << endl; + gLog << "Display title: " << displaytitle << endl; + gLog << endl; + gLog << "lvl1: " << lvl1 <150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut"); + // Filter to calculate further muon parameters +// MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&" +// "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut"); + MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&" + "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4", + "MuonSearchCut"); + + // Filter to fill the MHMuonPar + MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2", + "MuonFinalCut"); + // Filter to write Muons to Muon tree + //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut"); + + // --------------------------------------------------------- + + MStatusDisplay *d = new MStatusDisplay; + + MTaskList tlist; + + MParList plist2; + plist2.AddToList(&tlist); + plist2.AddToList(&bins1); + plist2.AddToList(&bins3); + plist2.AddToList(&bins4); + plist2.AddToList(&binsM1); + plist2.AddToList(&binsM2); + plist2.AddToList(&binsM3); + plist2.AddToList(&binsM4); + plist2.AddToList(&binsM5); + plist2.AddToList(&binsM6); + plist2.AddToList(&binsM7); + plist2.AddToList(&binsM8); + plist2.AddToList(&binsM9); + + MMuonSetup muonsetup; + plist2.AddToList(&muonsetup); + + MEvtLoop loop; + loop.SetDisplay(d); + loop.SetParList(&plist2); + + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + read.AddFile(callistofile); + + MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData"); + + MGeomApply apply; + + MImgCleanTime clean; + clean.SetMinCount(0); + clean.SetMinSize(25); + clean.SetDeltaT(1.25*17.5*0.1111); + + MHillasCalc hcalc; + hcalc.Disable(MHillasCalc::kCalcConc); + + MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/); + hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)"); + hrate.SetName("Rate"); + hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];"); + hrate.InitLabels(MH3::kLabelsX); + hrate.GetHist().SetMinimum(0); + + MFillH fillR(&hrate, "", "FillRate"); + + // ------------------ Setup histograms and fill tasks ---------------- + MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]"); + MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction"); + MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]"); + evtC1.SetErrorSpread(kFALSE); + evtC2.SetErrorSpread(kFALSE); + evtC2.SetThreshold(0); + + MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam"); + MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels"); + MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos"); + + MFillH fillD1("MHHillas", "MHillas", "FillHillas"); + MFillH fillD2("MHHillasExt", "", "FillHillasExt"); + MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc"); + MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar"); + MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar"); + MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize"); + + MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)"); + halpha.SetName("AvsD;Dist;Alpha"); + + MFillH filla(&halpha); + filla.SetNameTab("AvsD"); + filla.SetDrawOption("colz"); + + // ----------------------- Muon Analysis ---------------------- + //writem.SetFilter(&fmuon4); + + MHn mhn1; + mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg"); + mhn1.InitName("MuonRadius"); + mhn1.InitTitle("MuonRadius"); + + mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg"); + mhn1.InitName("MuonDeviation"); + mhn1.InitTitle("MuonDeviation"); + + mhn1.AddHist("MMuonSearchPar.fTime"); + mhn1.InitName("MuonTime"); + mhn1.InitTitle("MuonTime"); + + mhn1.AddHist("MMuonSearchPar.fTimeRms"); + mhn1.InitName("MuonTimeRms"); + mhn1.InitTitle("MuonTimeRms"); + + MHn mhn3; + mhn3.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg"); + mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation"); + mhn3.InitTitle("MuonRadius"); + mhn3.SetDrawOption("colz"); + + MHn mhn2; + mhn2.AddHist("MMuonCalibPar.fArcPhi"); + mhn2.InitName("MuonArcPhi"); + mhn2.InitTitle("MuonArcPhi"); + + mhn2.AddHist("MMuonCalibPar.fArcWidth"); + mhn2.InitName("MuonArcWidth"); + mhn2.InitTitle("MuonArcWidth"); + + mhn2.AddHist("MMuonCalibPar.fMuonSize"); + mhn2.InitName("MuonSize"); + mhn2.InitTitle("MuonSize"); + + mhn2.AddHist("MMuonCalibPar.fRelTimeSigma"); + mhn2.InitName("MuonRelTimeSigma"); + mhn2.InitTitle("MuonRelTimeSigma"); + + MHn mhn4; + mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth"); + mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth"); + mhn4.InitTitle("MuonArcPhi"); + mhn4.SetDrawOption("colz"); + + MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar"); + MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar"); + MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius"); + MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi"); + fillmhn1.SetNameTab("MuonS"); + fillmhn3.SetNameTab("DevVsRad"); + fillmhn2.SetNameTab("MuonC"); + fillmhn4.SetNameTab("WidthVsPhi"); + fillmhn2.SetFilter(&fmuon2); + fillmhn4.SetFilter(&fmuon2); + + MMuonSearchParCalc muscalc; + //muscalc.SetFilter(&fmuon1); + + MMuonCalibParCalc mcalc; + mcalc.SetFilter(&fmuon2); + + MFillH fillmuon("MHSingleMuon", "", "FillMuon"); + MFillH fillmpar("MHMuonPar", "", "FillMuonPar"); + fillmuon.SetFilter(&fmuon2); + fillmpar.SetFilter(&fmuon3); + fillmuon.SetBit(MFillH::kDoNotDisplay); + + // --------------------------------------------------------------- + + MWriteAsciiFile writeascii(txtfile, "MHillas"); + writeascii.AddColumns("MHillasExt"); + writeascii.AddColumns("MHillasSrc"); + + + MWriteRootFile write(outfile, "RECREATE", "Image parameters", 2); + write.AddContainer("MTime", "Events", kFALSE); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MImagePar", "Events"); + write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MRawEvtHeader", "Events"); + write.AddContainer("ThetaSquared", "Events"); + write.AddContainer("Disp", "Events"); + write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE); + write.AddContainer("MMcEvt", "Events", kFALSE); + write.AddContainer("IncidentAngle", "Events", kFALSE); + write.AddContainer("MPointingPos", "Events", kFALSE); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MGeomCam", "RunHeaders"); + write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); + + MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0", + "MuonHistCut"); + fillmhn2.SetFilter(&fmuonhn); + + TArrayD param(12); + // Parametrization of Disp + param[0] = 1.15136; // constant + param[8] = 0.0681437; // slope + param[9] = 2.62932; // leak + param[10] = 1.51279; // size-offset + param[11] = 0.0507821; // size-slope + // Parametrization for sign of disp (m3long, slope) + param[5] = -0.07; + param[6] = 7.2; + param[7] = 0.5; + // ThetaSq-Cut + param[1] = 0.11; // 0.215 + // Area-Cut + param[2] = 0.215468; + param[3] = 5.63973; + param[4] = 0.0836169; + MFMagicCuts cuts; + cuts.SetHadronnessCut(MFMagicCuts::kNoCut); + cuts.SetVariables(param); + + + tlist.AddToList(&read); + tlist.AddToList(&cont); + tlist.AddToList(&apply); + tlist.AddToList(&clean); + tlist.AddToList(&hcalc); + tlist.AddToList(&cuts); + tlist.AddToList(&fillC1); + tlist.AddToList(&fillC2); + tlist.AddToList(&fillC3); + tlist.AddToList(&fillR); + tlist.AddToList(&fillD1); + tlist.AddToList(&fillD2); + tlist.AddToList(&fillD3); + tlist.AddToList(&fillD4); + tlist.AddToList(&fillD5); + tlist.AddToList(&fillD6); + tlist.AddToList(&filla); + //tlist.AddToList(&fmuon1); + tlist.AddToList(&muscalc); + tlist.AddToList(&fillmhn1); + tlist.AddToList(&fillmhn3); + tlist.AddToList(&fmuon2); + tlist.AddToList(&fillmuon); + tlist.AddToList(&mcalc); + tlist.AddToList(&fmuonhn); + tlist.AddToList(&fillmhn2); + tlist.AddToList(&fillmhn4); + tlist.AddToList(&fmuon3); + tlist.AddToList(&fillmpar); + //tlist.AddToList(&fmuon4); + //tlist.AddToList(&writem); + tlist.AddToList(&write); + if (!txtfile.IsNull()) + tlist.AddToList(&writeascii); + + if (!loop.Eventloop()) + return 3; + + if (!loop.GetDisplay()) + return 4; + + d->SetTitle(displaytitle, kFALSE); + d->SaveAs(displayfile); + + return 0; +} + Index: branches/Mars_McMismatchStudy/fact/analysis/star.C =================================================================== --- branches/Mars_McMismatchStudy/fact/analysis/star.C (revision 17961) +++ branches/Mars_McMismatchStudy/fact/analysis/star.C (revision 17962) @@ -1,2 +1,3 @@ +//Star version from revision 17885 (used in the processing with the old image cleaning) #include "MLogManip.h" @@ -303,5 +304,5 @@ TString file = Form("/scratch/fact/sequences/%04d/%02d/%02d/%06d_%03d.seq", - night/10000, (night/100)%100, night%100, num); + night/10000, (night/100)%100, night%100, night, num); return star(file.Data(), lvl1, lvl2, inpath, outpath); Index: branches/Mars_McMismatchStudy/fact/analysis/star.C.ORIGINAL =================================================================== --- branches/Mars_McMismatchStudy/fact/analysis/star.C.ORIGINAL (revision 17962) +++ branches/Mars_McMismatchStudy/fact/analysis/star.C.ORIGINAL (revision 17962) @@ -0,0 +1,297 @@ +#include "MLogManip.h" + +int star(const char *seqfile="sequences/20111205_013.seq", + const char *inpath = "callisto_new", const char *outpath = "callisto_new") +{ + double deltat = 17.5; + + // The sequence file which defines the files for the analysis + MSequence seq(seqfile); + if (!seq.IsValid()) + { + gLog << err << "ERROR - Sequence invalid!" << endl; + return 1; + } + + // ------------------------------------------------------ + + gLog.Separator("Star"); + gLog << all << "Calculate image parameters of sequence "; + gLog << seq.GetFileName() << endl; + gLog << endl; + + // ------------------------------------------------------ + + gLog << "Inpath: " << inpath << endl; + gLog << "Outpath: " << outpath << endl; + + const TString rule(Form("s/([0-9]+_[0-9]+)_C[.]root?$/%s\\/$1_I.root/", + MJob::Esc(outpath).Data())); + gLog << "Rule: " << rule << endl; + + MDirIter iter; + if (seq.GetRuns(iter, MSequence::kFactCal, inpath)<=0) + { + gLog << err << "ERROR - Sequence valid but without files." << endl; + return 2; + } + + iter.Print(); + + gLog.Separator(); + + // -------------------------------------------------------- + + MBinning bins1( 36, -90, 90, "BinningAlpha"); + MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin"); + MBinning bins4( 25, 0, 2.5, "BinningDist"); + + MBinning binsM1(100, 0, 5, "BinningMuonRadius"); + MBinning binsM2( 60, 0, 0.3, "BinningMuonDeviation"); + MBinning binsM3(150, 0, 60, "BinningMuonTime"); + MBinning binsM4(300, 0, 30, "BinningMuonTimeRms"); + MBinning binsM5( 20, 0, 360, "BinningMuonArcPhi"); + MBinning binsM6( 50, 0, 0.5, "BinningMuonArcWidth"); + MBinning binsM7( 30, 5, 5000, "BinningMuonSize", "log"); + MBinning binsM8(100, 0, 5, "BinningMuonRelTimeSigma"); + + MBinning binsM9(100, 0, 5, "BinningRadius"); + + // --------------------------------------------------------- + + // Filter to start muon analysis + //MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut"); + // Filter to calculate further muon parameters +// MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&" +// "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut"); + MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&" + "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4", + "MuonSearchCut"); + + // Filter to fill the MHMuonPar + MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2", + "MuonFinalCut"); + // Filter to write Muons to Muon tree + //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut"); + + // --------------------------------------------------------- + + MStatusDisplay *d = new MStatusDisplay; + + MTaskList tlist; + + MParList plist2; + plist2.AddToList(&tlist); + plist2.AddToList(&bins1); + plist2.AddToList(&bins3); + plist2.AddToList(&bins4); + plist2.AddToList(&binsM1); + plist2.AddToList(&binsM2); + plist2.AddToList(&binsM3); + plist2.AddToList(&binsM4); + plist2.AddToList(&binsM5); + plist2.AddToList(&binsM6); + plist2.AddToList(&binsM7); + plist2.AddToList(&binsM8); + plist2.AddToList(&binsM9); + + MMuonSetup muonsetup; + plist2.AddToList(&muonsetup); + + MEvtLoop loop; + loop.SetDisplay(d); + loop.SetParList(&plist2); + + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + read.AddFiles(iter); + + MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData"); + //MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0", "SelectDat"); + //MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); + + MGeomApply apply; + + MImgCleanTime clean; + clean.SetMinCount(0); + clean.SetMinSize(25); + clean.SetDeltaT(1.25*17.5*0.1111); + + MHillasCalc hcalc; + hcalc.Disable(MHillasCalc::kCalcConc); + + MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/); + hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)"); + hrate.SetName("Rate"); + hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];"); + hrate.InitLabels(MH3::kLabelsX); + hrate.GetHist().SetMinimum(0); + + MFillH fillR(&hrate, "", "FillRate"); + + // ------------------ Setup histograms and fill tasks ---------------- + MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]"); + MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction"); + MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]"); + evtC1.SetErrorSpread(kFALSE); + evtC2.SetErrorSpread(kFALSE); + evtC2.SetThreshold(0); + + MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam"); + MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels"); + MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos"); + + MFillH fillD1("MHHillas", "MHillas", "FillHillas"); + MFillH fillD2("MHHillasExt", "", "FillHillasExt"); + MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc"); + MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar"); + MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar"); + MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize"); + + MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)"); + halpha.SetName("AvsD;Dist;Alpha"); + + MFillH filla(&halpha); + filla.SetNameTab("AvsD"); + filla.SetDrawOption("colz"); + + // ----------------------- Muon Analysis ---------------------- + //writem.SetFilter(&fmuon4); + + MHn mhn1; + mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg"); + mhn1.InitName("MuonRadius"); + mhn1.InitTitle("MuonRadius"); + + mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg"); + mhn1.InitName("MuonDeviation"); + mhn1.InitTitle("MuonDeviation"); + + mhn1.AddHist("MMuonSearchPar.fTime"); + mhn1.InitName("MuonTime"); + mhn1.InitTitle("MuonTime"); + + mhn1.AddHist("MMuonSearchPar.fTimeRms"); + mhn1.InitName("MuonTimeRms"); + mhn1.InitTitle("MuonTimeRms"); + + MHn mhn3; + mhn3.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg"); + mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation"); + mhn3.InitTitle("MuonRadius"); + mhn3.SetDrawOption("colz"); + + MHn mhn2; + mhn2.AddHist("MMuonCalibPar.fArcPhi"); + mhn2.InitName("MuonArcPhi"); + mhn2.InitTitle("MuonArcPhi"); + + mhn2.AddHist("MMuonCalibPar.fArcWidth"); + mhn2.InitName("MuonArcWidth"); + mhn2.InitTitle("MuonArcWidth"); + + mhn2.AddHist("MMuonCalibPar.fMuonSize"); + mhn2.InitName("MuonSize"); + mhn2.InitTitle("MuonSize"); + + mhn2.AddHist("MMuonCalibPar.fRelTimeSigma"); + mhn2.InitName("MuonRelTimeSigma"); + mhn2.InitTitle("MuonRelTimeSigma"); + + MHn mhn4; + mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth"); + mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth"); + mhn4.InitTitle("MuonArcPhi"); + mhn4.SetDrawOption("colz"); + + MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar"); + MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar"); + MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius"); + MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi"); + fillmhn1.SetNameTab("MuonS"); + fillmhn3.SetNameTab("DevVsRad"); + fillmhn2.SetNameTab("MuonC"); + fillmhn4.SetNameTab("WidthVsPhi"); + fillmhn2.SetFilter(&fmuon2); + fillmhn4.SetFilter(&fmuon2); + + MMuonSearchParCalc muscalc; + //muscalc.SetFilter(&fmuon1); + + MMuonCalibParCalc mcalc; + mcalc.SetFilter(&fmuon2); + + MFillH fillmuon("MHSingleMuon", "", "FillMuon"); + MFillH fillmpar("MHMuonPar", "", "FillMuonPar"); + fillmuon.SetFilter(&fmuon2); + fillmpar.SetFilter(&fmuon3); + fillmuon.SetBit(MFillH::kDoNotDisplay); + + // --------------------------------------------------------------- + + MWriteRootFile write(2, rule, "RECREATE", "Image parameters"); // EffectiveOnTime + write.AddContainer("MTime", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MImagePar", "Events"); + write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MRawEvtHeader", "Events"); + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MGeomCam", "RunHeaders"); + + MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0", + "MuonHistCut"); + fillmhn2.SetFilter(&fmuonhn); + + tlist.AddToList(&read); + tlist.AddToList(&cont); + tlist.AddToList(&apply); + tlist.AddToList(&clean); + tlist.AddToList(&hcalc); + tlist.AddToList(&fillC1); + tlist.AddToList(&fillC2); + tlist.AddToList(&fillC3); + tlist.AddToList(&fillR); + tlist.AddToList(&fillD1); + tlist.AddToList(&fillD2); + tlist.AddToList(&fillD3); + tlist.AddToList(&fillD4); + tlist.AddToList(&fillD5); + tlist.AddToList(&fillD6); + tlist.AddToList(&filla); + //tlist.AddToList(&fmuon1); + tlist.AddToList(&muscalc); + tlist.AddToList(&fillmhn1); + tlist.AddToList(&fillmhn3); + tlist.AddToList(&fmuon2); + tlist.AddToList(&fillmuon); + tlist.AddToList(&mcalc); + tlist.AddToList(&fmuonhn); + tlist.AddToList(&fillmhn2); + tlist.AddToList(&fillmhn4); + tlist.AddToList(&fmuon3); + tlist.AddToList(&fillmpar); + //tlist.AddToList(&fmuon4); + //tlist.AddToList(&writem); + tlist.AddToList(&write); + + if (!loop.Eventloop()) + return 3; + + if (!loop.GetDisplay()) + return 4; + + TString title = "-- Image parameters #"; + title += seq.GetSequence(); + title += " --"; + d->SetTitle(title, kFALSE); + + TString path; + path += Form("%s/20%06d_%03d-images.root", outpath, + seq.GetSequence()/1000, seq.GetSequence()%1000); + + d->SaveAs(path); + + return 0; +} Index: branches/Mars_McMismatchStudy/fact/processing/buildseqentries.C =================================================================== --- branches/Mars_McMismatchStudy/fact/processing/buildseqentries.C (revision 17961) +++ branches/Mars_McMismatchStudy/fact/processing/buildseqentries.C (revision 17962) @@ -584,5 +584,5 @@ TString query; - query = "SELECT fRunID FROM runinfo WHERE fExcludedFDAKEY=1 AND "; + query = "SELECT fRunID FROM RunInfo WHERE fExcludedFDAKEY=1 AND "; query += cond; query += " AND fRunTypeKEY=4 AND fRunID BETWEEN "; @@ -599,5 +599,5 @@ // Get DRS file with pedestal (roi<1024) - query = "SELECT MAX(fRunID) FROM runinfo WHERE "; + query = "SELECT MAX(fRunID) FROM RunInfo WHERE "; query += cond; query += " AND fHasDrsFile=1 AND (fDrsStep IS NULL OR fDrsStep=2) AND fRunTypeKEY=2"; Index: branches/Mars_McMismatchStudy/fact/processing/camhum.C =================================================================== --- branches/Mars_McMismatchStudy/fact/processing/camhum.C (revision 17962) +++ branches/Mars_McMismatchStudy/fact/processing/camhum.C (revision 17962) @@ -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/Mars_McMismatchStudy/fact/processing/drstemp.C =================================================================== --- branches/Mars_McMismatchStudy/fact/processing/drstemp.C (revision 17961) +++ branches/Mars_McMismatchStudy/fact/processing/drstemp.C (revision 17962) @@ -163,9 +163,9 @@ float temp[160]; - if (file.SetPtrAddress("temp", temp, 82)) - { - drstemp82(file, beg, end); - return; - } +// if (file.SetPtrAddress("temp", temp, 82)) +// { +// drstemp82(file, beg, end); +// return; +// } file.SetPtrAddress("temp", temp, 160); Index: branches/Mars_McMismatchStudy/fact/processing/lastth.C =================================================================== --- branches/Mars_McMismatchStudy/fact/processing/lastth.C (revision 17962) +++ branches/Mars_McMismatchStudy/fact/processing/lastth.C (revision 17962) @@ -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/Mars_McMismatchStudy/fact/processing/threshold.C =================================================================== --- branches/Mars_McMismatchStudy/fact/processing/threshold.C (revision 17962) +++ branches/Mars_McMismatchStudy/fact/processing/threshold.C (revision 17962) @@ -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/Mars_McMismatchStudy/macros/optim/optimenergy.C =================================================================== --- branches/Mars_McMismatchStudy/macros/optim/optimenergy.C (revision 17961) +++ branches/Mars_McMismatchStudy/macros/optim/optimenergy.C (revision 17962) @@ -46,5 +46,5 @@ */ - // opt.SetPathOut("optimenergy.root"); - opt.RunEnergy("ganymedmcpart.root", r); + // opt.SetPathOut("optimenergy.root"); + opt.RunEnergy("/unsaved_data/toscanos/FACT/CRAB/ganymedMC_Full_20140716-summary.root", r); } Index: branches/Mars_McMismatchStudy/mcore/PixelMap.h =================================================================== --- branches/Mars_McMismatchStudy/mcore/PixelMap.h (revision 17961) +++ branches/Mars_McMismatchStudy/mcore/PixelMap.h (revision 17962) @@ -27,5 +27,5 @@ int cbpx; /// Hardware index as CBPX int gapd; /// gAPD index -// float Vgapd; /// gAPD Bias voltage + float Vgapd; /// gAPD Bias voltage int hv_board; /// Bias suppply board int hv_channel; @@ -41,4 +41,5 @@ int count() const { return pixel()>3 ? 5 : 4; } int hv() const { return hv_channel+hv_board*32; } + bool isTM() const { return pixel()%9==8; } operator bool() const { return index>=0; } @@ -74,6 +75,5 @@ std::stringstream str(buf); - int idummy; - float fdummy; + int idummy; PixelMapEntry entry; @@ -84,5 +84,5 @@ str >> idummy; str >> entry.gapd; - str >> fdummy; //entry.Vgapd; + str >> entry.Vgapd; str >> entry.hv_board; str >> entry.hv_channel; @@ -228,5 +228,5 @@ std::ifstream fin(fname); - int l = 0; + long long int l = 0; std::string buf; Index: branches/Mars_McMismatchStudy/resmc/fact/fact-reflector-real_ceres.txt =================================================================== --- branches/Mars_McMismatchStudy/resmc/fact/fact-reflector-real_ceres.txt (revision 17962) +++ branches/Mars_McMismatchStudy/resmc/fact/fact-reflector-real_ceres.txt (revision 17962) @@ -0,0 +1,48 @@ +# x, y and z describe the x-, y- and z-coordinate of the center of the mirror (in cm) +# nx, ny and nz describe the x-, y- and z- component of the normal vector of the mirror +# which is normed +# F is the focal distance of the spherical mirror (in cm) +# Type: Hex is a hexagonal mirror +# fD = 60.6 is the distance of two parallel sides of the hexagon (in cm) +# the normalvectors point to (0,0,978.132) +# Ann-Kristin Overkemping, 03 2012 +# +# added information about PSF: theta1 is the angle according to the smearing radius +# of 1 sigma of the PSF area (in rad) +# Ann-Kristin Overkemping, 06 2013 +# +# recalculated all values according to real distances of mirrors to camera, +# real reflection positions and partly tilted mirror mountings +# Ann-Kristin Overkemping, 07 2013 +######################################################################################### +# x y z nx ny nz F Type fD +31 53.694 6.745 -0.032935 -0.055560 0.997912 487.624 Hex 60.6 +62 0 10.374 -0.064847 0.000127 0.997895 488.431 Hex 60.6 +31 -53.694 7.753 -0.031980 0.054888 0.997980 488.286 Hex 60.6 +-31 -53.694 8.055 0.031505 0.055296 0.997973 488.236 Hex 60.6 +-62 0 9.064 0.062825 -0.000215 0.998025 488.218 Hex 60.6 +-31 53.694 5.434 0.031224 -0.056422 0.997919 488.133 Hex 60.6 +93.072 52.915 13.993 -0.096734 -0.054764 0.993803 488.453 Hex 60.6 +92.516 -52.915 17.144 -0.096221 0.054676 0.993857 488.496 Hex 60.6 +0 -107 10.799 0.000122 0.109760 0.993958 488.506 Hex 60.6 +-92.829 -52.915 15.372 0.095753 0.054568 0.993908 488.637 Hex 60.6 +-93.367 52.915 12.319 0.095154 -0.055501 0.993914 488.696 Hex 60.6 +0 107 12.234 -0.000817 -0.111529 0.993761 488.441 Hex 60.6 +62 107.387 16.266 -0.064900 -0.111481 0.991645 488.783 Hex 60.6 +123.448 0 22.728 -0.129106 -0.000044 0.991631 488.938 Hex 60.6 +62 -107.387 14.922 -0.063895 0.110571 0.991812 488.791 Hex 60.6 +-62 -107.387 14.302 0.063484 0.110542 0.991842 489.000 Hex 60.6 +-123.969 0 19.773 0.127436 -0.000885 0.991846 489.075 Hex 60.6 +-62 107.387 16.163 0.063208 -0.111675 0.991733 488.767 Hex 60.6 +-31 161.168 27.391 0.031514 -0.168234 0.985243 490.212 Hex 60.6 +31 161.116 27.687 -0.032460 -0.168136 0.985229 490.199 Hex 60.6 +123.743 107.331 29.558 -0.129328 -0.112322 0.985220 489.448 Hex 60.6 +154.708 53.582 29.755 -0.161510 -0.056332 0.985262 490.166 Hex 60.6 +155.212 -53.582 26.899 -0.161225 0.055293 0.985368 490.109 Hex 60.6 +124.246 -107.331 26.702 -0.128675 0.110978 0.985458 489.448 Hex 60.6 +31 -161.481 25.619 -0.032036 0.167185 0.985405 490.071 Hex 60.6 +-31 -161.481 25.619 0.032047 0.166922 0.985449 489.978 Hex 60.6 +-124.420 -107.331 25.717 0.128309 0.111059 0.985496 489.479 Hex 60.6 +-154.795 -53.582 29.263 0.160523 0.055663 0.985461 489.975 Hex 60.6 +-155.229 53.582 26.801 0.160255 -0.056631 0.985450 489.956 Hex 60.6 +-124.038 107.331 27.884 0.127897 -0.112420 0.985395 489.439 Hex 60.6 Index: branches/Mars_McMismatchStudy/resmc/fact/night-sky-la-palma.txt =================================================================== --- branches/Mars_McMismatchStudy/resmc/fact/night-sky-la-palma.txt (revision 17962) +++ branches/Mars_McMismatchStudy/resmc/fact/night-sky-la-palma.txt (revision 17962) @@ -0,0 +1,1056 @@ +250.014 0.800477 +260.014 0.800477 +270.014 0.800477 +280.014 0.800477 +290.014 0.800477 +311.014 0.800477 +321.014 0.900477 +331.014 1.300477 +341.014 1.500477 +346.014 1.400477 +350.014 0.650477 +351.014 0.600477 +351.594 0.790559 +352.174 0.817864 +352.754 1.72687 +353.333 2.17965 +353.913 3.15677 +354.493 2.3622 +355.072 2.84384 +355.652 2.64883 +356.232 2.25675 +356.812 1.30812 +357.391 1.5174 +357.971 0.763332 +358.551 1.47381 +359.13 1.35968 +359.71 0.970784 +360.29 0.385688 +360.87 1.20294 +361.449 1.14437 +362.029 1.56325 +362.609 1.80576 +363.188 2.09278 +363.768 2.08063 +364.348 1.57732 +364.928 1.73802 +365.507 1.48751 +366.087 1.30638 +366.667 1.06186 +367.246 0.951596 +367.826 1.28396 +368.406 1.54215 +368.986 1.07401 +369.565 1.30309 +370.145 1.34888 +370.725 1.78687 +371.304 1.34982 +371.884 1.30309 +372.464 1.88863 +373.043 2.03467 +373.623 2.38421 +374.203 2.06903 +374.783 1.88457 +375.362 1.57732 +375.942 2.09278 +376.522 2.09278 +377.101 1.84677 +377.681 2.18803 +378.261 1.43243 +378.841 1.06186 +379.42 1.53294 +380 2.19191 +380.58 1.86504 +381.159 1.52089 +381.739 1.88482 +382.319 1.88409 +382.899 1.83505 +383.478 2.19812 +384.058 2.20885 +384.638 1.60353 +385.217 0.932006 +385.797 1.06186 +386.377 1.31959 +386.957 1.00041 +387.536 0.735206 +388.116 0.69447 +388.696 1.02703 +389.275 1.67438 +389.855 1.88482 +390.435 1.93853 +391.014 1.79855 +391.594 1.41249 +392.174 2.49735 +392.754 7.38119 +393.333 6.15931 +393.913 7.70165 +394.493 2.97863 +395.072 1.8909 +395.652 1.39714 +396.232 1.40765 +396.812 1.26982 +397.391 1.31959 +397.971 1.57732 +398.551 1.56621 +399.13 1.25935 +399.71 1.53544 +400.87 2.60825 +401.449 2.02782 +402.029 1.96907 +402.609 2.29094 +403.188 2.49029 +403.768 3.31409 +404.348 3.1445 +404.928 3.67122 +405.507 3.01399 +406.087 2.91575 +406.667 2.94285 +407.246 2.61936 +407.826 2.29248 +408.406 2.66079 +408.986 2.31192 +409.565 1.84651 +410.145 2.26986 +410.725 2.05305 +411.304 2.05653 +411.884 2.05703 +412.464 2.09278 +413.043 2.60825 +413.623 2.12761 +414.203 2.48597 +414.783 2.17986 +415.362 2.58449 +415.942 3.2583 +416.522 3.52439 +417.101 3.1126 +417.681 3.47758 +418.261 3.12371 +418.841 2.53188 +419.42 2.69254 +420 2.86794 +420.58 2.38121 +421.159 2.06385 +421.739 2.09278 +422.319 2.10494 +422.899 1.79426 +423.478 2.21439 +424.058 1.92404 +424.638 2.09278 +425.217 2.65481 +425.797 3.24531 +426.377 3.10964 +426.957 2.94544 +427.536 2.20757 +428.116 2.26693 +428.696 2.39956 +429.275 2.35052 +429.855 2.72267 +430.435 2.32886 +431.014 2.38403 +431.594 2.93163 +432.174 3.12371 +432.754 2.62232 +433.333 2.40666 +433.913 2.55193 +434.493 4.17461 +435.072 4.04616 +435.652 5.0667 +436.232 4.22379 +436.812 3.86105 +437.391 2.53719 +437.971 2.55944 +438.551 2.60825 +439.13 3.24115 +439.71 3.41008 +440.29 3.8619 +440.87 3.82629 +441.449 3.69268 +442.029 3.23343 +442.609 3.12371 +443.188 3.29338 +443.768 3.4923 +444.348 3.33068 +444.928 3.16149 +445.507 2.88247 +446.087 3.24531 +446.667 2.95135 +447.246 3.27059 +447.826 3.13532 +448.406 2.86598 +448.986 3.16149 +449.565 3.79291 +450.145 3.62429 +450.725 3.03391 +451.304 3.14992 +451.884 3.15825 +452.464 3.63918 +453.043 3.97445 +453.623 3.61542 +454.203 3.68355 +454.783 3.5948 +455.362 3.36266 +455.942 3.73911 +456.522 3.32631 +457.101 3.83126 +457.681 3.63918 +458.261 4.01804 +458.841 3.65133 +459.42 3.84172 +460 3.95073 +460.58 3.55489 +461.159 3.6014 +461.739 3.55149 +462.319 3.57773 +462.899 3.4923 +463.478 3.80635 +464.058 3.82799 +464.638 3.74889 +465.217 3.57773 +465.797 3.78406 +466.377 3.45841 +466.957 3.63918 +467.536 4.18085 +468.116 4.01804 +468.696 4.15464 +469.275 3.70296 +469.855 3.38144 +470.435 3.58102 +471.014 4.20584 +471.594 3.82446 +472.174 4.01001 +472.754 4.0153 +473.333 3.63918 +473.913 3.63918 +474.493 3.7738 +475.072 4.15464 +475.652 4.62907 +476.232 4.76668 +476.812 4.48129 +477.391 4.3343 +477.971 4.15464 +478.551 4.40447 +479.13 4.03304 +479.71 4.6701 +480.29 3.63918 +480.87 3.97357 +481.449 4.08505 +482.029 4.43678 +482.609 4.96319 +483.768 5.26193 +484.348 4.14128 +484.928 3.89457 +485.507 3.67146 +486.087 3.74943 +486.667 3.78844 +487.246 4.15464 +487.826 4.74566 +488.406 5.22458 +488.986 5.04468 +489.565 4.92783 +490.145 4.46774 +490.725 4.48954 +491.304 4.6701 +491.884 4.15464 +492.464 4.48992 +493.043 4.46141 +493.623 4.11361 +494.203 4.75747 +494.783 5.18557 +495.362 4.68171 +495.942 4.97688 +496.522 4.88704 +497.101 5.03132 +497.681 4.74983 +498.261 4.11686 +498.841 3.80077 +499.42 4.14303 +500 4.4704 +500.58 4.4704 +501.159 4.74983 +501.739 4.73507 +502.319 4.55532 +502.899 4.6701 +503.478 4.6701 +504.058 4.52202 +504.638 4.37158 +505.217 4.32878 +505.797 4.3626 +506.377 4.3355 +506.957 4.6474 +507.536 4.9776 +508.116 4.63527 +508.696 4.32878 +509.275 4.30889 +509.855 4.09821 +510.435 4.14303 +511.014 4.39948 +511.594 4.43692 +512.174 4.42526 +512.754 4.15464 +513.333 4.15464 +513.913 4.15464 +514.493 4.15464 +515.072 4.15464 +515.652 4.6474 +516.232 4.91495 +516.812 4.91858 +517.391 5.22604 +517.971 5.81129 +518.551 5.98805 +519.13 5.90899 +519.71 5.94227 +520.29 5.71319 +520.87 5.10842 +521.449 4.91134 +522.029 4.91495 +522.609 5.22994 +523.188 5.54678 +523.768 5.37438 +524.348 5.18557 +524.928 5.47259 +525.507 5.22905 +526.087 4.96863 +526.667 4.9776 +527.246 4.9776 +527.826 4.91495 +528.406 4.6701 +528.986 4.77982 +529.565 4.85892 +530.145 5.13676 +530.725 5.70103 +531.304 5.8819 +531.884 5.70103 +532.464 5.91797 +533.043 6.00853 +533.623 6.02768 +534.203 5.79761 +534.783 5.70103 +535.362 5.70103 +535.942 5.52046 +536.522 5.52046 +537.101 5.49307 +537.681 5.50835 +538.261 5.59985 +538.841 5.18557 +539.42 5.58015 +540 5.39353 +540.58 5.48409 +541.159 5.52046 +541.739 5.50133 +542.319 5.33982 +542.899 5.72174 +543.478 6.45252 +544.058 6.38324 +544.638 7.66951 +545.217 7.71453 +545.797 7.97983 +546.377 7.78141 +546.957 6.91994 +547.536 6.87478 +548.116 6.73196 +548.696 6.20689 +549.275 5.80221 +549.855 5.90899 +550.435 5.9797 +551.014 5.90899 +551.594 5.8386 +552.174 5.70103 +552.754 5.70103 +553.333 6.16443 +553.913 6.78316 +554.493 9.75775 +555.072 15.537 +555.652 34.2691 +556.232 39.5102 +556.812 59.2689 +557.391 64.412 +557.971 54.1549 +558.551 35.913 +559.13 28.9918 +559.71 12.2252 +560.29 8.45154 +560.87 6.73196 +561.449 6.93166 +562.029 7.00141 +562.609 6.67989 +563.188 6.37074 +563.768 6.50792 +564.348 6.50351 +564.928 6.75467 +565.507 7.06656 +566.087 7.00258 +566.667 7.24742 +567.246 7.54595 +567.826 7.45539 +568.406 7.49344 +568.986 7.60864 +569.565 7.27971 +570.145 6.73196 +570.725 6.88004 +571.304 7.01424 +571.884 6.98043 +572.464 7.08583 +573.043 6.73196 +573.623 7.3746 +574.203 7.19536 +574.783 7.44492 +575.362 7.74018 +575.942 8.02062 +576.522 8.25564 +577.101 8.56063 +577.681 8.51154 +578.261 8.59412 +578.841 8.79381 +579.42 7.88377 +580 7.72932 +580.58 7.52535 +581.159 7.642 +581.739 7.8212 +582.319 8.16863 +582.899 8.07865 +583.478 7.7306 +584.058 7.24742 +584.638 8.04125 +585.217 8.36343 +585.797 9.94521 +586.377 12.0639 +586.957 16.3176 +587.536 20.4668 +588.116 25.8064 +588.696 29.3807 +589.275 29.3102 +589.855 23.1 +590.435 21.595 +591.014 15.2585 +591.594 12.9361 +592.174 11.3938 +592.754 11.1909 +593.333 10.8764 +593.913 10.835 +594.493 10.4499 +595.072 10.6582 +595.652 10.3195 +596.232 9.97282 +596.812 10.0825 +597.391 9.76309 +597.971 9.95192 +598.551 9.29817 +599.13 8.94189 +599.71 9.14768 +600.29 8.8261 +600.87 8.38343 +601.449 8.51154 +602.029 8.58585 +602.609 8.50239 +603.188 8.605 +603.768 8.31318 +604.348 8.10868 +604.928 8.27835 +605.507 7.24742 +606.087 7.24742 +606.667 7.4401 +607.246 7.42829 +607.826 7.24742 +608.406 8.17327 +608.986 8.08954 +609.565 8.27835 +610.145 8.16863 +610.725 8.07039 +611.304 8.27835 +611.884 8.17327 +612.464 8.38186 +613.043 8.47586 +613.623 9.61242 +614.203 10.3052 +614.783 10.6869 +615.362 10.8557 +615.942 10.6134 +616.522 10.3186 +617.101 9.57928 +617.681 9.62084 +618.261 8.27835 +618.841 8.75033 +619.42 8.95541 +620 9.29817 +620.58 9.63593 +621.159 9.82474 +621.739 10.7321 +622.319 11.4407 +622.899 11.8267 +623.478 12.4137 +624.058 13.0447 +624.638 12.9175 +625.217 12.2959 +625.797 12.115 +626.377 12.2923 +626.957 14.9133 +627.536 21.1287 +628.116 36.6867 +628.696 41.3354 +629.275 57.6942 +629.855 59.1521 +630.435 51.2687 +631.014 34.9122 +631.594 27.8569 +632.174 14.5532 +632.754 12.1075 +633.333 11.2123 +633.913 13.7866 +634.493 19.0633 +635.072 20.1133 +635.652 25.8978 +636.232 26.0465 +636.812 20.1411 +637.391 19.3843 +637.971 13.1583 +638.551 10.2015 +639.13 8.32889 +639.71 8.13027 +640.29 8.27835 +640.87 8.27835 +641.449 8.27835 +642.029 8.17717 +642.609 7.8212 +643.188 7.71034 +643.768 8.63631 +644.348 9.31677 +644.928 11.0567 +645.507 10.9788 +646.087 11.8866 +646.667 11.3317 +647.246 11.1218 +647.826 10.6636 +648.406 10.9853 +648.986 11.3526 +649.565 11.8988 +650.145 11.5231 +650.725 11.1334 +651.304 10.6558 +651.884 10.3551 +652.464 10.9293 +653.043 11.3354 +653.623 11.9449 +654.203 12.3009 +654.783 12.2024 +655.362 12.5102 +655.942 12.3731 +656.522 11.8866 +657.101 11.0804 +657.681 10.654 +658.261 10.2695 +658.841 9.33821 +659.42 9.19449 +660 9.06443 +660.58 8.75806 +661.159 8.26107 +661.739 7.97085 +662.319 7.86796 +662.899 8.24606 +663.478 8.29906 +664.058 7.97085 +664.638 8.04516 +665.217 8.13027 +665.797 7.76289 +666.377 7.40167 +666.957 7.3525 +667.536 7.36831 +668.116 7.76289 +668.696 7.52535 +669.275 7.72511 +669.855 7.99133 +670.435 8.06141 +671.014 8.27835 +671.594 8.14078 +672.174 7.774 +672.754 7.43624 +673.333 7.45539 +673.913 7.4401 +674.493 7.45539 +675.072 7.65317 +675.652 7.50515 +676.232 7.35714 +676.812 7.55492 +677.391 7.56319 +677.971 7.53885 +678.551 7.55492 +679.13 7.55492 +679.71 7.81932 +680.29 9.65148 +680.87 11.3051 +681.449 14.5364 +682.029 16.5313 +682.609 18.8741 +683.188 18.4397 +683.768 17.8652 +684.348 16.3563 +684.928 17.4229 +685.507 17.5141 +686.087 18.0528 +686.667 16.3282 +687.246 15.7171 +687.826 13.8877 +688.406 12.725 +688.986 13.0041 +689.565 13.9203 +690.145 14.9206 +690.725 15.1415 +691.304 16.0103 +691.884 16.1584 +692.464 16.2809 +693.043 16.1719 +693.623 16.4814 +694.203 16.4533 +694.783 15.8885 +695.362 15.3143 +695.942 14.2308 +696.522 13.8519 +697.101 13.7062 +697.681 13.391 +698.261 12.6206 +698.841 12.5406 +699.42 12.178 +700 12.2478 +700.58 11.64 +701.159 10.9589 +701.739 9.8549 +702.319 9.38564 +702.899 9.14054 +703.478 9.77268 +704.058 9.98634 +704.638 10.2819 +705.217 9.97252 +705.797 9.92132 +706.377 9.67666 +706.957 9.64388 +707.536 9.51724 +708.116 9.52622 +708.696 9.24783 +709.275 8.69724 +709.855 8.63311 +710.435 8.24352 +711.014 8.48704 +711.594 8.75033 +712.174 8.95541 +712.754 9.25721 +713.333 9.55412 +713.913 9.30928 +714.493 9.30928 +715.072 9.5963 +715.652 9.82474 +716.232 9.82474 +716.812 8.79381 +717.391 9.05155 +717.971 8.86507 +718.551 8.22954 +719.13 8.29906 +719.71 8.63957 +720.29 8.87248 +720.87 9.3305 +721.449 15.637 +722.029 20.1401 +722.609 28.6894 +723.188 33.0478 +723.768 38.4977 +724.348 33.1741 +724.928 34.2699 +725.507 29.8928 +726.087 34.2373 +726.667 34.2889 +727.246 37.9426 +727.826 33.7338 +728.406 32.8606 +728.986 28.6109 +729.565 29.0202 +730.145 29.3213 +730.725 32.2093 +731.304 34.2558 +731.884 34.3695 +732.464 36.1769 +733.043 35.634 +733.623 33.8883 +734.203 32.9167 +734.783 31.33 +735.362 29.8979 +735.942 28.4604 +736.522 26.3706 +737.101 24.3426 +737.681 22.4484 +738.261 20.933 +738.841 21.3864 +739.42 21.8313 +740 19.1661 +740.58 18.5118 +741.159 16.6295 +741.739 15.6363 +742.319 15.4449 +742.899 15.8667 +743.478 15.4146 +744.058 15.6517 +744.638 15.9859 +745.217 18.4977 +745.797 20.7169 +746.377 24.0519 +746.957 25.8541 +747.536 23.9817 +748.116 21.1378 +748.696 17.8534 +749.275 16.2876 +749.855 17.0713 +750.435 21.7365 +751.014 23.9927 +751.594 28.8293 +752.174 27.5174 +752.754 25.3872 +753.333 21.0108 +753.913 18.8854 +754.493 18.5931 +755.072 20.9298 +755.652 22.2347 +756.232 25.2723 +756.812 25.7506 +757.391 23.2869 +757.971 21.6081 +758.551 18.7737 +759.13 16.1576 +759.71 14.8936 +760.29 14.111 +760.87 14.865 +761.449 15.7166 +762.029 14.9983 +762.609 14.3403 +763.188 12.9534 +763.768 12.9258 +764.348 14.4084 +764.928 14.856 +765.507 15.2049 +766.087 14.4824 +766.667 13.9169 +767.246 13.8262 +767.826 15.7917 +768.406 21.0755 +768.986 30.6862 +769.565 38.484 +770.145 45.6146 +770.725 53.9026 +771.304 48.2944 +771.884 48.9772 +772.464 40.7157 +773.043 42.2932 +773.623 44.9167 +774.203 53.3396 +774.783 63.4209 +775.362 67.1626 +775.942 75.1938 +776.522 76.6821 +777.101 71.6165 +777.681 61.3182 +778.261 59.5181 +778.841 49.1038 +779.42 46.4721 +780 45.7293 +780.58 46.372 +781.159 45.7281 +781.739 43.2748 +782.319 42.2241 +782.899 47.5645 +783.478 54.2501 +784.058 61.4242 +784.638 69.7246 +785.217 71.7094 +785.797 70.2685 +786.377 63.8373 +786.957 55.3266 +787.536 51.677 +788.116 45.8456 +788.696 49.3286 +789.275 54.1389 +789.855 65.8969 +790.435 71.46 +791.014 74.519 +791.594 61.2674 +792.174 57.2467 +792.754 44.4853 +793.333 38.2027 +793.913 41.1907 +794.493 43.6893 +795.072 49.8505 +795.652 54.6016 +796.232 57.742 +796.812 60.5992 +797.391 60.7104 +797.971 59.5874 +798.551 56.3547 +799.13 50.9674 +799.71 47.3366 +800.29 45.6159 +800.87 44.6832 +801.449 46.5385 +802.029 40.6538 +802.609 36.8755 +803.188 31.398 +803.768 27.8589 +804.348 27.4822 +804.928 29.3532 +805.507 27.1221 +806.087 25.1291 +806.667 21.7734 +807.246 19.052 +807.826 18.512 +808.406 19.7206 +808.986 20.5909 +809.565 18.743 +810.145 17.4667 +810.725 15.9447 +811.304 14.3856 +811.884 14.662 +812.464 15.2577 +813.043 15.8025 +813.623 15.4659 +814.203 13.8491 +814.783 14.0153 +815.362 12.3152 +815.942 13.0334 +816.522 14.2521 +817.101 16.1205 +817.681 16.899 +818.261 17.4205 +818.841 17.0197 +819.42 15.5314 +820 14.7387 +820.58 13.8677 +821.159 13.2426 +821.739 12.9175 +822.319 12.9765 +822.899 12.9175 +823.478 13.5005 +824.058 14.6433 +824.638 19.5446 +825.217 28.5358 +825.797 42.6672 +826.377 60.0186 +826.957 77.8219 +827.536 95.9089 +828.116 102.099 +828.696 94.9218 +829.275 78.75 +829.855 60.2716 +830.435 51.8195 +831.014 36.7166 +831.594 63.471 +832.174 70.6965 +832.754 102.346 +833.333 116.413 +833.913 115.91 +834.493 97.2797 +835.072 81.3286 +835.652 61.6605 +836.232 49.5733 +836.812 61.2955 +837.391 62.4773 +837.971 78.0635 +838.551 89.0819 +839.13 92.8436 +839.71 89.4163 +840.29 94.8599 +840.87 91.0374 +841.449 102.8 +842.029 105.877 +842.609 100.947 +843.188 91.777 +843.768 94.1602 +844.348 87.2941 +844.928 86.0332 +845.507 75.8283 +846.087 63.7297 +846.667 50.1956 +847.246 45.8875 +847.826 43.4738 +848.406 49.5997 +848.986 51.9194 +849.565 47.2826 +850.145 36.4516 +850.725 32.3637 +851.304 23.4339 +851.884 25.1516 +852.464 26.8526 +853.043 30.2968 +853.623 30.16 +854.203 25.9995 +854.783 23.1355 +855.362 18.692 +855.942 16.025 +856.522 17.0069 +857.101 19.5337 +857.681 20.0135 +858.261 23.1367 +858.841 25.7826 +859.42 32.514 +860 38.8799 +860.58 50.0294 +861.159 54.5849 +861.739 62.4711 +862.319 62.7262 +862.899 62.1884 +863.478 63.3225 +864.058 65.827 +864.638 65.7299 +865.217 63.0209 +865.797 54.0712 +866.377 48.4893 +866.957 38.977 +867.536 31.6818 +868.116 26.2303 +868.696 21.9726 +869.275 18.8798 +869.855 16.7428 +870.435 15.7515 +871.014 14.231 +871.594 14.9983 +872.174 20.936 +872.754 36.5435 +873.333 61.783 +873.913 91.5534 +874.493 131.882 +875.072 147.322 +875.652 179.038 +876.232 154.566 +876.812 136.851 +877.391 107.779 +877.971 75.3209 +878.551 59.6519 +879.13 54.7789 +879.71 101.333 +880.29 117.704 +880.87 163.486 +881.449 191.7 +882.029 178.108 +882.609 159.788 +883.188 119.377 +883.768 94.9878 +884.348 70.1019 +884.928 79.9134 +885.507 91.7946 +886.087 109.664 +886.667 124.143 +887.246 136.776 +887.826 138.387 +888.406 135.886 +888.986 141.823 +889.565 147.058 +890.145 152.041 +890.725 146.2 +891.304 132.171 +891.884 113.843 +892.464 108.223 +893.043 109.324 +893.623 120.206 +894.203 122.098 +894.783 109.378 +895.362 86.926 +895.942 81.9209 +896.522 59.7124 +897.101 75.2968 +897.681 71.3721 +898.261 84.7366 +898.841 72.9625 +899.42 62.1695 +900 47.2853 +900.58 38.2968 +901.159 33.6638 +901.739 40.6661 +902.319 46.503 +902.899 52.222 +903.478 50.4708 +904.058 42.7641 +904.638 36.7373 +905.217 27.1263 +905.797 22.6297 +906.377 21.1298 +906.957 23.8146 +907.536 26.2566 +908.116 28.9092 +908.696 27.8603 +909.275 24.2104 +909.855 20.8028 +910.435 15.8617 +911.014 14.0861 +911.594 13.2357 +912.174 15.1197 +912.754 17.5077 +913.333 19.4647 +913.913 21.4817 +914.493 20.4455 +915.072 18.9619 +915.652 16.1543 +916.232 14.1858 +916.812 11.485 +917.391 10.9527 +917.971 12.0169 +918.551 12.8542 +919.13 15.4948 +919.71 17.5039 +920.29 18.912 +920.87 19.3187 +921.449 17.9194 +922.029 17.1219 +922.609 19.6963 +923.188 25.7384 +923.768 29.4421 +924.348 35.7747 +924.928 37.1728 +925.507 35.0882 +926.087 29.8179 +926.667 36.1395 +927.246 57.3186 +927.826 106.151 +928.406 143.517 +928.986 189.79 +929.565 227.23 +930.145 218.653 +930.725 200.627 +931.304 149.621 +931.884 110.815 +932.464 75.7787 +933.043 52.4871 +933.623 67.6589 +934.203 97.4867 +934.783 162.024 +935.362 188.029 +935.942 240.881 +936.522 206.897 +937.101 173.414 +937.681 141.731 +938.261 103.435 +938.841 93.5804 +939.42 105.045 +940 126.395 +940.58 163.206 +941.159 183.921 +941.739 220.847 +942.319 241.179 +942.899 244.362 +943.478 214.486 +944.058 227.481 +944.638 216.387 +945.217 232.365 +945.797 231.24 +946.377 212.689 +946.957 181.278 +947.536 172.93 +948.116 145.411 +948.696 154.141 +949.275 147.864 +949.855 143.195 +950.435 124.966 +#951.014 105.634 +#951.594 91.6268 +#952.174 84.7906 +#952.754 104.656 +#953.333 121.646 +#953.913 147.06 +#954.493 161.154 +#955.072 153.611 +#955.652 117.388 +#956.232 98.9875 +#956.812 66.1366 +#957.391 47.8399 Index: branches/Mars_McMismatchStudy/sequence.def =================================================================== --- branches/Mars_McMismatchStudy/sequence.def (revision 17962) +++ branches/Mars_McMismatchStudy/sequence.def (revision 17962) @@ -0,0 +1,1 @@ +Telescope: 1 Index: branches/Mars_McMismatchStudy/sql.rc =================================================================== --- branches/Mars_McMismatchStudy/sql.rc (revision 17962) +++ branches/Mars_McMismatchStudy/sql.rc (revision 17962) @@ -0,0 +1,4 @@ +URL: mysql://lp-fact +User: operator +Password: Ics+eMM! +Database: factdata