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

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