source: trunk/MagicSoft/Mars/mtemp/mifae/programs/controlPlot.cc@ 6363

Last change on this file since 6363 was 5879, checked in by rico, 21 years ago
*** empty log message ***
File size: 9.4 KB
Line 
1////////////////////////////////////////////////////////////////////////////////////
2//
3// _____ Control Plot program_____
4//
5// Take as input root files with hillas parameters and generates several
6// control plots.
7// Output is a file with :
8// - Number of hits plot
9//
10// Ester Aliu <aliu@ifae.es>
11// Oscar Blanch <blanch@ifae.es>
12// Javier Rico <jrico@ifae.es>
13////////////////////////////////////////////////////////////////////////////////////
14
15#include <fstream>
16#include <iostream>
17
18#include "TString.h"
19
20#include "MHillasSrcCalc.h"
21#include "MHillasSrc.h"
22#include "MSrcRotate.h"
23#include "MSrcPosFromFile.h"
24#include "MSrcTranslate.h"
25#include "MParList.h"
26#include "MTaskList.h"
27#include "MHillas.h"
28#include "MReadTree.h"
29#include "MEvtLoop.h"
30#include "MLog.h"
31#include "MArgs.h"
32#include "MWriteRootFile.h"
33#include "MTime.h"
34#include "MControlPlots.h"
35#include "MIslands.h"
36#include "MContinue.h"
37#include "MF.h"
38
39using namespace std;
40
41Bool_t readDatacards(TString& filename);
42void controlPlot();
43
44//-----------------------------------------------------------------------------
45// declaration of variables read from datacards
46//-----------------------------------------------------------------------------
47
48TString inputFile;
49TString offFile;
50TString outputFile;
51TString filter;
52TString sizecut;
53TString distcut;
54TString widthcut;
55TString lengthcut;
56TString centercut;
57ULong_t nmaxevents=999999999;
58Float_t xsrcpos=0.;
59Float_t ysrcpos=0.;
60Double_t mjdpos=0;
61
62//-----------------------------------------------------------------------------
63// constants
64//-----------------------------------------------------------------------------
65
66const TString defaultcard="controlplot.datacard";
67const Float_t conver = 189./0.6; // conversion factor degrees to mm
68
69//-----------------------------------------------------------------------------
70
71static void Usage()
72{
73 gLog <<endl;
74 gLog << "Usage is:" << endl;
75 gLog << " controlPlot [-h] [-?] <datacards>" << endl << endl;
76 gLog << " <datacards>: datacards file name (dafault " << defaultcard <<")" << endl;
77 gLog << " -?/-h: This help" << endl << endl;
78}
79
80//-----------------------------------------------------------------------------
81int main(int argc, char **argv)
82{
83 // evaluate arguments
84 MArgs arg(argc, argv);
85 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
86 {
87 Usage();
88 return -1;
89 }
90
91 TString datacard = arg.GetArgumentStr(0);
92 if(!datacard.Length())
93 datacard = defaultcard;
94
95 if(!readDatacards(datacard))
96 {
97 cout << "Error reading datacards. Stoping" << endl;
98 return -1;
99 }
100
101 controlPlot();
102}
103
104//-----------------------------------------------------------------------------
105void controlPlot()
106{
107 // containers
108 MParList plist;
109 MTaskList tlist;
110 MIslands islands;
111
112 // include containers in parameter list
113 plist.AddToList(&tlist);
114
115 MF cut(filter);
116 MF cutS(sizecut);
117 MF cutD(distcut);
118 MF cutL(lengthcut);
119 MF cutW(widthcut);
120 MF cutC(centercut);
121
122 // tasks
123 MReadTree read("Parameters", inputFile);
124 read.DisableAutoScheme();
125
126 MControlPlots controlplots(outputFile);
127 controlplots.SetProduceFile(kFALSE);
128 controlplots.SetMode(MControlPlots::kOn);
129
130 MContinue applycutS(&cutS);
131 applycutS.SetInverted(kTRUE);
132
133 MContinue applycutD(&cutD);
134 applycutD.SetInverted(kTRUE);
135
136 MContinue applycutL(&cutL);
137 applycutL.SetInverted(kTRUE);
138
139 MContinue applycutW(&cutW);
140 applycutW.SetInverted(kTRUE);
141
142 MContinue applycutC(&cutC);
143 applycutC.SetInverted(kTRUE);
144
145 MContinue applycut(&cut);
146 applycut.SetInverted(kTRUE);
147
148
149 // include tasks in task list
150 tlist.AddToList(&read);
151 if(sizecut.Length())
152 tlist.AddToList(&applycutS);
153 if(distcut.Length())
154 tlist.AddToList(&applycutD);
155 if(lengthcut.Length())
156 tlist.AddToList(&applycutL);
157 if(widthcut.Length())
158 tlist.AddToList(&applycutW);
159 if(filter.Length())
160 tlist.AddToList(&applycut);
161 if(centercut.Length())
162 tlist.AddToList(&applycutC);
163 tlist.AddToList(&controlplots);
164
165 // Eventloop
166 MEvtLoop evtloop;
167 evtloop.SetParList(&plist);
168 if (!evtloop.Eventloop(nmaxevents))
169 return;
170
171 tlist.PrintStatistics();
172
173 // do off-data if input file was specified
174 if(!offFile.Length())
175 return;
176
177
178 // containers
179 MParList plist2;
180 MTaskList tlist2;
181 MIslands islands2;
182
183 // include containers in parameter list
184 plist2.AddToList(&tlist2);
185
186 // tasks
187 MReadTree read2("Parameters", offFile);
188 read2.DisableAutoScheme();
189 tlist2.AddToList(&read2);
190
191 controlplots.SetMode(MControlPlots::kOff);
192 controlplots.SetProduceFile(kTRUE);
193 controlplots.Reset();
194
195 // include tasks in task list
196 if(sizecut.Length())
197 tlist2.AddToList(&applycutS);
198 if(distcut.Length())
199 tlist2.AddToList(&applycutD);
200 if(lengthcut.Length())
201 tlist2.AddToList(&applycutL);
202 if(widthcut.Length())
203 tlist2.AddToList(&applycutW);
204 if(filter.Length())
205 tlist2.AddToList(&applycut);
206 if(centercut.Length())
207 tlist2.AddToList(&applycutC);
208 tlist2.AddToList(&controlplots);
209
210 // Eventloop
211 MEvtLoop evtloop2;
212 evtloop2.SetParList(&plist2);
213 if (!evtloop2.Eventloop(nmaxevents))
214 return;
215
216 tlist2.PrintStatistics();
217
218 // do off-data if input file was specified
219 if(!offFile.Length())
220 return;
221
222}
223
224//-----------------------------------------------------------------------
225
226Bool_t readDatacards(TString& filename)
227{
228 ifstream ifun(filename.Data());
229 if(!ifun)
230 {
231 cout << "File " << filename << " not found" << endl;
232 return kFALSE;
233 }
234
235 TString word;
236
237 while(ifun >> word)
238 {
239 // skip comments
240 if(word[0]=='/' && word[1]=='/')
241 {
242 while(ifun.get()!='\n'); // skip line
243 continue;
244 }
245
246 // number of events
247 if(strcmp(word.Data(),"NEVENTS")==0)
248 ifun >> nmaxevents;
249
250 // input file name
251 if(strcmp(word.Data(),"ONFILES")==0)
252 {
253 if(inputFile.Length())
254 cout << "readDataCards Warning: overriding on-data file name" << endl;
255 ifun >> inputFile;
256 }
257
258 // off-data file name
259 if(strcmp(word.Data(),"OFFFILES")==0)
260 {
261 if(offFile.Length())
262 cout << "readDataCards Warning: overriding off-data file name" << endl;
263 ifun >> offFile;
264 }
265
266 // output file name
267 if(strcmp(word.Data(),"HITSFILE")==0)
268 {
269 if(outputFile.Length())
270 cout << "readDataCards Warning: overriding output file name" << endl;
271 ifun >> outputFile;
272 }
273 // exclusion cut
274 if(strcmp(word.Data(),"FILTER")==0)
275 {
276 if(filter.Length())
277 cout << "readDataCards Warning: overriding existing cut" << endl;
278
279 char ch;
280 while((ch=ifun.get())!='\n')
281 filter.Append(ch);
282 }
283 if(strcmp(word.Data(),"SIZECUT")==0)
284 {
285 if(sizecut.Length())
286 cout << "readDataCards Warning: overriding existing size cut" << endl;
287 char ch[10];
288 sizecut="(MHillas.fSize>";
289 ifun >> ch;
290 sizecut+=ch;
291 sizecut+=") && (MHillas.fSize<";
292 ifun >> ch;
293 sizecut+=ch;
294 sizecut+=")";
295 }
296 if(strcmp(word.Data(),"DISTCUT")==0)
297 {
298 if(distcut.Length())
299 cout << "readDataCards Warning: overriding existing DIST cut" << endl;
300 char ch[10];
301 distcut="(MHillasSrc.fDist>";
302 ifun >> ch;
303 distcut+=ch;
304 distcut+="*";
305 distcut+=conver;
306 distcut+=") && (MHillasSrc.fDist<";
307 ifun >> ch;
308 distcut+=ch;
309 distcut+="*";
310 distcut+=conver;
311 distcut+=")";
312 }
313 if(strcmp(word.Data(),"WIDTHCUT")==0)
314 {
315 if(widthcut.Length())
316 cout << "readDataCards Warning: overriding existing width cut" << endl;
317 char ch[10];
318 widthcut="(MHillas.fWidth>";
319 ifun >> ch;
320 widthcut+=ch;
321 widthcut+="*";
322 widthcut+=conver;
323 widthcut+=") && (MHillas.fWidth<";
324 ifun >> ch;
325 widthcut+=ch;
326 widthcut+="*";
327 widthcut+=conver;
328 widthcut+=")";
329 }
330 if(strcmp(word.Data(),"LENGTHCUT")==0)
331 {
332 if(lengthcut.Length())
333 cout << "readDataCards Warning: overriding existing length cut" << endl;
334 char ch[10];
335 lengthcut="(MHillas.fLength>";
336 ifun >> ch;
337 lengthcut+=ch;
338 lengthcut+="*";
339 lengthcut+=conver;
340 lengthcut+=") && (MHillas.fLength<";
341 ifun >> ch;
342 lengthcut+=ch;
343 lengthcut+="*";
344 lengthcut+=conver;
345 lengthcut+=")";
346 }
347 if(strcmp(word.Data(),"CENTERCUT")==0)
348 {
349 if(centercut.Length())
350 cout << "readDataCards Warning: overriding existing center cut" << endl;
351 char ch[10];
352 centercut="sqrt(MHillas.fMeanX*MHillas.fMeanX+MHillas.fMeanY*MHillas.fMeanY)<";
353 ifun >> ch;
354 centercut+=ch;
355 centercut+="*";
356 centercut+=conver;
357 }
358 }
359
360 // check compulsory values
361 if(!inputFile.Length())
362 {
363 cout << "No on-data file name specified" << endl;
364 return kFALSE;
365 }
366
367 if(!offFile.Length())
368 {
369 cout << "No off-data file name specified" << endl;
370 return kFALSE;
371 }
372
373 if(!outputFile.Length())
374 {
375 cout << "No output file name specified" << endl;
376 return kFALSE;
377 }
378
379 MTime thetime(mjdpos);
380
381 // Dump read values
382 cout << "************************************************" << endl;
383 cout << "* Datacards read from file " << filename << endl;
384 cout << "************************************************" << endl;
385 cout << "Maximum number of input events: " << nmaxevents << endl;
386 cout << "Input file name(s): " << inputFile << endl;
387 cout << "Input file name(s) for off data: " << offFile << endl;
388 cout << "Output file name for # hits plot : " << outputFile << endl;
389 if(filter.Length())
390 cout << "Applying rejection cut: " << filter << endl;
391 cout << "***********" << endl << endl;
392
393 return kTRUE;
394}
Note: See TracBrowser for help on using the repository browser.