source: trunk/Mars/fact/analysis/star.C@ 17667

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