Index: /releases/Mars.2014.05.26/fact/analysis/callisto_data.C =================================================================== --- /releases/Mars.2014.05.26/fact/analysis/callisto_data.C (revision 18028) +++ /releases/Mars.2014.05.26/fact/analysis/callisto_data.C (revision 18028) @@ -0,0 +1,442 @@ +#include "MLogManip.h" + +int callisto_data(const char *datafile, const char *drsfile, const char *drstime, const char *delays, const char *outpath=".") +{ + // ====================================================== + + // 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; + + //const char *pulse_template = "template-pulse.root"; + + // ------------------------------------------------------ + + 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 + + // ------------------------------------------------------ + + // nominal gain as extracted with gain extraction + // That is the number which describes + // the pulse amplitude in the file! + double gain = 241.; + + // integrate around maximum or from leading edge on? + const bool maximum = true; + + // ------------------------------------------------------ + + // 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; + const int fall_time_cal = maximum ? 120 : 160; + + const int rise_time_dat = maximum ? 5 : 0; + const int fall_time_dat = maximum ? 25 : 30; + + // Extraction type: Extract integral and half leading edge + + const int 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 + + // ====================================================== + + if (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 11; + } + + // ------------------------------------------------------- + + TFile file(drstime); + if (file.IsZombie()) + { + gLog << err << "ERROR - Could not open " << drstime << endl; + return 21; + } + + MDrsCalibrationTime *timecam = 0; + file.GetObject("MDrsCalibrationTime", timecam); + if (!timecam) + { + gLog << err << "ERROR - Could not get MDrsCalibrationTime from " << drstime << endl; + return 22; + } + + // ------------------------------------------------------ + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile)) + return 31; + + // ------------------------------------------------------- + + //if (use_delays) + // timecam.SetDelays(*evt1h.GetHist()); + + // ------------------------------------------------------- + + // Range in which function/spline is defined + int minx = -10; + int maxx = 125; + + // Number of points + int N = (maxx-minx)*2; // two points per nanosecond (2GHz) + + TF1 f("pulse", "[0]*(1-exp(-[1]*x*x))*exp(-[2]*x)*(x>0)", float(minx), float(maxx)); + f.SetNpx(N); + + // Normalize pulse + f.SetParameters(1.9*gain/30, 0.10, 0.053); + + // Convert to graph + TGraph g(&f); + + // Convert to float + MArrayF data(g.GetN()); + for (int i=0; iSetOptFit(kTRUE); + + // =================================================================== + + gLog << endl; + gLog.Separator("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(datafile); + read5.LoadMap(map); + + 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); + + // --- + 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(); + + TString fname = gSystem->ConcatFileName(outpath, gSystem->BaseName(datafile)); + fname.ReplaceAll(".fits.gz", "_C.root"); + fname.ReplaceAll(".fits.fz", "_C.root"); + fname.ReplaceAll(".fits", "_C.root"); + + // The second rule is for the case reading raw-files! + MWriteRootFile write5(fname, "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 ---------------- + + //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(&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(max)) + return 18; + + if (!loop5.GetDisplay()) + return 19; + + // ============================================================ + + TFile *ofile = gROOT->GetListOfFiles()->FindObject(fname); + if (!ofile || !ofile->IsOpen() || ofile->IsZombie()) + { + gLog << err << "File " << fname << " not found" << endl; + return 20; + } + + TString title = "-- Calibrated signal ["; + title += datafile; + title += "] --"; + d->SetTitle(title, kFALSE); + + ofile->cd(); + d->Write(); + + return 0; +} Index: /releases/Mars.2014.05.26/fact/analysis/callisto_drstime.C =================================================================== --- /releases/Mars.2014.05.26/fact/analysis/callisto_drstime.C (revision 18028) +++ /releases/Mars.2014.05.26/fact/analysis/callisto_drstime.C (revision 18028) @@ -0,0 +1,91 @@ +#include "MLogManip.h" + +int callisto_drstime(const char *timfile, const char *drs1024, const char *outfile="callisto-drstime.root") +{ + // ====================================================== + + // Maximum number of events to be procesed + Long_t max = 0; + + // ====================================================== + + MDrsCalibration drscalib1024; + if (!drscalib1024.ReadFits(drs1024)) + return 52; + + // ====================================================== + + MStatusDisplay *d = new MStatusDisplay; + + // ====================================================== + + gLog.Separator("Callisto"); + gLog << all; + gLog << "DRS Timing " << timfile << '\n'; + gLog << "DRS 1024 " << drs1024 << '\n'; + gLog << endl; + + + MDrsCalibrationTime timecam; + + gStyle->SetOptFit(kTRUE); + + // ====================================================== + + gLog << endl; + gLog.Separator("Processing DRS timing calibration run"); + + MTaskList tlist0; + + MParList plist0; + plist0.AddToList(&tlist0); + plist0.AddToList(&drscalib1024); + 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(max)) + return 8; + + if (!loop0.GetDisplay()) + return 9; + + //MHDrsCalibrationTime *t = (MHDrsCalibrationTime*)plist0.FindObject("MHDrsCalibrationTime"); + //t->SetDisplay(d); + //t->PlotAll(); + + // ====================================================== + + TString title = "-- DRS timing "; + title += " ["; + title += timfile; + title += "] --"; + d->SetTitle(title, kFALSE); + d->SaveAs(outfile); + + TFile file(outfile, "UPDATE"); + timecam.Write(); + + return 0; +} Index: /releases/Mars.2014.05.26/fact/analysis/merpp_file.C =================================================================== --- /releases/Mars.2014.05.26/fact/analysis/merpp_file.C (revision 18028) +++ /releases/Mars.2014.05.26/fact/analysis/merpp_file.C (revision 18028) @@ -0,0 +1,94 @@ +#include "MLogManip.h" + +void merpp(TString &froot, const char *id, const char *faux) +{ + gLog.Separator(froot+" - "+TString(id)); + + TFile file(froot); + if (file.IsZombie()) + { + gLog << err << "merpp: Invalid file '" << froot << "'" << endl; + return; + } + + MRawRunHeader *h = NULL; + TTree *tree = file.Get("RunHeaders"); + if (!tree) + { + gLog << err << "merpp: Tree 'RunHeaders' not found... skipped." << endl; + return; + } + + if (tree->GetEntries()!=1) + { + gLog << err << "merpp: Number of RunHeaders do not match 1... skipped." << endl; + return; + } + + tree->SetBranchAddress("MRawRunHeader.", &h); + tree->GetEntry(0); + + if (file.Get(id)) + { + gLog << warn << "WARNING - Tree '" << id << "' already existing... skipped." << endl; + return; + } + + file.Close(); + + UInt_t night = h->GetRunNumber(); + + TString ffits = Form("/fact/aux/%4d/%02d/%02d/%6d.%s.fits", + night/10000, (night/100)%100, night%100, night, faux); + + TString report = Form("MReport%s", id); + TString time = Form("MTime%s", id); + + gLog << all; + gLog << " --- Fits file: " << ffits << endl; + gLog << " --- Root file: " << froot << endl; + gLog << " --- Tree: " << id << " (" << faux << ")" << endl; + + MReportFitsRead read(ffits); + read.SetReportName(report); + read.SetTimeStart(h->GetRunStart()); + read.SetTimeStop(h->GetRunEnd()); + + // FIXME: Write also last event BEFORE start of run + + MWriteRootFile write(froot, "UPDATE"); + write.AddContainer(report, id); + write.AddContainer(time, id); + + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + + tlist.AddToList(&read); + tlist.AddToList(&write); + + MEvtLoop loop; + loop.SetParList(&plist); + + if (!loop.Eventloop()) + return; +} + +int merpp_file(const char *datafile) +{ + // ------------------------------------------------------ + + gLog.Separator("Merpp"); + gLog << all; + gLog << "Merge slow control data of "; + gLog << datafile << endl; + gLog << endl; + + merpp(datafile, "Weather", "MAGIC_WEATHER_DATA"); + merpp(datafile, "Drive", "DRIVE_CONTROL_TRACKING_POSITION"); + merpp(datafile, "Rates", "FTM_CONTROL_TRIGGER_RATES"); + merpp(datafile, "Temperatures", "FSC_CONTROL_TEMPERATURE"); + merpp(datafile, "Humidity", "FSC_CONTROL_HUMIDITY"); + + return 0; +} Index: /releases/Mars.2014.05.26/fact/analysis/star_file.C =================================================================== --- /releases/Mars.2014.05.26/fact/analysis/star_file.C (revision 18028) +++ /releases/Mars.2014.05.26/fact/analysis/star_file.C (revision 18028) @@ -0,0 +1,284 @@ +#include "MLogManip.h" + +int star_file(const char *datafile, Double_t lvl1=7.8, Double_t lvl2=3.9, const char *outpath = ".") +{ + double deltat = 17.5; // ns/deg + // lvl1 = 2.5; + // lvl2 = 0.5; + + // ------------------------------------------------------ + + gLog.Separator("Star"); + gLog << all << "Calculate image parameters of sequence "; + gLog << datafile << endl; + gLog << endl; + + // ------------------------------------------------------ + + gLog << "Outpath: " << outpath << endl; + + TString fname = gSystem->ConcatFileName(outpath, gSystem->BaseName(datafile)); + fname.ReplaceAll("_C.root", "_I.root"); + + gLog.Separator(); + + // -------------------------------------------------------- + + MBinning bins1( 36, -90, 90, "BinningAlpha"); + MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin"); + MBinning bins4( 25, 0, 2.5, "BinningDist"); + + MBinning binsM1(100, 0, 5, "BinningMuonRadius"); + MBinning binsM2( 60, 0, 0.3, "BinningMuonDeviation"); + MBinning binsM3(150, 0, 60, "BinningMuonTime"); + MBinning binsM4(300, 0, 30, "BinningMuonTimeRms"); + MBinning binsM5( 20, 0, 360, "BinningMuonArcPhi"); + MBinning binsM6( 50, 0, 0.5, "BinningMuonArcWidth"); + MBinning binsM7( 30, 5, 5000, "BinningMuonSize", "log"); + MBinning binsM8(100, 0, 5, "BinningMuonRelTimeSigma"); + + MBinning binsM9(100, 0, 5, "BinningRadius"); + + // --------------------------------------------------------- + + // Filter to start muon analysis + //MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut"); + // Filter to calculate further muon parameters +// MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&" +// "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut"); + MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&" + "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4", + "MuonSearchCut"); + + // Filter to fill the MHMuonPar + MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2", + "MuonFinalCut"); + // Filter to write Muons to Muon tree + //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut"); + + // --------------------------------------------------------- + + MStatusDisplay *d = new MStatusDisplay; + + MTaskList tlist; + + MParList plist2; + plist2.AddToList(&tlist); + plist2.AddToList(&bins1); + plist2.AddToList(&bins3); + plist2.AddToList(&bins4); + plist2.AddToList(&binsM1); + plist2.AddToList(&binsM2); + plist2.AddToList(&binsM3); + plist2.AddToList(&binsM4); + plist2.AddToList(&binsM5); + plist2.AddToList(&binsM6); + plist2.AddToList(&binsM7); + plist2.AddToList(&binsM8); + plist2.AddToList(&binsM9); + + MMuonSetup muonsetup; + plist2.AddToList(&muonsetup); + + MEvtLoop loop; + loop.SetDisplay(d); + loop.SetParList(&plist2); + + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + read.AddFile(datafile); + + 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); + + // --------------------------------------------------------------- + + MWriteRootFile write(fname, "RECREATE", "Image parameters", 2); // EffectiveOnTime + write.AddContainer("MTime", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MImagePar", "Events"); + write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MRawEvtHeader", "Events"); + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MGeomCam", "RunHeaders"); + + MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0", + "MuonHistCut"); + fillmhn2.SetFilter(&fmuonhn); + + tlist.AddToList(&read); + tlist.AddToList(&cont); + tlist.AddToList(&apply); + tlist.AddToList(&clean); + tlist.AddToList(&hcalc); + tlist.AddToList(&fillC1); + tlist.AddToList(&fillC2); + tlist.AddToList(&fillC3); + tlist.AddToList(&fillR); + tlist.AddToList(&fillD1); + tlist.AddToList(&fillD2); + tlist.AddToList(&fillD3); + tlist.AddToList(&fillD4); + tlist.AddToList(&fillD5); + tlist.AddToList(&fillD6); + tlist.AddToList(&filla); + //tlist.AddToList(&fmuon1); + tlist.AddToList(&muscalc); + tlist.AddToList(&fillmhn1); + tlist.AddToList(&fillmhn3); + tlist.AddToList(&fmuon2); + tlist.AddToList(&fillmuon); + tlist.AddToList(&mcalc); + tlist.AddToList(&fmuonhn); + tlist.AddToList(&fillmhn2); + tlist.AddToList(&fillmhn4); + tlist.AddToList(&fmuon3); + tlist.AddToList(&fillmpar); + //tlist.AddToList(&fmuon4); + //tlist.AddToList(&writem); + tlist.AddToList(&write); + + if (!loop.Eventloop()) + return 3; + + if (!loop.GetDisplay()) + return 4; + + // ============================================================ + + TFile *ofile = gROOT->GetListOfFiles()->FindObject(fname); + if (!ofile || !ofile->IsOpen() || ofile->IsZombie()) + { + gLog << err << "File " << fname << " not found" << endl; + return 20; + } + + TString title = "-- Image parameters ["; + title += datafile; + title += "] --"; + d->SetTitle(title, kFALSE); + + ofile->cd(); + d->Write(); + + return 0; +}