source: releases/Mars.2014.05.26/fact/analysis/star_file.C

Last change on this file was 18028, checked in by dneise, 10 years ago
Dedicated analysis macros for QLA. They have been written by Thomas Bretz for the QLA.
File size: 9.6 KB
Line 
1#include "MLogManip.h"
2
3int star_file(const char *datafile, Double_t lvl1=7.8, Double_t lvl2=3.9, const char *outpath = ".")
4{
5 double deltat = 17.5; // ns/deg
6 // lvl1 = 2.5;
7 // lvl2 = 0.5;
8
9 // ------------------------------------------------------
10
11 gLog.Separator("Star");
12 gLog << all << "Calculate image parameters of sequence ";
13 gLog << datafile << endl;
14 gLog << endl;
15
16 // ------------------------------------------------------
17
18 gLog << "Outpath: " << outpath << endl;
19
20 TString fname = gSystem->ConcatFileName(outpath, gSystem->BaseName(datafile));
21 fname.ReplaceAll("_C.root", "_I.root");
22
23 gLog.Separator();
24
25 // --------------------------------------------------------
26
27 MBinning bins1( 36, -90, 90, "BinningAlpha");
28 MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
29 MBinning bins4( 25, 0, 2.5, "BinningDist");
30
31 MBinning binsM1(100, 0, 5, "BinningMuonRadius");
32 MBinning binsM2( 60, 0, 0.3, "BinningMuonDeviation");
33 MBinning binsM3(150, 0, 60, "BinningMuonTime");
34 MBinning binsM4(300, 0, 30, "BinningMuonTimeRms");
35 MBinning binsM5( 20, 0, 360, "BinningMuonArcPhi");
36 MBinning binsM6( 50, 0, 0.5, "BinningMuonArcWidth");
37 MBinning binsM7( 30, 5, 5000, "BinningMuonSize", "log");
38 MBinning binsM8(100, 0, 5, "BinningMuonRelTimeSigma");
39
40 MBinning binsM9(100, 0, 5, "BinningRadius");
41
42 // ---------------------------------------------------------
43
44 // Filter to start muon analysis
45 //MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut");
46 // Filter to calculate further muon parameters
47// MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&"
48// "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut");
49 MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&"
50 "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4",
51 "MuonSearchCut");
52
53 // Filter to fill the MHMuonPar
54 MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2",
55 "MuonFinalCut");
56 // Filter to write Muons to Muon tree
57 //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
58
59 // ---------------------------------------------------------
60
61 MStatusDisplay *d = new MStatusDisplay;
62
63 MTaskList tlist;
64
65 MParList plist2;
66 plist2.AddToList(&tlist);
67 plist2.AddToList(&bins1);
68 plist2.AddToList(&bins3);
69 plist2.AddToList(&bins4);
70 plist2.AddToList(&binsM1);
71 plist2.AddToList(&binsM2);
72 plist2.AddToList(&binsM3);
73 plist2.AddToList(&binsM4);
74 plist2.AddToList(&binsM5);
75 plist2.AddToList(&binsM6);
76 plist2.AddToList(&binsM7);
77 plist2.AddToList(&binsM8);
78 plist2.AddToList(&binsM9);
79
80 MMuonSetup muonsetup;
81 plist2.AddToList(&muonsetup);
82
83 MEvtLoop loop;
84 loop.SetDisplay(d);
85 loop.SetParList(&plist2);
86
87 MReadMarsFile read("Events");
88 read.DisableAutoScheme();
89 read.AddFile(datafile);
90
91 MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData");
92
93 MGeomApply apply;
94
95 MImgCleanStd clean(lvl1, lvl2);
96 clean.SetMethod(MImgCleanStd::kAbsolute);
97 clean.SetTimeLvl1(deltat);//17.5/*1.3*/);
98 clean.SetTimeLvl2(deltat);//17.5/*1.3*/);
99 clean.SetPostCleanType(3);
100
101 MHillasCalc hcalc;
102 hcalc.Disable(MHillasCalc::kCalcConc);
103
104 MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
105 hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
106 hrate.SetName("Rate");
107 hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];");
108 hrate.InitLabels(MH3::kLabelsX);
109 hrate.GetHist().SetMinimum(0);
110
111 MFillH fillR(&hrate, "", "FillRate");
112
113 // ------------------ Setup histograms and fill tasks ----------------
114 MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]");
115 MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
116 MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]");
117 evtC1.SetErrorSpread(kFALSE);
118 evtC2.SetErrorSpread(kFALSE);
119 evtC2.SetThreshold(0);
120
121 MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam");
122 MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels");
123 MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos");
124
125 MFillH fillD1("MHHillas", "MHillas", "FillHillas");
126 MFillH fillD2("MHHillasExt", "", "FillHillasExt");
127 MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
128 MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar");
129 MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
130 MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize");
131
132 MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)");
133 halpha.SetName("AvsD;Dist;Alpha");
134
135 MFillH filla(&halpha);
136 filla.SetNameTab("AvsD");
137 filla.SetDrawOption("colz");
138
139 // ----------------------- Muon Analysis ----------------------
140 //writem.SetFilter(&fmuon4);
141
142 MHn mhn1;
143 mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg");
144 mhn1.InitName("MuonRadius");
145 mhn1.InitTitle("MuonRadius");
146
147 mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
148 mhn1.InitName("MuonDeviation");
149 mhn1.InitTitle("MuonDeviation");
150
151 mhn1.AddHist("MMuonSearchPar.fTime");
152 mhn1.InitName("MuonTime");
153 mhn1.InitTitle("MuonTime");
154
155 mhn1.AddHist("MMuonSearchPar.fTimeRms");
156 mhn1.InitName("MuonTimeRms");
157 mhn1.InitTitle("MuonTimeRms");
158
159 MHn mhn3;
160 mhn3.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
161 mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation");
162 mhn3.InitTitle("MuonRadius");
163 mhn3.SetDrawOption("colz");
164
165 MHn mhn2;
166 mhn2.AddHist("MMuonCalibPar.fArcPhi");
167 mhn2.InitName("MuonArcPhi");
168 mhn2.InitTitle("MuonArcPhi");
169
170 mhn2.AddHist("MMuonCalibPar.fArcWidth");
171 mhn2.InitName("MuonArcWidth");
172 mhn2.InitTitle("MuonArcWidth");
173
174 mhn2.AddHist("MMuonCalibPar.fMuonSize");
175 mhn2.InitName("MuonSize");
176 mhn2.InitTitle("MuonSize");
177
178 mhn2.AddHist("MMuonCalibPar.fRelTimeSigma");
179 mhn2.InitName("MuonRelTimeSigma");
180 mhn2.InitTitle("MuonRelTimeSigma");
181
182 MHn mhn4;
183 mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth");
184 mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth");
185 mhn4.InitTitle("MuonArcPhi");
186 mhn4.SetDrawOption("colz");
187
188 MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar");
189 MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar");
190 MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius");
191 MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi");
192 fillmhn1.SetNameTab("MuonS");
193 fillmhn3.SetNameTab("DevVsRad");
194 fillmhn2.SetNameTab("MuonC");
195 fillmhn4.SetNameTab("WidthVsPhi");
196 fillmhn2.SetFilter(&fmuon2);
197 fillmhn4.SetFilter(&fmuon2);
198
199 MMuonSearchParCalc muscalc;
200 //muscalc.SetFilter(&fmuon1);
201
202 MMuonCalibParCalc mcalc;
203 mcalc.SetFilter(&fmuon2);
204
205 MFillH fillmuon("MHSingleMuon", "", "FillMuon");
206 MFillH fillmpar("MHMuonPar", "", "FillMuonPar");
207 fillmuon.SetFilter(&fmuon2);
208 fillmpar.SetFilter(&fmuon3);
209 fillmuon.SetBit(MFillH::kDoNotDisplay);
210
211 // ---------------------------------------------------------------
212
213 MWriteRootFile write(fname, "RECREATE", "Image parameters", 2); // EffectiveOnTime
214 write.AddContainer("MTime", "Events");
215 write.AddContainer("MHillas", "Events");
216 write.AddContainer("MHillasExt", "Events");
217 write.AddContainer("MHillasSrc", "Events");
218 write.AddContainer("MImagePar", "Events");
219 write.AddContainer("MNewImagePar", "Events");
220 write.AddContainer("MRawEvtHeader", "Events");
221 write.AddContainer("MRawRunHeader", "RunHeaders");
222 write.AddContainer("MGeomCam", "RunHeaders");
223
224 MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0",
225 "MuonHistCut");
226 fillmhn2.SetFilter(&fmuonhn);
227
228 tlist.AddToList(&read);
229 tlist.AddToList(&cont);
230 tlist.AddToList(&apply);
231 tlist.AddToList(&clean);
232 tlist.AddToList(&hcalc);
233 tlist.AddToList(&fillC1);
234 tlist.AddToList(&fillC2);
235 tlist.AddToList(&fillC3);
236 tlist.AddToList(&fillR);
237 tlist.AddToList(&fillD1);
238 tlist.AddToList(&fillD2);
239 tlist.AddToList(&fillD3);
240 tlist.AddToList(&fillD4);
241 tlist.AddToList(&fillD5);
242 tlist.AddToList(&fillD6);
243 tlist.AddToList(&filla);
244 //tlist.AddToList(&fmuon1);
245 tlist.AddToList(&muscalc);
246 tlist.AddToList(&fillmhn1);
247 tlist.AddToList(&fillmhn3);
248 tlist.AddToList(&fmuon2);
249 tlist.AddToList(&fillmuon);
250 tlist.AddToList(&mcalc);
251 tlist.AddToList(&fmuonhn);
252 tlist.AddToList(&fillmhn2);
253 tlist.AddToList(&fillmhn4);
254 tlist.AddToList(&fmuon3);
255 tlist.AddToList(&fillmpar);
256 //tlist.AddToList(&fmuon4);
257 //tlist.AddToList(&writem);
258 tlist.AddToList(&write);
259
260 if (!loop.Eventloop())
261 return 3;
262
263 if (!loop.GetDisplay())
264 return 4;
265
266 // ============================================================
267
268 TFile *ofile = gROOT->GetListOfFiles()->FindObject(fname);
269 if (!ofile || !ofile->IsOpen() || ofile->IsZombie())
270 {
271 gLog << err << "File " << fname << " not found" << endl;
272 return 20;
273 }
274
275 TString title = "-- Image parameters [";
276 title += datafile;
277 title += "] --";
278 d->SetTitle(title, kFALSE);
279
280 ofile->cd();
281 d->Write();
282
283 return 0;
284}
Note: See TracBrowser for help on using the repository browser.