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

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