source: trunk/Mars/fact/analysis/mc/star.C@ 17309

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