source: branches/removing_cpp11_features/fact/analysis/ganymed.C@ 18086

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