source: branches/Mars_MC/fact/analysis/mc/star.C@ 17839

Last change on this file since 17839 was 17056, checked in by dneise, 11 years ago
edited both scripts, so we can use them on phido for processing of simulated data
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.