Context Navigation

source: branches/MarsMoreSimulationTruth/fact/analysis/star_file.C@ 19829

Visit:

Last change on this file since 19829 was 18578, checked in by tbretz, 10 years ago
MSoftwareTrigger is not required, added missing containers to the output produced by the MC chain
File size: 10.3 KB

Line
1	#include "MLogManip.h"
2
3	int star_file(const char datafile, Double_t lvl1=7.8, Double_t lvl2=3.9, const char outpath = ".")
4	{
5	double deltat = 17.5; // ns/deg
6	// lvl1 = 2.5;
7	// lvl2 = 0.5;
8
9	// ------------------------------------------------------
10
11	gLog.Separator("Star");
12	gLog << all << "Calculate image parameters of sequence ";
13	gLog << datafile << endl;
14	gLog << endl;
15
16	// ------------------------------------------------------
17
18	gLog << "Outpath: " << outpath << endl;
19
20	TString fname = gSystem->ConcatFileName(outpath, gSystem->BaseName(datafile));
21	fname.ReplaceAll("_C.root", "_I.root");
22
23	gLog.Separator();
24
25	// --------------------------------------------------------
26
27	MBinning bins1( 36, -90, 90, "BinningAlpha");
28	MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
29	MBinning bins4( 25, 0, 2.5, "BinningDist");
30
31	MBinning binsM1(100, 0, 5, "BinningMuonRadius");
32	MBinning binsM2( 60, 0, 0.3, "BinningMuonDeviation");
33	MBinning binsM3(150, 0, 60, "BinningMuonTime");
34	MBinning binsM4(300, 0, 30, "BinningMuonTimeRms");
35	MBinning binsM5( 20, 0, 360, "BinningMuonArcPhi");
36	MBinning binsM6( 50, 0, 0.5, "BinningMuonArcWidth");
37	MBinning binsM7( 30, 5, 5000, "BinningMuonSize", "log");
38	MBinning binsM8(100, 0, 5, "BinningMuonRelTimeSigma");
39
40	MBinning binsM9(100, 0, 5, "BinningRadius");
41
42	// ---------------------------------------------------------
43
44	// Filter to start muon analysis
45	//MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut");
46	// Filter to calculate further muon parameters
47	// MFDataPhrase fmuon2("(MMuonSearchPar.fRadiusMGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadiusMGeomCam.fConvMm2Deg<1.35) &&"
48	// "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut");
49	MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&"
50	"MMuonSearchPar.fDeviationMGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadiusMGeomCam.fConvMm2Deg-0.1)*0.15/0.4",
51	"MuonSearchCut");
52
53	// Filter to fill the MHMuonPar
54	MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2",
55	"MuonFinalCut");
56	// Filter to write Muons to Muon tree
57	//MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
58
59	// ---------------------------------------------------------
60
61	MStatusDisplay *d = new MStatusDisplay;
62
63	MTaskList tlist;
64
65	MParList plist2;
66	plist2.AddToList(&tlist);
67	plist2.AddToList(&bins1);
68	plist2.AddToList(&bins3);
69	plist2.AddToList(&bins4);
70	plist2.AddToList(&binsM1);
71	plist2.AddToList(&binsM2);
72	plist2.AddToList(&binsM3);
73	plist2.AddToList(&binsM4);
74	plist2.AddToList(&binsM5);
75	plist2.AddToList(&binsM6);
76	plist2.AddToList(&binsM7);
77	plist2.AddToList(&binsM8);
78	plist2.AddToList(&binsM9);
79
80	MMuonSetup muonsetup;
81	plist2.AddToList(&muonsetup);
82
83	MEvtLoop loop;
84	loop.SetDisplay(d);
85	loop.SetParList(&plist2);
86
87	MReadMarsFile read("Events");
88	read.DisableAutoScheme();
89	read.AddFile(datafile);
90
91	MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData");
92
93	MGeomApply apply;
94
95	MImgCleanStd clean(lvl1, lvl2);
96	clean.SetMethod(MImgCleanStd::kAbsolute);
97	clean.SetTimeLvl1(deltat);//17.5/1.3/);
98	clean.SetTimeLvl2(deltat);//17.5/1.3/);
99	clean.SetPostCleanType(3);
100
101	MHillasCalc hcalc;
102	hcalc.Disable(MHillasCalc::kCalcConc);
103
104	MH3 hrate("MRawRunHeader.GetFileID"/, "MRawEvtHeader.GetTriggerID"/);
105	hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
106	hrate.SetName("Rate");
107	hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];");
108	hrate.InitLabels(MH3::kLabelsX);
109	hrate.GetHist().SetMinimum(0);
110
111	MFillH fillR(&hrate, "", "FillRate");
112
113	// ------------------ Setup histograms and fill tasks ----------------
114	MHCamEvent evtC1(0, "Cleaned", "Average signal after Cleaning;;S [phe]");
115	MHCamEvent evtC2(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
116	MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]");
117	evtC1.SetErrorSpread(kFALSE);
118	evtC2.SetErrorSpread(kFALSE);
119	evtC2.SetThreshold(0);
120
121	MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam");
122	MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels");
123	MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos");
124
125	MFillH fillD1("MHHillas", "MHillas", "FillHillas");
126	MFillH fillD2("MHHillasExt", "", "FillHillasExt");
127	MFillH fillD3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
128	MFillH fillD4("MHImagePar", "MImagePar", "FillImagePar");
129	MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
130	MFillH fillD6("MHVsSize", "MHillasSrc", "FillVsSize");
131
132	MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)");
133	halpha.SetName("AvsD;Dist;Alpha");
134
135	MFillH filla(&halpha);
136	filla.SetNameTab("AvsD");
137	filla.SetDrawOption("colz");
138
139	// ----------------------- Muon Analysis ----------------------
140	//writem.SetFilter(&fmuon4);
141
142	MHn mhn1;
143	mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg");
144	mhn1.InitName("MuonRadius");
145	mhn1.InitTitle("MuonRadius");
146
147	mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
148	mhn1.InitName("MuonDeviation");
149	mhn1.InitTitle("MuonDeviation");
150
151	mhn1.AddHist("MMuonSearchPar.fTime");
152	mhn1.InitName("MuonTime");
153	mhn1.InitTitle("MuonTime");
154
155	mhn1.AddHist("MMuonSearchPar.fTimeRms");
156	mhn1.InitName("MuonTimeRms");
157	mhn1.InitTitle("MuonTimeRms");
158
159	MHn mhn3;
160	mhn3.AddHist("MMuonSearchPar.fRadiusMGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviationMGeomCam.fConvMm2Deg");
161	mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation");
162	mhn3.InitTitle("MuonRadius");
163	mhn3.SetDrawOption("colz");
164
165	MHn mhn2;
166	mhn2.AddHist("MMuonCalibPar.fArcPhi");
167	mhn2.InitName("MuonArcPhi");
168	mhn2.InitTitle("MuonArcPhi");
169
170	mhn2.AddHist("MMuonCalibPar.fArcWidth");
171	mhn2.InitName("MuonArcWidth");
172	mhn2.InitTitle("MuonArcWidth");
173
174	mhn2.AddHist("MMuonCalibPar.fMuonSize");
175	mhn2.InitName("MuonSize");
176	mhn2.InitTitle("MuonSize");
177
178	mhn2.AddHist("MMuonCalibPar.fRelTimeSigma");
179	mhn2.InitName("MuonRelTimeSigma");
180	mhn2.InitTitle("MuonRelTimeSigma");
181
182	MHn mhn4;
183	mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth");
184	mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth");
185	mhn4.InitTitle("MuonArcPhi");
186	mhn4.SetDrawOption("colz");
187
188	MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar");
189	MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar");
190	MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius");
191	MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi");
192	fillmhn1.SetNameTab("MuonS");
193	fillmhn3.SetNameTab("DevVsRad");
194	fillmhn2.SetNameTab("MuonC");
195	fillmhn4.SetNameTab("WidthVsPhi");
196	fillmhn2.SetFilter(&fmuon2);
197	fillmhn4.SetFilter(&fmuon2);
198
199	MMuonSearchParCalc muscalc;
200	//muscalc.SetFilter(&fmuon1);
201
202	MMuonCalibParCalc mcalc;
203	mcalc.SetFilter(&fmuon2);
204
205	MFillH fillmuon("MHSingleMuon", "", "FillMuon");
206	MFillH fillmpar("MHMuonPar", "", "FillMuonPar");
207	fillmuon.SetFilter(&fmuon2);
208	fillmpar.SetFilter(&fmuon3);
209	fillmuon.SetBit(MFillH::kDoNotDisplay);
210
211	// ---------------------------------------------------------------
212
213	MWriteRootFile write(fname, "RECREATE", "Image parameters", 2); // EffectiveOnTime
214	write.AddContainer("MTime", "Events");
215	write.AddContainer("MHillas", "Events");
216	write.AddContainer("MHillasExt", "Events");
217	write.AddContainer("MHillasSrc", "Events");
218	write.AddContainer("MImagePar", "Events");
219	write.AddContainer("MNewImagePar", "Events");
220	write.AddContainer("MRawEvtHeader", "Events");
221	write.AddContainer("MSoftwareTrigger","Events", kFALSE);
222	write.AddContainer("MRawRunHeader", "RunHeaders");
223	write.AddContainer("MGeomCam", "RunHeaders");
224
225	write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
226	write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
227	write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
228	write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
229	write.AddContainer("MMcEvt", "Events", kFALSE);
230	write.AddContainer("IncidentAngle", "Events", kFALSE);
231	write.AddContainer("MPointingPos", "Events", kFALSE);
232	write.AddContainer("MTime", "Events", kFALSE);
233	write.AddContainer("MSrcPosCam", "Events", kFALSE);
234
235	MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0",
236	"MuonHistCut");
237	fillmhn2.SetFilter(&fmuonhn);
238
239	tlist.AddToList(&read);
240	tlist.AddToList(&cont);
241	tlist.AddToList(&apply);
242	tlist.AddToList(&clean);
243	tlist.AddToList(&hcalc);
244	tlist.AddToList(&fillC1);
245	tlist.AddToList(&fillC2);
246	tlist.AddToList(&fillC3);
247	tlist.AddToList(&fillR);
248	tlist.AddToList(&fillD1);
249	tlist.AddToList(&fillD2);
250	tlist.AddToList(&fillD3);
251	tlist.AddToList(&fillD4);
252	tlist.AddToList(&fillD5);
253	tlist.AddToList(&fillD6);
254	tlist.AddToList(&filla);
255	//tlist.AddToList(&fmuon1);
256	tlist.AddToList(&muscalc);
257	tlist.AddToList(&fillmhn1);
258	tlist.AddToList(&fillmhn3);
259	tlist.AddToList(&fmuon2);
260	tlist.AddToList(&fillmuon);
261	tlist.AddToList(&mcalc);
262	tlist.AddToList(&fmuonhn);
263	tlist.AddToList(&fillmhn2);
264	tlist.AddToList(&fillmhn4);
265	tlist.AddToList(&fmuon3);
266	tlist.AddToList(&fillmpar);
267	//tlist.AddToList(&fmuon4);
268	//tlist.AddToList(&writem);
269	tlist.AddToList(&write);
270
271	if (!loop.Eventloop())
272	return 3;
273
274	if (!loop.GetDisplay())
275	return 4;
276
277	// ============================================================
278
279	TFile *ofile = gROOT->GetListOfFiles()->FindObject(fname);
280	if (!ofile \|\| !ofile->IsOpen() \|\| ofile->IsZombie())
281	{
282	gLog << err << "File " << fname << " not found" << endl;
283	return 20;
284	}
285
286	TString title = "-- Image parameters [";
287	title += datafile;
288	title += "] --";
289	d->SetTitle(title, kFALSE);
290
291	ofile->cd();
292	d->Write();
293
294	return 0;
295	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: