| 1 | #include "MLogManip.h"
|
|---|
| 2 |
|
|---|
| 3 | /*****************************************************************
|
|---|
| 4 |
|
|---|
| 5 | STAR -- STandard Analysis and Reconstruction
|
|---|
| 6 |
|
|---|
| 7 | datafile:
|
|---|
| 8 | A calibarted root file which came out of callist,
|
|---|
| 9 | e.g. 20170727_006_C.root
|
|---|
| 10 |
|
|---|
| 11 | lvl1, lvl2:
|
|---|
| 12 | Lower (1) and higher (2) cleaning level
|
|---|
| 13 |
|
|---|
| 14 | outpath:
|
|---|
| 15 | The path to which the output files are written
|
|---|
| 16 |
|
|---|
| 17 |
|
|---|
| 18 | To run the macro from the command line (assuming you are in a directory
|
|---|
| 19 | Mars/build where you have built your Mars environment) ou can do
|
|---|
| 20 |
|
|---|
| 21 | root ../hawc/star.C\(\"filename_C.root\"\)
|
|---|
| 22 |
|
|---|
| 23 | or from within root
|
|---|
| 24 |
|
|---|
| 25 | [0] .x ../hawc/star.C("filename_C.root")
|
|---|
| 26 |
|
|---|
| 27 | ******************************************************************/
|
|---|
| 28 | int star(const char *datafile, const char *outpath=".")
|
|---|
| 29 | {
|
|---|
| 30 | // The delta t [ns/deg] for the image cleaning
|
|---|
| 31 | double deltat = 17.5;
|
|---|
| 32 |
|
|---|
| 33 | // --------------------------------------------------------
|
|---|
| 34 |
|
|---|
| 35 | // Define the binning for the plots
|
|---|
| 36 | MBinning binse( 40, 10 , 1e4, "BinningSize", "log");
|
|---|
| 37 | MBinning bins1( 36, -90, 90, "BinningAlpha");
|
|---|
| 38 | MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
|
|---|
| 39 | MBinning bins4( 40, 0, 8, "BinningDist");
|
|---|
| 40 | MBinning bins5( 40, 0, 2, "BinningWidth");
|
|---|
| 41 | MBinning bins6( 40, 0, 2, "BinningLength");
|
|---|
| 42 | MBinning bins7( 41, -2, 2, "BinningM3Long");
|
|---|
| 43 | MBinning bins8( 41, -2, 2, "BinningM3Trans");
|
|---|
| 44 | MBinning bins9( 61, -3, 3, "BinningAsym");
|
|---|
| 45 | MBinning binsA( 100, 0, 10, "BinningArea");
|
|---|
| 46 |
|
|---|
| 47 | MBinning binsM9(100, 0, 5, "BinningRadius");
|
|---|
| 48 |
|
|---|
| 49 | // ------------------------------------------------------
|
|---|
| 50 |
|
|---|
| 51 | gLog.Separator("Star");
|
|---|
| 52 | gLog << all << "Calculate image parameters of sequence ";
|
|---|
| 53 | gLog << datafile << endl;
|
|---|
| 54 | gLog << endl;
|
|---|
| 55 |
|
|---|
| 56 | // ------------------------------------------------------
|
|---|
| 57 |
|
|---|
| 58 | gLog << "Outpath: " << outpath << endl;
|
|---|
| 59 |
|
|---|
| 60 | TString fname = gSystem->ConcatFileName(outpath, gSystem->BaseName(datafile));
|
|---|
| 61 | fname.ReplaceAll("_C_", "_I_");
|
|---|
| 62 |
|
|---|
| 63 | gSystem->ExpandPathName(fname);
|
|---|
| 64 |
|
|---|
| 65 | gLog.Separator();
|
|---|
| 66 |
|
|---|
| 67 | // ---------------------------------------------------------
|
|---|
| 68 |
|
|---|
| 69 | // Allocate/open the status display
|
|---|
| 70 | MStatusDisplay *d = new MStatusDisplay;
|
|---|
| 71 |
|
|---|
| 72 | // Instantiate the list of tasks
|
|---|
| 73 | MTaskList tlist;
|
|---|
| 74 |
|
|---|
| 75 | // Instantiate the list of parameter containers
|
|---|
| 76 | MParList plist2;
|
|---|
| 77 | plist2.AddToList(&tlist);
|
|---|
| 78 |
|
|---|
| 79 | // Fill it with the binnings
|
|---|
| 80 | plist2.AddToList(&binse);
|
|---|
| 81 | plist2.AddToList(&bins1);
|
|---|
| 82 | plist2.AddToList(&bins3);
|
|---|
| 83 | plist2.AddToList(&bins4);
|
|---|
| 84 | plist2.AddToList(&bins5);
|
|---|
| 85 | plist2.AddToList(&bins6);
|
|---|
| 86 | plist2.AddToList(&bins7);
|
|---|
| 87 | plist2.AddToList(&bins8);
|
|---|
| 88 | plist2.AddToList(&bins9);
|
|---|
| 89 | plist2.AddToList(&binsA);
|
|---|
| 90 |
|
|---|
| 91 | // Instatiate the event loop
|
|---|
| 92 | MEvtLoop loop;
|
|---|
| 93 | loop.SetDisplay(d);
|
|---|
| 94 | loop.SetParList(&plist2);
|
|---|
| 95 |
|
|---|
| 96 | MDirIter files(datafile);
|
|---|
| 97 |
|
|---|
| 98 | // Instantiate the reading task
|
|---|
| 99 | // You can use
|
|---|
| 100 | // read.AddFiles("*_C.root")
|
|---|
| 101 | // for example to read more than one file at once
|
|---|
| 102 | MReadMarsFile read("Events");
|
|---|
| 103 | read.DisableAutoScheme();
|
|---|
| 104 | read.AddFiles(files);
|
|---|
| 105 |
|
|---|
| 106 | // Instantiate the task which takes care of the size of all containers
|
|---|
| 107 | MGeomApply apply;
|
|---|
| 108 |
|
|---|
| 109 | // Instantiate the image cleaning
|
|---|
| 110 |
|
|---|
| 111 | // This is the most simple image cleaning possible, standard
|
|---|
| 112 | // tail cut based on absolute amplitudes
|
|---|
| 113 | // MImgCleanStd clean(lvl1, lvl2);
|
|---|
| 114 | // clean.SetMethod(MImgCleanStd::kAbsolute);
|
|---|
| 115 | // clean.SetPostCleanType(3);
|
|---|
| 116 |
|
|---|
| 117 | // Instantiate the image cleaning as described in the TRAC
|
|---|
| 118 | MImgCleanTime clean;
|
|---|
| 119 | clean.SetMinCount(0);
|
|---|
| 120 | clean.SetMinSize(55);
|
|---|
| 121 | clean.SetDeltaT(17.5);
|
|---|
| 122 |
|
|---|
| 123 | // Instantiate the calculation of the image parameters
|
|---|
| 124 | MHillasCalc hcalc;
|
|---|
| 125 | hcalc.Disable(MHillasCalc::kCalcConc);
|
|---|
| 126 |
|
|---|
| 127 | // Instantiate a plot for the rate after cleaning
|
|---|
| 128 | MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
|
|---|
| 129 | hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1.)");
|
|---|
| 130 | hrate.SetName("Rate");
|
|---|
| 131 | hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];");
|
|---|
| 132 | hrate.InitLabels(MH3::kLabelsX);
|
|---|
| 133 | hrate.GetHist().SetMinimum(0);
|
|---|
| 134 |
|
|---|
| 135 | // Instantiate the task to fill the rate plot
|
|---|
| 136 | MFillH fillR(&hrate, "", "FillRate");
|
|---|
| 137 |
|
|---|
| 138 | // Instantiate the camera displays to be filled
|
|---|
| 139 | MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]");
|
|---|
| 140 | MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
|
|---|
| 141 | MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]");
|
|---|
| 142 | evtC1.SetErrorSpread(kFALSE);
|
|---|
| 143 | evtC2.SetErrorSpread(kFALSE);
|
|---|
| 144 | evtC2.SetThreshold(0);
|
|---|
| 145 |
|
|---|
| 146 | // Instantiate tasks to fill the camera displays
|
|---|
| 147 | MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam");
|
|---|
| 148 | MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels");
|
|---|
| 149 | MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos");
|
|---|
| 150 |
|
|---|
| 151 | // Instantiate tasks to fill image parameter histograms
|
|---|
| 152 | MFillH fillD1("MHHillas", "MHillas", "FillHillas");
|
|---|
| 153 | MFillH fillD2("MHHillasExt", "", "FillHillasExt");
|
|---|
| 154 | MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
|
|---|
| 155 | MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar");
|
|---|
| 156 | MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
|
|---|
| 157 | MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize");
|
|---|
| 158 |
|
|---|
| 159 | // Instantiate an example for a user defined plot
|
|---|
| 160 | MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)");
|
|---|
| 161 | halpha.SetName("AvsD;Dist;Alpha");
|
|---|
| 162 |
|
|---|
| 163 | // Instantiate a task to fill that plot
|
|---|
| 164 | MFillH filla(&halpha);
|
|---|
| 165 | filla.SetNameTab("AvsD");
|
|---|
| 166 | filla.SetDrawOption("colz");
|
|---|
| 167 |
|
|---|
| 168 | // This is an alternative if more than one file should be read in a single loop
|
|---|
| 169 | // In this case the output file name of each file is created with this rule
|
|---|
| 170 | // from the input file name. Note that the filename of the output file
|
|---|
| 171 | // with the status display must then be explicitly given.
|
|---|
| 172 | //
|
|---|
| 173 | //const TString fname(Form("s/([0-9]+_[0-9]+)[.]_C.root$/%s\\/$1_I.root/",
|
|---|
| 174 | // MJob::Esc(outpath).Data()));
|
|---|
| 175 | //MWriteRootFile write5(2, fname, "RECREATE", "Image parameters");
|
|---|
| 176 |
|
|---|
| 177 | const TString rule(Form("s/(([0-9]+_)?[0-9.]+)_Y_(.*)([.]root)$/%s\\/$1_Y_$3.root/",
|
|---|
| 178 | MJob:: Esc(outpath).Data()));
|
|---|
| 179 |
|
|---|
| 180 | // Instantiate writing the file
|
|---|
| 181 | MWriteRootFile write(2, rule, "RECREATE", "Image parameters");
|
|---|
| 182 | write.AddContainer("MTime", "Events");
|
|---|
| 183 | write.AddContainer("MHillas", "Events");
|
|---|
| 184 | write.AddContainer("MHillasExt", "Events");
|
|---|
| 185 | write.AddContainer("MHillasSrc", "Events");
|
|---|
| 186 | write.AddContainer("MImagePar", "Events");
|
|---|
| 187 | write.AddContainer("MNewImagePar", "Events");
|
|---|
| 188 | write.AddContainer("MRawEvtHeader", "Events");
|
|---|
| 189 | write.AddContainer("MSoftwareTrigger","Events");
|
|---|
| 190 | write.AddContainer("MRawRunHeader", "RunHeaders");
|
|---|
| 191 | write.AddContainer("MGeomCam", "RunHeaders");
|
|---|
| 192 |
|
|---|
| 193 | write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
|
|---|
| 194 | write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
|
|---|
| 195 | write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
|
|---|
| 196 | write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
|
|---|
| 197 | write.AddContainer("MMcEvt", "Events", kFALSE);
|
|---|
| 198 | write.AddContainer("IncidentAngle", "Events", kFALSE);
|
|---|
| 199 | write.AddContainer("MPointingPos", "Events", kFALSE);
|
|---|
| 200 | write.AddContainer("MTime", "Events", kFALSE);
|
|---|
| 201 | write.AddContainer("MSrcPosCam", "Events", kFALSE);
|
|---|
| 202 |
|
|---|
| 203 | // Setup the task list
|
|---|
| 204 | tlist.AddToList(&read);
|
|---|
| 205 | tlist.AddToList(&apply);
|
|---|
| 206 | tlist.AddToList(&clean);
|
|---|
| 207 | tlist.AddToList(&hcalc);
|
|---|
| 208 | tlist.AddToList(&fillC1);
|
|---|
| 209 | tlist.AddToList(&fillC2);
|
|---|
| 210 | tlist.AddToList(&fillC3);
|
|---|
| 211 | tlist.AddToList(&fillR);
|
|---|
| 212 | tlist.AddToList(&fillD1);
|
|---|
| 213 | tlist.AddToList(&fillD2);
|
|---|
| 214 | tlist.AddToList(&fillD3);
|
|---|
| 215 | tlist.AddToList(&fillD4);
|
|---|
| 216 | tlist.AddToList(&fillD5);
|
|---|
| 217 | tlist.AddToList(&fillD6);
|
|---|
| 218 | tlist.AddToList(&filla);
|
|---|
| 219 | tlist.AddToList(&write);
|
|---|
| 220 |
|
|---|
| 221 | // Run the eventloop
|
|---|
| 222 | if (!loop.Eventloop())
|
|---|
| 223 | return 3;
|
|---|
| 224 |
|
|---|
| 225 | // Check if the display was closed (deleted) by the user
|
|---|
| 226 | if (!loop.GetDisplay())
|
|---|
| 227 | return 4;
|
|---|
| 228 |
|
|---|
| 229 | // ============================================================
|
|---|
| 230 |
|
|---|
| 231 | TString fname = write.GetFileName();
|
|---|
| 232 | fname.ReplaceAll(".root", "-display.root");
|
|---|
| 233 |
|
|---|
| 234 | // Write the status display to the file
|
|---|
| 235 | TString title = "-- Image parameters [";
|
|---|
| 236 | title += datafile;
|
|---|
| 237 | title += "] --";
|
|---|
| 238 | d->SetTitle(title, kFALSE);
|
|---|
| 239 |
|
|---|
| 240 | d->SaveAsRoot(fname);
|
|---|
| 241 |
|
|---|
| 242 | return 0;
|
|---|
| 243 | }
|
|---|
