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

Last change on this file since 17885 was 17885, checked in by tbretz, 10 years ago
Added calls which allow to use the sequence number instead of the sequence file.
File size: 10.4 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}
298
299int star(const ULong64_t seqnum, Double_t lvl1=4.0, Double_t lvl2=2.5, const char *inpath = "output", const char *outpath = "output")
300{
301 UInt_t night = seqnum/1000;
302 UInt_t num = seqnum%1000;
303
304 TString file = Form("/scratch/fact/sequences/%04d/%02d/%02d/%06d_%03d.seq",
305 night/10000, (night/100)%100, night%100, num);
306
307 return star(file.Data(), lvl1, lvl2, inpath, outpath);
308}
Note: See TracBrowser for help on using the repository browser.