| 1 | #include <sstream>
|
|---|
| 2 | #include <iostream>
|
|---|
| 3 |
|
|---|
| 4 | #include "MLog.h"
|
|---|
| 5 | #include "MLogManip.h"
|
|---|
| 6 |
|
|---|
| 7 | #if (!defined(__CINT__) && !defined(__CLING__)) || defined(__MAKECINT__)
|
|---|
| 8 |
|
|---|
| 9 | #include "TH1F.h"
|
|---|
| 10 | #include "TFile.h"
|
|---|
| 11 | #include "TStyle.h"
|
|---|
| 12 | #include "TGraph.h"
|
|---|
| 13 | #include "TLine.h"
|
|---|
| 14 |
|
|---|
| 15 | #include "../mcore/DrsCalib.h"
|
|---|
| 16 | #include "MDrsCalibration.h"
|
|---|
| 17 | #include "MExtralgoSpline.h"
|
|---|
| 18 | #include "MSequence.h"
|
|---|
| 19 | #include "MStatusArray.h"
|
|---|
| 20 | #include "MHCamera.h"
|
|---|
| 21 | #include "MJob.h"
|
|---|
| 22 | #include "MWriteRootFile.h"
|
|---|
| 23 | #include "MHCamera.h"
|
|---|
| 24 | #include "MBadPixelsCam.h"
|
|---|
| 25 | #include "MBadPixelsPix.h"
|
|---|
| 26 | #include "MDirIter.h"
|
|---|
| 27 | #include "MTaskList.h"
|
|---|
| 28 | #include "MFDataPhrase.h"
|
|---|
| 29 | #include "MArrayF.h"
|
|---|
| 30 | #include "MBadPixelsTreat.h"
|
|---|
| 31 | #include "MCalibrateDrsTimes.h"
|
|---|
| 32 | #include "MHSectorVsTime.h"
|
|---|
| 33 | #include "MHCamEvent.h"
|
|---|
| 34 | #include "MExtractFACT.h"
|
|---|
| 35 | #include "MFillH.h"
|
|---|
| 36 | #include "MDrsCalibApply.h"
|
|---|
| 37 | #include "MGeomApply.h"
|
|---|
| 38 | #include "MContinue.h"
|
|---|
| 39 | #include "MRawFitsRead.h"
|
|---|
| 40 | #include "MReadMarsFile.h"
|
|---|
| 41 | #include "MEvtLoop.h"
|
|---|
| 42 | #include "MParList.h"
|
|---|
| 43 | #include "MStatusDisplay.h"
|
|---|
| 44 | #include "MDrsCalibrationTime.h"
|
|---|
| 45 | #include "MH3.h"
|
|---|
| 46 | #include "MGeomCamFACT.h"
|
|---|
| 47 | #include "MCalibrateFact.h"
|
|---|
| 48 | #include "MParameters.h"
|
|---|
| 49 | #include "MWriteAsciiFile.h"
|
|---|
| 50 |
|
|---|
| 51 | #endif
|
|---|
| 52 |
|
|---|
| 53 | using namespace std;
|
|---|
| 54 |
|
|---|
| 55 | /* Maybe you wanna use this macro like this:
|
|---|
| 56 | *
|
|---|
| 57 | * 0.) ---- call root ----
|
|---|
| 58 | * root -b
|
|---|
| 59 | *
|
|---|
| 60 | * 1.) ---- compile the stuff ----
|
|---|
| 61 | * .L fact/analysis/callisto_buildable_no_sequence_file.C++
|
|---|
| 62 | * <read a lot of warnings>
|
|---|
| 63 | *
|
|---|
| 64 | * 2.) ---- you can call it then ----
|
|---|
| 65 | * Therefore you need to specify all the paths ... see below.
|
|---|
| 66 | *
|
|---|
| 67 | * When you wanna call the stuff directly from the bash make sure to
|
|---|
| 68 | * escape the bracets and quotes correctly.
|
|---|
| 69 | *
|
|---|
| 70 | * your can do:
|
|---|
| 71 | * root -b -q callisto_buildable_no_sequence_file.C++'("path1","path2",...)'
|
|---|
| 72 | * or:
|
|---|
| 73 | * root -b -q callisto_buildable_no_sequence_file.C++(\"path1\",\"$HOME\",...)
|
|---|
| 74 | * using bash enviroment variables like $HOME is not possible in the upper variant.
|
|---|
| 75 | */
|
|---|
| 76 |
|
|---|
| 77 | int callisto_data(const TString drsfile="test300samples2.drs.fits.gz",
|
|---|
| 78 | const TString datfile="00000003.387_D_MonteCarlo010_Events.fits",
|
|---|
| 79 | TString outpath = "",
|
|---|
| 80 | TString displayfile = "", TString displaytitle = "")
|
|---|
| 81 | {
|
|---|
| 82 |
|
|---|
| 83 | // ======================================================
|
|---|
| 84 |
|
|---|
| 85 | if (displaytitle.IsNull())
|
|---|
| 86 | displaytitle = gSystem->BaseName(datfile);
|
|---|
| 87 |
|
|---|
| 88 | FileStat_t fstat;
|
|---|
| 89 | int rc = gSystem->GetPathInfo(outpath, fstat);
|
|---|
| 90 | bool isdir = !rc || R_ISDIR(fstat.fMode);
|
|---|
| 91 |
|
|---|
| 92 | TString filename = datfile + "_callisto.root";
|
|---|
| 93 | filename.Replace(0, filename.Last('/')+1, "");
|
|---|
| 94 | const char *buf = gSystem->ConcatFileName(outpath, filename);
|
|---|
| 95 | TString outfile = buf;
|
|---|
| 96 | delete [] buf;
|
|---|
| 97 |
|
|---|
| 98 | if (displayfile.IsNull())
|
|---|
| 99 | {
|
|---|
| 100 | displayfile = outfile;
|
|---|
| 101 | displayfile.Insert(displayfile.Last('.'), "-display");
|
|---|
| 102 | }
|
|---|
| 103 | else
|
|---|
| 104 | {
|
|---|
| 105 | if (isdir && gSystem->DirName(displayfile)==TString("."))
|
|---|
| 106 | {
|
|---|
| 107 | buf = gSystem->ConcatFileName(outfile, displayfile);
|
|---|
| 108 | displayfile = buf;
|
|---|
| 109 | delete [] buf;
|
|---|
| 110 | }
|
|---|
| 111 | }
|
|---|
| 112 |
|
|---|
| 113 | // ======================================================
|
|---|
| 114 |
|
|---|
| 115 | // true: Display correctly mapped pixels in the camera displays
|
|---|
| 116 | // but the value-vs-index plot is in software/spiral indices
|
|---|
| 117 | // false: Display pixels in hardware/linear indices,
|
|---|
| 118 | // but the order is the camera display is distorted
|
|---|
| 119 | bool usemap = false;
|
|---|
| 120 |
|
|---|
| 121 | // map file to use (get that from La Palma!)
|
|---|
| 122 | const char *pmap = usemap ? "FACTmap111030.txt" : NULL;
|
|---|
| 123 |
|
|---|
| 124 | // ------------------------------------------------------
|
|---|
| 125 |
|
|---|
| 126 | MStatusDisplay *d = new MStatusDisplay;
|
|---|
| 127 |
|
|---|
| 128 | MBadPixelsCam badpixels;
|
|---|
| 129 | badpixels.InitSize(1440);
|
|---|
| 130 | badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
|---|
| 131 | badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
|---|
| 132 | badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
|---|
| 133 | badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
|---|
| 134 | badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
|---|
| 135 | badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
|---|
| 136 | // Twin pixel
|
|---|
| 137 | // 113
|
|---|
| 138 | // 115
|
|---|
| 139 | // 354
|
|---|
| 140 | // 423
|
|---|
| 141 | // 1195
|
|---|
| 142 | // 1393
|
|---|
| 143 |
|
|---|
| 144 | // ------------------------------------------------------
|
|---|
| 145 |
|
|---|
| 146 | // Calib: 51 / 90 / 197 (20% TH)
|
|---|
| 147 | // Data: 52 / 64 / 104 (20% TH)
|
|---|
| 148 |
|
|---|
| 149 | // Extraction range in slices. It will always(!) contain the full range
|
|---|
| 150 | // of integration
|
|---|
| 151 | const int first_slice = 25; // 10ns
|
|---|
| 152 | const int last_slice = 225; // 125ns
|
|---|
| 153 |
|
|---|
| 154 | Long_t max = 0; // All
|
|---|
| 155 | Long_t max3 = max; // Pedestal Rndm
|
|---|
| 156 | Long_t max4 = max; // Pedestal Ext
|
|---|
| 157 |
|
|---|
| 158 | // ======================================================
|
|---|
| 159 |
|
|---|
| 160 | //double scale = 0.1;
|
|---|
| 161 | double scale = 0.1024; // 0.00389429
|
|---|
| 162 |
|
|---|
| 163 | // ======================================================
|
|---|
| 164 |
|
|---|
| 165 | if (pmap && gSystem->AccessPathName(pmap, kFileExists))
|
|---|
| 166 | {
|
|---|
| 167 | gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl;
|
|---|
| 168 | return 1;
|
|---|
| 169 | }
|
|---|
| 170 |
|
|---|
| 171 | gLog.Separator("Callisto");
|
|---|
| 172 | gLog << all;
|
|---|
| 173 | gLog << "Data File: " << datfile << '\n';
|
|---|
| 174 | gLog << "DRS calib 300: " << drsfile << endl;;
|
|---|
| 175 |
|
|---|
| 176 | MDrsCalibration drscalib300;
|
|---|
| 177 | if (!drscalib300.ReadFits(drsfile.Data())) {
|
|---|
| 178 | gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl;
|
|---|
| 179 | return 5;
|
|---|
| 180 | }
|
|---|
| 181 | gLog << all;
|
|---|
| 182 | gLog << "Output file: " << outfile << '\n';
|
|---|
| 183 | gLog << "Display file: " << displayfile << '\n';
|
|---|
| 184 | gLog << "Display title: " << displaytitle << endl;
|
|---|
| 185 |
|
|---|
| 186 | // ------------------------------------------------------
|
|---|
| 187 |
|
|---|
| 188 | // ======================================================
|
|---|
| 189 |
|
|---|
| 190 | // Plot the trigger pattern rates vs. run-number
|
|---|
| 191 | MH3 hrate("MRawRunHeader.GetFileID", "MRawEvtHeader.GetTriggerID&0xff00");
|
|---|
| 192 | hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
|
|---|
| 193 | hrate.SetName("Rate");
|
|---|
| 194 | hrate.SetTitle("Event rate [Hz];File Id;Trigger Type;");
|
|---|
| 195 | hrate.InitLabels(MH3::kLabelsXY);
|
|---|
| 196 | hrate.DefineLabelY( 0, "Data"); // What if TriggerID==0 already???
|
|---|
| 197 | hrate.DefineLabelY(0x100, "Cal");
|
|---|
| 198 | hrate.DefineLabelY(0x400, "Ped");
|
|---|
| 199 | // hrate.DefaultLabelY("ERROR");
|
|---|
| 200 | gStyle->SetOptFit(kTRUE);
|
|---|
| 201 |
|
|---|
| 202 |
|
|---|
| 203 | // ========================= Result ==================================
|
|---|
| 204 |
|
|---|
| 205 | gLog << endl;
|
|---|
| 206 | gLog.Separator("Extracting and calibrating data");
|
|---|
| 207 |
|
|---|
| 208 | MTaskList tlist5;
|
|---|
| 209 |
|
|---|
| 210 | MParList plist5;
|
|---|
| 211 | plist5.AddToList(&tlist5);
|
|---|
| 212 | plist5.AddToList(&drscalib300);
|
|---|
| 213 | plist5.AddToList(&badpixels);
|
|---|
| 214 |
|
|---|
| 215 | MEvtLoop loop5("CalibratingData");
|
|---|
| 216 | loop5.SetDisplay(d);
|
|---|
| 217 | loop5.SetParList(&plist5);
|
|---|
| 218 |
|
|---|
| 219 | // ------------------ Setup the tasks ---------------
|
|---|
| 220 |
|
|---|
| 221 | MRawFitsRead read5a;
|
|---|
| 222 | MReadMarsFile read5b("Events");
|
|---|
| 223 | if (pmap)
|
|---|
| 224 | read5a.LoadMap(pmap);
|
|---|
| 225 | read5a.AddFile(datfile);
|
|---|
| 226 | read5b.DisableAutoScheme();
|
|---|
| 227 | read5b.AddFile(datfile);
|
|---|
| 228 |
|
|---|
| 229 | MRead &read5 = datfile.EndsWith(".root") ? static_cast<MRead&>(read5b) : static_cast<MRead&>(read5a);
|
|---|
| 230 |
|
|---|
| 231 | MFillH fill5a(&hrate);
|
|---|
| 232 |
|
|---|
| 233 | MGeomApply apply5;
|
|---|
| 234 |
|
|---|
| 235 | MDrsCalibApply drsapply5;
|
|---|
| 236 |
|
|---|
| 237 | MTreatSaturation treatsat5;
|
|---|
| 238 |
|
|---|
| 239 | MFilterData filterdata5;
|
|---|
| 240 |
|
|---|
| 241 | MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0", "SelectDat");
|
|---|
| 242 | MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal");
|
|---|
| 243 | MFDataPhrase filterped("(MRawEvtHeader.GetTriggerID&0xff00)==0x400", "SelectPed");
|
|---|
| 244 | MFDataPhrase filterncl("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectNonCal");
|
|---|
| 245 |
|
|---|
| 246 | //MContinue cont4("MRawEvtHeader.GetTriggerID!=4", "SelectData");
|
|---|
| 247 |
|
|---|
| 248 | // ---
|
|---|
| 249 |
|
|---|
| 250 | MExtractFACT extractor5dat;
|
|---|
| 251 | extractor5dat.SetRange(first_slice, last_slice);
|
|---|
| 252 | extractor5dat.SetNoiseCalculation(kFALSE);
|
|---|
| 253 |
|
|---|
| 254 | MExtractFACT extractor5cal;
|
|---|
| 255 | extractor5cal.SetRange(first_slice, last_slice);
|
|---|
| 256 | extractor5cal.SetNoiseCalculation(kFALSE);
|
|---|
| 257 |
|
|---|
| 258 | MExtractFACT extractor5tm("ExtractTM");
|
|---|
| 259 | extractor5tm.SetRange(last_slice, 294);
|
|---|
| 260 | extractor5tm.SetNoiseCalculation(kFALSE);
|
|---|
| 261 | extractor5tm.SetNameSignalCam("TimeMarkerAmplitude");
|
|---|
| 262 | extractor5tm.SetNameTimeCam("TimeMarkerTime");
|
|---|
| 263 |
|
|---|
| 264 | extractor5dat.SetFilter(&filterncl);
|
|---|
| 265 | extractor5cal.SetFilter(&filtercal);
|
|---|
| 266 | //extractor4tm.SetFilter(&filtercal);
|
|---|
| 267 |
|
|---|
| 268 | // ---
|
|---|
| 269 | MCalibrateFact conv5;
|
|---|
| 270 | conv5.SetScale(scale);
|
|---|
| 271 | //conv5.SetCalibConst(calib);
|
|---|
| 272 |
|
|---|
| 273 | MCalibrateDrsTimes calctm5;
|
|---|
| 274 | calctm5.SetNameUncalibrated("UncalibratedTimes");
|
|---|
| 275 |
|
|---|
| 276 | MCalibrateDrsTimes calctm5tm("CalibrateTimeMarker");
|
|---|
| 277 | calctm5tm.SetNameArrivalTime("TimeMarkerTime");
|
|---|
| 278 | calctm5tm.SetNameUncalibrated("UncalTimeMarker");
|
|---|
| 279 | calctm5tm.SetNameCalibrated("TimeMarker");
|
|---|
| 280 | calctm5tm.SetTimeMarker();
|
|---|
| 281 | //calctm4tm.SetFilter(&filtercal);
|
|---|
| 282 |
|
|---|
| 283 | MBadPixelsTreat treat5;
|
|---|
| 284 | treat5.SetProcessPedestalRun(kFALSE);
|
|---|
| 285 | treat5.SetProcessPedestalEvt(kFALSE);
|
|---|
| 286 |
|
|---|
| 287 | MHCamEvent evt5b(0, "ExtSig", "Extracted signal;;S [mV#dot sl]");
|
|---|
| 288 | MHCamEvent evt5c(0, "CalSig", "Calibrated and interpolated signal;;S [~phe]");
|
|---|
| 289 | MHCamEvent evt5d(4, "ExtSigTm", "Extracted time;;T [sl]");
|
|---|
| 290 | MHCamEvent evt5e(6, "CalSigTm", "Calibrated and interpolated time;;T [ns]");
|
|---|
| 291 |
|
|---|
| 292 | MFillH fill5b(&evt5b, "MExtractedSignalCam", "FillExtSig");
|
|---|
| 293 | MFillH fill5c(&evt5c, "MSignalCam", "FillCalSig");
|
|---|
| 294 | MFillH fill5d(&evt5d, "MArrivalTimeCam", "FillExtTm");
|
|---|
| 295 | MFillH fill5e(&evt5e, "MSignalCam", "FillCalTm");
|
|---|
| 296 |
|
|---|
| 297 | fill5c.SetDrawOption("gaus");
|
|---|
| 298 | fill5d.SetDrawOption("gaus");
|
|---|
| 299 | fill5e.SetDrawOption("gaus");
|
|---|
| 300 |
|
|---|
| 301 | /*
|
|---|
| 302 | fill4b.SetFilter(&filterdat);
|
|---|
| 303 | fill4c.SetFilter(&filterdat);
|
|---|
| 304 | fill4d.SetFilter(&filterdat);
|
|---|
| 305 | fill4e.SetFilter(&filterdat);
|
|---|
| 306 | */
|
|---|
| 307 |
|
|---|
| 308 | //MFSoftwareTrigger swtrig;
|
|---|
| 309 | //MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger");
|
|---|
| 310 | //contsw.SetInverted();
|
|---|
| 311 |
|
|---|
| 312 | // The second rule is for the case reading raw-files!
|
|---|
| 313 |
|
|---|
| 314 | MWriteRootFile write5(outfile, "RECREATE", "Calibrated Data", 2);
|
|---|
| 315 | write5.AddContainer("MRawRunHeader", "RunHeaders");
|
|---|
| 316 | write5.AddContainer("MGeomCam", "RunHeaders");
|
|---|
| 317 | write5.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
|
|---|
| 318 | write5.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
|
|---|
| 319 | write5.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
|
|---|
| 320 |
|
|---|
| 321 | // Common events
|
|---|
| 322 | write5.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
|
|---|
| 323 | write5.AddContainer("MMcEvt", "Events", kFALSE);
|
|---|
| 324 | write5.AddContainer("IncidentAngle", "Events", kFALSE);
|
|---|
| 325 | write5.AddContainer("MPointingPos", "Events", kFALSE);
|
|---|
| 326 | write5.AddContainer("MSignalCam", "Events");
|
|---|
| 327 | write5.AddContainer("MTime", "Events", kFALSE);
|
|---|
| 328 | write5.AddContainer("MRawEvtHeader", "Events");
|
|---|
| 329 | //write.AddContainer("MTriggerPattern", "Events");
|
|---|
| 330 |
|
|---|
| 331 | // ------------------ Setup histograms and fill tasks ----------------
|
|---|
| 332 |
|
|---|
| 333 | MContinue test;
|
|---|
| 334 | test.SetFilter(&filterncl);
|
|---|
| 335 |
|
|---|
| 336 | MTaskList tlist5tm;
|
|---|
| 337 | tlist5tm.AddToList(&extractor5tm);
|
|---|
| 338 | tlist5tm.AddToList(&calctm5tm);
|
|---|
| 339 | tlist5tm.SetFilter(&filtercal);
|
|---|
| 340 |
|
|---|
| 341 | MTaskList tlist5dat;
|
|---|
| 342 | tlist5dat.AddToList(&fill5b);
|
|---|
| 343 | tlist5dat.AddToList(&fill5c);
|
|---|
| 344 | tlist5dat.AddToList(&fill5d);
|
|---|
| 345 | tlist5dat.AddToList(&fill5e);
|
|---|
| 346 | tlist5dat.SetFilter(&filterdat);
|
|---|
| 347 |
|
|---|
| 348 | tlist5.AddToList(&read5);
|
|---|
| 349 | tlist5.AddToList(&apply5);
|
|---|
| 350 | tlist5.AddToList(&drsapply5);
|
|---|
| 351 | tlist5.AddToList(&filterncl);
|
|---|
| 352 | //tlist5.AddToList(&test);
|
|---|
| 353 | tlist5.AddToList(&filterdat);
|
|---|
| 354 | tlist5.AddToList(&filtercal);
|
|---|
| 355 | tlist5.AddToList(&filterped);
|
|---|
| 356 | tlist5.AddToList(&fill5a);
|
|---|
| 357 | tlist5.AddToList(&treatsat5);
|
|---|
| 358 | tlist5.AddToList(&filterdata5);
|
|---|
| 359 | tlist5.AddToList(&extractor5dat);
|
|---|
| 360 | tlist5.AddToList(&extractor5cal);
|
|---|
| 361 | tlist5.AddToList(&calctm5);
|
|---|
| 362 | tlist5.AddToList(&tlist5tm);
|
|---|
| 363 | tlist5.AddToList(&conv5);
|
|---|
| 364 | tlist5.AddToList(&treat5);
|
|---|
| 365 | tlist5.AddToList(&tlist5dat);
|
|---|
| 366 | tlist5.AddToList(&write5);
|
|---|
| 367 |
|
|---|
| 368 | if (!loop5.Eventloop(max4))
|
|---|
| 369 | return 18;
|
|---|
| 370 |
|
|---|
| 371 | if (!loop5.GetDisplay())
|
|---|
| 372 | return 19;
|
|---|
| 373 |
|
|---|
| 374 | d->SetTitle(displaytitle, kFALSE);
|
|---|
| 375 | d->SaveAs(displayfile);
|
|---|
| 376 |
|
|---|
| 377 | return 0;
|
|---|
| 378 | }
|
|---|
