| 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, Double_t lvl1=7.8, Double_t lvl2=3.9, 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.root", "_I.root");
|
|---|
| 62 |
|
|---|
| 63 | gLog.Separator();
|
|---|
| 64 |
|
|---|
| 65 | // ---------------------------------------------------------
|
|---|
| 66 |
|
|---|
| 67 | // Allocate/open the status display
|
|---|
| 68 | MStatusDisplay *d = new MStatusDisplay;
|
|---|
| 69 |
|
|---|
| 70 | // Instantiate the list of tasks
|
|---|
| 71 | MTaskList tlist;
|
|---|
| 72 |
|
|---|
| 73 | // Instantiate the list of parameter containers
|
|---|
| 74 | MParList plist2;
|
|---|
| 75 | plist2.AddToList(&tlist);
|
|---|
| 76 |
|
|---|
| 77 | // Fill it with the binnings
|
|---|
| 78 | plist2.AddToList(&binse);
|
|---|
| 79 | plist2.AddToList(&bins1);
|
|---|
| 80 | plist2.AddToList(&bins3);
|
|---|
| 81 | plist2.AddToList(&bins4);
|
|---|
| 82 | plist2.AddToList(&bins5);
|
|---|
| 83 | plist2.AddToList(&bins6);
|
|---|
| 84 | plist2.AddToList(&bins7);
|
|---|
| 85 | plist2.AddToList(&bins8);
|
|---|
| 86 | plist2.AddToList(&bins9);
|
|---|
| 87 | plist2.AddToList(&binsA);
|
|---|
| 88 |
|
|---|
| 89 | // Instatiate the event loop
|
|---|
| 90 | MEvtLoop loop;
|
|---|
| 91 | loop.SetDisplay(d);
|
|---|
| 92 | loop.SetParList(&plist2);
|
|---|
| 93 |
|
|---|
| 94 | // Instantiate the reading task
|
|---|
| 95 | // You can use
|
|---|
| 96 | // read.AddFiles("*_C.root")
|
|---|
| 97 | // for example to read more than one file at once
|
|---|
| 98 | MReadMarsFile read("Events");
|
|---|
| 99 | read.DisableAutoScheme();
|
|---|
| 100 | read.AddFile(datafile);
|
|---|
| 101 |
|
|---|
| 102 | // Instantiate the task which takes care of the size of all containers
|
|---|
| 103 | MGeomApply apply;
|
|---|
| 104 |
|
|---|
| 105 | // Instantiate the image cleaning
|
|---|
| 106 |
|
|---|
| 107 | // This is the most simple image cleaning possible, standard
|
|---|
| 108 | // tail cut based on absolute amplitudes
|
|---|
| 109 | // MImgCleanStd clean(lvl1, lvl2);
|
|---|
| 110 | // clean.SetMethod(MImgCleanStd::kAbsolute);
|
|---|
| 111 | // clean.SetPostCleanType(3);
|
|---|
| 112 |
|
|---|
| 113 | // Instantiate the image cleaning as described in the TRAC
|
|---|
| 114 | MImgCleanTime clean;
|
|---|
| 115 | clean.SetMinCount(0);
|
|---|
| 116 | clean.SetMinSize(50);
|
|---|
| 117 | clean.SetDeltaT(17.5);
|
|---|
| 118 |
|
|---|
| 119 | // Instantiate the calculation of the image parameters
|
|---|
| 120 | MHillasCalc hcalc;
|
|---|
| 121 | hcalc.Disable(MHillasCalc::kCalcConc);
|
|---|
| 122 |
|
|---|
| 123 | // Instantiate a plot for the rate after cleaning
|
|---|
| 124 | MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
|
|---|
| 125 | hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1.)");
|
|---|
| 126 | hrate.SetName("Rate");
|
|---|
| 127 | hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];");
|
|---|
| 128 | hrate.InitLabels(MH3::kLabelsX);
|
|---|
| 129 | hrate.GetHist().SetMinimum(0);
|
|---|
| 130 |
|
|---|
| 131 | // Instantiate the task to fill the rate plot
|
|---|
| 132 | MFillH fillR(&hrate, "", "FillRate");
|
|---|
| 133 |
|
|---|
| 134 | // Instantiate the camera displays to be filled
|
|---|
| 135 | MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]");
|
|---|
| 136 | MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
|
|---|
| 137 | MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]");
|
|---|
| 138 | evtC1.SetErrorSpread(kFALSE);
|
|---|
| 139 | evtC2.SetErrorSpread(kFALSE);
|
|---|
| 140 | evtC2.SetThreshold(0);
|
|---|
| 141 |
|
|---|
| 142 | // Instantiate tasks to fill the camera displays
|
|---|
| 143 | MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam");
|
|---|
| 144 | MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels");
|
|---|
| 145 | MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos");
|
|---|
| 146 |
|
|---|
| 147 | // Instantiate tasks to fill image parameter histograms
|
|---|
| 148 | MFillH fillD1("MHHillas", "MHillas", "FillHillas");
|
|---|
| 149 | MFillH fillD2("MHHillasExt", "", "FillHillasExt");
|
|---|
| 150 | MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
|
|---|
| 151 | MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar");
|
|---|
| 152 | MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
|
|---|
| 153 | MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize");
|
|---|
| 154 |
|
|---|
| 155 | // Instantiate an example for a user defined plot
|
|---|
| 156 | MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)");
|
|---|
| 157 | halpha.SetName("AvsD;Dist;Alpha");
|
|---|
| 158 |
|
|---|
| 159 | // Instantiate a task to fill that plot
|
|---|
| 160 | MFillH filla(&halpha);
|
|---|
| 161 | filla.SetNameTab("AvsD");
|
|---|
| 162 | filla.SetDrawOption("colz");
|
|---|
| 163 |
|
|---|
| 164 | // This is an alternative if more than one file should be read in a single loop
|
|---|
| 165 | // In this case the output file name of each file is created with this rule
|
|---|
| 166 | // from the input file name. Note that the filename of the output file
|
|---|
| 167 | // with the status display must then be explicitly given.
|
|---|
| 168 | //
|
|---|
| 169 | //const TString fname(Form("s/([0-9]+_[0-9]+)[.]_C.root$/%s\\/$1_I.root/",
|
|---|
| 170 | // MJob::Esc(outpath).Data()));
|
|---|
| 171 | //MWriteRootFile write5(2, fname, "RECREATE", "Image parameters");
|
|---|
| 172 |
|
|---|
| 173 | // Instantiate writing the file
|
|---|
| 174 | MWriteRootFile write(fname, "RECREATE", "Image parameters", 2);
|
|---|
| 175 | write.AddContainer("MTime", "Events");
|
|---|
| 176 | write.AddContainer("MHillas", "Events");
|
|---|
| 177 | write.AddContainer("MHillasExt", "Events");
|
|---|
| 178 | write.AddContainer("MHillasSrc", "Events");
|
|---|
| 179 | write.AddContainer("MImagePar", "Events");
|
|---|
| 180 | write.AddContainer("MNewImagePar", "Events");
|
|---|
| 181 | write.AddContainer("MRawEvtHeader", "Events");
|
|---|
| 182 | write.AddContainer("MSoftwareTrigger","Events");
|
|---|
| 183 | write.AddContainer("MRawRunHeader", "RunHeaders");
|
|---|
| 184 | write.AddContainer("MGeomCam", "RunHeaders");
|
|---|
| 185 |
|
|---|
| 186 | write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
|
|---|
| 187 | write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
|
|---|
| 188 | write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
|
|---|
| 189 | write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
|
|---|
| 190 | write.AddContainer("MMcEvt", "Events", kFALSE);
|
|---|
| 191 | write.AddContainer("IncidentAngle", "Events", kFALSE);
|
|---|
| 192 | write.AddContainer("MPointingPos", "Events", kFALSE);
|
|---|
| 193 | write.AddContainer("MTime", "Events", kFALSE);
|
|---|
| 194 | write.AddContainer("MSrcPosCam", "Events", kFALSE);
|
|---|
| 195 |
|
|---|
| 196 | // Setup the task list
|
|---|
| 197 | tlist.AddToList(&read);
|
|---|
| 198 | tlist.AddToList(&apply);
|
|---|
| 199 | tlist.AddToList(&clean);
|
|---|
| 200 | tlist.AddToList(&hcalc);
|
|---|
| 201 | tlist.AddToList(&fillC1);
|
|---|
| 202 | tlist.AddToList(&fillC2);
|
|---|
| 203 | tlist.AddToList(&fillC3);
|
|---|
| 204 | tlist.AddToList(&fillR);
|
|---|
| 205 | tlist.AddToList(&fillD1);
|
|---|
| 206 | tlist.AddToList(&fillD2);
|
|---|
| 207 | tlist.AddToList(&fillD3);
|
|---|
| 208 | tlist.AddToList(&fillD4);
|
|---|
| 209 | tlist.AddToList(&fillD5);
|
|---|
| 210 | tlist.AddToList(&fillD6);
|
|---|
| 211 | tlist.AddToList(&filla);
|
|---|
| 212 | tlist.AddToList(&write);
|
|---|
| 213 |
|
|---|
| 214 | // Run the eventloop
|
|---|
| 215 | if (!loop.Eventloop())
|
|---|
| 216 | return 3;
|
|---|
| 217 |
|
|---|
| 218 | // Check if the display was closed (deleted) by the user
|
|---|
| 219 | if (!loop.GetDisplay())
|
|---|
| 220 | return 4;
|
|---|
| 221 |
|
|---|
| 222 | // ============================================================
|
|---|
| 223 |
|
|---|
| 224 | // Check if the output file is still accessible from root
|
|---|
| 225 | TFile *ofile = (TFile*)gROOT->GetListOfFiles()->FindObject(fname);
|
|---|
| 226 | if (!ofile || !ofile->IsOpen() || ofile->IsZombie())
|
|---|
| 227 | {
|
|---|
| 228 | gLog << err << "File " << fname << " not found" << endl;
|
|---|
| 229 | return 20;
|
|---|
| 230 | }
|
|---|
| 231 |
|
|---|
| 232 | // Write the status display to the file
|
|---|
| 233 | TString title = "-- Image parameters [";
|
|---|
| 234 | title += datafile;
|
|---|
| 235 | title += "] --";
|
|---|
| 236 | d->SetTitle(title, kFALSE);
|
|---|
| 237 |
|
|---|
| 238 | ofile->cd();
|
|---|
| 239 | d->Write();
|
|---|
| 240 |
|
|---|
| 241 | return 0;
|
|---|
| 242 | }
|
|---|
