source: branches/Corsika7500Compatibility/fact/analysis/ganymed.C@ 18459

Last change on this file since 18459 was 18105, checked in by Daniela Dorner, 10 years ago
new coefficients for new FACT disp parametrization
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( 67, -0.005, 0.665, "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.SetFlags(MHillasCalc::kCalcHillasSrc);
243 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
244 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
245 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
246
247 MSrcPosRndm srcrndm;
248
249 tlist2.AddToList(&scalc);
250 //tlist2.AddToList(&scor);
251 tlist2.AddToList(&srcrndm);
252 tlist2.AddToList(&hcalc);
253 tlist2.AddToList(&hcalc2);
254 tlist2.AddToList(&setrunnum);
255 tlist2.AddToList(&setevtnum);
256 tlist2.AddToList(&cuts);
257 tlist2.AddToList(write0);
258 tlist2.AddToList(&cont0); // Post/Pre shown cuts
259 tlist2.AddToList(&cont1); // Calculate ThetaSq
260 tlist2.AddToList(&fill2a);
261 tlist2.AddToList(&ftheta); // Calculate theta Cut
262 tlist2.AddToList(&falpha);
263 tlist2.AddToList(&fill3c); // With theta cut
264 tlist2.AddToList(&cont2);
265 tlist2.AddToList(&fill3a);
266 tlist2.AddToList(&fill4a);
267 tlist2.AddToList(&fill5a);
268 tlist2.AddToList(&fill6a);
269 tlist2.AddToList(&fill7a);
270 tlist2.AddToList(&cont3);
271 tlist2.AddToList(write1);
272
273 MPointingPosCalc pcalc;
274 MPointingDevCalc devcalc;
275
276 tlist.AddToList(&readoffdata);
277 tlist.AddToList(&pcalc, "Drive");
278 tlist.AddToList(&contq, "Events");
279 tlist.AddToList(&fill1a, "Events");
280 tlist.AddToList(&tlist2, "Events");
281
282 // by setting it here it is distributed to all consecutive tasks
283 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
284
285 par.SetVal(0);
286
287 MStatusDisplay *d = new MStatusDisplay;
288
289 // Create and setup the eventloop
290 MEvtLoop evtloop("Off");
291 evtloop.SetParList(&plist);
292 evtloop.SetDisplay(d);
293
294 // Execute first analysis
295 if (!evtloop.Eventloop())
296 return 2;
297
298 if (!evtloop.GetDisplay())
299 return 2;
300
301 // ------------- Loop On Data --------------------
302 MReadReports readondata;
303 // readoffdata.DisableAutoScheme();
304 readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
305 readondata.AddTree("Drive", MReadReports::kRequired);
306 readondata.AddTree("Rates", MReadReports::kRequired);
307 readondata.AddTree("Humidity", MReadReports::kRequired);
308 readondata.AddTree("Temperatures", MReadReports::kRequired);
309
310 readondata.AddFiles(files);
311
312 scalc.SetMode(MSrcPosCalc::kDefault);
313 scalc.SetNumRandomOffPositions(fNumOffSourcePos);
314
315 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
316 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
317 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
318 MFillH fill3d("MHVsSizeOnTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta");
319 MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
320 MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
321 MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
322 MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
323 //MFillH fill9b("MHEffOnTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
324 fill1b.SetNameTab("PreCut");
325 fill2b.SetNameTab("PostCut");
326 fill3b.SetNameTab("VsSize");
327 fill3d.SetNameTab("CutT");
328 fill4b.SetNameTab("HilExt");
329 fill5b.SetNameTab("HilSrc");
330 fill6b.SetNameTab("ImgPar");
331 fill7b.SetNameTab("NewPar");
332 //fill9b.SetNameTab("EffOnT");
333 fill1b.SetDrawOption("same");
334 fill2b.SetDrawOption("same");
335 fill3b.SetDrawOption("same");
336 fill3d.SetDrawOption("same");
337 fill4b.SetDrawOption("same");
338 fill5b.SetDrawOption("same");
339 fill6b.SetDrawOption("same");
340 fill7b.SetDrawOption("same");
341 //fill9b.SetFilter(&fbin);
342
343 fill3d.SetFilter(&ftheta);
344
345 MHThetaSq *halphaon = plist.FindCreateObj("MHThetaSq", "Hist");
346 halphaon->ForceUsingSize();
347
348 MFillH falpha2(halphaon, "", "FillHist");
349
350 tlist2.SetNumPasses();
351
352 tlist.Replace(&readondata);
353
354 MHSrcPosCam hsrcpos(kTRUE);
355
356 MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
357 fillsrc.SetNameTab("SrcPos");
358
359 tlist2.AddToListBefore(&fillsrc, &hcalc);
360
361 MH3 hvs1("MPointingPos.fZd");
362 MH3 hvs2("MPointingPos.fZd");
363 MHVsTime hvstM("MTimeDrive.GetMjd");
364 MHVsTime hvst0("MPointingPos.fZd");
365 MHVsTime hvst1("MReportRates.fTriggerRate");
366 MH3 hvst2("MPointingPos.fZd", "MReportRates.fTriggerRate");
367 MHVsTime hvst("MReportTemperatures.GetTempMean");
368 MHVsTime hvsh("MReportHumidity.GetMean");
369 hvs1.SetName("ObsTime;Theta");
370 hvs2.SetName("OnTime;Theta");
371 hvstM.SetName("MJD");
372 hvst0.SetName("ZdTime");
373 hvst1.SetName("Rate");
374 hvst2.SetName("RateZd;Theta;Rate");
375 hvst.SetName("CamTemp");
376 hvsh.SetName("CamHum");
377 hvs1.SetTitle("Observation time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
378 hvs2.SetTitle("On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
379 hvstM.SetTitle("MJD vs time;Time;MJD");
380 hvst0.SetTitle("Zenith distance vs time;Time;Zd [\\circ]");
381 hvst1.SetTitle("Camera trigger rate vs time;Time;Rate [Hz]");
382 hvst2.SetTitle("Camera trigger rate vs Zd; Zd [\\circ];Rate [Hz]");
383 hvst.SetTitle("Average sensor compartment temperature;Time;T [\\circ C]");
384 hvsh.SetTitle("Average camera humidity;Time;Hum [%]");
385
386 MFillH fillvs1(&hvs1, "", "FillObsTime");
387 MFillH fillvs2(&hvs2, "", "FillOnTime");
388 MFillH fillvstM(&hvstM, "MTime", "FillMjdVsTime");
389 MFillH fillvst0(&hvst0, "MTimeDrive", "FillZdVsTime");
390 MFillH fillvst1(&hvst1, "MTimeRates", "FillRate");
391 MFillH fillvst2(&hvst2, "", "FillRateVsZd");
392 MFillH fillvst(&hvst, "MTimeTemperatures", "FillCamTemp");
393 MFillH fillvsh(&hvsh, "MTimeHumidity", "FillCamHum");
394 fillvs1.SetNameTab("ObsTime");
395 fillvs2.SetNameTab("OnTime");
396 fillvst1.SetNameTab("Rate");
397 fillvst2.SetNameTab("RateZd");
398 fillvst.SetNameTab("CamTemp");
399 fillvsh.SetNameTab("CamHum");
400
401 hvs1.SetWeight("MReportRates.fElapsedTime");
402 hvs2.SetWeight("MReportRates.fElapsedOnTime");
403
404 hvst1.SetMinimum(0);
405 hvst2.SetDrawOption("colz");
406
407 MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
408 hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
409 hrate.SetName("RateVsRun");
410 hrate.SetTitle("Event rate after quality cuts;File Id;Event Rate [Hz];");
411 hrate.InitLabels(MH3::kLabelsX);
412 hrate.GetHist().SetMinimum(0);
413
414 MFillH fillr(&hrate, "", "FillRateVsRun");
415 tlist.AddToListAfter(&fillr, &contq, "Events");
416
417 tlist.Replace(&fill1b);
418
419 tlist2.Replace(&fill2b);
420 tlist2.Replace(&fill3b);
421 tlist2.Replace(&fill3d);
422 tlist2.Replace(&fill4b);
423 tlist2.Replace(&fill5b);
424 tlist2.Replace(&fill6b);
425 tlist2.Replace(&fill7b);
426 tlist2.Replace(&falpha2);
427 tlist.AddToList(&fillvs1, "Rates");
428 tlist.AddToList(&fillvs2, "Rates");
429 tlist.AddToList(&fillvst1, "Rates");
430 tlist.AddToList(&fillvst2, "Rates");
431 tlist.AddToList(&fillvst, "Temperatures");
432 tlist.AddToList(&fillvsh, "Humidity");
433 tlist.AddToList(&fillvst0, "Drive");
434 tlist.AddToList(&fillvstM, "Drive");
435
436 // by setting it here it is distributed to all consecutive tasks
437 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
438
439 par.SetVal(1);
440
441 // Execute first analysis
442 if (!evtloop.Eventloop())
443 return 2;
444
445 if (!evtloop.GetDisplay())
446 return 2;
447
448 TObjArray cont;
449 //cont.Add(const_cast<TEnv*>(GetEnv()));
450 //cont.Add(const_cast<MSequence*>(&fSequence));
451
452 //TNamed cmdline("CommandLine", fCommandLine.Data());
453 //cont.Add(&cmdline);
454
455 if (d)
456 {
457 TString title = "-- Analysis #";
458 title += gSystem->BaseName(outfile);
459 title += " --";
460 d->SetTitle(title, kFALSE);
461
462 cont.Add(d);
463 TString path;
464 path += fname1;
465
466 d->SaveAs(path);
467 }
468
469// if (!WriteContainer(cont, fname1, "RECREATE"))
470// return;
471
472 return 0;
473}
474
475int ganymed(const char *dataset, const char *outfile, Float_t ra, Float_t dec)
476{
477 MDirIter files;
478 if (files.ReadList(dataset)<=0)
479 return 2;
480 files.Sort();
481 return process(files, outfile, ra, dec);
482}
483
484
485int ganymed(Float_t ra, Float_t dec, const char *starfile, const char *outfile)
486{
487 MDirIter files;
488 files.AddFile(starfile);
489 return process(files, outfile, ra, dec);
490}
491
492
Note: See TracBrowser for help on using the repository browser.