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

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