#include "MLogManip.h" /***************************************************************** STAR -- STandard Analysis and Reconstruction datafile: A calibarted root file which came out of callist, e.g. 20170727_006_Y.root or 00000015.003_Y_MonteCarlo003_Events.root lvl1, lvl2: Lower (1) and higher (2) cleaning level outpath: The path to which the output files are written To run the macro from the command line (assuming you are in a directory Mars/build where you have built your Mars environment) ou can do root ../hawc/star.C\(\"00000015.003_Y_MonteCarlo003_Events.root\"\) or from within root [0] .x ../hawc/star.C("00000015.003_Y_MonteCarlo003_Events.root") ******************************************************************/ int star(const char *datafile, const char *outpath=".") { // minimum number of pixels that can make an island int mincount = 0; // minimum size for a valid island double minsize = 25; // The delta t [ns/deg] for the image cleaning double deltat = 2.5; // -------------------------------------------------------- // Define the binning for the plots MBinning binse( 40, 10 , 1e4, "BinningSize", "log"); MBinning bins1( 36, -90, 90, "BinningAlpha"); MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin"); MBinning bins4( 40, 0, 8, "BinningDist"); MBinning bins5( 30, 0, 6, "BinningWidth"); MBinning bins6( 30, 0, 6, "BinningLength"); MBinning bins7( 31, -6, 6, "BinningM3Long"); MBinning bins8( 31, -6, 6, "BinningM3Trans"); MBinning bins9( 31, -6, 6, "BinningAsym"); MBinning binsA( 100, 0, 10, "BinningArea"); MBinning binsT( 30, 0, 6, "BinningTimeSpread"); // ------------------------------------------------------ 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_", "_I_"); // gSystem->ExpandPathName(fname); gLog.Separator(); // --------------------------------------------------------- // Allocate/open the status display MStatusDisplay *d = new MStatusDisplay; // Instantiate the list of tasks MTaskList tlist; // Instantiate the list of parameter containers MParList plist2; plist2.AddToList(&tlist); // Fill it with the binnings plist2.AddToList(&binse); plist2.AddToList(&bins1); plist2.AddToList(&bins3); plist2.AddToList(&bins4); plist2.AddToList(&bins5); plist2.AddToList(&bins6); plist2.AddToList(&bins7); plist2.AddToList(&bins8); plist2.AddToList(&bins9); plist2.AddToList(&binsA); plist2.AddToList(&binsT); // Instatiate the event loop MEvtLoop loop; loop.SetDisplay(d); loop.SetParList(&plist2); MDirIter files(datafile); // Instantiate the reading task // You can use // read.AddFiles("*_C.root") // for example to read more than one file at once MReadMarsFile read("Events"); read.DisableAutoScheme(); read.AddFiles(files); // Instantiate the task which takes care of the size of all containers MGeomApply apply; // Instantiate the image cleaning // This is the most simple image cleaning possible, standard // tail cut based on absolute amplitudes // MImgCleanStd clean(lvl1, lvl2); // clean.SetMethod(MImgCleanStd::kAbsolute); // clean.SetPostCleanType(3); // Instantiate the image cleaning as described in the TRAC MImgCleanTime clean; clean.SetMinCount(mincount); clean.SetMinSize(minsize); clean.SetDeltaT(deltat); // Instantiate the calculation of the image parameters MHillasCalc hcalc; hcalc.Disable(MHillasCalc::kCalcConc); // Instantiate a plot for the rate after cleaning 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); // Instantiate the task to fill the rate plot MFillH fillR(&hrate, "", "FillRate"); // Instantiate the camera displays to be filled 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); // Instantiate tasks to fill the camera displays MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam"); MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels"); MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos"); // Instantiate tasks to fill image parameter histograms 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"); // Instantiate an example for a user defined plot MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)"); halpha.SetName("AvsD;Dist;Alpha"); // Instantiate a task to fill that plot MFillH filla(&halpha); filla.SetNameTab("AvsD"); filla.SetDrawOption("colz"); // This is an alternative if more than one file should be read in a single loop // In this case the output file name of each file is created with this rule // from the input file name. Note that the filename of the output file // with the status display must then be explicitly given. // //const TString fname(Form("s/([0-9]+_[0-9]+)[.]_C.root$/%s\\/$1_I.root/", // MJob::Esc(outpath).Data())); //MWriteRootFile write5(2, fname, "RECREATE", "Image parameters"); const TString rule(Form("s/([0-9]+)([._])([0-9]+)(_Y)?([_].*)?[.]root$/%s\\/$1$2$3_I$5.root/", MJob:: Esc(outpath).Data())); // Instantiate writing the file MWriteRootFile write(2, rule, "RECREATE", "Image parameters"); 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("MSoftwareTrigger","Events"); write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MGeomCam", "RunHeaders"); write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE); write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE); write.AddContainer("MMcEvt", "Events", kFALSE); write.AddContainer("IncidentAngle", "Events", kFALSE); write.AddContainer("MPointingPos", "Events", kFALSE); write.AddContainer("MTime", "Events", kFALSE); write.AddContainer("MSrcPosCam", "Events", kFALSE); // Setup the task list tlist.AddToList(&read); 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(&write); // Run the eventloop if (!loop.Eventloop()) return 3; // Check if the display was closed (deleted) by the user if (!loop.GetDisplay()) return 4; // ============================================================ TString fname = write.GetFileName(); fname.ReplaceAll(".root", "-display.root"); // Write the status display to the file TString title = "-- Image parameters ["; title += datafile; title += "] --"; d->SetTitle(title, kFALSE); d->SaveAsRoot(fname); return 0; }