source: branches/Mars_McMismatchStudy/fact/analysis/star.C

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