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

Last change on this file since 4800 was 4395, checked in by blanch, 20 years ago
*** empty log message ***
File size: 8.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="srcpos.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 plist.AddToList(&islands);
115
116 MF cut(filter);
117 MF cutS(sizecut);
118 MF cutD(distcut);
119 MF cutL(lengthcut);
120 MF cutW(widthcut);
121 MF cutC(centercut);
122
123 // tasks
124 MReadTree read("Parameters", inputFile);
125 read.DisableAutoScheme();
126
127 MControlPlots controlplots(outputFile);
128 controlplots.SetProduceFile(kFALSE);
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(centercut.Length())
160 tlist.AddToList(&applycutC);
161 if(filter.Length())
162 tlist.AddToList(&applycut);
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 MReadTree read2("Parameters", offFile);
178 read2.DisableAutoScheme();
179 tlist.AddToListBefore(&read2, &read, "All");
180 tlist.RemoveFromList(&read);
181
182 controlplots.SetMode(MControlPlots::kOff);
183 controlplots.SetProduceFile(kTRUE);
184
185 if (!evtloop.Eventloop(nmaxevents))
186 return;
187
188 tlist.PrintStatistics();
189}
190//-----------------------------------------------------------------------
191
192Bool_t readDatacards(TString& filename)
193{
194 ifstream ifun(filename.Data());
195 if(!ifun)
196 {
197 cout << "File " << filename << " not found" << endl;
198 return kFALSE;
199 }
200
201 TString word;
202
203 while(ifun >> word)
204 {
205 // skip comments
206 if(word[0]=='/' && word[1]=='/')
207 {
208 while(ifun.get()!='\n'); // skip line
209 continue;
210 }
211
212 // number of events
213 if(strcmp(word.Data(),"NEVENTS")==0)
214 ifun >> nmaxevents;
215
216 // input file name
217 if(strcmp(word.Data(),"ONFILES")==0)
218 {
219 if(inputFile.Length())
220 cout << "readDataCards Warning: overriding on-data file name" << endl;
221 ifun >> inputFile;
222 }
223
224 // off-data file name
225 if(strcmp(word.Data(),"OFFFILES")==0)
226 {
227 if(offFile.Length())
228 cout << "readDataCards Warning: overriding off-data file name" << endl;
229 ifun >> offFile;
230 }
231
232 // output file name
233 if(strcmp(word.Data(),"HITSFILE")==0)
234 {
235 if(outputFile.Length())
236 cout << "readDataCards Warning: overriding output file name" << endl;
237 ifun >> outputFile;
238 }
239 // exclusion cut
240 if(strcmp(word.Data(),"FILTER")==0)
241 {
242 if(filter.Length())
243 cout << "readDataCards Warning: overriding existing cut" << endl;
244
245 char ch;
246 while((ch=ifun.get())!='\n')
247 filter.Append(ch);
248 }
249 if(strcmp(word.Data(),"SIZECUT")==0)
250 {
251 if(sizecut.Length())
252 cout << "readDataCards Warning: overriding existing size cut" << endl;
253 char ch[10];
254 sizecut="(MHillas.fSize>";
255 ifun >> ch;
256 sizecut+=ch;
257 sizecut+=") && (MHillas.fSize<";
258 ifun >> ch;
259 sizecut+=ch;
260 sizecut+=")";
261 }
262 if(strcmp(word.Data(),"DISTCUT")==0)
263 {
264 if(distcut.Length())
265 cout << "readDataCards Warning: overriding existing DIST cut" << endl;
266 char ch[10];
267 distcut="(MHillasSrc.fDist>";
268 ifun >> ch;
269 distcut+=ch;
270 distcut+=") && (MHillasSrc.fDist<";
271 ifun >> ch;
272 distcut+=ch;
273 distcut+=")";
274 }
275 if(strcmp(word.Data(),"WIDTHCUT")==0)
276 {
277 if(widthcut.Length())
278 cout << "readDataCards Warning: overriding existing width cut" << endl;
279 char ch[10];
280 widthcut="(MHillas.fWidth>";
281 ifun >> ch;
282 widthcut+=ch;
283 widthcut+=") && (MHillas.fWidth<";
284 ifun >> ch;
285 widthcut+=ch;
286 widthcut+=")";
287 }
288 if(strcmp(word.Data(),"LENGTHCUT")==0)
289 {
290 if(lengthcut.Length())
291 cout << "readDataCards Warning: overriding existing length cut" << endl;
292 char ch[10];
293 lengthcut="(MHillas.fLength>";
294 ifun >> ch;
295 lengthcut+=ch;
296 lengthcut+=") && (MHillas.fLength<";
297 ifun >> ch;
298 lengthcut+=ch;
299 lengthcut+=")";
300 }
301 if(strcmp(word.Data(),"CENTERCUT")==0)
302 {
303 if(centercut.Length())
304 cout << "readDataCards Warning: overriding existing center cut" << endl;
305 char ch[10];
306 centercut="sqrt(MHillas.fMeanX*MHillas.fMeanX+MHillas.fMeanY*MHillas.fMeanY) < 0.00317*";
307 ifun >> ch;
308 centercut+=ch;
309 }
310 }
311
312 // check compulsory values
313 if(!inputFile.Length())
314 {
315 cout << "No on-data file name specified" << endl;
316 return kFALSE;
317 }
318
319 if(!offFile.Length())
320 {
321 cout << "No off-data file name specified" << endl;
322 return kFALSE;
323 }
324
325 if(!outputFile.Length())
326 {
327 cout << "No output file name specified" << endl;
328 return kFALSE;
329 }
330
331 MTime thetime(mjdpos);
332
333 // Dump read values
334 cout << "************************************************" << endl;
335 cout << "* Datacards read from file " << filename << endl;
336 cout << "************************************************" << endl;
337 cout << "Maximum number of input events: " << nmaxevents << endl;
338 cout << "Input file name(s): " << inputFile << endl;
339 cout << "Input file name(s) for off data: " << offFile << endl;
340 cout << "Output file name for # hits plot : " << outputFile << endl;
341 if(filter.Length())
342 cout << "Applying rejection cut: " << filter << endl;
343 cout << "***********" << endl << endl;
344
345 return kTRUE;
346}
Note: See TracBrowser for help on using the repository browser.