source: branches/Mars_MC/fact/analysis/sandbox_dneise/star.C@ 17430

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