source: trunk/Mars/fact/analysis/ganymed.C@ 20112

Last change on this file since 20112 was 19160, checked in by Daniela Dorner, 6 years ago
new disp parametrization for analysis + cut that had not been committed yet
File size: 16.9 KB
Line 
1void SetupWriter(MWriteRootFile *write, const char *name) const
2{
3 if (!write)
4 return;
5
6 write->SetName(name);
7 write->AddContainer("MHillas", "Events");
8 write->AddContainer("MHillasSrc", "Events");
9 write->AddContainer("MHillasExt", "Events");
10 write->AddContainer("MPointingPos", "Events");
11 write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
12 write->AddContainer("MImagePar", "Events", kFALSE);
13 write->AddContainer("MNewImagePar", "Events", kFALSE);
14 write->AddContainer("MNewImagePar2", "Events", kFALSE);
15 write->AddContainer("Hadronness", "Events", kFALSE);
16 write->AddContainer("MSrcPosCam", "Events", kFALSE);
17 write->AddContainer("MSrcPosAnti", "Events", kFALSE);
18 write->AddContainer("ThetaSquared", "Events", kFALSE);
19 write->AddContainer("OpticalAxis", "Events", kFALSE);
20 write->AddContainer("Disp", "Events", kFALSE);
21 write->AddContainer("Ghostbuster", "Events", kFALSE);
22 write->AddContainer("MEnergyEst", "Events", kFALSE);
23 write->AddContainer("MTime", "Events", kFALSE);
24 write->AddContainer("MMcEvt", "Events", kFALSE);
25 write->AddContainer("DataType", "Events");
26 write->AddContainer("FileId", "Events");
27 write->AddContainer("EvtNumber", "Events");
28 // write->AddContainer("MMuonSearchPar", "Events", kFALSE);
29 // write->AddContainer("MMuonCalibPar", "Events", kFALSE);
30}
31
32int process(MDirIter &files, const char *outfile, Float_t ra, Float_t dec)
33{
34 /*
35 MDirIter files;
36 if (files.ReadList(dataset)<=0)
37 return 2;
38 */
39 files.Sort();
40
41 // =============================================================
42
43 // Crab (05h34'32" 22d00'52")
44 MPointingPos source("MSourcePos");
45 //source.SetSkyPosition(5.5755, 22.0144); // ra[h], dec[deg] // Crab
46 //source.SetSkyPosition(12.356089, 30.176897); // ra[h], dec[deg]//1218
47 //source.SetSkyPosition(16.897867,39.760201);// mrk501
48 //source.SetSkyPosition(11.0742, 38.2089); // ra[h], dec[deg]//mrk421
49 source.SetSkyPosition(ra, dec); // ra[h], dec[deg]//mrk421
50
51 UInt_t fNumOffSourcePos = 5;
52
53 // =============================================================
54
55 // FIXME: If fPathIn read cuts and energy estimator from file!
56 MContinue contq("MNewImagePar.fNumUsedPixels<5.5"
57 // "|| MImagePar.fNumIslands>3.5"
58 "|| MNewImagePar.fLeakage1>0.1",
59 "CutQ");
60
61 MContinue cont0("MHillas.fSize<0"
62 "|| MHillas.GetArea > log10(MHillas.fSize)*898-1535 "
63 , "Cut0");
64
65 // Cuts before plots
66 MContinue cont2("0"
67 "|| MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg<-0.07"
68 "|| log10(MNewImagePar.fConc1)<-log10(MHillas.fSize)*0.57+0.13"
69 //"|| log10(MNewImagePar.fConc1)>-log10(MHillas.fSize)*0.57+0.5"
70 //"|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>log10(MHillas.fSize)*0.42+0.29"
71 "|| MHillas.GetArea*MGeomCam.fConvMm2Deg^2<0.01"
72 "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>log10(MHillas.fSize)*0.5+0.095"
73 "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg<0.4"
74 "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>1.2",
75 "Cut2");
76
77 TArrayD param(11);
78 // Parametrization of Disp
79 // for new cleaning
80 //param[0] = 1.39252; // constant
81 param[0] = 1.34; // constant
82 //param[8] = 0.154247; // slope
83 param[8] = 0.0755; // slope
84 param[9] = 1.67972; // leak
85 param[10] = 4.86232; // leak
86 // for old cleaning
87 //param[0] = 1.34723; // constant
88 //param[8] = 0.15214; // slope
89 //param[9] = 0.970704; // leak
90 //param[10] = 8.89826; // leak
91 // Parametrization for sign of disp (m3long, slope)
92 param[5] = -0.07;
93 param[6] = 7.2;
94 param[7] = 0.5;
95 // ThetaSq-Cut
96 param[1] = 0.156; // 0.215
97 // Area-Cut
98 param[2] = 0.215468;
99 param[3] = 5.63973;
100 param[4] = 0.0836169;
101
102 // --------------------------------------------------------------------------------
103
104 // Setup fitter and histograms
105 MAlphaFitter fit;
106 fit.SetBackgroundFitMin(0.09);
107 fit.SetBackgroundFitMax(0.25);
108 fit.SetPolynomOrder(1);
109 fit.SetSignalFunction(MAlphaFitter::kThetaSq);
110 fit.SetScaleUser(1./fNumOffSourcePos); // includes fit.SetScaleMode(MAlphaFitter::kUserScale);
111
112 // --------------------------------------------------------------------------------
113
114 //MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
115 MBinning bins3(100, -0.005, 0.995, "BinningTheta", "asin");
116 MBinning binsD( 25, 0, 2.5, "BinningDist");
117 MBinning binsC( 25, 0, 2.5, "BinningDistC");
118 MBinning binsR( 60, 0, 300, "BinningRate");
119
120 // --------------------------------------------------------------------------------
121
122 MParList plist;
123 plist.AddToList(&source);
124 plist.AddToList(&fit);
125
126 // Initialize default binnings
127 plist.AddToList(&bins3);
128 plist.AddToList(&binsD);
129 plist.AddToList(&binsC);
130 plist.AddToList(&binsR);
131
132 MParameterI par("DataType");
133 plist.AddToList(&par);
134
135 // Setup Tasklist
136 MTaskList tlist;
137 plist.AddToList(&tlist);
138
139 // La Palma Magic1
140 MObservatory obs;
141 plist.AddToList(&obs);
142
143 // --------------------------------------------------------------------------------
144
145 MParameterD scale;
146 scale.SetVal(1./fNumOffSourcePos);
147
148 MHThetaSq *halphaoff = plist.FindCreateObj("MHThetaSq", "HistOff");
149 halphaoff->ForceUsingSize();
150
151 MFillH falpha(halphaoff, "", "FillHist");
152
153 MFMagicCuts cuts;
154 cuts.SetHadronnessCut(MFMagicCuts::kNoCut);
155 cuts.SetVariables(param);
156
157 // Cuts before thetasq
158 MContinue cont1(&cuts, "Cut1");
159
160 // Cuts before writing
161 MContinue cont3("", "Cut3");
162 contq.SetAllowEmpty();
163 cont0.SetAllowEmpty();
164 cont1.SetAllowEmpty();
165 cont2.SetAllowEmpty();
166 cont3.SetAllowEmpty();
167
168 cont1.SetInverted();
169
170 // Filter for VsSize
171 MFDataPhrase ftheta("ThetaSquared.fVal < ThetaSquaredCut.fVal", "CutT");
172
173 // ------------- Loop Off Data --------------------
174 MReadReports readoffdata;
175 // readoffdata.DisableAutoScheme();
176 readoffdata.AddTree("Events", "MTime.", MReadReports::kMaster);
177 readoffdata.AddTree("Drive", MReadReports::kRequired);
178
179 readoffdata.AddFiles(files);
180
181 TString fname0 = Form("%s-summary.root", outfile);
182 TString fname1 = Form("%s-analysis.root", outfile);
183
184 MWriteRootFile dummy0(fname0.Data(), "RECREATE");
185 MWriteRootFile dummy1(fname1.Data(), "RECREATE");
186
187 MWriteRootFile *write0 = &dummy0;
188 MWriteRootFile *write1 = &dummy1;
189 SetupWriter(write0, "WriteAfterCut0");
190 SetupWriter(write1, "WriteAfterCut3");
191
192 MParameterCalc setevtnum("MRawEvtHeader.fDAQEvtNumber", "SetEvtNumber");
193 setevtnum.SetNameParameter("EvtNumber");
194
195 MParameterCalc setrunnum("MRawRunHeader.GetFileID", "SetFileId");
196 setrunnum.SetNameParameter("FileId");
197
198 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
199 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
200 MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
201 MFillH fill3c("MHVsSizeOffTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta");
202 MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
203 MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
204 MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
205 MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
206 //MFillH fill9a("MHEffOffTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
207 fill1a.SetNameTab("PreCut");
208 fill2a.SetNameTab("PostCut");
209 fill3a.SetNameTab("VsSize");
210 fill3c.SetNameTab("CutT");
211 fill4a.SetNameTab("HilExt");
212 fill5a.SetNameTab("HilSrc");
213 fill6a.SetNameTab("ImgPar");
214 fill7a.SetNameTab("NewPar");
215 //fill9a.SetNameTab("EffOffT");
216
217 fill3c.SetFilter(&ftheta);
218
219 //MFDataMember fbin("Bin.fVal", '>', 0);
220 //fill9a.SetFilter(&fbin);
221
222 MTaskList tlist2;
223 tlist2.SetNumPasses(fNumOffSourcePos);
224 fill2a.SetWeight(&scale);
225 fill3a.SetWeight(&scale);
226 fill3c.SetWeight(&scale);
227 fill4a.SetWeight(&scale);
228 fill5a.SetWeight(&scale);
229 fill6a.SetWeight(&scale);
230 fill7a.SetWeight(&scale);
231
232 // How to get source position from off- and on-data?
233 MSrcPosCalc scalc;
234 scalc.SetMode(MSrcPosCalc::kWobble);
235 scalc.SetCallback(&tlist2);
236
237 MSrcPosCorrect scor;
238
239 MHillasCalc hcalc;
240 MHillasCalc hcalc2("MHillasCalcAnti");
241 hcalc.SetAbberation(1.02);
242 hcalc2.SetAbberation(1.02);
243 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
244 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
245 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
246 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
247
248 MSrcPosRndm srcrndm;
249
250 tlist2.AddToList(&scalc);
251 //tlist2.AddToList(&scor);
252 tlist2.AddToList(&srcrndm);
253 tlist2.AddToList(&hcalc);
254 tlist2.AddToList(&hcalc2);
255 tlist2.AddToList(&setrunnum);
256 tlist2.AddToList(&setevtnum);
257 tlist2.AddToList(&cuts);
258 tlist2.AddToList(write0);
259 tlist2.AddToList(&cont0); // Post/Pre shown cuts
260 tlist2.AddToList(&cont1); // Calculate ThetaSq
261 tlist2.AddToList(&fill2a);
262 tlist2.AddToList(&ftheta); // Calculate theta Cut
263 tlist2.AddToList(&falpha);
264 tlist2.AddToList(&fill3c); // With theta cut
265 tlist2.AddToList(&cont2);
266 tlist2.AddToList(&fill3a);
267 tlist2.AddToList(&fill4a);
268 tlist2.AddToList(&fill5a);
269 tlist2.AddToList(&fill6a);
270 tlist2.AddToList(&fill7a);
271 tlist2.AddToList(&cont3);
272 tlist2.AddToList(write1);
273
274 MPointingPosCalc pcalc;
275 MPointingDevCalc devcalc;
276
277 tlist.AddToList(&readoffdata);
278 tlist.AddToList(&pcalc, "Drive");
279 tlist.AddToList(&contq, "Events");
280 tlist.AddToList(&fill1a, "Events");
281 tlist.AddToList(&tlist2, "Events");
282
283 // by setting it here it is distributed to all consecutive tasks
284 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
285
286 par.SetVal(0);
287
288 MStatusDisplay *d = new MStatusDisplay;
289
290 // Create and setup the eventloop
291 MEvtLoop evtloop("Off");
292 evtloop.SetParList(&plist);
293 evtloop.SetDisplay(d);
294
295 // Execute first analysis
296 if (!evtloop.Eventloop())
297 return 2;
298
299 if (!evtloop.GetDisplay())
300 return 2;
301
302 // ------------- Loop On Data --------------------
303 MReadReports readondata;
304 // readoffdata.DisableAutoScheme();
305 readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
306 readondata.AddTree("Drive", MReadReports::kRequired);
307 readondata.AddTree("Rates", MReadReports::kRequired);
308 readondata.AddTree("Humidity", MReadReports::kRequired);
309 readondata.AddTree("Temperatures", MReadReports::kRequired);
310
311 readondata.AddFiles(files);
312
313 scalc.SetMode(MSrcPosCalc::kDefault);
314 scalc.SetNumRandomOffPositions(fNumOffSourcePos);
315
316 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
317 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
318 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
319 MFillH fill3d("MHVsSizeOnTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta");
320 MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
321 MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
322 MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
323 MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
324 //MFillH fill9b("MHEffOnTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
325 fill1b.SetNameTab("PreCut");
326 fill2b.SetNameTab("PostCut");
327 fill3b.SetNameTab("VsSize");
328 fill3d.SetNameTab("CutT");
329 fill4b.SetNameTab("HilExt");
330 fill5b.SetNameTab("HilSrc");
331 fill6b.SetNameTab("ImgPar");
332 fill7b.SetNameTab("NewPar");
333 //fill9b.SetNameTab("EffOnT");
334 fill1b.SetDrawOption("same");
335 fill2b.SetDrawOption("same");
336 fill3b.SetDrawOption("same");
337 fill3d.SetDrawOption("same");
338 fill4b.SetDrawOption("same");
339 fill5b.SetDrawOption("same");
340 fill6b.SetDrawOption("same");
341 fill7b.SetDrawOption("same");
342 //fill9b.SetFilter(&fbin);
343
344 fill3d.SetFilter(&ftheta);
345
346 MHThetaSq *halphaon = plist.FindCreateObj("MHThetaSq", "Hist");
347 halphaon->ForceUsingSize();
348
349 MFillH falpha2(halphaon, "", "FillHist");
350
351 tlist2.SetNumPasses();
352
353 tlist.Replace(&readondata);
354
355 MHSrcPosCam hsrcpos(kTRUE);
356
357 MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
358 fillsrc.SetNameTab("SrcPos");
359
360 tlist2.AddToListBefore(&fillsrc, &hcalc);
361
362 MH3 hvs1("MPointingPos.fZd");
363 MH3 hvs2("MPointingPos.fZd");
364 MHVsTime hvstM("MTimeDrive.GetMjd");
365 MHVsTime hvst0("MPointingPos.fZd");
366 MHVsTime hvst1("MReportRates.fTriggerRate");
367 MH3 hvst2("MPointingPos.fZd", "MReportRates.fTriggerRate");
368 MHVsTime hvst("MReportTemperatures.GetTempMean");
369 MHVsTime hvsh("MReportHumidity.GetMean");
370 hvs1.SetName("ObsTime;Theta");
371 hvs2.SetName("OnTime;Theta");
372 hvstM.SetName("MJD");
373 hvst0.SetName("ZdTime");
374 hvst1.SetName("Rate");
375 hvst2.SetName("RateZd;Theta;Rate");
376 hvst.SetName("CamTemp");
377 hvsh.SetName("CamHum");
378 hvs1.SetTitle("Observation time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
379 hvs2.SetTitle("On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
380 hvstM.SetTitle("MJD vs time;Time;MJD");
381 hvst0.SetTitle("Zenith distance vs time;Time;Zd [\\circ]");
382 hvst1.SetTitle("Camera trigger rate vs time;Time;Rate [Hz]");
383 hvst2.SetTitle("Camera trigger rate vs Zd; Zd [\\circ];Rate [Hz]");
384 hvst.SetTitle("Average sensor compartment temperature;Time;T [\\circ C]");
385 hvsh.SetTitle("Average camera humidity;Time;Hum [%]");
386
387 MFillH fillvs1(&hvs1, "", "FillObsTime");
388 MFillH fillvs2(&hvs2, "", "FillOnTime");
389 MFillH fillvstM(&hvstM, "MTime", "FillMjdVsTime");
390 MFillH fillvst0(&hvst0, "MTimeDrive", "FillZdVsTime");
391 MFillH fillvst1(&hvst1, "MTimeRates", "FillRate");
392 MFillH fillvst2(&hvst2, "", "FillRateVsZd");
393 MFillH fillvst(&hvst, "MTimeTemperatures", "FillCamTemp");
394 MFillH fillvsh(&hvsh, "MTimeHumidity", "FillCamHum");
395 fillvs1.SetNameTab("ObsTime");
396 fillvs2.SetNameTab("OnTime");
397 fillvst1.SetNameTab("Rate");
398 fillvst2.SetNameTab("RateZd");
399 fillvst.SetNameTab("CamTemp");
400 fillvsh.SetNameTab("CamHum");
401
402 hvs1.SetWeight("MReportRates.fElapsedTime");
403 hvs2.SetWeight("MReportRates.fElapsedOnTime");
404
405 hvst1.SetMinimum(0);
406 hvst2.SetDrawOption("colz");
407
408 MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
409 hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
410 hrate.SetName("RateVsRun");
411 hrate.SetTitle("Event rate after quality cuts;File Id;Event Rate [Hz];");
412 hrate.InitLabels(MH3::kLabelsX);
413 hrate.GetHist().SetMinimum(0);
414
415 MFillH fillr(&hrate, "", "FillRateVsRun");
416 tlist.AddToListAfter(&fillr, &contq, "Events");
417
418 tlist.Replace(&fill1b);
419
420 tlist2.Replace(&fill2b);
421 tlist2.Replace(&fill3b);
422 tlist2.Replace(&fill3d);
423 tlist2.Replace(&fill4b);
424 tlist2.Replace(&fill5b);
425 tlist2.Replace(&fill6b);
426 tlist2.Replace(&fill7b);
427 tlist2.Replace(&falpha2);
428 tlist.AddToList(&fillvs1, "Rates");
429 tlist.AddToList(&fillvs2, "Rates");
430 tlist.AddToList(&fillvst1, "Rates");
431 tlist.AddToList(&fillvst2, "Rates");
432 tlist.AddToList(&fillvst, "Temperatures");
433 tlist.AddToList(&fillvsh, "Humidity");
434 tlist.AddToList(&fillvst0, "Drive");
435 tlist.AddToList(&fillvstM, "Drive");
436
437 // by setting it here it is distributed to all consecutive tasks
438 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
439
440 par.SetVal(1);
441
442 // Execute first analysis
443 if (!evtloop.Eventloop())
444 return 2;
445
446 if (!evtloop.GetDisplay())
447 return 2;
448
449 TObjArray cont;
450 //cont.Add(const_cast<TEnv*>(GetEnv()));
451 //cont.Add(const_cast<MSequence*>(&fSequence));
452
453 //TNamed cmdline("CommandLine", fCommandLine.Data());
454 //cont.Add(&cmdline);
455
456 if (d)
457 {
458 TString title = "-- Analysis #";
459 title += gSystem->BaseName(outfile);
460 title += " --";
461 d->SetTitle(title, kFALSE);
462
463 cont.Add(d);
464 TString path;
465 path += fname1;
466
467 //d->SaveAs(path);
468 d->Write();
469 }
470
471// if (!WriteContainer(cont, fname1, "RECREATE"))
472// return;
473
474 return 0;
475}
476
477int ganymed(const char *dataset, const char *outfile, Float_t ra, Float_t dec)
478{
479 MDirIter files;
480 if (files.ReadList(dataset)<=0)
481 return 2;
482 files.Sort();
483 return process(files, outfile, ra, dec);
484}
485
486
487int ganymed(Float_t ra, Float_t dec, const char *starfile, const char *outfile)
488{
489 MDirIter files;
490 files.AddFile(starfile);
491 return process(files, outfile, ra, dec);
492}
493
494
Note: See TracBrowser for help on using the repository browser.