source: branches/MarsMoreSimulationTruth/fact/analysis/ganymed.C@ 18764

Last change on this file since 18764 was 18468, checked in by tbretz, 10 years ago
The theta binning was too short... we have data until almost 90deg
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.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.