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

Last change on this file since 19971 was 19266, checked in by tbretz, 6 years ago
Fixed some minor issues with the cling based interpreter - this should not have any effect at runtime maybe except stability
File size: 13.5 KB
Line 
1#include <sstream>
2#include <iostream>
3
4#include "MLog.h"
5#include "MLogManip.h"
6
7#if (!defined(__CINT__) && !defined(__CLING__)) || defined(__MAKECINT__)
8
9#error test1
10
11#include "TH1F.h"
12#include "TFile.h"
13#include "TStyle.h"
14#include "TGraph.h"
15#include "TLine.h"
16
17#include "../mcore/DrsCalib.h"
18#include "MDrsCalibration.h"
19//#include "MLogManip.h"
20#include "MExtralgoSpline.h"
21#include "MSequence.h"
22#include "MStatusArray.h"
23#include "MHCamera.h"
24#include "MJob.h"
25#include "MWriteRootFile.h"
26#include "MHCamera.h"
27#include "MBadPixelsCam.h"
28#include "MBadPixelsPix.h"
29#include "MDirIter.h"
30#include "MTaskList.h"
31#include "MFDataPhrase.h"
32#include "MArrayF.h"
33#include "MBadPixelsTreat.h"
34#include "MCalibrateDrsTimes.h"
35#include "MHSectorVsTime.h"
36#include "MHCamEvent.h"
37#include "MExtractTimeAndChargeSpline.h"
38#include "MFillH.h"
39#include "MDrsCalibApply.h"
40#include "MGeomApply.h"
41#include "MContinue.h"
42#include "MRawFitsRead.h"
43#include "MEvtLoop.h"
44#include "MParList.h"
45#include "MStatusDisplay.h"
46#include "MDrsCalibrationTime.h"
47#include "MH3.h"
48#include "MGeomCamFACT.h"
49#include "MCalibrateFact.h"
50#include "MParameters.h"
51#include "MWriteAsciiFile.h"
52#include "MFMagicCuts.h"
53
54#include "MMuonSetup.h"
55#include "MReadMarsFile.h"
56#include "MHillasCalc.h"
57#include "MHn.h"
58#include "MMuonSearchParCalc.h"
59#include "MMuonCalibParCalc.h"
60#include "MBinning.h"
61#include "MImgCleanStd.h"
62
63#endif
64
65using namespace std;
66
67int star(const TString callistofile, TString outfile,
68 TString displayfile, TString displaytitle,
69 const TString txtfile,
70 Double_t lvl1=4.0, Double_t lvl2=2.5, Double_t deltat = 17.5);
71
72
73int star(const TString infile = "callisto.root", TString outfile = "",
74 Double_t lvl1=4.0, Double_t lvl2=2.5, Double_t deltat = 17.5)
75{
76 return star(infile, outfile, "", "", "", lvl1, lvl2, deltat);
77}
78
79int star(const TString callistofile, TString outfile,
80 TString displayfile, TString displaytitle,
81 const TString txtfile,
82 Double_t lvl1, Double_t lvl2, Double_t deltat)
83{
84
85 //Double_t lvl1=7.8,
86 //Double_t lvl2=3.9,
87 // lvl1 = 2.5;
88 // lvl2 = 0.5;
89 // ------------------------------------------------------
90
91 if (displaytitle.IsNull())
92 displaytitle = gSystem->BaseName(callistofile);
93
94 FileStat_t fstat;
95 int rc = gSystem->GetPathInfo(outfile, fstat);
96 bool isdir = !rc || R_ISDIR(fstat.fMode);
97
98 char *buf = gSystem->ConcatFileName(outfile, "star.root");
99 outfile = buf;
100 delete [] buf;
101
102 if (displayfile.IsNull())
103 {
104 displayfile = outfile;
105 displayfile.Insert(displayfile.Last('.'), "-display");
106 }
107 else
108 {
109 if (isdir && gSystem->DirName(displayfile)==TString("."))
110 {
111 buf = gSystem->ConcatFileName(outfile, displayfile);
112 displayfile = buf;
113 delete [] buf;
114 }
115 }
116
117 // ------------------------------------------------------
118
119 gLog.Separator("Star");
120 gLog << all << "Calculate image parameters of sequence ";
121 gLog << endl;
122 gLog << "Callisto file: " << callistofile << endl;
123 gLog << "Output file: " << outfile << endl;
124 gLog << "Display file: " << displayfile << endl;
125 gLog << "Display title: " << displaytitle << endl;
126 gLog << endl;
127 gLog << "lvl1: " << lvl1 <<endl;
128 gLog << "lvl2: " << lvl2 <<endl;
129 gLog << "deltat: " << deltat <<endl;
130
131
132 gLog << endl;
133
134 // ------------------------------------------------------
135
136 gLog.Separator();
137
138 // --------------------------------------------------------
139
140 MBinning bins1( 36, -90, 90, "BinningAlpha");
141 MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
142 MBinning bins4( 25, 0, 2.5, "BinningDist");
143 MBinning binsM1(100, 0, 5, "BinningMuonRadius");
144 MBinning binsM2( 60, 0, 0.3, "BinningMuonDeviation");
145 MBinning binsM3(150, 0, 60, "BinningMuonTime");
146 MBinning binsM4(300, 0, 30, "BinningMuonTimeRms");
147 MBinning binsM5( 20, 0, 360, "BinningMuonArcPhi");
148 MBinning binsM6( 50, 0, 0.5, "BinningMuonArcWidth");
149 MBinning binsM7( 30, 5, 5000, "BinningMuonSize", "log");
150 MBinning binsM8(100, 0, 5, "BinningMuonRelTimeSigma");
151 MBinning binsM9(100, 0, 5, "BinningRadius");
152
153 // ---------------------------------------------------------
154
155 // Filter to start muon analysis
156 //MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut");
157 // Filter to calculate further muon parameters
158// MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&"
159// "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut");
160 MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&"
161 "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4",
162 "MuonSearchCut");
163
164 // Filter to fill the MHMuonPar
165 MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2",
166 "MuonFinalCut");
167 // Filter to write Muons to Muon tree
168 //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
169
170 // ---------------------------------------------------------
171
172 MStatusDisplay *d = new MStatusDisplay;
173
174 MTaskList tlist;
175
176 MParList plist2;
177 plist2.AddToList(&tlist);
178 plist2.AddToList(&bins1);
179 plist2.AddToList(&bins3);
180 plist2.AddToList(&bins4);
181 plist2.AddToList(&binsM1);
182 plist2.AddToList(&binsM2);
183 plist2.AddToList(&binsM3);
184 plist2.AddToList(&binsM4);
185 plist2.AddToList(&binsM5);
186 plist2.AddToList(&binsM6);
187 plist2.AddToList(&binsM7);
188 plist2.AddToList(&binsM8);
189 plist2.AddToList(&binsM9);
190
191 MMuonSetup muonsetup;
192 plist2.AddToList(&muonsetup);
193
194 MEvtLoop loop;
195 loop.SetDisplay(d);
196 loop.SetParList(&plist2);
197
198 MReadMarsFile read("Events");
199 read.DisableAutoScheme();
200 read.AddFile(callistofile);
201
202 MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData");
203
204 MGeomApply apply;
205
206 MImgCleanStd clean(lvl1, lvl2);
207 clean.SetMethod(MImgCleanStd::kAbsolute);
208 clean.SetTimeLvl1(deltat);//17.5/*1.3*/);
209 clean.SetTimeLvl2(deltat);//17.5/*1.3*/);
210 clean.SetPostCleanType(3);
211
212 MHillasCalc hcalc;
213 hcalc.Disable(MHillasCalc::kCalcConc);
214
215 MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
216 hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1.)");
217 hrate.SetName("Rate");
218 hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];");
219 hrate.InitLabels(MH3::kLabelsX);
220 hrate.GetHist().SetMinimum(0);
221
222 MFillH fillR(&hrate, "", "FillRate");
223
224 // ------------------ Setup histograms and fill tasks ----------------
225 MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]");
226 MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
227 MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]");
228 evtC1.SetErrorSpread(kFALSE);
229 evtC2.SetErrorSpread(kFALSE);
230 evtC2.SetThreshold(0);
231
232 MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam");
233 MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels");
234 MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos");
235
236 MFillH fillD1("MHHillas", "MHillas", "FillHillas");
237 MFillH fillD2("MHHillasExt", "", "FillHillasExt");
238 MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
239 MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar");
240 MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
241 MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize");
242
243 MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)");
244 halpha.SetName("AvsD;Dist;Alpha");
245
246 MFillH filla(&halpha);
247 filla.SetNameTab("AvsD");
248 filla.SetDrawOption("colz");
249
250 // ----------------------- Muon Analysis ----------------------
251 //writem.SetFilter(&fmuon4);
252
253 MHn mhn1;
254 mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg");
255 mhn1.InitName("MuonRadius");
256 mhn1.InitTitle("MuonRadius");
257
258 mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
259 mhn1.InitName("MuonDeviation");
260 mhn1.InitTitle("MuonDeviation");
261
262 mhn1.AddHist("MMuonSearchPar.fTime");
263 mhn1.InitName("MuonTime");
264 mhn1.InitTitle("MuonTime");
265
266 mhn1.AddHist("MMuonSearchPar.fTimeRms");
267 mhn1.InitName("MuonTimeRms");
268 mhn1.InitTitle("MuonTimeRms");
269
270 MHn mhn3;
271 mhn3.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
272 mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation");
273 mhn3.InitTitle("MuonRadius");
274 mhn3.SetDrawOption("colz");
275
276 MHn mhn2;
277 mhn2.AddHist("MMuonCalibPar.fArcPhi");
278 mhn2.InitName("MuonArcPhi");
279 mhn2.InitTitle("MuonArcPhi");
280
281 mhn2.AddHist("MMuonCalibPar.fArcWidth");
282 mhn2.InitName("MuonArcWidth");
283 mhn2.InitTitle("MuonArcWidth");
284
285 mhn2.AddHist("MMuonCalibPar.fMuonSize");
286 mhn2.InitName("MuonSize");
287 mhn2.InitTitle("MuonSize");
288
289 mhn2.AddHist("MMuonCalibPar.fRelTimeSigma");
290 mhn2.InitName("MuonRelTimeSigma");
291 mhn2.InitTitle("MuonRelTimeSigma");
292
293 MHn mhn4;
294 mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth");
295 mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth");
296 mhn4.InitTitle("MuonArcPhi");
297 mhn4.SetDrawOption("colz");
298
299 MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar");
300 MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar");
301 MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius");
302 MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi");
303 fillmhn1.SetNameTab("MuonS");
304 fillmhn3.SetNameTab("DevVsRad");
305 fillmhn2.SetNameTab("MuonC");
306 fillmhn4.SetNameTab("WidthVsPhi");
307 fillmhn2.SetFilter(&fmuon2);
308 fillmhn4.SetFilter(&fmuon2);
309
310 MMuonSearchParCalc muscalc;
311 //muscalc.SetFilter(&fmuon1);
312
313 MMuonCalibParCalc mcalc;
314 mcalc.SetFilter(&fmuon2);
315
316 MFillH fillmuon("MHSingleMuon", "", "FillMuon");
317 MFillH fillmpar("MHMuonPar", "", "FillMuonPar");
318 fillmuon.SetFilter(&fmuon2);
319 fillmpar.SetFilter(&fmuon3);
320 fillmuon.SetBit(MFillH::kDoNotDisplay);
321
322 // ---------------------------------------------------------------
323
324 MWriteAsciiFile writeascii(txtfile, "MHillas");
325 writeascii.AddColumns("MHillasExt");
326 writeascii.AddColumns("MHillasSrc");
327
328
329 MWriteRootFile write(outfile, "RECREATE", "Image parameters", 2);
330 write.AddContainer("MTime", "Events", kFALSE);
331 write.AddContainer("MHillas", "Events");
332 write.AddContainer("MHillasExt", "Events");
333 write.AddContainer("MHillasSrc", "Events");
334 write.AddContainer("MImagePar", "Events");
335 write.AddContainer("MNewImagePar", "Events");
336 write.AddContainer("MRawEvtHeader", "Events");
337 write.AddContainer("ThetaSquared", "Events");
338 write.AddContainer("Disp", "Events");
339 write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
340 write.AddContainer("MMcEvt", "Events", kFALSE);
341 write.AddContainer("IncidentAngle", "Events", kFALSE);
342 write.AddContainer("MPointingPos", "Events", kFALSE);
343
344 write.AddContainer("MRawRunHeader", "RunHeaders");
345 write.AddContainer("MGeomCam", "RunHeaders");
346 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
347 write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
348 write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
349
350 MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0",
351 "MuonHistCut");
352 fillmhn2.SetFilter(&fmuonhn);
353
354 TArrayD param(12);
355 // Parametrization of Disp
356 param[0] = 1.15136; // constant
357 param[8] = 0.0681437; // slope
358 param[9] = 2.62932; // leak
359 param[10] = 1.51279; // size-offset
360 param[11] = 0.0507821; // size-slope
361 // Parametrization for sign of disp (m3long, slope)
362 param[5] = -0.07;
363 param[6] = 7.2;
364 param[7] = 0.5;
365 // ThetaSq-Cut
366 param[1] = 0.11; // 0.215
367 // Area-Cut
368 param[2] = 0.215468;
369 param[3] = 5.63973;
370 param[4] = 0.0836169;
371 MFMagicCuts cuts;
372 cuts.SetHadronnessCut(MFMagicCuts::kNoCut);
373 cuts.SetVariables(param);
374
375
376 tlist.AddToList(&read);
377 tlist.AddToList(&cont);
378 tlist.AddToList(&apply);
379 tlist.AddToList(&clean);
380 tlist.AddToList(&hcalc);
381 tlist.AddToList(&cuts);
382 tlist.AddToList(&fillC1);
383 tlist.AddToList(&fillC2);
384 tlist.AddToList(&fillC3);
385 tlist.AddToList(&fillR);
386 tlist.AddToList(&fillD1);
387 tlist.AddToList(&fillD2);
388 tlist.AddToList(&fillD3);
389 tlist.AddToList(&fillD4);
390 tlist.AddToList(&fillD5);
391 tlist.AddToList(&fillD6);
392 tlist.AddToList(&filla);
393 //tlist.AddToList(&fmuon1);
394 tlist.AddToList(&muscalc);
395 tlist.AddToList(&fillmhn1);
396 tlist.AddToList(&fillmhn3);
397 tlist.AddToList(&fmuon2);
398 tlist.AddToList(&fillmuon);
399 tlist.AddToList(&mcalc);
400 tlist.AddToList(&fmuonhn);
401 tlist.AddToList(&fillmhn2);
402 tlist.AddToList(&fillmhn4);
403 tlist.AddToList(&fmuon3);
404 tlist.AddToList(&fillmpar);
405 //tlist.AddToList(&fmuon4);
406 //tlist.AddToList(&writem);
407 tlist.AddToList(&write);
408 if (!txtfile.IsNull())
409 tlist.AddToList(&writeascii);
410
411 if (!loop.Eventloop())
412 return 3;
413
414 if (!loop.GetDisplay())
415 return 4;
416
417 d->SetTitle(displaytitle, kFALSE);
418 d->SaveAs(displayfile);
419
420 return 0;
421}
422
Note: See TracBrowser for help on using the repository browser.