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

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