source: trunk/MagicSoft/Mars/mtemp/mifae/programs/srcPos.cc@ 5868

Last change on this file since 5868 was 5191, checked in by rico, 21 years ago
*** empty log message ***
File size: 9.9 KB
Line 
1////////////////////////////////////////////////////////////////////////////////////
2//
3// _____ Source Position macro_____
4//
5// Take as input root files with hillas parameters and recompute the ones depending
6// on the source position for a new (input) position, optionally rotating it in an
7// event by event basis. Output is a file with recomputed hillas parameters
8//
9// Ester Aliu <aliu@ifae.es>
10// Oscar Blanch <blanch@ifae.es>
11// Javier Rico <jrico@ifae.es>
12////////////////////////////////////////////////////////////////////////////////////
13
14#include <fstream>
15#include <iostream>
16
17#include "TString.h"
18
19#include "MHillasSrcCalc.h"
20#include "MHillasSrc.h"
21#include "MSrcRotate.h"
22#include "MSrcPosFromFile.h"
23#include "MSrcTranslate.h"
24#include "MParList.h"
25#include "MTaskList.h"
26#include "MHillas.h"
27#include "MReadTree.h"
28#include "MEvtLoop.h"
29#include "MLog.h"
30#include "MArgs.h"
31#include "MWriteRootFile.h"
32#include "MTime.h"
33#include "MControlPlots.h"
34#include "MIslands.h"
35#include "MArgs.h"
36
37using namespace std;
38
39Bool_t readDatacards(TString& filename);
40void srcPos();
41
42//-----------------------------------------------------------------------------
43// declaration of variables read from datacards
44//-----------------------------------------------------------------------------
45
46TString inputFile;
47TString offFile;
48TString outputFile;
49TString offOutputFile;
50ULong_t nmaxevents=999999999;
51Float_t xsrcpos=0.;
52Float_t ysrcpos=0.;
53Double_t mjdpos=0;
54Bool_t kRotating=0;
55Bool_t kSrcPolicy=kFALSE;
56Double_t fRA= -1.;
57Double_t fDEC= -1.;
58TString srcFile;
59TString controlplotsfilename="controlplots.ps";
60
61//-----------------------------------------------------------------------------
62// constants
63//-----------------------------------------------------------------------------
64
65const TString defaultcard="srcpos.datacard";
66const Float_t conver = 189./0.6; // conversion factor degrees to mm
67
68//-----------------------------------------------------------------------------
69
70static void Usage()
71{
72 gLog <<endl;
73 gLog << "Usage is:" << endl;
74 gLog << " srcPos [-h] [-?] <datacards>" << endl << endl;
75 gLog << " <datacards>: datacards file name (dafault " << defaultcard <<")" << endl;
76 gLog << " -?/-h: This help" << endl << endl;
77}
78
79//-----------------------------------------------------------------------------
80int main(int argc, char **argv)
81{
82 // evaluate arguments
83 MArgs arg(argc, argv);
84 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
85 {
86 Usage();
87 return -1;
88 }
89
90 TString datacard = arg.GetArgumentStr(0);
91 if(!datacard.Length())
92 datacard = defaultcard;
93
94 if(!readDatacards(datacard))
95 {
96 cout << "Error reading datacards. Stoping" << endl;
97 return -1;
98 }
99
100 srcPos();
101}
102
103//-----------------------------------------------------------------------------
104void srcPos()
105{
106 // variable declaration
107 Float_t xpos=xsrcpos*conver; // [mm]
108 Float_t ypos=ysrcpos*conver; // [mm]
109
110 // containers
111 MParList plist;
112 MTaskList tlist;
113 MIslands islands;
114
115 // include containers in parameter list
116 plist.AddToList(&tlist);
117 plist.AddToList(&islands);
118
119 // tasks
120 MReadTree read("Parameters", inputFile);
121 read.DisableAutoScheme();
122
123 MSrcTranslate srctranslate;
124 MSrcRotate srcrotate;
125 MSrcPosFromFile srcfromfile;
126
127 if(srcFile.Length())
128 {
129 srcfromfile.SetInputFileName(srcFile);
130 srcfromfile.SetMode(MSrcPlace::kOn);
131 }
132 else
133 {
134 srctranslate.SetTranslation(xpos,ypos);
135 srctranslate.SetRelativeTranslation(kSrcPolicy);
136 srctranslate.SetCreateHisto(kFALSE);
137 srcrotate.SetRAandDECandRefMJD(fRA,fDEC,mjdpos);
138 srcrotate.SetMode(MSrcPlace::kOn);
139 }
140
141 MHillasSrcCalc csrc1;
142
143 MWriteRootFile write(outputFile,"RECREATE");
144 write.AddContainer("MHillas" , "Parameters");
145 write.AddContainer("MHillasSrc" , "Parameters");
146 write.AddContainer("MHillasExt" , "Parameters");
147 write.AddContainer("MNewImagePar" , "Parameters");
148 write.AddContainer("MRawEvtHeader" , "Parameters");
149 write.AddContainer("MRawRunHeader" , "Parameters");
150 write.AddContainer("MTime" , "Parameters");
151 write.AddContainer("MConcentration" , "Parameters");
152 write.AddContainer("MSrcPosCam" , "Parameters");
153 write.AddContainer("MIslands" , "Parameters");
154 //write.AddContainer("MIslands2" , "Parameters");
155
156 MControlPlots controlplots(controlplotsfilename);
157 controlplots.SetProduceFile(kFALSE);
158
159 // include tasks in task list
160 tlist.AddToList(&read);
161 if(srcFile.Length())
162 tlist.AddToList(&srcfromfile);
163 else
164 {
165 tlist.AddToList(&srctranslate);
166 if(kRotating)
167 tlist.AddToList(&srcrotate);
168 }
169 tlist.AddToList(&csrc1);
170 tlist.AddToList(&controlplots);
171 tlist.AddToList(&write);
172
173 // Eventloop
174 MEvtLoop evtloop;
175 evtloop.SetParList(&plist);
176 if (!evtloop.Eventloop(nmaxevents))
177 return;
178
179 tlist.PrintStatistics();
180
181 // do off-data if input file was specified
182 if(!offFile.Length())
183 return;
184
185 MReadTree read2("Parameters", offFile);
186 read2.DisableAutoScheme();
187 tlist.AddToListBefore(&read2, &read, "All");
188 tlist.RemoveFromList(&read);
189
190 MWriteRootFile write2(offOutputFile,"RECREATE");
191 write2.AddContainer("MHillas" , "Parameters");
192 write2.AddContainer("MHillasSrc" , "Parameters");
193 write2.AddContainer("MHillasExt" , "Parameters");
194 write2.AddContainer("MNewImagePar" , "Parameters");
195 write2.AddContainer("MRawEvtHeader" , "Parameters");
196 write2.AddContainer("MRawRunHeader" , "Parameters");
197 write2.AddContainer("MTime" , "Parameters");
198 write2.AddContainer("MConcentration" , "Parameters");
199 write2.AddContainer("MSrcPosCam" , "Parameters");
200 write2.AddContainer("MIslands" , "Parameters");
201 // write2.AddContainer("MIslands2" , "Parameters");
202 tlist.AddToListBefore(&write2,&write,"All");
203 tlist.RemoveFromList(&write);
204
205 if(srcFile.Length())
206 srcfromfile.SetMode(MSrcPlace::kOff);
207 else
208 srcrotate.SetMode(MSrcPlace::kOff);
209
210 controlplots.SetMode(MControlPlots::kOff);
211 controlplots.SetProduceFile(kTRUE);
212
213 if (!evtloop.Eventloop(nmaxevents))
214 return;
215
216 tlist.PrintStatistics();
217}
218//-----------------------------------------------------------------------
219
220Bool_t readDatacards(TString& filename)
221{
222 ifstream ifun(filename.Data());
223 if(!ifun)
224 {
225 cout << "File " << filename << " not found" << endl;
226 return kFALSE;
227 }
228
229 TString word;
230
231 while(ifun >> word)
232 {
233 // skip comments
234 if(word[0]=='/' && word[1]=='/')
235 {
236 while(ifun.get()!='\n'); // skip line
237 continue;
238 }
239
240 // number of events
241 if(strcmp(word.Data(),"NEVENTS")==0)
242 ifun >> nmaxevents;
243
244 // input file name
245 if(strcmp(word.Data(),"INPUTFILES")==0)
246 {
247 if(inputFile.Length())
248 cout << "readDataCards Warning: overriding on-data file name" << endl;
249 ifun >> inputFile;
250 }
251
252 // off-data file name
253 if(strcmp(word.Data(),"OFFFILES")==0)
254 {
255 if(offFile.Length())
256 cout << "readDataCards Warning: overriding off-data file name" << endl;
257 ifun >> offFile;
258 }
259
260 // output file name
261 if(strcmp(word.Data(),"OUTFILE")==0)
262 {
263 if(outputFile.Length())
264 cout << "readDataCards Warning: overriding output file name" << endl;
265 ifun >> outputFile;
266 }
267
268 // output file name (off data)
269 if(strcmp(word.Data(),"OFFOUTFILE")==0)
270 {
271 if(offOutputFile.Length())
272 cout << "readDataCards Warning: overriding output file name for off data" << endl;
273 ifun >> offOutputFile;
274 }
275
276 // source position input file
277 if(strcmp(word.Data(),"SRCFILE")==0)
278 {
279 if(srcFile.Length())
280 cout << "readDataCards Warning: overriding source-position file name" << endl;
281 ifun >> srcFile;
282 }
283
284 // source position
285 if(strcmp(word.Data(),"SRCPOS")==0)
286 {
287 ifun >> xsrcpos;
288 ifun >> ysrcpos;
289 ifun >> mjdpos;
290 }
291
292 // source celestial coordinates
293 if(strcmp(word.Data(),"SRCCOORDS")==0)
294 {
295 ifun >> fRA;
296 ifun >> fDEC;
297 }
298
299 // source movement policy flag
300 if(strcmp(word.Data(),"SRCABS")==0)
301 ifun >> kSrcPolicy;
302
303 // rotation flag
304 if(strcmp(word.Data(),"ROTFLAG")==0)
305 ifun >> kRotating;
306 }
307
308 // check compulsory values
309 if(!inputFile.Length())
310 {
311 cout << "No on-data file name specified" << endl;
312 return kFALSE;
313 }
314
315 if(!outputFile.Length())
316 {
317 cout << "No output file name specified" << endl;
318 return kFALSE;
319 }
320
321 if(offFile.Length() && !offOutputFile.Length())
322 {
323 cout << "No output file name specified for off data" << endl;
324 return kFALSE;
325 }
326
327 if(xsrcpos==0 && ysrcpos==0)
328 {
329 cout << "Source position is center of the camera (as in input file)" << endl;
330 return kFALSE;
331 }
332
333 MTime thetime(mjdpos);
334
335 // Dump read values
336 cout << "************************************************" << endl;
337 cout << "* Datacards read from file " << filename << endl;
338 cout << "************************************************" << endl;
339 cout << "Maximum number of input events: " << nmaxevents << endl;
340 cout << "Input file name(s): " << inputFile << endl;
341 if(offFile.Length())
342 cout << "Input file name(s) for off data: " << offFile << endl;
343 cout << "Output file name: " << outputFile << endl;
344 if(offFile.Length())
345 cout << "Output file name for off data: " << offOutputFile << endl;
346 if(srcFile.Length())
347 cout << "Reading source position from file " << srcFile << endl;
348 else
349 {
350 cout << "Source position (degrees) X=" << xsrcpos << ", Y="<<ysrcpos;
351 if(kSrcPolicy)
352 cout << " (RELATIVE TO INITIAL SOURCE POSITION)";
353 else
354 cout << " (ABSOLUTE POSITION IN THE CAMERA)";
355 cout << ", at " << thetime.GetSqlDateTime() << endl;
356 cout << "De-rotation flag " << kRotating << endl;
357 cout << "Source celestial coordiantes (rad): RA = " << fRA << ", DEC = " << fDEC << endl;
358 }
359 cout << "***********" << endl << endl;
360
361 return kTRUE;
362}
Note: See TracBrowser for help on using the repository browser.