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