source: branches/Mars_IncreaseNsb/fact/analysis/mc/callisto_data.C

Last change on this file was 18579, checked in by tbretz, 8 years ago
A copy of mc/callisto.C, removed the pedestal processing, adapted the mapping, removed the timecam from the tasklist to avoid processing the drs time calibration --> all this should go to callisto.C one day where MC files should be processed with the same code than data.
File size: 11.6 KB
Line 
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 "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
53using 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
77int 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·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}
Note: See TracBrowser for help on using the repository browser.