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

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