source: branches/MarsISDCBranchBasedOn17887/fact/analysis/mc/star_NEWCLEANING.C@ 19084

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