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

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