source: branches/Mars_McMismatchStudy/fact/analysis/star.C.ORIGINAL@ 20029

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