source: branches/Corsika7500Compatibility/fact/analysis/mc/ganymed.C@ 19366

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