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

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