1 | #include <sstream>
|
---|
2 | #include <iostream>
|
---|
3 |
|
---|
4 | #include "MLog.h"
|
---|
5 | #include "MLogManip.h"
|
---|
6 |
|
---|
7 | #if !defined(__CINT__) || 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 "MExtractTimeAndChargeSpline.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(const TString drsfile="test300samples.drs.fits",
|
---|
78 | const TString pedfile="00000001.001_P_MonteCarlo000_Events.fits",
|
---|
79 | const TString datfile="00000003.387_D_MonteCarlo010_Events.fits",
|
---|
80 | TString outfile = "",
|
---|
81 | TString displayfile = "", TString displaytitle = "")
|
---|
82 | {
|
---|
83 |
|
---|
84 | // ======================================================
|
---|
85 |
|
---|
86 | if (displaytitle.IsNull())
|
---|
87 | displaytitle = gSystem->BaseName(datfile);
|
---|
88 |
|
---|
89 | FileStat_t fstat;
|
---|
90 | int rc = gSystem->GetPathInfo(outfile, fstat);
|
---|
91 | bool isdir = !rc || R_ISDIR(fstat.fMode);
|
---|
92 |
|
---|
93 | const char *buf = gSystem->ConcatFileName(outfile, "callisto.root");
|
---|
94 | outfile = buf;
|
---|
95 | delete [] buf;
|
---|
96 |
|
---|
97 | if (displayfile.IsNull())
|
---|
98 | {
|
---|
99 | displayfile = outfile;
|
---|
100 | displayfile.Insert(displayfile.Last('.'), "-display");
|
---|
101 | }
|
---|
102 | else
|
---|
103 | {
|
---|
104 | if (isdir && gSystem->DirName(displayfile)==TString("."))
|
---|
105 | {
|
---|
106 | buf = gSystem->ConcatFileName(outfile, displayfile);
|
---|
107 | displayfile = buf;
|
---|
108 | delete [] buf;
|
---|
109 | }
|
---|
110 | }
|
---|
111 |
|
---|
112 | // ======================================================
|
---|
113 |
|
---|
114 | // true: Display correctly mapped pixels in the camera displays
|
---|
115 | // but the value-vs-index plot is in software/spiral indices
|
---|
116 | // false: Display pixels in hardware/linear indices,
|
---|
117 | // but the order is the camera display is distorted
|
---|
118 | bool usemap = true;
|
---|
119 |
|
---|
120 | // map file to use (get that from La Palma!)
|
---|
121 | const char *pmap = usemap ? "TestForThomas/FACT/FACTmap111030.txt" : NULL;
|
---|
122 |
|
---|
123 | Bool_t maximum = kTRUE;
|
---|
124 |
|
---|
125 | //const char *lp_template = maximum ?
|
---|
126 | // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-maximum.root" :
|
---|
127 | // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-leading-edge.root";
|
---|
128 |
|
---|
129 | const char *pulse_template = "TestForThomas/FACT/template-pulse.root";
|
---|
130 |
|
---|
131 | // ------------------------------------------------------
|
---|
132 |
|
---|
133 | // Calib: 51 / 90 / 197 (20% TH)
|
---|
134 | // Data: 52 / 64 / 104 (20% TH)
|
---|
135 |
|
---|
136 | // Extraction range in slices. It will always(!) contain the full range
|
---|
137 | // of integration
|
---|
138 | const int first_slice = 20; // 10ns
|
---|
139 | const int last_slice = 250; // 125ns
|
---|
140 |
|
---|
141 | // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel:
|
---|
142 | //
|
---|
143 | // IntegralFixed:
|
---|
144 | // * fRiseTime: Number of slices left from arrival time
|
---|
145 | // * fFallTime: Number of slices right from arrival time
|
---|
146 | // IntegralRel:
|
---|
147 | // * fRiseTime: Number of slices left from maximum time
|
---|
148 | // * fFallTime: Number of slices right from maximum time
|
---|
149 | //
|
---|
150 | const int rise_time_cal = maximum ? 40 : 10; // was 13; 5ns
|
---|
151 | const int fall_time_cal = maximum ? 120 : 160; // was 23; 80ns
|
---|
152 |
|
---|
153 | const int rise_time_dat = maximum ? 10 : 2; // was 13; was 10; 1ns
|
---|
154 | const int fall_time_dat = maximum ? 40 : 48; // was 23; was 40; 24ns
|
---|
155 |
|
---|
156 | // Extraction type: Extract integral and half leading edge
|
---|
157 |
|
---|
158 | const MExtralgoSpline::ExtractionType_t type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed);
|
---|
159 | //const int type = MExtralgoSpline::kIntegralFixed;
|
---|
160 |
|
---|
161 |
|
---|
162 | const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 }
|
---|
163 |
|
---|
164 | Long_t max = 0; // All
|
---|
165 | Long_t max3 = max; // Pedestal Rndm
|
---|
166 | Long_t max4 = max; // Pedestal Ext
|
---|
167 |
|
---|
168 | // ======================================================
|
---|
169 |
|
---|
170 | if (pmap && gSystem->AccessPathName(pmap, kFileExists))
|
---|
171 | {
|
---|
172 | gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl;
|
---|
173 | return 1;
|
---|
174 | }
|
---|
175 |
|
---|
176 | gLog.Separator("Callisto");
|
---|
177 | gLog << all;
|
---|
178 | gLog << "Data File: " << datfile << '\n';
|
---|
179 | gLog << "DRS calib 300: " << drsfile << endl;;
|
---|
180 |
|
---|
181 | MDrsCalibration drscalib300;
|
---|
182 | if (!drscalib300.ReadFits(drsfile.Data())) {
|
---|
183 | gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl;
|
---|
184 | return 5;
|
---|
185 | }
|
---|
186 | gLog << all;
|
---|
187 | gLog << "Pedestal file: " << pedfile << '\n';
|
---|
188 | gLog << "Output file: " << outfile << '\n';
|
---|
189 | gLog << "Display file: " << displayfile << '\n';
|
---|
190 | gLog << "Display title: " << displaytitle << endl;
|
---|
191 |
|
---|
192 | // ------------------------------------------------------
|
---|
193 | MStatusArray arrt, arrp;
|
---|
194 |
|
---|
195 | // TFile ft(lp_template);
|
---|
196 | // if (arrt.Read()<=0)
|
---|
197 | // {
|
---|
198 | // gLog << err << "ERROR - Reading LP template from " << lp_template << endl;
|
---|
199 | // return 100;
|
---|
200 | // }
|
---|
201 |
|
---|
202 | // MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam");
|
---|
203 | // if (!lpref)
|
---|
204 | // {
|
---|
205 | // gLog << err << "ERROR - LP Template not found in " << lp_template << endl;
|
---|
206 | // return 101;
|
---|
207 | // }
|
---|
208 | // lpref->SetDirectory(0);
|
---|
209 |
|
---|
210 | // MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain");
|
---|
211 | // if (!gain)
|
---|
212 | // {
|
---|
213 | // gLog << err << "ERROR - Gain not found in " << lp_template << endl;
|
---|
214 | // return 101;
|
---|
215 | // }
|
---|
216 | // gain->SetDirectory(0);
|
---|
217 |
|
---|
218 | TFile fp(pulse_template);
|
---|
219 | if (arrp.Read()<=0)
|
---|
220 | {
|
---|
221 | gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;
|
---|
222 | return 102;
|
---|
223 | }
|
---|
224 |
|
---|
225 | TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0");
|
---|
226 | if (!hpulse)
|
---|
227 | {
|
---|
228 | gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;
|
---|
229 | return 103;
|
---|
230 | }
|
---|
231 | hpulse->SetDirectory(0);
|
---|
232 | // ======================================================
|
---|
233 |
|
---|
234 | MStatusDisplay *d = new MStatusDisplay;
|
---|
235 |
|
---|
236 | MBadPixelsCam badpixels;
|
---|
237 | badpixels.InitSize(1440);
|
---|
238 | /*
|
---|
239 | badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
240 | badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
241 | badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
242 | badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
243 | badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
244 | badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
245 | */
|
---|
246 | // Twin pixel
|
---|
247 | // 113
|
---|
248 | // 115
|
---|
249 | // 354
|
---|
250 | // 423
|
---|
251 | // 1195
|
---|
252 | // 1393
|
---|
253 |
|
---|
254 | MDrsCalibrationTime timecam;
|
---|
255 |
|
---|
256 | // Plot the trigger pattern rates vs. run-number
|
---|
257 | MH3 hrate("MRawRunHeader.GetFileID", "MRawEvtHeader.GetTriggerID&0xff00");
|
---|
258 | hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
|
---|
259 | hrate.SetName("Rate");
|
---|
260 | hrate.SetTitle("Event rate [Hz];File Id;Trigger Type;");
|
---|
261 | hrate.InitLabels(MH3::kLabelsXY);
|
---|
262 | hrate.DefineLabelY( 0, "Data"); // What if TriggerID==0 already???
|
---|
263 | hrate.DefineLabelY(0x100, "Cal");
|
---|
264 | hrate.DefineLabelY(0x400, "Ped");
|
---|
265 | // hrate.DefaultLabelY("ERROR");
|
---|
266 | gStyle->SetOptFit(kTRUE);
|
---|
267 |
|
---|
268 |
|
---|
269 | // ========================= Result ==================================
|
---|
270 |
|
---|
271 | //~ Double_t avgS = evt1f.GetHist()->GetMean();
|
---|
272 | //~ Double_t medS = evt1f.GetHist()->GetMedian();
|
---|
273 | //~ Double_t rmsS = evt1f.GetHist()->GetRMS();
|
---|
274 | //~ Double_t maxS = evt1f.GetHist()->GetMaximum();
|
---|
275 |
|
---|
276 | MArrayF der1(hpulse->GetNbinsX());
|
---|
277 | MArrayF der2(hpulse->GetNbinsX());
|
---|
278 |
|
---|
279 | MExtralgoSpline spline(hpulse->GetArray()+1, hpulse->GetNbinsX(),
|
---|
280 | der1.GetArray(), der2.GetArray());
|
---|
281 | spline.SetRiseFallTime(rise_time_dat, fall_time_dat);
|
---|
282 | spline.SetExtractionType(type);
|
---|
283 | spline.SetHeightTm(heighttm);
|
---|
284 |
|
---|
285 | spline.Extract(hpulse->GetMaximumBin()-1);
|
---|
286 |
|
---|
287 | // The pulser signal is most probably around 400mV/9.5mV
|
---|
288 | // IntegraFixed 2/48 corresponds to roughly 215mV*50slices
|
---|
289 | Double_t scale = 1./spline.GetSignal();
|
---|
290 |
|
---|
291 | MArrayD calib(1440);
|
---|
292 | for (int i=0; i<1440; i++)
|
---|
293 | calib[i] =1.;
|
---|
294 |
|
---|
295 | gROOT->SetSelectedPad(0);
|
---|
296 | d->AddTab("PulseTemp");
|
---|
297 | gPad->SetGrid();
|
---|
298 | hpulse->SetNameTitle("Pulse", "Single p.e. pulse template");
|
---|
299 | hpulse->SetDirectory(0);
|
---|
300 | hpulse->SetLineColor(kBlack);
|
---|
301 | hpulse->DrawCopy();
|
---|
302 |
|
---|
303 | TAxis *ax = hpulse->GetXaxis();
|
---|
304 |
|
---|
305 | Double_t w = hpulse->GetBinWidth(1);
|
---|
306 | Double_t T = w*(spline.GetTime()+0.5) +ax->GetXmin();
|
---|
307 | //~ Double_t H = w*(hpulse->GetMaximumBin()+0.5)+ax->GetXmin();
|
---|
308 |
|
---|
309 | TLine line;
|
---|
310 | line.SetLineColor(kRed);
|
---|
311 | line.DrawLine(T-rise_time_dat*w, spline.GetHeight(),
|
---|
312 | T+fall_time_dat*w, spline.GetHeight());
|
---|
313 | line.DrawLine(T, spline.GetHeight()/4, T, 3*spline.GetHeight()/4);
|
---|
314 | line.DrawLine(T-rise_time_dat*w, 0,
|
---|
315 | T-rise_time_dat*w, spline.GetHeight());
|
---|
316 | line.DrawLine(T+fall_time_dat*w, 0,
|
---|
317 | T+fall_time_dat*w, spline.GetHeight());
|
---|
318 |
|
---|
319 | TGraph gg;
|
---|
320 | for (int ix=1; ix<=hpulse->GetNbinsX(); ix++)
|
---|
321 | for (int i=0; i<10; i++)
|
---|
322 | {
|
---|
323 | Double_t x = hpulse->GetBinLowEdge(ix)+i*hpulse->GetBinWidth(ix)/10.;
|
---|
324 | gg.SetPoint(gg.GetN(), x+w/2, spline.EvalAt(ix-1+i/10.));
|
---|
325 | }
|
---|
326 |
|
---|
327 | gg.SetLineColor(kBlue);
|
---|
328 | gg.SetMarkerColor(kBlue);
|
---|
329 | gg.SetMarkerStyle(kFullDotMedium);
|
---|
330 | gg.DrawClone("L");
|
---|
331 |
|
---|
332 | gROOT->SetSelectedPad(0);
|
---|
333 | d->AddTab("CalConst");
|
---|
334 | MGeomCamFACT fact;
|
---|
335 | MHCamera hcalco(fact);
|
---|
336 | hcalco.SetName("CalConst");
|
---|
337 | hcalco.SetTitle(Form("Relative calibration constant [%.0f/pe]", 1./scale));
|
---|
338 | hcalco.SetCamContent(calib);
|
---|
339 | hcalco.SetAllUsed();
|
---|
340 | //hcalco.Scale(scale);
|
---|
341 | hcalco.DrawCopy();
|
---|
342 |
|
---|
343 | // ======================================================
|
---|
344 |
|
---|
345 | gLog << endl;
|
---|
346 | gLog.Separator("Extracting random pedestal");
|
---|
347 |
|
---|
348 | MTaskList tlist3;
|
---|
349 |
|
---|
350 | MParList plist3;
|
---|
351 | plist3.AddToList(&tlist3);
|
---|
352 | plist3.AddToList(&drscalib300);
|
---|
353 | plist3.AddToList(&badpixels);
|
---|
354 | plist3.AddToList(&timecam);
|
---|
355 |
|
---|
356 | MEvtLoop loop3("DetermineRndmPed");
|
---|
357 | loop3.SetDisplay(d);
|
---|
358 | loop3.SetParList(&plist3);
|
---|
359 |
|
---|
360 | // ------------------ Setup the tasks ---------------
|
---|
361 |
|
---|
362 | MRawFitsRead read3;
|
---|
363 | read3.LoadMap(pmap);
|
---|
364 | read3.AddFile(pedfile);
|
---|
365 |
|
---|
366 | MFillH fill3a(&hrate);
|
---|
367 |
|
---|
368 | MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed");
|
---|
369 |
|
---|
370 | MGeomApply apply3;
|
---|
371 |
|
---|
372 | MDrsCalibApply drsapply3;
|
---|
373 |
|
---|
374 | //---
|
---|
375 |
|
---|
376 | MExtractTimeAndChargeSpline extractor3;
|
---|
377 | extractor3.SetRange(first_slice, last_slice);
|
---|
378 | extractor3.SetRiseTimeHiGain(rise_time_dat);
|
---|
379 | extractor3.SetFallTimeHiGain(fall_time_dat);
|
---|
380 | extractor3.SetHeightTm(heighttm);
|
---|
381 | extractor3.SetChargeType(type);
|
---|
382 | extractor3.SetSaturationLimit(600000);
|
---|
383 | extractor3.SetNoiseCalculation(kTRUE);
|
---|
384 |
|
---|
385 | // MHCamEvent evt2a(0, "PedRdm", "Extracted Pedestal Signal;;S");
|
---|
386 |
|
---|
387 | // MFillH fill2a(&evt2a, "MExtractedSignalCam", "FillPedRndm");
|
---|
388 |
|
---|
389 | // Use this for data, but not for calibration events
|
---|
390 | // evt2a.SetErrorSpread(kFALSE);
|
---|
391 |
|
---|
392 | /*
|
---|
393 | MCalibrateData conv3;
|
---|
394 | conv3.SetCalibrationMode(MCalibrateData::kNone);
|
---|
395 | conv3.SetPedestalFlag(MCalibrateData::kNo);
|
---|
396 | conv3.SetCalibConvMinLimit(0);
|
---|
397 | conv3.SetCalibConvMaxLimit(10000);
|
---|
398 | conv3.SetScaleFactor(scale);
|
---|
399 | */
|
---|
400 |
|
---|
401 | MCalibrateFact conv3;
|
---|
402 | conv3.SetScale(scale);
|
---|
403 | conv3.SetCalibConst(calib);
|
---|
404 |
|
---|
405 | MBadPixelsTreat treat3;
|
---|
406 | treat3.SetProcessPedestalRun(kFALSE);
|
---|
407 | treat3.SetProcessPedestalEvt(kFALSE);
|
---|
408 | treat3.SetProcessTimes(kFALSE);
|
---|
409 |
|
---|
410 | MHCamEvent evt3b(0, "PedRdm","Interpolated random pedestal;;Signal [~phe]");
|
---|
411 | //evt2b.SetErrorSpread(kFALSE);
|
---|
412 |
|
---|
413 | MFillH fill3b(&evt3b, "MSignalCam", "FillPedRdm");
|
---|
414 | fill3b.SetDrawOption("gaus");
|
---|
415 |
|
---|
416 | // ------------------ Setup eventloop and run analysis ---------------
|
---|
417 |
|
---|
418 | tlist3.AddToList(&read3);
|
---|
419 | tlist3.AddToList(&apply3);
|
---|
420 | tlist3.AddToList(&drsapply3);
|
---|
421 | tlist3.AddToList(&cont3);
|
---|
422 | tlist3.AddToList(&extractor3);
|
---|
423 | // tlist3.AddToList(&fill3a);
|
---|
424 | tlist3.AddToList(&conv3);
|
---|
425 | tlist3.AddToList(&treat3);
|
---|
426 | tlist3.AddToList(&fill3b);
|
---|
427 |
|
---|
428 | if (!loop3.Eventloop(max3))
|
---|
429 | return 14;
|
---|
430 |
|
---|
431 | if (!loop3.GetDisplay())
|
---|
432 | return 15;
|
---|
433 |
|
---|
434 | // ======================================================
|
---|
435 |
|
---|
436 | gLog << endl;
|
---|
437 | gLog.Separator("Extracting pedestal");
|
---|
438 |
|
---|
439 | MTaskList tlist4;
|
---|
440 |
|
---|
441 | MParList plist4;
|
---|
442 | plist4.AddToList(&tlist4);
|
---|
443 | plist4.AddToList(&drscalib300);
|
---|
444 | plist4.AddToList(&badpixels);
|
---|
445 | plist4.AddToList(&timecam);
|
---|
446 |
|
---|
447 | MEvtLoop loop4("DetermineExtractedPed");
|
---|
448 | loop4.SetDisplay(d);
|
---|
449 | loop4.SetParList(&plist4);
|
---|
450 |
|
---|
451 | // ------------------ Setup the tasks ---------------
|
---|
452 |
|
---|
453 | MRawFitsRead read4;
|
---|
454 | read4.LoadMap(pmap);
|
---|
455 | read4.AddFile(pedfile);
|
---|
456 |
|
---|
457 | MContinue cont4("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed");
|
---|
458 |
|
---|
459 | MGeomApply apply4;
|
---|
460 |
|
---|
461 | MDrsCalibApply drsapply4;
|
---|
462 |
|
---|
463 | MExtractTimeAndChargeSpline extractor4;
|
---|
464 | extractor4.SetRange(first_slice, last_slice);
|
---|
465 | extractor4.SetRiseTimeHiGain(rise_time_dat);
|
---|
466 | extractor4.SetFallTimeHiGain(fall_time_dat);
|
---|
467 | extractor4.SetHeightTm(heighttm);
|
---|
468 | extractor4.SetChargeType(type);
|
---|
469 | extractor4.SetSaturationLimit(600000);
|
---|
470 | extractor4.SetNoiseCalculation(kFALSE);
|
---|
471 |
|
---|
472 | // MHCamEvent evt3a(0, "PedExt", "Extracted Pedestal Signal;;S");
|
---|
473 |
|
---|
474 | // MFillH fill3a(&evt3a, "MExtractedSignalCam", "FillPedExt");
|
---|
475 |
|
---|
476 | // Use this for data, but not for calibration events
|
---|
477 | // evt3a.SetErrorSpread(kFALSE);
|
---|
478 | /*
|
---|
479 | MCalibrateData conv4;
|
---|
480 | conv4.SetCalibrationMode(MCalibrateData::kNone);
|
---|
481 | conv4.SetPedestalFlag(MCalibrateData::kNo);
|
---|
482 | conv4.SetCalibConvMinLimit(0);
|
---|
483 | conv4.SetCalibConvMaxLimit(10000);
|
---|
484 | conv4.SetScaleFactor(scale);
|
---|
485 | */
|
---|
486 | MCalibrateFact conv4;
|
---|
487 | conv4.SetScale(scale);
|
---|
488 | conv4.SetCalibConst(calib);
|
---|
489 |
|
---|
490 | MBadPixelsTreat treat4;
|
---|
491 | treat4.SetProcessPedestalRun(kFALSE);
|
---|
492 | treat4.SetProcessPedestalEvt(kFALSE);
|
---|
493 |
|
---|
494 | MHCamEvent evt4b(0, "PedExt","Interpolated extracted pedestal;;Signal [~phe]");
|
---|
495 | //evt4b.SetErrorSpread(kFALSE);
|
---|
496 |
|
---|
497 | MFillH fill4b(&evt4b, "MSignalCam", "FillPedExt");
|
---|
498 | fill4b.SetDrawOption("gaus");
|
---|
499 |
|
---|
500 | // ------------------ Setup eventloop and run analysis ---------------
|
---|
501 |
|
---|
502 | tlist4.AddToList(&read4);
|
---|
503 | tlist4.AddToList(&apply4);
|
---|
504 | tlist4.AddToList(&drsapply4);
|
---|
505 | tlist4.AddToList(&cont4);
|
---|
506 | tlist4.AddToList(&extractor4);
|
---|
507 | // tlist4.AddToList(&fill4a);
|
---|
508 | tlist4.AddToList(&conv4);
|
---|
509 | tlist4.AddToList(&treat4);
|
---|
510 | tlist4.AddToList(&fill4b);
|
---|
511 |
|
---|
512 | if (!loop4.Eventloop(max4))
|
---|
513 | return 15;
|
---|
514 |
|
---|
515 | if (!loop4.GetDisplay())
|
---|
516 | return 16;
|
---|
517 |
|
---|
518 | // ===================================================================
|
---|
519 |
|
---|
520 | gLog << endl;
|
---|
521 | gLog.Separator("Extracting and calibration data");
|
---|
522 |
|
---|
523 | MTaskList tlist5;
|
---|
524 |
|
---|
525 | MParList plist5;
|
---|
526 | plist5.AddToList(&tlist5);
|
---|
527 | plist5.AddToList(&drscalib300);
|
---|
528 | plist5.AddToList(&badpixels);
|
---|
529 | plist5.AddToList(&timecam);
|
---|
530 |
|
---|
531 | MEvtLoop loop5("CalibratingData");
|
---|
532 | loop5.SetDisplay(d);
|
---|
533 | loop5.SetParList(&plist5);
|
---|
534 |
|
---|
535 | // ------------------ Setup the tasks ---------------
|
---|
536 |
|
---|
537 | MRawFitsRead read5a;
|
---|
538 | MReadMarsFile read5b("Events");
|
---|
539 | read5a.LoadMap(pmap);
|
---|
540 | read5a.AddFile(datfile);
|
---|
541 | read5b.DisableAutoScheme();
|
---|
542 | read5b.AddFile(datfile);
|
---|
543 |
|
---|
544 | MRead &read5 = datfile.EndsWith(".root") ? static_cast<MRead&>(read5b) : static_cast<MRead&>(read5a);
|
---|
545 |
|
---|
546 | MFillH fill5a(&hrate);
|
---|
547 |
|
---|
548 | MGeomApply apply5;
|
---|
549 |
|
---|
550 | MDrsCalibApply drsapply5;
|
---|
551 |
|
---|
552 | MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0", "SelectDat");
|
---|
553 | MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal");
|
---|
554 | MFDataPhrase filterped("(MRawEvtHeader.GetTriggerID&0xff00)==0x400", "SelectPed");
|
---|
555 | MFDataPhrase filterncl("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectNonCal");
|
---|
556 |
|
---|
557 | //MContinue cont4("MRawEvtHeader.GetTriggerID!=4", "SelectData");
|
---|
558 |
|
---|
559 | // ---
|
---|
560 |
|
---|
561 | MExtractTimeAndChargeSpline extractor5dat;
|
---|
562 | extractor5dat.SetRange(first_slice, last_slice);
|
---|
563 | extractor5dat.SetRiseTimeHiGain(rise_time_dat);
|
---|
564 | extractor5dat.SetFallTimeHiGain(fall_time_dat);
|
---|
565 | extractor5dat.SetHeightTm(heighttm);
|
---|
566 | extractor5dat.SetChargeType(type);
|
---|
567 | extractor5dat.SetSaturationLimit(600000);
|
---|
568 | extractor5dat.SetNoiseCalculation(kFALSE);
|
---|
569 |
|
---|
570 | MExtractTimeAndChargeSpline extractor5cal;
|
---|
571 | extractor5cal.SetRange(first_slice, last_slice);
|
---|
572 | extractor5cal.SetRiseTimeHiGain(rise_time_cal);
|
---|
573 | extractor5cal.SetFallTimeHiGain(fall_time_cal);
|
---|
574 | extractor5cal.SetHeightTm(heighttm);
|
---|
575 | extractor5cal.SetChargeType(type);
|
---|
576 | extractor5cal.SetSaturationLimit(600000);
|
---|
577 | extractor5cal.SetNoiseCalculation(kFALSE);
|
---|
578 |
|
---|
579 | MExtractTimeAndChargeSpline extractor5tm("ExtractTM");
|
---|
580 | extractor5tm.SetRange(last_slice, 294);
|
---|
581 | extractor5tm.SetRiseTimeHiGain(1);
|
---|
582 | extractor5tm.SetFallTimeHiGain(1);
|
---|
583 | extractor5tm.SetHeightTm(0.5);
|
---|
584 | extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel);
|
---|
585 | extractor5tm.SetSaturationLimit(600000);
|
---|
586 | extractor5tm.SetNoiseCalculation(kFALSE);
|
---|
587 | extractor5tm.SetNameSignalCam("TimeMarkerAmplitude");
|
---|
588 | extractor5tm.SetNameTimeCam("TimeMarkerTime");
|
---|
589 |
|
---|
590 | extractor5dat.SetFilter(&filterncl);
|
---|
591 | extractor5cal.SetFilter(&filtercal);
|
---|
592 | //extractor4tm.SetFilter(&filtercal);
|
---|
593 |
|
---|
594 | // ---
|
---|
595 | /*
|
---|
596 | MCalibrateData conv5;
|
---|
597 | conv5.SetCalibrationMode(MCalibrateData::kNone);
|
---|
598 | conv5.SetPedestalFlag(MCalibrateData::kNo);
|
---|
599 | conv5.SetCalibConvMinLimit(0);
|
---|
600 | conv5.SetCalibConvMaxLimit(10000);
|
---|
601 | conv5.SetScaleFactor(scale);
|
---|
602 | */
|
---|
603 | MCalibrateFact conv5;
|
---|
604 | conv5.SetScale(scale);
|
---|
605 | conv5.SetCalibConst(calib);
|
---|
606 |
|
---|
607 | MCalibrateDrsTimes calctm5;
|
---|
608 | calctm5.SetNameUncalibrated("UncalibratedTimes");
|
---|
609 |
|
---|
610 | MCalibrateDrsTimes calctm5tm("CalibrateTimeMarker");
|
---|
611 | calctm5tm.SetNameArrivalTime("TimeMarkerTime");
|
---|
612 | calctm5tm.SetNameUncalibrated("UncalTimeMarker");
|
---|
613 | calctm5tm.SetNameCalibrated("TimeMarker");
|
---|
614 | calctm5tm.SetTimeMarker();
|
---|
615 | //calctm4tm.SetFilter(&filtercal);
|
---|
616 |
|
---|
617 | MBadPixelsTreat treat5;
|
---|
618 | treat5.SetProcessPedestalRun(kFALSE);
|
---|
619 | treat5.SetProcessPedestalEvt(kFALSE);
|
---|
620 |
|
---|
621 | MHCamEvent evt5b(0, "ExtSig", "Extracted signal;;S [mV·sl]");
|
---|
622 | MHCamEvent evt5c(0, "CalSig", "Calibrated and interpolated signal;;S [~phe]");
|
---|
623 | MHCamEvent evt5d(4, "ExtSigTm", "Extracted time;;T [sl]");
|
---|
624 | MHCamEvent evt5e(6, "CalSigTm", "Calibrated and interpolated time;;T [ns]");
|
---|
625 |
|
---|
626 | MFillH fill5b(&evt5b, "MExtractedSignalCam", "FillExtSig");
|
---|
627 | MFillH fill5c(&evt5c, "MSignalCam", "FillCalSig");
|
---|
628 | MFillH fill5d(&evt5d, "MArrivalTimeCam", "FillExtTm");
|
---|
629 | MFillH fill5e(&evt5e, "MSignalCam", "FillCalTm");
|
---|
630 |
|
---|
631 | fill5c.SetDrawOption("gaus");
|
---|
632 | fill5d.SetDrawOption("gaus");
|
---|
633 | fill5e.SetDrawOption("gaus");
|
---|
634 |
|
---|
635 | /*
|
---|
636 | fill4b.SetFilter(&filterdat);
|
---|
637 | fill4c.SetFilter(&filterdat);
|
---|
638 | fill4d.SetFilter(&filterdat);
|
---|
639 | fill4e.SetFilter(&filterdat);
|
---|
640 | */
|
---|
641 |
|
---|
642 | //MFSoftwareTrigger swtrig;
|
---|
643 | //MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger");
|
---|
644 | //contsw.SetInverted();
|
---|
645 |
|
---|
646 | // The second rule is for the case reading raw-files!
|
---|
647 |
|
---|
648 | MWriteRootFile write5(outfile, "RECREATE", "Calibrated Data", 2);
|
---|
649 | write5.AddContainer("MRawRunHeader", "RunHeaders");
|
---|
650 | write5.AddContainer("MGeomCam", "RunHeaders");
|
---|
651 | write5.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
|
---|
652 | write5.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
|
---|
653 | write5.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
|
---|
654 |
|
---|
655 | // Common events
|
---|
656 | write5.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
|
---|
657 | write5.AddContainer("MMcEvt", "Events", kFALSE);
|
---|
658 | write5.AddContainer("IncidentAngle", "Events", kFALSE);
|
---|
659 | write5.AddContainer("MPointingPos", "Events", kFALSE);
|
---|
660 | write5.AddContainer("MSignalCam", "Events");
|
---|
661 | write5.AddContainer("MTime", "Events", kFALSE);
|
---|
662 | write5.AddContainer("MRawEvtHeader", "Events");
|
---|
663 | //write.AddContainer("MTriggerPattern", "Events");
|
---|
664 |
|
---|
665 | // ------------------ Setup histograms and fill tasks ----------------
|
---|
666 |
|
---|
667 | MContinue test;
|
---|
668 | test.SetFilter(&filterncl);
|
---|
669 |
|
---|
670 | MTaskList tlist5tm;
|
---|
671 | tlist5tm.AddToList(&extractor5tm);
|
---|
672 | tlist5tm.AddToList(&calctm5tm);
|
---|
673 | tlist5tm.SetFilter(&filtercal);
|
---|
674 |
|
---|
675 | MTaskList tlist5dat;
|
---|
676 | tlist5dat.AddToList(&fill5b);
|
---|
677 | tlist5dat.AddToList(&fill5c);
|
---|
678 | tlist5dat.AddToList(&fill5d);
|
---|
679 | tlist5dat.AddToList(&fill5e);
|
---|
680 | tlist5dat.SetFilter(&filterdat);
|
---|
681 |
|
---|
682 | tlist5.AddToList(&read5);
|
---|
683 | tlist5.AddToList(&apply5);
|
---|
684 | tlist5.AddToList(&drsapply5);
|
---|
685 | tlist5.AddToList(&filterncl);
|
---|
686 | //tlist5.AddToList(&test);
|
---|
687 | tlist5.AddToList(&filterdat);
|
---|
688 | tlist5.AddToList(&filtercal);
|
---|
689 | tlist5.AddToList(&filterped);
|
---|
690 | tlist5.AddToList(&fill5a);
|
---|
691 | tlist5.AddToList(&extractor5dat);
|
---|
692 | tlist5.AddToList(&extractor5cal);
|
---|
693 | tlist5.AddToList(&calctm5);
|
---|
694 | tlist5.AddToList(&tlist5tm);
|
---|
695 | tlist5.AddToList(&conv5);
|
---|
696 | tlist5.AddToList(&treat5);
|
---|
697 | tlist5.AddToList(&tlist5dat);
|
---|
698 | tlist5.AddToList(&write5);
|
---|
699 |
|
---|
700 | if (!loop5.Eventloop(max4))
|
---|
701 | return 18;
|
---|
702 |
|
---|
703 | if (!loop5.GetDisplay())
|
---|
704 | return 19;
|
---|
705 |
|
---|
706 | d->SetTitle(displaytitle, kFALSE);
|
---|
707 | d->SaveAs(displayfile);
|
---|
708 |
|
---|
709 | return 0;
|
---|
710 | }
|
---|