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

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