source: branches/Mars_IncreaseNsb/fact/analysis/star_file.C@ 20041

Last change on this file since 20041 was 18578, checked in by tbretz, 9 years ago
MSoftwareTrigger is not required, added missing containers to the output produced by the MC chain
File size: 10.3 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("MSoftwareTrigger","Events", kFALSE);
222 write.AddContainer("MRawRunHeader", "RunHeaders");
223 write.AddContainer("MGeomCam", "RunHeaders");
224
225 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
226 write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
227 write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
228 write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
229 write.AddContainer("MMcEvt", "Events", kFALSE);
230 write.AddContainer("IncidentAngle", "Events", kFALSE);
231 write.AddContainer("MPointingPos", "Events", kFALSE);
232 write.AddContainer("MTime", "Events", kFALSE);
233 write.AddContainer("MSrcPosCam", "Events", kFALSE);
234
235 MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0",
236 "MuonHistCut");
237 fillmhn2.SetFilter(&fmuonhn);
238
239 tlist.AddToList(&read);
240 tlist.AddToList(&cont);
241 tlist.AddToList(&apply);
242 tlist.AddToList(&clean);
243 tlist.AddToList(&hcalc);
244 tlist.AddToList(&fillC1);
245 tlist.AddToList(&fillC2);
246 tlist.AddToList(&fillC3);
247 tlist.AddToList(&fillR);
248 tlist.AddToList(&fillD1);
249 tlist.AddToList(&fillD2);
250 tlist.AddToList(&fillD3);
251 tlist.AddToList(&fillD4);
252 tlist.AddToList(&fillD5);
253 tlist.AddToList(&fillD6);
254 tlist.AddToList(&filla);
255 //tlist.AddToList(&fmuon1);
256 tlist.AddToList(&muscalc);
257 tlist.AddToList(&fillmhn1);
258 tlist.AddToList(&fillmhn3);
259 tlist.AddToList(&fmuon2);
260 tlist.AddToList(&fillmuon);
261 tlist.AddToList(&mcalc);
262 tlist.AddToList(&fmuonhn);
263 tlist.AddToList(&fillmhn2);
264 tlist.AddToList(&fillmhn4);
265 tlist.AddToList(&fmuon3);
266 tlist.AddToList(&fillmpar);
267 //tlist.AddToList(&fmuon4);
268 //tlist.AddToList(&writem);
269 tlist.AddToList(&write);
270
271 if (!loop.Eventloop())
272 return 3;
273
274 if (!loop.GetDisplay())
275 return 4;
276
277 // ============================================================
278
279 TFile *ofile = gROOT->GetListOfFiles()->FindObject(fname);
280 if (!ofile || !ofile->IsOpen() || ofile->IsZombie())
281 {
282 gLog << err << "File " << fname << " not found" << endl;
283 return 20;
284 }
285
286 TString title = "-- Image parameters [";
287 title += datafile;
288 title += "] --";
289 d->SetTitle(title, kFALSE);
290
291 ofile->cd();
292 d->Write();
293
294 return 0;
295}
Note: See TracBrowser for help on using the repository browser.