Index: /trunk/Mars/fact/analysis/mc/callisto.C =================================================================== --- /trunk/Mars/fact/analysis/mc/callisto.C (revision 17148) +++ /trunk/Mars/fact/analysis/mc/callisto.C (revision 17148) @@ -0,0 +1,670 @@ + +#include +#include + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.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" + +/* 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. + */ + +using namespace std; + +int callisto_for_monte_carlo_simulated_data( + const char *drs_calib_300_path="/fhgfs/groups/app/fact/mc_test/testingdrsfile/test300samples.drs.fits", + const char *pedestal_file_path="/fhgfs/groups/app/fact/mc_test/ceresfitstest/mcFilesForTests/mcNsbPedestal/00000001.001_P_MonteCarlo000_Events.fits", + const char *data_file_path="/fhgfs/groups/app/fact/mc_test/ceresfitstest/mcFilesForTests/mcProton/00000003.387_D_MonteCarlo010_Events.fits", + + const char *root_file_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto.root", + const char *status_display_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto_status_display.root", + const char *status_display_title = "callisto_status_display") +{ + + // ====================================================== + + // 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 *map = usemap ? "/fhgfs/groups/app/fact/resources/monte_carlo_FACTmap.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 = "/cm/shared/apps/fact/Mars_svn_LP/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 (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 1; + } + + TString datfile = TString(data_file_path); + TString drsfile = TString(drs_calib_300_path); + TString pedfile = TString(pedestal_file_path); + + gLog.Separator("Callisto"); + gLog << all; + gLog << "Data File : " << datfile << "\n"; + gLog << "DRS calib 300: " << drsfile << '\n'; + + 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 << "root_file_output_path: " << root_file_output_path << endl; + gLog << "status_display_output_path: " << status_display_output_path << endl; + gLog << "status_display_title: " << status_display_title << 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(map); + read3.AddFile(pedfile); + + MFillH fill3a(&hrate); + + MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); + + MGeomApply apply3; + + MDrsCalibApply drsapply3; + + //--- + + 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(map); + 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 read5; + read5.LoadMap(map); + read5.AddFile(datfile); + + MFillH fill5a(&hrate); + + MGeomApply apply5; + + MDrsCalibApply drsapply5; + + 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(root_file_output_path, "RECREATE", "Calibrated Data", 2); + write5.AddContainer("MRawRunHeader", "RunHeaders"); + write5.AddContainer("MGeomCam", "RunHeaders"); + write5.AddContainer("MSignalCam", "Events"); + write5.AddContainer("MTime", "Events"); + 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(&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(status_display_title, kFALSE); + d->SaveAs(status_display_output_path); + + return 0; +} Index: /trunk/Mars/fact/analysis/mc/star.C =================================================================== --- /trunk/Mars/fact/analysis/mc/star.C (revision 17148) +++ /trunk/Mars/fact/analysis/mc/star.C (revision 17148) @@ -0,0 +1,339 @@ +#include +#include + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.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 "MMuonSetup.h" +#include "MReadMarsFile.h" +#include "MHillasCalc.h" +#include "MHn.h" +#include "MMuonSearchParCalc.h" +#include "MMuonCalibParCalc.h" +#include "MBinning.h" +#include "MImgCleanStd.h" + + +using namespace std; + +int star_for_monte_carlo_simulated_data( + const char *mars_data_file_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto.root", + const char *root_output_file_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/star.root", + const char *status_display_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/star_status_display.root", + const char *status_display_title = "star_status_display_title", + const char *hillas_txt_file_path = "/fhgfs/groups/app/callisto_star_test/star_hillas.txt", + Double_t lvl1=4.0, + Double_t lvl2=2.5, + Double_t deltat = 17.5) +{ + //Double_t lvl1=7.8, + //Double_t lvl2=3.9, + // lvl1 = 2.5; + // lvl2 = 0.5; + // ------------------------------------------------------ + + gLog.Separator("Star"); + gLog << all << "Calculate image parameters of sequence "; + gLog << endl; + gLog << "mars_data_file_path: " << mars_data_file_path << endl; + gLog << "root_output_file_path: " << root_output_file_path << endl; + gLog << "status_display_output_path: " << status_display_output_path << endl; + gLog << "status_display_title: " << status_display_title << 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(mars_data_file_path); + + MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData"); + + MGeomApply apply; + + MImgCleanStd clean(lvl1, lvl2); + clean.SetMethod(MImgCleanStd::kAbsolute); + clean.SetTimeLvl1(deltat);//17.5/*1.3*/); + clean.SetTimeLvl2(deltat);//17.5/*1.3*/); + clean.SetPostCleanType(3); + + 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 write_ascii_hillas(hillas_txt_file_path, "MHillas"); + write_ascii_hillas.AddColumns("MHillasExt"); + write_ascii_hillas.AddColumns("MHillasSrc"); + + + MWriteRootFile write(root_output_file_path, "RECREATE", "Image parameters", 2); + 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); + tlist.AddToList(&write_ascii_hillas); + + if (!loop.Eventloop()) + return 3; + + if (!loop.GetDisplay()) + return 4; + + d->SetTitle(status_display_title, kFALSE); + d->SaveAs(status_display_output_path); + + return 0; +} Index: /trunk/Mars/fact/analysis/sandbox_dneise/callisto_buildable_no_sequence_file.C =================================================================== --- /trunk/Mars/fact/analysis/sandbox_dneise/callisto_buildable_no_sequence_file.C (revision 17148) +++ /trunk/Mars/fact/analysis/sandbox_dneise/callisto_buildable_no_sequence_file.C (revision 17148) @@ -0,0 +1,1069 @@ + +#include +#include + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.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" + +/* 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. + */ + +using namespace std; + +int callisto_buildable_no_sequence_file( + const char *drs_calib_300_path="/path/to/drs_calibration_300.drs.fits.gz", + const char *drs_calib_1024_path="/path/to/drs_calibration_1024.drs.fits.gz", + const char *time_calib_path="/path/to/time_calibration.fits.gz", + const char *pedestal_file_path="/path/to/pedestal_file.fits.gz", + const char *ext_lighpulser_file_path="/path/to/ext_lightpulser_file.fits.gz", + const char *data_file_path="/path/to/data_file.fits.gz", + const char *title_for_output="use yyyymmdd_rrr this seems to be default", + const char *outpath = "outpath", + const char *output_fname = "output_fname", + bool use_delays=true) +{ + + gLog << endl; + gLog.Separator("Callisto"); + + /* const char *seqfile: + * relativer Pfad zu einem SeqFile, welches über die zu analysierenden Daten informiert + * + * ostringstream drsname: (104) + * relativer Pfad zu einem sogenannten "Drs Sequence File" mit Endung ".drs.seq" + * Er wird aus dem Pfad zum vorigen SeqFile sowie: + * seq.GetNight().GetNightAsInt() und + * seq.GetDrsSequence() + * gebildet. D.h. in dem vorigen SeqFile muss es eine Info names "DrsSequence" geben + * sowie "neben" dem vorigen SeqFile muss ein DrsSeqFile liegen. + * + * Die beiden MSequence instancen heissen: + * seq(97) und + * drs(117) + * + * Wofür werden diese beiden MSequence Instanzen denn nun eigentlich + * benötigt? + * ----------------------------------------------------------------------- + * TString drsfile = seq.GetFileName(0, MSequence::kRawDrs); + * TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); + * TString calfile = seq.GetFileName(0, MSequence::kFitsCal); + * TString timfile = drs.GetFileName(0, MSequence::kFitsDat); + * TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs); + * Ausserdem wird offenbar ein MDirIter mit Hilfe einer MSequence instanz + * erzeugt: + * MDirIter iter; + * seq.GetRuns(iter, MSequence::kFitsDat); + * Ich vermute das liesst sich so: + * Mache einen leeren MDirIter namens iter + * Hey MSequence, teile dem iter mal deine "FitsDat" artigen files + * mit. Oder so...... + * ----------------------------------------------------------------------- + * + * Der MDirIter iter, der schon in Zeile 179 erzeugt wird, wird dann in + * Zeile 756 folgendermaßen benutzt: + * MRawFitsRead read5; + * read5.LoadMap(map); ---> map=/home/fact/FACT++/FACTmap111030.txt" + * read5.AddFiles(iter); + * + * Zur Info: Fast alle anderen MRawFitsRead Instanzen, die benutzt werden, + * werden so ähnlich verwendet. Nur wird statt AddFiles(MDirIter iter) + * AddFile(TString file_path) gerufen. Alle ausser MRawFitsRead read0. + * Das wird so erzeugt: + * MRawFitsRead read0(timfile); + * Da wird also *scheinbar* map=/home/fact/FACT++/FACTmap111030.txt" + * nicht benutzt. Bedeutet das, dass die DRS4 Zeitkalibrationskonstanten + * nicht für die selbe Pixel ID bestimmt werden, wie die, die im Rest + * dieses macros?? + * + * Eine Frage ist, ob eine MRawFitsRead Instanz einen vollständigen Pfad + * (relativ oder absolut) braucht, oder nur einen Filename? Die Frage ergibt + * sich, da die TStrings, die "geadded" werden so erzeugt werden: + * + * TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); + * MRawFitsRead read3; + * read3.LoadMap(map); + * read3.AddFile(pedfile); + * Da wird also explizit GetFile*Name* und nicht GetFile*Path* gerufen. + * Ich hab mich mal ein wenig in MRawFitsRead umgesehen und keine Stelle + * gefunden an der Pfade zusammengesetzt werden, also returned GetFileName + * bereits vollständige Pfade, gut so! + * ----------------------------------------------------------------------- + * ----------------------------------------------------------------------- + * Da wir Callisto ja immer nur für 1 daten file rufen wollen benötigen + * wir den MDirIter garnicht. + * Und da wir gerne explizite file names angeben, die verarbeitet werden + * sollen und unser Callisto.C sowieso per (Python-)Skript rufen macht + * es uns auch garnix aus alle benötigten Pfade sofort explizit beim + * macro call zu übergeben. + * + * Daher finde ich nun alle Pfade und übergebe sie explizit beim macro + * call. + */ + + + + // ====================================================== + + // 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 *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"; + + // ------------------------------------------------------ + + // 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 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; + } + + TString datfile = TString(data_file_path); + TString drsfile = TString(drs_calib_300_path); + TString drs1024 = TString(drs_calib_1024_path); + TString timfile = TString(time_calib_path); + TString pedfile = TString(pedestal_file_path); + TString calfile = TString(ext_lighpulser_file_path); + TString title_part = TString(title_for_output); + + gLog << all; + gLog << "Data File : " << datfile << "\n"; + gLog << "DRS calib 300: " << drsfile << '\n'; + gLog << "DRS calib 1024: " << drs1024 << "\n\n"; + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile.Data())) { + gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl; + return 5; + } + MDrsCalibration drscalib1024; + if (!drscalib1024.ReadFits(drs1024.Data())) { + gLog << err << "ERROR - Cannot access drscallib1024 file '" << drs1024 << "'" << endl; + return 6; + } + gLog << all; + gLog << "Time calibration : " << timfile << '\n'; + gLog << "Pedestal file: " << pedfile << '\n'; + gLog << "Light Pulser file: " << calfile << '\n' << endl; + gLog << "Output will be witten to:" << outpath << 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"); + + Bool_t isinteg = + (type&MExtralgoSpline::kIntegral) || + (type&MExtralgoSpline::kFixedWidth) || + (type&MExtralgoSpline::kDynWidth) + ? kTRUE : kFALSE; + + gStyle->SetOptFit(kTRUE); + + // ====================================================== + + gLog << endl; + gLog.Separator("Processing DRS timing calibration run"); + + MTaskList tlist0; + + MParList plist0; + plist0.AddToList(&tlist0); + plist0.AddToList(&drscalib1024); + plist0.AddToList(&badpixels); + plist0.AddToList(&timecam); + + MEvtLoop loop0("DetermineTimeCal"); + loop0.SetDisplay(d); + loop0.SetParList(&plist0); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read0(timfile); + + MContinue cont0("MRawEvtHeader.GetTriggerID!=33792", "SelectTim"); + + MGeomApply apply0; + + MDrsCalibApply drsapply0; + + MFillH fill0("MHDrsCalibrationTime"); + fill0.SetNameTab("DeltaT"); + + tlist0.AddToList(&read0); + tlist0.AddToList(&apply0); + tlist0.AddToList(&drsapply0); + tlist0.AddToList(&cont0); + tlist0.AddToList(&fill0); + + if (!loop0.Eventloop(max0)) + return 8; + + if (!loop0.GetDisplay()) + return 9; + + /* + MHDrsCalibrationT *t = (MHDrsCalibrationT*)plist4.FindObject("MHDrsCalibrationT"); + t->SetDisplay(d); + t->PlotAll(); + */ + + // ====================================================== + + gLog << endl; + gLog.Separator("Processing external light pulser run"); + + MTaskList tlist1; + + MParList plist1; + plist1.AddToList(&tlist1); + plist1.AddToList(&drscalib300); + plist1.AddToList(&badpixels); + plist1.AddToList(&timecam); + + MEvtLoop loop1("DetermineCalConst"); + loop1.SetDisplay(d); + loop1.SetParList(&plist1); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read1; + read1.LoadMap(map); + read1.AddFile(calfile); + + MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal"); + + MGeomApply apply1; + + MDrsCalibApply drsapply1; + + /* + MPedestalCam fPedestalCamOut1a; + MPedestalCam fPedestalCamOut1b; + MPedestalCam fPedestalCamOut1c; + + MPedCalcPedRun pedcalc1a; + MPedCalcPedRun pedcalc1b; + MPedCalcPedRun pedcalc1c; + pedcalc1a.SetPedestalsOut(&fPedestalCamOut1a); + pedcalc1b.SetPedestalsOut(&fPedestalCamOut1b); + pedcalc1c.SetPedestalsOut(&fPedestalCamOut1c); + + MExtractTimeAndChargeSpline extractor1ab; + extractor1a.SetRange(first_slice, last_slice); + extractor1a.SetRiseTimeHiGain(rise_time); + extractor1a.SetFallTimeHiGain(fall_time); + extractor1a.SetChargeType(type); + extractor1a.SetSaturationLimit(600000); + extractor1a.SetNoiseCalculation(kFALSE); + + pedcalc1a.SetRandomCalculation(kTRUE); + pedcalc1b.SetRandomCalculation(kFALSE); + pedcalc1a.SetExtractor(&extractor1a); + pedcalc1b.SetExtractor($extractor1a); + pedcalc1c.SetRangeFromExtractor(&extractor1a); + */ + + // --- + + MExtractTimeAndChargeSpline 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); + + // --- + + MHCamEvent evt1f(0, "ExtCalSig", "Extracted calibration signal;;S [mV·sl]"); + MHCamEvent evt1g(4, "ExtCalTm", "Extracted arrival times;;T [sl]"); + MHCamEvent evt1h(6, "CalCalTm", "Calibrated arrival times;;T [sl]"); + + MHSectorVsTime hist1rmsb("ExtSigVsTm"); + MHSectorVsTime hist1tmb("CalTmVsTm"); + hist1rmsb.SetTitle("Extracted calibration vs event number;;S [mV·sl]"); + hist1rmsb.SetType(0); + hist1tmb.SetTitle("Extracted arrival time vs event number;;T [sl]"); + //hist1tmb.SetType(4); + hist1tmb.SetType(6); + + MFillH fill1f(&evt1f, "MExtractedSignalCam", "FillExtSig"); + MFillH fill1g(&evt1g, "MArrivalTimeCam", "FillExtTm"); + MFillH fill1h(&evt1h, "MSignalCam", "FillCalTm"); + MFillH fill1r(&hist1rmsb, "MExtractedSignalCam", "FillExtSigVsTm"); + //MFillH fill1j(&hist1tmb, "MArrivalTimeCam", "FillExtTmVsTm"); + MFillH fill1j(&hist1tmb, "MSignalCam", "FillCalTmVsTm"); + + fill1f.SetDrawOption("gaus"); + fill1h.SetDrawOption("gaus"); + + // --- + + MCalibrateDrsTimes calctm1a("CalibrateCalEvents"); + calctm1a.SetNameUncalibrated("UncalibratedTimes"); + + MBadPixelsTreat treat1; + treat1.SetProcessPedestalRun(kFALSE); + treat1.SetProcessPedestalEvt(kFALSE); + + // --- + + MHCamEvent evt1c(6, "ExtCalTmShift", "Relative extracted arrival time of calibration pulse (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt1d(6, "CalCalTmShift", "Relative calibrated arrival time of calibration pulse (w.r.t. event-median);;\\Delta T [ns]"); + + evt1c.SetMedianShift(); + evt1d.SetMedianShift(); + + MFillH fill1c(&evt1c, "UncalibratedTimes", "FillExtCalTm"); + MFillH fill1d(&evt1d, "MSignalCam", "FillCalCalTm"); + fill1d.SetDrawOption("gaus"); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist1.AddToList(&read1); + tlist1.AddToList(&apply1); + tlist1.AddToList(&drsapply1); + tlist1.AddToList(&cont1); + tlist1.AddToList(&extractor1b); + if (isinteg) + { + tlist1.AddToList(&extractor1c); + tlist1.AddToList(&fill1a); + } + tlist1.AddToList(&calctm1a); + tlist1.AddToList(&treat1); + tlist1.AddToList(&fill1f); + tlist1.AddToList(&fill1g); + tlist1.AddToList(&fill1h); + tlist1.AddToList(&fill1r); + tlist1.AddToList(&fill1j); + tlist1.AddToList(&fill1c); + tlist1.AddToList(&fill1d); + + if (!loop1.Eventloop(max1)) + return 10; + + if (!loop1.GetDisplay()) + return 11; + + if (use_delays) + 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(); + + // ====================================================== + + 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(map); + read3.AddFile(pedfile); + + MFillH fill3a(&hrate); + + MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); + + MGeomApply apply3; + + MDrsCalibApply drsapply3; + + //--- + + 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(map); + 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 read5; + read5.LoadMap(map); + read5.AddFile(datfile); + + MFillH fill5a(&hrate); + + MGeomApply apply5; + + MDrsCalibApply drsapply5; + + 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); + + MHCamEvent evt5m(6, "ExtTm", "Extracted arrival times of calibration pulse;;\\Delta T [ns]"); + MHCamEvent evt5n(6, "CalTm", "Calibrated arrival times of calibration pulse;;\\Delta T [ns]"); + MHCamEvent evt5q(6, "ExtTmShift", "Relative extracted arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt5r(6, "CalTmShift", "Relative calibrated arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt5s(6, "ExtTM", "Extracted absolute time marker position;;T [sl]"); + MHCamEvent evt5t(6, "CalTM", "Calibrated absolute time marker position;;T [ns]"); + MHCamEvent evt5u(6, "ExtTMshift", "Relative extracted time marker position (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt5v(6, "CalTMshift", "Relative calibrated time marker position (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt5w(6, "ExtDiff", "Difference between extracted arrival time of time marker and calibration pulse;;\\Delta T [ns]"); + MHCamEvent evt5x(6, "CalDiff", "Difference between calibrated arrival time of time marker and calibration pulse;;\\Delta T [ns]"); + + evt5w.SetNameSub("UncalibratedTimes"); + evt5x.SetNameSub("MSignalCam"); + + evt5q.SetMedianShift(); + evt5r.SetMedianShift(); + evt5u.SetMedianShift(); + evt5v.SetMedianShift(); + //evt4w.SetMedianShift(); + //evt4x.SetMedianShift(); + + MFillH fill5m(&evt5m, "UncalibratedTimes", "FillExtTm"); + MFillH fill5n(&evt5n, "MSignalCam", "FillCalTm"); + MFillH fill5q(&evt5q, "UncalibratedTimes", "FillExtTmShift"); + MFillH fill5r(&evt5r, "MSignalCam" , "FillCalTmShift"); + MFillH fill5s(&evt5s, "UncalTimeMarker", "FillExtTM"); + MFillH fill5t(&evt5t, "TimeMarker", "FillCalTM"); + MFillH fill5u(&evt5u, "UncalTimeMarker", "FillExtTMshift"); + MFillH fill5v(&evt5v, "TimeMarker", "FillCalTMshift"); + MFillH fill5w(&evt5w, "UncalTimeMarker", "FillExtDiff"); + MFillH fill5x(&evt5x, "TimeMarker", "FillCalDiff"); + + fill5m.SetDrawOption("gaus"); + fill5n.SetDrawOption("gaus"); + fill5q.SetDrawOption("gaus"); + fill5r.SetDrawOption("gaus"); + //fill5s.SetDrawOption("gaus"); + //fill5t.SetDrawOption("gaus"); + //fill5u.SetDrawOption("gaus"); + //fill5v.SetDrawOption("gaus"); + //fill5w.SetDrawOption("gaus"); + //fill5x.SetDrawOption("gaus"); + + + MBadPixelsTreat treat5; + treat5.SetProcessPedestalRun(kFALSE); + treat5.SetProcessPedestalEvt(kFALSE); + + MHSectorVsTime hist5cal("CalVsTm"); + MHSectorVsTime hist5ped("PedVsTm"); + hist5cal.SetTitle("Median calibrated calibration signal vs event number;;Signal [~phe]"); + hist5ped.SetTitle("Median calibrated pedestal signal vs event number;;Signal [~phe]"); + hist5cal.SetType(0); + hist5ped.SetType(0); + hist5cal.SetMinimum(0); + hist5ped.SetMinimum(0); + hist5cal.SetUseMedian(); + hist5ped.SetUseMedian(); + hist5cal.SetNameTime("MTime"); + hist5ped.SetNameTime("MTime"); + + MFillH fill5cal(&hist5cal, "MSignalCam", "FillCalVsTm"); + MFillH fill5ped(&hist5ped, "MSignalCam", "FillPedVsTm"); + fill5cal.SetFilter(&filtercal); + fill5ped.SetFilter(&filterped); + + 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(); + + const TString fname(Form("s/([0-9]+_[0-9]+)[.]fits([.]gz)?$/%s\\/$1_C.root/", + MJob::Esc(outpath).Data())); + + // The second rule is for the case reading raw-files! + MWriteRootFile write5(2, fname, "RECREATE", "Calibrated Data"); + write5.AddContainer("MRawRunHeader", "RunHeaders"); + write5.AddContainer("MGeomCam", "RunHeaders"); + write5.AddContainer("MSignalCam", "Events"); + write5.AddContainer("MTime", "Events"); + 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.AddToList(&fill5m); + tlist5tm.AddToList(&fill5n); + tlist5tm.AddToList(&fill5q); + tlist5tm.AddToList(&fill5r); + //tlist5tm.AddToList(&fill5s); + //tlist5tm.AddToList(&fill5t); + tlist5tm.AddToList(&fill5u); + tlist5tm.AddToList(&fill5v); + tlist5tm.AddToList(&fill5w); + tlist5tm.AddToList(&fill5x); + 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(&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(&fill5ped); + tlist5.AddToList(&fill5cal); + tlist5.AddToList(&tlist5dat); + tlist5.AddToList(&write5); + + if (!loop5.Eventloop(max4)) + return 18; + + if (!loop5.GetDisplay()) + return 19; + + TString title = "-- Calibrated signal #"; + title += title_part; + title += " ("; + title += drsfile; + title += ") --"; + d->SetTitle(title, kFALSE); + + TString path; + path += Form("%s/%s", outpath, output_fname); + + d->SaveAs(path); + + return 0; +} Index: /trunk/Mars/fact/analysis/sandbox_dneise/callisto_buildable_no_sequence_file_orig.C =================================================================== --- /trunk/Mars/fact/analysis/sandbox_dneise/callisto_buildable_no_sequence_file_orig.C (revision 17148) +++ /trunk/Mars/fact/analysis/sandbox_dneise/callisto_buildable_no_sequence_file_orig.C (revision 17148) @@ -0,0 +1,1081 @@ + +#include +#include + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.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" + +/* 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. + */ + +using namespace std; + +int callisto_buildable_no_sequence_file( + const char *drs_calib_300_path="/fhgfs/groups/app/fact/raw/2013/08/01/20130801_007.drs.fits.gz", + const char *drs_calib_1024_path="/fhgfs/groups/app/fact/raw/2013/08/01/20130801_004.drs.fits.gz", + const char *time_calib_path="/fhgfs/groups/app/fact/raw/2013/08/01/20130801_006.fits.gz", + const char *pedestal_file_path="/fhgfs/groups/app/fact/raw/2013/08/01/20130801_010.fits.gz", + const char *ext_lighpulser_file_path="/fhgfs/groups/app/fact/raw/2013/08/01/20130801_011.fits.gz", + const char *data_file_path="/fhgfs/groups/app/fact/raw/2013/08/01/20130801_012.fits.gz", + + const char *root_file_output_path = "/fhgfs/groups/app/callisto_star_test/callisto.root", + const char *raw_run_header_txt_file_path = "/fhgfs/groups/app/callisto_star_test/raw_run_header.txt", + const char *raw_event_header_txt_file_path = "/fhgfs/groups/app/callisto_star_test/raw_event_header.txt", + const char *geom_cam_txt_file_path = "/fhgfs/groups/app/callisto_star_test/geom_cam.txt", + const char *signal_cam_txt_file_path = "/fhgfs/groups/app/callisto_star_test/signal_cam.txt", + const char *time_txt_file_path = "/fhgfs/groups/app/callisto_star_test/time.txt", + const char *status_display_output_path= "/fhgfs/groups/app/callisto_star_test/callisto_status_display.root", + const char *status_display_title = "callisto_status_display", + + bool use_delays=true) +{ + + /* const char *seqfile: + * relativer Pfad zu einem SeqFile, welches über die zu analysierenden Daten informiert + * + * ostringstream drsname: (104) + * relativer Pfad zu einem sogenannten "Drs Sequence File" mit Endung ".drs.seq" + * Er wird aus dem Pfad zum vorigen SeqFile sowie: + * seq.GetNight().GetNightAsInt() und + * seq.GetDrsSequence() + * gebildet. D.h. in dem vorigen SeqFile muss es eine Info names "DrsSequence" geben + * sowie "neben" dem vorigen SeqFile muss ein DrsSeqFile liegen. + * + * Die beiden MSequence instancen heissen: + * seq(97) und + * drs(117) + * + * Wofür werden diese beiden MSequence Instanzen denn nun eigentlich + * benötigt? + * ----------------------------------------------------------------------- + * TString drsfile = seq.GetFileName(0, MSequence::kRawDrs); + * TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); + * TString calfile = seq.GetFileName(0, MSequence::kFitsCal); + * TString timfile = drs.GetFileName(0, MSequence::kFitsDat); + * TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs); + * Ausserdem wird offenbar ein MDirIter mit Hilfe einer MSequence instanz + * erzeugt: + * MDirIter iter; + * seq.GetRuns(iter, MSequence::kFitsDat); + * Ich vermute das liesst sich so: + * Mache einen leeren MDirIter namens iter + * Hey MSequence, teile dem iter mal deine "FitsDat" artigen files + * mit. Oder so...... + * ----------------------------------------------------------------------- + * + * Der MDirIter iter, der schon in Zeile 179 erzeugt wird, wird dann in + * Zeile 756 folgendermaßen benutzt: + * MRawFitsRead read5; + * read5.LoadMap(map); ---> map=/home/fact/FACT++/FACTmap111030.txt" + * read5.AddFiles(iter); + * + * Zur Info: Fast alle anderen MRawFitsRead Instanzen, die benutzt werden, + * werden so ähnlich verwendet. Nur wird statt AddFiles(MDirIter iter) + * AddFile(TString file_path) gerufen. Alle ausser MRawFitsRead read0. + * Das wird so erzeugt: + * MRawFitsRead read0(timfile); + * Da wird also *scheinbar* map=/home/fact/FACT++/FACTmap111030.txt" + * nicht benutzt. Bedeutet das, dass die DRS4 Zeitkalibrationskonstanten + * nicht für die selbe Pixel ID bestimmt werden, wie die, die im Rest + * dieses macros?? + * + * Eine Frage ist, ob eine MRawFitsRead Instanz einen vollständigen Pfad + * (relativ oder absolut) braucht, oder nur einen Filename? Die Frage ergibt + * sich, da die TStrings, die "geadded" werden so erzeugt werden: + * + * TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); + * MRawFitsRead read3; + * read3.LoadMap(map); + * read3.AddFile(pedfile); + * Da wird also explizit GetFile*Name* und nicht GetFile*Path* gerufen. + * Ich hab mich mal ein wenig in MRawFitsRead umgesehen und keine Stelle + * gefunden an der Pfade zusammengesetzt werden, also returned GetFileName + * bereits vollständige Pfade, gut so! + * ----------------------------------------------------------------------- + * ----------------------------------------------------------------------- + * Da wir Callisto ja immer nur für 1 daten file rufen wollen benötigen + * wir den MDirIter garnicht. + * Und da wir gerne explizite file names angeben, die verarbeitet werden + * sollen und unser Callisto.C sowieso per (Python-)Skript rufen macht + * es uns auch garnix aus alle benötigten Pfade sofort explizit beim + * macro call zu übergeben. + * + * Daher finde ich nun alle Pfade und übergebe sie explizit beim macro + * call. + */ + + + + // ====================================================== + + // 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 *map = usemap ? "/fhgfs/groups/app/fact/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 *pulse_template = "/cm/shared/apps/fact/Mars_svn_LP/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 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; + } + + TString datfile = TString(data_file_path); + TString drsfile = TString(drs_calib_300_path); + TString drs1024 = TString(drs_calib_1024_path); + TString timfile = TString(time_calib_path); + TString pedfile = TString(pedestal_file_path); + TString calfile = TString(ext_lighpulser_file_path); + + gLog.Separator("Callisto"); + gLog << all; + gLog << "Data File : " << datfile << "\n"; + gLog << "DRS calib 300: " << drsfile << '\n'; + gLog << "DRS calib 1024: " << drs1024 << "\n\n"; + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile.Data())) { + gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl; + return 5; + } + MDrsCalibration drscalib1024; + if (!drscalib1024.ReadFits(drs1024.Data())) { + gLog << err << "ERROR - Cannot access drscallib1024 file '" << drs1024 << "'" << endl; + return 6; + } + gLog << all; + gLog << "Time calibration : " << timfile << '\n'; + gLog << "Pedestal file: " << pedfile << '\n'; + gLog << "Light Pulser file: " << calfile << '\n' << endl; + + + gLog << "root_file_output_path: " << root_file_output_path << endl; + gLog << "raw_run_header_txt_file_path: " << raw_run_header_txt_file_path << endl; + gLog << "raw_event_header_txt_file_path: " << raw_event_header_txt_file_path << endl; + gLog << "geom_cam_txt_file_path: " << geom_cam_txt_file_path << endl; + gLog << "signal_cam_txt_file_path: " << signal_cam_txt_file_path << endl; + gLog << "time_txt_file_path: " << time_txt_file_path << endl; + gLog << "status_display_output_path: " << status_display_output_path << endl; + gLog << "status_display_title: " << status_display_title << 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"); + + Bool_t isinteg = + (type&MExtralgoSpline::kIntegral) || + (type&MExtralgoSpline::kFixedWidth) || + (type&MExtralgoSpline::kDynWidth) + ? kTRUE : kFALSE; + + gStyle->SetOptFit(kTRUE); + + // ====================================================== + + gLog << endl; + gLog.Separator("Processing DRS timing calibration run"); + + MTaskList tlist0; + + MParList plist0; + plist0.AddToList(&tlist0); + plist0.AddToList(&drscalib1024); + plist0.AddToList(&badpixels); + plist0.AddToList(&timecam); + + MEvtLoop loop0("DetermineTimeCal"); + loop0.SetDisplay(d); + loop0.SetParList(&plist0); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read0(timfile); + + MContinue cont0("MRawEvtHeader.GetTriggerID!=33792", "SelectTim"); + + MGeomApply apply0; + + MDrsCalibApply drsapply0; + + MFillH fill0("MHDrsCalibrationTime"); + fill0.SetNameTab("DeltaT"); + + tlist0.AddToList(&read0); + tlist0.AddToList(&apply0); + tlist0.AddToList(&drsapply0); + tlist0.AddToList(&cont0); + tlist0.AddToList(&fill0); + + if (!loop0.Eventloop(max0)) + return 8; + + if (!loop0.GetDisplay()) + return 9; + + /* + MHDrsCalibrationT *t = (MHDrsCalibrationT*)plist4.FindObject("MHDrsCalibrationT"); + t->SetDisplay(d); + t->PlotAll(); + */ + + // ====================================================== + + gLog << endl; + gLog.Separator("Processing external light pulser run"); + + MTaskList tlist1; + + MParList plist1; + plist1.AddToList(&tlist1); + plist1.AddToList(&drscalib300); + plist1.AddToList(&badpixels); + plist1.AddToList(&timecam); + + MEvtLoop loop1("DetermineCalConst"); + loop1.SetDisplay(d); + loop1.SetParList(&plist1); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read1; + read1.LoadMap(map); + read1.AddFile(calfile); + + MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal"); + + MGeomApply apply1; + + MDrsCalibApply drsapply1; + + /* + MPedestalCam fPedestalCamOut1a; + MPedestalCam fPedestalCamOut1b; + MPedestalCam fPedestalCamOut1c; + + MPedCalcPedRun pedcalc1a; + MPedCalcPedRun pedcalc1b; + MPedCalcPedRun pedcalc1c; + pedcalc1a.SetPedestalsOut(&fPedestalCamOut1a); + pedcalc1b.SetPedestalsOut(&fPedestalCamOut1b); + pedcalc1c.SetPedestalsOut(&fPedestalCamOut1c); + + MExtractTimeAndChargeSpline extractor1ab; + extractor1a.SetRange(first_slice, last_slice); + extractor1a.SetRiseTimeHiGain(rise_time); + extractor1a.SetFallTimeHiGain(fall_time); + extractor1a.SetChargeType(type); + extractor1a.SetSaturationLimit(600000); + extractor1a.SetNoiseCalculation(kFALSE); + + pedcalc1a.SetRandomCalculation(kTRUE); + pedcalc1b.SetRandomCalculation(kFALSE); + pedcalc1a.SetExtractor(&extractor1a); + pedcalc1b.SetExtractor($extractor1a); + pedcalc1c.SetRangeFromExtractor(&extractor1a); + */ + + // --- + + MExtractTimeAndChargeSpline 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); + + // --- + + MHCamEvent evt1f(0, "ExtCalSig", "Extracted calibration signal;;S [mV·sl]"); + MHCamEvent evt1g(4, "ExtCalTm", "Extracted arrival times;;T [sl]"); + MHCamEvent evt1h(6, "CalCalTm", "Calibrated arrival times;;T [sl]"); + + MHSectorVsTime hist1rmsb("ExtSigVsTm"); + MHSectorVsTime hist1tmb("CalTmVsTm"); + hist1rmsb.SetTitle("Extracted calibration vs event number;;S [mV·sl]"); + hist1rmsb.SetType(0); + hist1tmb.SetTitle("Extracted arrival time vs event number;;T [sl]"); + //hist1tmb.SetType(4); + hist1tmb.SetType(6); + + MFillH fill1f(&evt1f, "MExtractedSignalCam", "FillExtSig"); + MFillH fill1g(&evt1g, "MArrivalTimeCam", "FillExtTm"); + MFillH fill1h(&evt1h, "MSignalCam", "FillCalTm"); + MFillH fill1r(&hist1rmsb, "MExtractedSignalCam", "FillExtSigVsTm"); + //MFillH fill1j(&hist1tmb, "MArrivalTimeCam", "FillExtTmVsTm"); + MFillH fill1j(&hist1tmb, "MSignalCam", "FillCalTmVsTm"); + + fill1f.SetDrawOption("gaus"); + fill1h.SetDrawOption("gaus"); + + // --- + + MCalibrateDrsTimes calctm1a("CalibrateCalEvents"); + calctm1a.SetNameUncalibrated("UncalibratedTimes"); + + MBadPixelsTreat treat1; + treat1.SetProcessPedestalRun(kFALSE); + treat1.SetProcessPedestalEvt(kFALSE); + + // --- + + MHCamEvent evt1c(6, "ExtCalTmShift", "Relative extracted arrival time of calibration pulse (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt1d(6, "CalCalTmShift", "Relative calibrated arrival time of calibration pulse (w.r.t. event-median);;\\Delta T [ns]"); + + evt1c.SetMedianShift(); + evt1d.SetMedianShift(); + + MFillH fill1c(&evt1c, "UncalibratedTimes", "FillExtCalTm"); + MFillH fill1d(&evt1d, "MSignalCam", "FillCalCalTm"); + fill1d.SetDrawOption("gaus"); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist1.AddToList(&read1); + tlist1.AddToList(&apply1); + tlist1.AddToList(&drsapply1); + tlist1.AddToList(&cont1); + tlist1.AddToList(&extractor1b); + if (isinteg) + { + tlist1.AddToList(&extractor1c); + tlist1.AddToList(&fill1a); + } + tlist1.AddToList(&calctm1a); + tlist1.AddToList(&treat1); + tlist1.AddToList(&fill1f); + tlist1.AddToList(&fill1g); + tlist1.AddToList(&fill1h); + tlist1.AddToList(&fill1r); + tlist1.AddToList(&fill1j); + tlist1.AddToList(&fill1c); + tlist1.AddToList(&fill1d); + + if (!loop1.Eventloop(max1)) + return 10; + + if (!loop1.GetDisplay()) + return 11; + + if (use_delays) + 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(); + + // ====================================================== + + 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(map); + read3.AddFile(pedfile); + + MFillH fill3a(&hrate); + + MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); + + MGeomApply apply3; + + MDrsCalibApply drsapply3; + + //--- + + 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(map); + 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 read5; + read5.LoadMap(map); + read5.AddFile(datfile); + + MFillH fill5a(&hrate); + + MGeomApply apply5; + + MDrsCalibApply drsapply5; + + 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); + + MHCamEvent evt5m(6, "ExtTm", "Extracted arrival times of calibration pulse;;\\Delta T [ns]"); + MHCamEvent evt5n(6, "CalTm", "Calibrated arrival times of calibration pulse;;\\Delta T [ns]"); + MHCamEvent evt5q(6, "ExtTmShift", "Relative extracted arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt5r(6, "CalTmShift", "Relative calibrated arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt5s(6, "ExtTM", "Extracted absolute time marker position;;T [sl]"); + MHCamEvent evt5t(6, "CalTM", "Calibrated absolute time marker position;;T [ns]"); + MHCamEvent evt5u(6, "ExtTMshift", "Relative extracted time marker position (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt5v(6, "CalTMshift", "Relative calibrated time marker position (w.r.t. event-median);;\\Delta T [ns]"); + MHCamEvent evt5w(6, "ExtDiff", "Difference between extracted arrival time of time marker and calibration pulse;;\\Delta T [ns]"); + MHCamEvent evt5x(6, "CalDiff", "Difference between calibrated arrival time of time marker and calibration pulse;;\\Delta T [ns]"); + + evt5w.SetNameSub("UncalibratedTimes"); + evt5x.SetNameSub("MSignalCam"); + + evt5q.SetMedianShift(); + evt5r.SetMedianShift(); + evt5u.SetMedianShift(); + evt5v.SetMedianShift(); + //evt4w.SetMedianShift(); + //evt4x.SetMedianShift(); + + MFillH fill5m(&evt5m, "UncalibratedTimes", "FillExtTm"); + MFillH fill5n(&evt5n, "MSignalCam", "FillCalTm"); + MFillH fill5q(&evt5q, "UncalibratedTimes", "FillExtTmShift"); + MFillH fill5r(&evt5r, "MSignalCam" , "FillCalTmShift"); + MFillH fill5s(&evt5s, "UncalTimeMarker", "FillExtTM"); + MFillH fill5t(&evt5t, "TimeMarker", "FillCalTM"); + MFillH fill5u(&evt5u, "UncalTimeMarker", "FillExtTMshift"); + MFillH fill5v(&evt5v, "TimeMarker", "FillCalTMshift"); + MFillH fill5w(&evt5w, "UncalTimeMarker", "FillExtDiff"); + MFillH fill5x(&evt5x, "TimeMarker", "FillCalDiff"); + + fill5m.SetDrawOption("gaus"); + fill5n.SetDrawOption("gaus"); + fill5q.SetDrawOption("gaus"); + fill5r.SetDrawOption("gaus"); + //fill5s.SetDrawOption("gaus"); + //fill5t.SetDrawOption("gaus"); + //fill5u.SetDrawOption("gaus"); + //fill5v.SetDrawOption("gaus"); + //fill5w.SetDrawOption("gaus"); + //fill5x.SetDrawOption("gaus"); + + + MBadPixelsTreat treat5; + treat5.SetProcessPedestalRun(kFALSE); + treat5.SetProcessPedestalEvt(kFALSE); + + MHSectorVsTime hist5cal("CalVsTm"); + MHSectorVsTime hist5ped("PedVsTm"); + hist5cal.SetTitle("Median calibrated calibration signal vs event number;;Signal [~phe]"); + hist5ped.SetTitle("Median calibrated pedestal signal vs event number;;Signal [~phe]"); + hist5cal.SetType(0); + hist5ped.SetType(0); + hist5cal.SetMinimum(0); + hist5ped.SetMinimum(0); + hist5cal.SetUseMedian(); + hist5ped.SetUseMedian(); + hist5cal.SetNameTime("MTime"); + hist5ped.SetNameTime("MTime"); + + MFillH fill5cal(&hist5cal, "MSignalCam", "FillCalVsTm"); + MFillH fill5ped(&hist5ped, "MSignalCam", "FillPedVsTm"); + fill5cal.SetFilter(&filtercal); + fill5ped.SetFilter(&filterped); + + 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(root_file_output_path, "RECREATE", "Calibrated Data", 2); + write5.AddContainer("MRawRunHeader", "RunHeaders"); + write5.AddContainer("MGeomCam", "RunHeaders"); + write5.AddContainer("MSignalCam", "Events"); + write5.AddContainer("MTime", "Events"); + write5.AddContainer("MRawEvtHeader", "Events"); + //write.AddContainer("MTriggerPattern", "Events"); + + MWriteAsciiFile raw_run_header_writer_task(raw_run_header_txt_file_path, "MRawRunHeader"); + MWriteAsciiFile raw_event_header_writer_task(raw_event_header_txt_file_path, "MRawEvtHeader"); + MWriteAsciiFile geom_cam_writer_task(geom_cam_txt_file_path, "MGeomCam"); + MWriteAsciiFile signal_cam_writer_task(signal_cam_txt_file_path, "MSignalCam"); + MWriteAsciiFile time_writer_task(time_txt_file_path, "MTime"); + + // ------------------ Setup histograms and fill tasks ---------------- + + MContinue test; + test.SetFilter(&filterncl); + + MTaskList tlist5tm; + tlist5tm.AddToList(&extractor5tm); + tlist5tm.AddToList(&calctm5tm); + tlist5tm.AddToList(&fill5m); + tlist5tm.AddToList(&fill5n); + tlist5tm.AddToList(&fill5q); + tlist5tm.AddToList(&fill5r); + //tlist5tm.AddToList(&fill5s); + //tlist5tm.AddToList(&fill5t); + tlist5tm.AddToList(&fill5u); + tlist5tm.AddToList(&fill5v); + tlist5tm.AddToList(&fill5w); + tlist5tm.AddToList(&fill5x); + 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(&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(&fill5ped); + tlist5.AddToList(&fill5cal); + tlist5.AddToList(&tlist5dat); + tlist5.AddToList(&write5); + + tlist5.AddToList(&raw_run_header_writer_task); + tlist5.AddToList(&raw_event_header_writer_task); + tlist5.AddToList(&geom_cam_writer_task); + tlist5.AddToList(&signal_cam_writer_task); + tlist5.AddToList(&time_writer_task); + + if (!loop5.Eventloop(max4)) + return 18; + + if (!loop5.GetDisplay()) + return 19; + + d->SetTitle(status_display_title, kFALSE); + d->SaveAs(status_display_output_path); + + return 0; +} Index: /trunk/Mars/fact/analysis/sandbox_dneise/callisto_for_monte_carlo_simulated_data.C =================================================================== --- /trunk/Mars/fact/analysis/sandbox_dneise/callisto_for_monte_carlo_simulated_data.C (revision 17148) +++ /trunk/Mars/fact/analysis/sandbox_dneise/callisto_for_monte_carlo_simulated_data.C (revision 17148) @@ -0,0 +1,746 @@ + +#include +#include + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.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" + +/* 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. + */ + +using namespace std; + +int callisto_for_monte_carlo_simulated_data( + const char *drs_calib_300_path="/fhgfs/groups/app/fact/mc_test/testingdrsfile/test300samples.drs.fits", + const char *pedestal_file_path="/fhgfs/groups/app/fact/mc_test/ceresfitstest/mcFilesForTests/mcNsbPedestal/00000001.001_P_MonteCarlo000_Events.fits", + const char *data_file_path="/fhgfs/groups/app/fact/mc_test/ceresfitstest/mcFilesForTests/mcProton/00000003.387_D_MonteCarlo010_Events.fits", + + const char *root_file_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto.root", + const char *status_display_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto_status_display.root", + const char *status_display_title = "callisto_status_display") +{ + + /* const char *seqfile: + * relativer Pfad zu einem SeqFile, welches über die zu analysierenden Daten informiert + * + * ostringstream drsname: (104) + * relativer Pfad zu einem sogenannten "Drs Sequence File" mit Endung ".drs.seq" + * Er wird aus dem Pfad zum vorigen SeqFile sowie: + * seq.GetNight().GetNightAsInt() und + * seq.GetDrsSequence() + * gebildet. D.h. in dem vorigen SeqFile muss es eine Info names "DrsSequence" geben + * sowie "neben" dem vorigen SeqFile muss ein DrsSeqFile liegen. + * + * Die beiden MSequence instancen heissen: + * seq(97) und + * drs(117) + * + * Wofür werden diese beiden MSequence Instanzen denn nun eigentlich + * benötigt? + * ----------------------------------------------------------------------- + * TString drsfile = seq.GetFileName(0, MSequence::kRawDrs); + * TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); + * TString calfile = seq.GetFileName(0, MSequence::kFitsCal); + * TString timfile = drs.GetFileName(0, MSequence::kFitsDat); + * TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs); + * Ausserdem wird offenbar ein MDirIter mit Hilfe einer MSequence instanz + * erzeugt: + * MDirIter iter; + * seq.GetRuns(iter, MSequence::kFitsDat); + * Ich vermute das liesst sich so: + * Mache einen leeren MDirIter namens iter + * Hey MSequence, teile dem iter mal deine "FitsDat" artigen files + * mit. Oder so...... + * ----------------------------------------------------------------------- + * + * Der MDirIter iter, der schon in Zeile 179 erzeugt wird, wird dann in + * Zeile 756 folgendermaßen benutzt: + * MRawFitsRead read5; + * read5.LoadMap(map); ---> map=/home/fact/FACT++/FACTmap111030.txt" + * read5.AddFiles(iter); + * + * Zur Info: Fast alle anderen MRawFitsRead Instanzen, die benutzt werden, + * werden so ähnlich verwendet. Nur wird statt AddFiles(MDirIter iter) + * AddFile(TString file_path) gerufen. Alle ausser MRawFitsRead read0. + * Das wird so erzeugt: + * MRawFitsRead read0(timfile); + * Da wird also *scheinbar* map=/home/fact/FACT++/FACTmap111030.txt" + * nicht benutzt. Bedeutet das, dass die DRS4 Zeitkalibrationskonstanten + * nicht für die selbe Pixel ID bestimmt werden, wie die, die im Rest + * dieses macros?? + * + * Eine Frage ist, ob eine MRawFitsRead Instanz einen vollständigen Pfad + * (relativ oder absolut) braucht, oder nur einen Filename? Die Frage ergibt + * sich, da die TStrings, die "geadded" werden so erzeugt werden: + * + * TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); + * MRawFitsRead read3; + * read3.LoadMap(map); + * read3.AddFile(pedfile); + * Da wird also explizit GetFile*Name* und nicht GetFile*Path* gerufen. + * Ich hab mich mal ein wenig in MRawFitsRead umgesehen und keine Stelle + * gefunden an der Pfade zusammengesetzt werden, also returned GetFileName + * bereits vollständige Pfade, gut so! + * ----------------------------------------------------------------------- + * ----------------------------------------------------------------------- + * Da wir Callisto ja immer nur für 1 daten file rufen wollen benötigen + * wir den MDirIter garnicht. + * Und da wir gerne explizite file names angeben, die verarbeitet werden + * sollen und unser Callisto.C sowieso per (Python-)Skript rufen macht + * es uns auch garnix aus alle benötigten Pfade sofort explizit beim + * macro call zu übergeben. + * + * Daher finde ich nun alle Pfade und übergebe sie explizit beim macro + * call. + */ + + + + // ====================================================== + + // 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 *map = usemap ? "/fhgfs/groups/app/fact/resources/monte_carlo_FACTmap.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 = "/cm/shared/apps/fact/Mars_svn_LP/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 (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 1; + } + + TString datfile = TString(data_file_path); + TString drsfile = TString(drs_calib_300_path); + TString pedfile = TString(pedestal_file_path); + + gLog.Separator("Callisto"); + gLog << all; + gLog << "Data File : " << datfile << "\n"; + gLog << "DRS calib 300: " << drsfile << '\n'; + + 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 << "root_file_output_path: " << root_file_output_path << endl; + gLog << "status_display_output_path: " << status_display_output_path << endl; + gLog << "status_display_title: " << status_display_title << 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(map); + read3.AddFile(pedfile); + + MFillH fill3a(&hrate); + + MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); + + MGeomApply apply3; + + MDrsCalibApply drsapply3; + + //--- + + 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(map); + 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 read5; + read5.LoadMap(map); + read5.AddFile(datfile); + + MFillH fill5a(&hrate); + + MGeomApply apply5; + + MDrsCalibApply drsapply5; + + 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(root_file_output_path, "RECREATE", "Calibrated Data", 2); + write5.AddContainer("MRawRunHeader", "RunHeaders"); + write5.AddContainer("MGeomCam", "RunHeaders"); + write5.AddContainer("MSignalCam", "Events"); + write5.AddContainer("MTime", "Events"); + 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(&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(status_display_title, kFALSE); + d->SaveAs(status_display_output_path); + + return 0; +} Index: /trunk/Mars/fact/analysis/sandbox_dneise/star.C =================================================================== --- /trunk/Mars/fact/analysis/sandbox_dneise/star.C (revision 17148) +++ /trunk/Mars/fact/analysis/sandbox_dneise/star.C (revision 17148) @@ -0,0 +1,341 @@ +#include +#include + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.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 "MMuonSetup.h" +#include "MReadMarsFile.h" +#include "MHillasCalc.h" +#include "MHn.h" +#include "MMuonSearchParCalc.h" +#include "MMuonCalibParCalc.h" +#include "MBinning.h" +#include "MImgCleanStd.h" + + +using namespace std; + +int star( + Double_t lvl1=4.0, + Double_t lvl2=2.5, + //Double_t lvl1=7.8, + //Double_t lvl2=3.9, + const char *mars_data_file_path = "/fhgfs/groups/app/callisto_star_test/callisto.root", + const char *root_output_file_path = "/fhgfs/groups/app/callisto_star_test/star.root", + const char *status_display_output_path = "/fhgfs/groups/app/callisto_star_test/star_status_display.root", + const char *status_display_title = "star_status_display_title", + + const char *hillas_txt_file_path = "/fhgfs/groups/app/callisto_star_test/star_hillas.txt" + ) +{ + double deltat = 17.5; + // lvl1 = 2.5; + // lvl2 = 0.5; + // ------------------------------------------------------ + + gLog.Separator("Star"); + gLog << all << "Calculate image parameters of sequence "; + gLog << endl; + gLog << "mars_data_file_path: " << mars_data_file_path << endl; + gLog << "root_output_file_path: " << root_output_file_path << endl; + gLog << "status_display_output_path: " << status_display_output_path << endl; + gLog << "status_display_title: " << status_display_title << 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(mars_data_file_path); + + MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData"); + + MGeomApply apply; + + MImgCleanStd clean(lvl1, lvl2); + clean.SetMethod(MImgCleanStd::kAbsolute); + clean.SetTimeLvl1(deltat);//17.5/*1.3*/); + clean.SetTimeLvl2(deltat);//17.5/*1.3*/); + clean.SetPostCleanType(3); + + 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 write_ascii_hillas(hillas_txt_file_path, "MHillas"); + write_ascii_hillas.AddColumns("MHillasExt"); + write_ascii_hillas.AddColumns("MHillasSrc"); + + + MWriteRootFile write(root_output_file_path, "RECREATE", "Image parameters", 2); + 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); + tlist.AddToList(&write_ascii_hillas); + + if (!loop.Eventloop()) + return 3; + + if (!loop.GetDisplay()) + return 4; + + d->SetTitle(status_display_title, kFALSE); + d->SaveAs(status_display_output_path); + + return 0; +} Index: /trunk/Mars/fact/analysis/sandbox_dneise/star_for_monte_carlo_simulated_data.C =================================================================== --- /trunk/Mars/fact/analysis/sandbox_dneise/star_for_monte_carlo_simulated_data.C (revision 17148) +++ /trunk/Mars/fact/analysis/sandbox_dneise/star_for_monte_carlo_simulated_data.C (revision 17148) @@ -0,0 +1,339 @@ +#include +#include + +#include "TH1F.h" +#include "TFile.h" +#include "TStyle.h" +#include "TGraph.h" +#include "TLine.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 "MMuonSetup.h" +#include "MReadMarsFile.h" +#include "MHillasCalc.h" +#include "MHn.h" +#include "MMuonSearchParCalc.h" +#include "MMuonCalibParCalc.h" +#include "MBinning.h" +#include "MImgCleanStd.h" + + +using namespace std; + +int star_for_monte_carlo_simulated_data( + const char *mars_data_file_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto.root", + const char *root_output_file_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/star.root", + const char *status_display_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/star_status_display.root", + const char *status_display_title = "star_status_display_title", + const char *hillas_txt_file_path = "/fhgfs/groups/app/callisto_star_test/star_hillas.txt", + Double_t lvl1=4.0, + Double_t lvl2=2.5, + Double_t deltat = 17.5) +{ + //Double_t lvl1=7.8, + //Double_t lvl2=3.9, + // lvl1 = 2.5; + // lvl2 = 0.5; + // ------------------------------------------------------ + + gLog.Separator("Star"); + gLog << all << "Calculate image parameters of sequence "; + gLog << endl; + gLog << "mars_data_file_path: " << mars_data_file_path << endl; + gLog << "root_output_file_path: " << root_output_file_path << endl; + gLog << "status_display_output_path: " << status_display_output_path << endl; + gLog << "status_display_title: " << status_display_title << 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(mars_data_file_path); + + MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData"); + + MGeomApply apply; + + MImgCleanStd clean(lvl1, lvl2); + clean.SetMethod(MImgCleanStd::kAbsolute); + clean.SetTimeLvl1(deltat);//17.5/*1.3*/); + clean.SetTimeLvl2(deltat);//17.5/*1.3*/); + clean.SetPostCleanType(3); + + 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 write_ascii_hillas(hillas_txt_file_path, "MHillas"); + write_ascii_hillas.AddColumns("MHillasExt"); + write_ascii_hillas.AddColumns("MHillasSrc"); + + + MWriteRootFile write(root_output_file_path, "RECREATE", "Image parameters", 2); + 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); + tlist.AddToList(&write_ascii_hillas); + + if (!loop.Eventloop()) + return 3; + + if (!loop.GetDisplay()) + return 4; + + d->SetTitle(status_display_title, kFALSE); + d->SaveAs(status_display_output_path); + + return 0; +} Index: /trunk/Mars/macros/fact/fake_drs_calibration.C =================================================================== --- /trunk/Mars/macros/fact/fake_drs_calibration.C (revision 17148) +++ /trunk/Mars/macros/fact/fake_drs_calibration.C (revision 17148) @@ -0,0 +1,62 @@ +#include "DrsCalib.h" + + +using namespace std; + +// -------------------------------------------------------------------------- + +int fake_drs_calibration() +{ + DrsCalibration drscalib; + uint16_t fRoi = 1024; + const string filename = "/fhgfs/groups/app/fact/mc_test/testingdrsfile/test1024samples.drs.fits"; + uint16_t fNumTm = 0; + float baselineMean = 3720 * 2000.0/4096.0; + float baselineRms = 0; + float gainMean = 1907.35 * 16; + float gainRms = 0; + float triggerOffsetMean = 0; + float triggerOffsetRms = 0; + float triggerOffsetTmMean = 0; + float triggerOffsetTmRms = 0; + + const size_t n = 1440*1024*4 + 1440*fRoi*2 + fNumTm*fRoi*2 + 3; + vector vec; + vec.resize(n); + + reinterpret_cast(vec.data())[0] = fRoi; + reinterpret_cast(vec.data())[1] = 1; + reinterpret_cast(vec.data())[2] = 2; + reinterpret_cast(vec.data())[3] = 3; + + for (int i = 0 ; i < 1024*1440 ; i++) + { + vec.data()[i+4] = baselineMean; + vec.data()[i+4+1*1024*1440] = baselineRms; + vec.data()[i+4+2*1024*1440] = gainMean; + vec.data()[i+4+3*1024*1440] = gainRms; + } + for (int i = 0 ; i < fRoi*1440 ; i++) + { + vec.data()[i+4+4*1024*1440] = triggerOffsetMean; + vec.data()[i+4+4*1024*1440+1*fRoi*1440] = triggerOffsetRms; + } + for (int i = 0 ; i < fRoi*fNumTm ; i++) + { + vec.data()[i+4+2*1024*1440+2*fRoi*1440] = triggerOffsetTmMean; + vec.data()[i+4+3*1024*1440+2*fRoi*1440+fRoi*fNumTm] = triggerOffsetTmRms; + } + + drscalib.fRoi = fRoi; + drscalib.fNumTm = fNumTm; + drscalib.fStep = 2; + + drscalib.fNumOffset = 1000; + drscalib.fNumGain = 1000*1953125; + drscalib.fNumTrgOff = 1000; + + drscalib.WriteFitsImp(filename,vec); + + return 0; +} + Index: /trunk/Mars/mextralgo/MExtralgoSpline.cc =================================================================== --- /trunk/Mars/mextralgo/MExtralgoSpline.cc (revision 17147) +++ /trunk/Mars/mextralgo/MExtralgoSpline.cc (revision 17148) @@ -207,5 +207,5 @@ // at x, until x==0. If y is not found or out of range -2 is returned. // -Double_t MExtralgoSpline::SearchYdn(Float_t x, Float_t y) const +Double_t MExtralgoSpline::SearchYdn(Double_t x, Double_t y) const { if (x>=fNum-1) @@ -228,5 +228,5 @@ // at x, until x==fNum-1. If y is not found or out of range -2 is returned. // -Double_t MExtralgoSpline::SearchYup(Float_t x, Float_t y) const +Double_t MExtralgoSpline::SearchYup(Double_t x, Double_t y) const { if (x<0) Index: /trunk/Mars/mextralgo/MExtralgoSpline.h =================================================================== --- /trunk/Mars/mextralgo/MExtralgoSpline.h (revision 17147) +++ /trunk/Mars/mextralgo/MExtralgoSpline.h (revision 17148) @@ -347,6 +347,6 @@ Float_t Deriv1(const Float_t x) const; - Double_t SearchYdn(Float_t maxpos, Float_t y) const; - Double_t SearchYup(Float_t maxpos, Float_t y) const; + Double_t SearchYdn(Double_t maxpos, Double_t y) const; + Double_t SearchYup(Double_t maxpos, Double_t y) const; Double_t SearchYdn(Float_t y) const { return SearchYdn(fNum, y); } Index: /trunk/Mars/mfileio/MWriteFitsFile.cc =================================================================== --- /trunk/Mars/mfileio/MWriteFitsFile.cc (revision 17147) +++ /trunk/Mars/mfileio/MWriteFitsFile.cc (revision 17148) @@ -15,4 +15,6 @@ #include "MWriteFitsFile.h" #include "MFitsArray.h" + +#include "MTime.h" //for find @@ -337,4 +339,81 @@ return kTRUE; } + +void MWriteFitsFile::SetupHeaderKeys(MRawRunHeader &header) +{ + const MTime now(-1); + SetHeaderKey("ISMC",true,"Bool if File is Montecarlo File"); + SetHeaderKey("TELESCOP", "FACT", ""); + SetHeaderKey("PACKAGE", "MARS Cheobs", ""); + SetHeaderKey("VERSION", "1.0", ""); + SetHeaderKey("CREATOR", "Ceres", ""); + SetHeaderKey("EXTREL", 1., ""); + SetHeaderKey("COMPILED", __DATE__" "__TIME__, ""); + //SetHeaderKey("REVISION", "0", ""); + SetHeaderKey("ORIGIN", "FACT", ""); + SetHeaderKey("DATE", now.GetStringFmt("%Y-%m-%dT%H:%M:%S").Data(), ""); + SetHeaderKey("NIGHT", now.GetNightAsInt(), ""); + SetHeaderKey("TIMESYS", "UTC", ""); + SetHeaderKey("TIMEUNIT", "d", ""); + SetHeaderKey("MJDREF", 40587, ""); + //SetHeaderKey("BLDVER", 1, ""); + SetHeaderKey("RUNID", header.GetRunNumber(), ""); + SetHeaderKey("NBOARD", 40, ""); + SetHeaderKey("NPIX", header.GetNumPixel(), ""); + SetHeaderKey("NROI", header.GetNumSamplesHiGain(), ""); + SetHeaderKey("NROITM", 0, ""); + SetHeaderKey("TMSHIFT", 0, ""); + SetHeaderKey("CAMERA", "MGeomCamFACT", "Montecarlo File"); + SetHeaderKey("DAQ", "DRS4", "Montecarlo File"); + + // FTemme: ADCRANGE and ADC have to be calculated, using the values for + // the fadctype. +// SetHeaderKey("ADCRANGE", 2000, "Dynamic range in mV"); +// SetHeaderKey("ADC", 12, "Resolution in bits"); + + switch(header.GetRunType()) + { + case MRawRunHeader::kRTData|MRawRunHeader::kRTMonteCarlo: + SetHeaderKey("RUNTYPE", "data", ""); + break; + case MRawRunHeader::kRTPedestal|MRawRunHeader::kRTMonteCarlo: + SetHeaderKey("RUNTYPE", "pedestal", ""); + break; + case MRawRunHeader::kRTCalibration|MRawRunHeader::kRTMonteCarlo: + SetHeaderKey("RUNTYPE", "calibration", ""); + break; + } +// SetHeaderKey("ID", 777, "Board 0: Board ID"); +// SetHeaderKey("FMVER", 532, "Board 0: Firmware Version"); +// SetHeaderKey("DNA", "0", ""); +// SetHeaderKey("BOARD", 0, ""); +// SetHeaderKey("PRESC", 40, ""); +// SetHeaderKey("PHASE", 0, ""); +// SetHeaderKey("DAC0", 26500, ""); +// SetHeaderKey("DAC1", 0, ""); +// SetHeaderKey("DAC2", 0, ""); +// SetHeaderKey("DAC3", 0, ""); +// SetHeaderKey("DAC4", 28800, ""); +// SetHeaderKey("DAC5", 28800, ""); +// SetHeaderKey("DAC6", 28800, ""); +// SetHeaderKey("DAC7", 28800, ""); + SetHeaderKey("REFCLK", header.GetFreqSampling(), ""); + SetHeaderKey("DRSCALIB", false, ""); +// SetHeaderKey("TSTARTI", 0, ""); +// SetHeaderKey("TSTARTF", 0., ""); +// SetHeaderKey("TSTOPI", 0, ""); +// SetHeaderKey("TSTOPF", 0., ""); +// SetHeaderKey("DATE-OBS", "1970-01-01T00:00:00", ""); +// SetHeaderKey("DATE-END", "1970-01-01T00:00:00", ""); +// SetHeaderKey("NTRG", 0, ""); +// SetHeaderKey("NTRGPED", 0, ""); +// SetHeaderKey("NTRGLPE", 0, ""); +// SetHeaderKey("NTRGTIM", 0, ""); +// SetHeaderKey("NTRGLPI", 0, ""); +// SetHeaderKey("NTRGEXT1", 0, ""); +// SetHeaderKey("NTRGEXT2", 0, ""); +// SetHeaderKey("NTRGMISC", 0, ""); +} + template<> std::string MWriteFitsFile::GetFitsString(const double& value) @@ -1127,4 +1206,8 @@ return kFALSE; + + MRawRunHeader* header = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + SetupHeaderKeys(*header); + if (GetContainer(pList) == kFALSE) return kFALSE; Index: /trunk/Mars/mfileio/MWriteFitsFile.h =================================================================== --- /trunk/Mars/mfileio/MWriteFitsFile.h (revision 17147) +++ /trunk/Mars/mfileio/MWriteFitsFile.h (revision 17148) @@ -19,4 +19,10 @@ #include "MFitsArray.h" #endif + +#ifndef MARS_MRawRunHeader +#include "MRawRunHeader.h" +#endif + +class MRawRunHeader; /////////////////////////////////////////////////////////////////////////////// @@ -137,4 +143,6 @@ Int_t PreProcess(MParList *pList); Int_t PostProcess(); + + void SetupHeaderKeys(MRawRunHeader &header); //Header keys related stuff Index: /trunk/Mars/mjobs/MJSimulation.cc =================================================================== --- /trunk/Mars/mjobs/MJSimulation.cc (revision 17147) +++ /trunk/Mars/mjobs/MJSimulation.cc (revision 17148) @@ -260,4 +260,92 @@ write.AddContainer("IncidentAngle", "Events", kFALSE); } + +void MJSimulation::SetupHeaderKeys(MWriteFitsFile &write,MRawRunHeader &header) const +{ + const MTime now(-1); + write.SetHeaderKey("ISMC",true,"Bool if File is Montecarlo File"); + write.SetHeaderKey("TELESCOP", "FACT", ""); + write.SetHeaderKey("PACKAGE", "MARS Cheobs", ""); + write.SetHeaderKey("VERSION", "1.0", ""); + write.SetHeaderKey("CREATOR", "Ceres", ""); + write.SetHeaderKey("EXTREL", 1., ""); + write.SetHeaderKey("COMPILED", __DATE__" "__TIME__, ""); + write.SetHeaderKey("REVISION", "0", ""); + write.SetHeaderKey("ORIGIN", "FACT", ""); + write.SetHeaderKey("DATE", now.GetStringFmt("%Y-%m-%dT%H:%M:%S").Data(), ""); + write.SetHeaderKey("NIGHT", now.GetNightAsInt(), ""); + write.SetHeaderKey("TIMESYS", "UTC", ""); + write.SetHeaderKey("TIMEUNIT", "d", ""); + write.SetHeaderKey("MJDREF", 40587, ""); + //write.SetHeaderKey("BLDVER", 1, ""); + write.SetHeaderKey("RUNID", header.GetRunNumber(), ""); + write.SetHeaderKey("NBOARD", 40, ""); + write.SetHeaderKey("NPIX", header.GetNumPixel(), ""); + write.SetHeaderKey("NROI", header.GetNumSamplesHiGain(), ""); + write.SetHeaderKey("NROITM", 0, ""); + write.SetHeaderKey("TMSHIFT", 0, ""); + write.SetHeaderKey("CAMERA", "MGeomCamFACT", ""); + write.SetHeaderKey("DAQ", "DRS4", ""); + + // FTemme: ADCRANGE and ADC have to be calculated, using the values for + // the fadctype. +// write.SetHeaderKey("ADCRANGE", 2000, "Dynamic range in mV"); +// write.SetHeaderKey("ADC", 12, "Resolution in bits"); + + switch(header.GetRunType()) + { + case MRawRunHeader::kRTData|MRawRunHeader::kRTMonteCarlo: + write.SetHeaderKey("RUNTYPE", "data", ""); + break; + case MRawRunHeader::kRTPedestal|MRawRunHeader::kRTMonteCarlo: + write.SetHeaderKey("RUNTYPE", "pedestal", ""); + break; + case MRawRunHeader::kRTCalibration|MRawRunHeader::kRTMonteCarlo: + write.SetHeaderKey("RUNTYPE", "calibration", ""); + break; + } +// write.SetHeaderKey("ID", 777, "Board 0: Board ID"); +// write.SetHeaderKey("FMVER", 532, "Board 0: Firmware Version"); +// write.SetHeaderKey("DNA", "0", ""); +// write.SetHeaderKey("BOARD", 0, ""); +// write.SetHeaderKey("PRESC", 40, ""); +// write.SetHeaderKey("PHASE", 0, ""); +// write.SetHeaderKey("DAC0", 26500, ""); +// write.SetHeaderKey("DAC1", 0, ""); +// write.SetHeaderKey("DAC2", 0, ""); +// write.SetHeaderKey("DAC3", 0, ""); +// write.SetHeaderKey("DAC4", 28800, ""); +// write.SetHeaderKey("DAC5", 28800, ""); +// write.SetHeaderKey("DAC6", 28800, ""); +// write.SetHeaderKey("DAC7", 28800, ""); + write.SetHeaderKey("REFCLK", header.GetFreqSampling(), ""); + write.SetHeaderKey("DRSCALIB", false, ""); +// write.SetHeaderKey("TSTARTI", 0, ""); +// write.SetHeaderKey("TSTARTF", 0., ""); +// write.SetHeaderKey("TSTOPI", 0, ""); +// write.SetHeaderKey("TSTOPF", 0., ""); +// write.SetHeaderKey("DATE-OBS", "1970-01-01T00:00:00", ""); +// write.SetHeaderKey("DATE-END", "1970-01-01T00:00:00", ""); +// write.SetHeaderKey("NTRG", 0, ""); +// write.SetHeaderKey("NTRGPED", 0, ""); +// write.SetHeaderKey("NTRGLPE", 0, ""); +// write.SetHeaderKey("NTRGTIM", 0, ""); +// write.SetHeaderKey("NTRGLPI", 0, ""); +// write.SetHeaderKey("NTRGEXT1", 0, ""); +// write.SetHeaderKey("NTRGEXT2", 0, ""); +// write.SetHeaderKey("NTRGMISC", 0, ""); +} + +void MJSimulation::SetupVetoColumns(MWriteFitsFile &write) const +{ + write.VetoColumn("MParameterD.fVal"); + write.VetoColumn("MRawEvtData.fLoGainPixId"); + write.VetoColumn("MRawEvtData.fLoGainFadcSamples"); + write.VetoColumn("MRawEvtData.fABFlags"); + write.VetoColumn("MRawEvtHeader.fNumTrigLvl2"); + //write.VetoColumn("MRawEvtHeader.fTrigPattern"); + write.VetoColumn("MRawEvtHeader.fNumLoGainOn"); +} + Bool_t MJSimulation::Process(const MArgs &args, const MSequence &seq) { @@ -558,66 +646,6 @@ MTask &write3a = fWriteFitsFile ? static_cast(write3af) : static_cast(write3ar); - write3af.SetHeaderKey("DUMMY0", 3, "First dummy"); - write3af.SetHeaderKey("DUMMY1", 3.14159265358979323846, "Second dummy"); - write3af.SetHeaderKey("DUMMY2", true, "Third dummy"); - write3af.SetHeaderKey("DUMMY3", "one value", "Fourth dummy"); - - write3af.VetoColumn("MParameterD.fVal"); - write3af.VetoColumn("MCorsikaEvtHeader.fEvtNumber"); - write3af.VetoColumn("MCorsikaEvtHeader.fNumReuse"); - write3af.VetoColumn("MCorsikaEvtHeader.fTotalEnergy"); - write3af.VetoColumn("MCorsikaEvtHeader.fStartAltitude"); - write3af.VetoColumn("MCorsikaEvtHeader.fFirstTargetNum"); - write3af.VetoColumn("MCorsikaEvtHeader.fFirstInteractionHeight"); - write3af.VetoColumn("MCorsikaEvtHeader.fMomentumX"); - write3af.VetoColumn("MCorsikaEvtHeader.fMomentumY"); - write3af.VetoColumn("MCorsikaEvtHeader.fMomentumZ"); - write3af.VetoColumn("MCorsikaEvtHeader.fAz"); - write3af.VetoColumn("MCorsikaEvtHeader.fWeightedNumPhotons"); - write3af.VetoColumn("MCorsikaEvtHeader.fZd"); - write3af.VetoColumn("MCorsikaEvtHeader.fAd"); - write3af.VetoColumn("MCorsikaEvtHeader.fX"); - write3af.VetoColumn("MCorsikaEvtHeader.fY"); - write3af.VetoColumn("MCorsikaEvtHeader.fWeightNumPhotons"); - write3af.VetoColumn("MMcEvt.fEvtNumber"); - write3af.VetoColumn("MMcEvt.fThick0"); - write3af.VetoColumn("MMcEvt.fFirstTarget"); - write3af.VetoColumn("MMcEvt.fZFirstInteraction"); - write3af.VetoColumn("MMcEvt.fCoreD"); - write3af.VetoColumn("MMcEvt.fCoreX"); - write3af.VetoColumn("MMcEvt.fCoreY"); - write3af.VetoColumn("MMcEvt.fTimeFirst"); - write3af.VetoColumn("MMcEvt.fTimeLast"); - write3af.VetoColumn("MMcEvt.fLongiNmax"); - write3af.VetoColumn("MMcEvt.fLongit0"); - write3af.VetoColumn("MMcEvt.fLongitmax"); - write3af.VetoColumn("MMcEvt.fLongia"); - write3af.VetoColumn("MMcEvt.fLongib"); - write3af.VetoColumn("MMcEvt.fLongic"); - write3af.VetoColumn("MMcEvt.fLongichi2"); - write3af.VetoColumn("MMcEvt.fPhotIni"); - write3af.VetoColumn("MMcEvt.fPassPhotAtm"); - write3af.VetoColumn("MMcEvt.fPassPhotRef"); - write3af.VetoColumn("MMcEvt.fPassPhotCone"); - write3af.VetoColumn("MMcEvt.fPhotElfromShower"); - write3af.VetoColumn("MMcEvt.fPhotElinCamera"); - write3af.VetoColumn("MMcEvt.fElecCphFraction"); - write3af.VetoColumn("MMcEvt.fMuonCphFraction"); - write3af.VetoColumn("MMcEvt.fOtherCphFraction"); - write3af.VetoColumn("MMcEvt.fFadcTimeJitter"); - write3af.VetoColumn("MMcEvt.fEventReuse"); - - write3af.VetoColumn("MRawEvtData.fHiGainPixId"); - write3af.VetoColumn("MRawEvtData.fLoGainPixId"); - write3af.VetoColumn("MRawEvtData.fLoGainFadcSamples"); - write3af.VetoColumn("MRawEvtData.fABFlags"); - write3af.VetoColumn("MRawEvtData.fStartCells"); - write3af.VetoColumn("MRawEvtData.fNumBytesPerSample"); - write3af.VetoColumn("MRawEvtData.fIsSigned"); - write3af.VetoColumn("MRawEvtHeader.fDAQEvtNumber"); //EventNum ? - write3af.VetoColumn("MRawEvtHeader.fNumTrigLvl1"); - write3af.VetoColumn("MRawEvtHeader.fNumTrigLvl2"); - write3af.VetoColumn("MRawEvtHeader.fTrigPattern"); - write3af.VetoColumn("MRawEvtHeader.fNumLoGainOn"); + //SetupHeaderKeys(write3af,header); + SetupVetoColumns(write3af); write3af.SetBytesPerSample("Data", 2); @@ -649,4 +677,6 @@ write3af.AddContainer("IntendedPulsePos", "RunHeaders"); write3af.AddContainer("MRawEvtData", "Events"); + write3af.AddContainer("MMcEvtBasic", "Events"); + write3af.AddContainer("MPhotonEvent","Events"); write3ar.AddContainer("ElectronicNoise", "RunHeaders"); Index: /trunk/Mars/mjobs/MJSimulation.h =================================================================== --- /trunk/Mars/mjobs/MJSimulation.h (revision 17147) +++ /trunk/Mars/mjobs/MJSimulation.h (revision 17148) @@ -12,4 +12,5 @@ class MWriteRootFile; class MWriteFitsFile; +class MRawRunHeader; class MJSimulation : public MJob @@ -34,4 +35,6 @@ void SetupCommonFileStructure(MWriteRootFile &write) const; void SetupCommonFileStructure(MWriteFitsFile& write) const; + void SetupHeaderKeys(MWriteFitsFile& write, MRawRunHeader &header) const; + void SetupVetoColumns(MWriteFitsFile& write) const; public: Index: /trunk/Mars/mmc/MMcEvt.hxx =================================================================== --- /trunk/Mars/mmc/MMcEvt.hxx (revision 17147) +++ /trunk/Mars/mmc/MMcEvt.hxx (revision 17148) @@ -10,5 +10,5 @@ { private: - UInt_t fEvtNumber; + UInt_t fEvtNumber; // [fits: name=EventNum ; unit=uint32] Float_t fThick0; // [g/cm2] Float_t fFirstTarget; // [] Index: /trunk/Mars/mraw/MRawEvtData.cc =================================================================== --- /trunk/Mars/mraw/MRawEvtData.cc (revision 17147) +++ /trunk/Mars/mraw/MRawEvtData.cc (revision 17148) @@ -114,5 +114,5 @@ // MRawEvtData::MRawEvtData(const char *name, const char *title) - : fRunHeader(0), fNumBytesPerSample(1), fIsSigned(false) + : fRunHeader(0), fNumBytesPerSample(1), fTriggerType(0), fIsSigned(false) { fName = name ? name : "MRawEvtData"; @@ -359,4 +359,24 @@ fNumBytesPerSample = fRunHeader ? fRunHeader->GetNumBytesPerSample() : 1; + + + if (fRunHeader == 0) + { + fTriggerType = 0x0000; + } + else + { + switch(fRunHeader->GetRunType()) + { + case MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTData: + fTriggerType = 0x0004; break; + case MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTPedestal: + fTriggerType = 0x0400; break; + case MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTCalibration: + fTriggerType = 0x0164; break; + default: + fTriggerType = 0x0000; break; + } + } fHiGainPixId = new MArrayS(numconnected); Index: /trunk/Mars/mraw/MRawEvtData.h =================================================================== --- /trunk/Mars/mraw/MRawEvtData.h (revision 17147) +++ /trunk/Mars/mraw/MRawEvtData.h (revision 17148) @@ -12,4 +12,5 @@ //class ifstream; #include +#include class MRawRunHeader; @@ -36,5 +37,7 @@ MArrayB *fABFlags; //-> A information about the exact trigger position - MArrayS *fStartCells; // + MArrayS *fStartCells; // {fits: unit=uint16 ; name=StartCellData } + + UShort_t fTriggerType; // {fits: unit=uint16 ; name=TriggerType } UShort_t fNumBytesPerSample; Index: /trunk/Mars/mraw/MRawEvtHeader.h =================================================================== --- /trunk/Mars/mraw/MRawEvtHeader.h (revision 17147) +++ /trunk/Mars/mraw/MRawEvtHeader.h (revision 17148) @@ -53,5 +53,5 @@ MTime *fTime; //! object to store the time in (ReadEvt) - UInt_t fDAQEvtNumber; // Number of Event [fits: name=EventNum] + UInt_t fDAQEvtNumber; // Number of Event {fits: name=EventNum ; unit=uint32} UInt_t fNumTrigLvl1; // Number of 1st level tiggers between 2 events Index: /trunk/Mars/msimcamera/MSimCamera.cc =================================================================== --- /trunk/Mars/msimcamera/MSimCamera.cc (revision 17147) +++ /trunk/Mars/msimcamera/MSimCamera.cc (revision 17148) @@ -43,4 +43,5 @@ #include #include +#include #include "MLog.h" @@ -63,4 +64,6 @@ #include "MMcEvt.hxx" // To be replaced by a CheObs class #include "MRawRunHeader.h" + +#include "math.h" ClassImp(MSimCamera); @@ -75,5 +78,6 @@ : fEvt(0), fStat(0), fRunHeader(0), fElectronicNoise(0), fGain(0), fCamera(0), fMcEvt(0), fSpline(0), fBaselineGain(kFALSE), - fDefaultOffset(-1), fDefaultNoise(-1), fDefaultGain(-1) + fDefaultOffset(-1), fDefaultNoise(-1), fDefaultGain(-1), fACFudgeFactor(0), + fACTimeConstant(0) { @@ -135,4 +139,11 @@ return kFALSE; + fAccidentalPhotons = (MPedestalCam*)pList->FindObject("AccidentalPhotonRates","MPedestalCam"); + if(!fAccidentalPhotons) + { + *fLog << err << "AccidentalPhotonRates [MPedestalCam] not found... aborting." << endl; + return kFALSE; + } + MParSpline *pulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline"); if (!pulse) @@ -142,10 +153,10 @@ } - if (fRunHeader->GetFreqSampling()!=1000) - { - *fLog << err << "ERROR - Sampling frequencies others than 1GHz are not yet supported." << endl; - *fLog << warn << "FIXME - SCALE MPulsShape WITH THE SAMPLING FREQUENCY." << endl; - return kFALSE; - } +// if (fRunHeader->GetFreqSampling()!=1000) +// { +// *fLog << err << "ERROR - Sampling frequencies others than 1GHz are not yet supported." << endl; +// *fLog << warn << "FIXME - SCALE MPulsShape WITH THE SAMPLING FREQUENCY." << endl; +// return kFALSE; +// } fSpline = pulse->GetSpline(); @@ -238,7 +249,57 @@ const Double_t rms = pix.GetPedestalRms(); + // FTemme: Implementation of AC-coupling: + // to calculate the value of the accoupling per slice I use the + // following equation: + // accouplingPerSlice = accidentalPhotonRate * (1 + crossTalkProb) + // * areaOfOnePulse / samplingRate; + // Therefore I need the following variables + Double_t accidentalPhotonRate = 0; // [MHz] + Float_t crossTalkProb = 0; // [1] + Double_t areaOfOnePulse = 0; // [ADC-Counts * s] + Double_t samplingRate = 0; // [slices * MHz] + + // The accidental photon rate is stored in GHz, so we have to multiply + // with 1E3 to get MHz: + MPedestalPix &accPhoPix = (*fAccidentalPhotons)[i]; + accidentalPhotonRate = accPhoPix.GetPedestal() * 1E3; + Double_t currentAccidentalPhotonRate = accidentalPhotonRate; + + if(fACTimeConstant!=0) + { + Double_t accidentalPhotons = fACTimeConstant * accidentalPhotonRate; + Double_t sigmaAccidentalPhotons = sqrt(accidentalPhotons); + + TRandom3 *random = new TRandom3(0); + Double_t gaus = random->Gaus(accidentalPhotons,sigmaAccidentalPhotons); + currentAccidentalPhotonRate = gaus / fACTimeConstant; + } + + // I don't know how to get the variable fCrosstalkProb from + // the class APD (see MAvalanchePhotoDiode.h), because there is no + // getter for the APD array(fAPDs) in MSimAPD. + // So I set the crossTalkProb hardcoded to the value 0.15, which is + // equal to the value of the apd of type 4 + crossTalkProb = 0.15; + + // To get the area of one Pulse, I only need to calculate the Integral + // of the Pulse Shape, which is stored in fSpline. Because the spline is + // normalized to a maximal amplitude of 1.0, I had to multiply it with + // the Default gain: + areaOfOnePulse = fSpline->Integral() * fDefaultGain; + + // The sampling rate I get from the RunHeader: + samplingRate = fRunHeader->GetFreqSampling(); + + Double_t accouplingPerSlice = currentAccidentalPhotonRate + * (1 + crossTalkProb + fACFudgeFactor) + * areaOfOnePulse / samplingRate; + + // The accoupling is substracted from the timeline by decreasing the + // mean of the gaussian noise which is added + if (!fBaselineGain) { - (*fCamera)[i].AddGaussianNoise(rms, val); + (*fCamera)[i].AddGaussianNoise(rms, val - accouplingPerSlice); continue; } @@ -250,5 +311,6 @@ // FIXME: How stable is the offset? // FIXME: Should we write a container AppliedGain for MSImTrigger? - (*fCamera)[i].AddGaussianNoise(rms*gain, val*gain); + + (*fCamera)[i].AddGaussianNoise(rms*gain, (val - accouplingPerSlice)*gain); } @@ -323,4 +385,14 @@ fDefaultGain = GetEnvValue(env, prefix, "DefaultGain", fDefaultGain); } + if (IsEnvDefined(env, prefix, "fACFudgeFactor", print)) + { + rc = kTRUE; + fACFudgeFactor = GetEnvValue(env, prefix, "fACFudgeFactor", fACFudgeFactor); + } + if (IsEnvDefined(env, prefix, "fACTimeConstant", print)) + { + rc = kTRUE; + fACTimeConstant = GetEnvValue(env, prefix, "fACTimeConstant", fACTimeConstant); + } return rc; Index: /trunk/Mars/msimcamera/MSimCamera.h =================================================================== --- /trunk/Mars/msimcamera/MSimCamera.h (revision 17147) +++ /trunk/Mars/msimcamera/MSimCamera.h (revision 17148) @@ -24,9 +24,10 @@ MPedestalCam *fElectronicNoise; //! Electronic noise (baseline and rms) MPedestalCam *fGain; //! Electronic noise (baseline and rms) + MPedestalCam *fAccidentalPhotons;//! Accidental photon rates MAnalogChannels *fCamera; //! Output of the analog signals MMcEvt *fMcEvt; //! For information stored in MMcEvt - const MSpline3 *fSpline; // Pulse Shape + const MSpline3 *fSpline; // Pulse Shape Bool_t fBaselineGain; // Should the gain be applied to baseline and electronic noise? @@ -35,4 +36,8 @@ Double_t fDefaultNoise; // Default noise (RMS) added to all samples Double_t fDefaultGain; // Default gain (multiplication factor to the given/used pulse shape) + + Double_t fACFudgeFactor; + + Double_t fACTimeConstant; // MParContainer Index: /trunk/Mars/msimcamera/MSimReadout.cc =================================================================== --- /trunk/Mars/msimcamera/MSimReadout.cc (revision 17147) +++ /trunk/Mars/msimcamera/MSimReadout.cc (revision 17148) @@ -145,4 +145,5 @@ fData->InitRead(fRunHeader); fData->ResetPixels(); + fData->InitStartCells(); fData->SetIndices(); @@ -202,5 +203,5 @@ // Check if the whole requested signal can be digitized - if (trig+nslices>=nsamp) + if (trig+nslices>nsamp) { *fLog << err << "ERROR - Trigger position beyond valid analog signal range." << endl;