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

Last change on this file since 4450 was 4424, checked in by rico, 20 years ago
*** empty log message ***
File size: 9.8 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("MTime" , "Parameters");
150 write.AddContainer("MConcentration" , "Parameters");
151 write.AddContainer("MSrcPosCam" , "Parameters");
152 write.AddContainer("MIslands" , "Parameters");
153 //write.AddContainer("MIslands2" , "Parameters");
154
155 MControlPlots controlplots(controlplotsfilename);
156 controlplots.SetProduceFile(kFALSE);
157
158 // include tasks in task list
159 tlist.AddToList(&read);
160 if(srcFile.Length())
161 tlist.AddToList(&srcfromfile);
162 else
163 {
164 tlist.AddToList(&srctranslate);
165 if(kRotate)
166 tlist.AddToList(&srcrotate);
167 }
168 tlist.AddToList(&csrc1);
169 tlist.AddToList(&controlplots);
170 tlist.AddToList(&write);
171
172 // Eventloop
173 MEvtLoop evtloop;
174 evtloop.SetParList(&plist);
175 if (!evtloop.Eventloop(nmaxevents))
176 return;
177
178 tlist.PrintStatistics();
179
180 // do off-data if input file was specified
181 if(!offFile.Length())
182 return;
183
184 MReadTree read2("Parameters", offFile);
185 read2.DisableAutoScheme();
186 tlist.AddToListBefore(&read2, &read, "All");
187 tlist.RemoveFromList(&read);
188
189 MWriteRootFile write2(offOutputFile,"RECREATE");
190 write2.AddContainer("MHillas" , "Parameters");
191 write2.AddContainer("MHillasSrc" , "Parameters");
192 write2.AddContainer("MHillasExt" , "Parameters");
193 write2.AddContainer("MNewImagePar" , "Parameters");
194 write2.AddContainer("MRawEvtHeader" , "Parameters");
195 write2.AddContainer("MRawRunHeader" , "Parameters");
196 write2.AddContainer("MTime" , "Parameters");
197 write2.AddContainer("MConcentration" , "Parameters");
198 write2.AddContainer("MSrcPosCam" , "Parameters");
199 write2.AddContainer("MIslands" , "Parameters");
200 // write2.AddContainer("MIslands2" , "Parameters");
201 tlist.AddToListBefore(&write2,&write,"All");
202 tlist.RemoveFromList(&write);
203
204 if(srcFile.Length())
205 srcfromfile.SetMode(MSrcPlace::kOff);
206 else
207 srcrotate.SetMode(MSrcPlace::kOff);
208
209 controlplots.SetMode(MControlPlots::kOff);
210 controlplots.SetProduceFile(kTRUE);
211
212 if (!evtloop.Eventloop(nmaxevents))
213 return;
214
215 tlist.PrintStatistics();
216}
217//-----------------------------------------------------------------------
218
219Bool_t readDatacards(TString& filename)
220{
221 ifstream ifun(filename.Data());
222 if(!ifun)
223 {
224 cout << "File " << filename << " not found" << endl;
225 return kFALSE;
226 }
227
228 TString word;
229
230 while(ifun >> word)
231 {
232 // skip comments
233 if(word[0]=='/' && word[1]=='/')
234 {
235 while(ifun.get()!='\n'); // skip line
236 continue;
237 }
238
239 // number of events
240 if(strcmp(word.Data(),"NEVENTS")==0)
241 ifun >> nmaxevents;
242
243 // input file name
244 if(strcmp(word.Data(),"INPUTFILES")==0)
245 {
246 if(inputFile.Length())
247 cout << "readDataCards Warning: overriding on-data file name" << endl;
248 ifun >> inputFile;
249 }
250
251 // off-data file name
252 if(strcmp(word.Data(),"OFFFILES")==0)
253 {
254 if(offFile.Length())
255 cout << "readDataCards Warning: overriding off-data file name" << endl;
256 ifun >> offFile;
257 }
258
259 // output file name
260 if(strcmp(word.Data(),"OUTFILE")==0)
261 {
262 if(outputFile.Length())
263 cout << "readDataCards Warning: overriding output file name" << endl;
264 ifun >> outputFile;
265 }
266
267 // output file name (off data)
268 if(strcmp(word.Data(),"OFFOUTFILE")==0)
269 {
270 if(offOutputFile.Length())
271 cout << "readDataCards Warning: overriding output file name for off data" << endl;
272 ifun >> offOutputFile;
273 }
274
275 // source position input file
276 if(strcmp(word.Data(),"SRCFILE")==0)
277 {
278 if(srcFile.Length())
279 cout << "readDataCards Warning: overriding source-position file name" << endl;
280 ifun >> srcFile;
281 }
282
283 // source position
284 if(strcmp(word.Data(),"SRCPOS")==0)
285 {
286 ifun >> xsrcpos;
287 ifun >> ysrcpos;
288 ifun >> mjdpos;
289 }
290
291 // source celestial coordinates
292 if(strcmp(word.Data(),"SRCCOORDS")==0)
293 {
294 ifun >> fRA;
295 ifun >> fDEC;
296 }
297
298 // source movement policy flag
299 if(strcmp(word.Data(),"SRCABS")==0)
300 ifun >> kSrcPolicy;
301
302 // rotation flag
303 if(strcmp(word.Data(),"ROTFLAG")==0)
304 ifun >> kRotate;
305 }
306
307 // check compulsory values
308 if(!inputFile.Length())
309 {
310 cout << "No on-data file name specified" << endl;
311 return kFALSE;
312 }
313
314 if(!outputFile.Length())
315 {
316 cout << "No output file name specified" << endl;
317 return kFALSE;
318 }
319
320 if(offFile.Length() && !offOutputFile.Length())
321 {
322 cout << "No output file name specified for off data" << endl;
323 return kFALSE;
324 }
325
326 if(xsrcpos==0 && ysrcpos==0)
327 {
328 cout << "Source position is center of the camera (as in input file)" << endl;
329 return kFALSE;
330 }
331
332 MTime thetime(mjdpos);
333
334 // Dump read values
335 cout << "************************************************" << endl;
336 cout << "* Datacards read from file " << filename << endl;
337 cout << "************************************************" << endl;
338 cout << "Maximum number of input events: " << nmaxevents << endl;
339 cout << "Input file name(s): " << inputFile << endl;
340 if(offFile.Length())
341 cout << "Input file name(s) for off data: " << offFile << endl;
342 cout << "Output file name: " << outputFile << endl;
343 if(offFile.Length())
344 cout << "Output file name for off data: " << offOutputFile << endl;
345 if(srcFile.Length())
346 cout << "Reading source position from file " << srcFile << endl;
347 else
348 {
349 cout << "Source position (degrees) X=" << xsrcpos << ", Y="<<ysrcpos;
350 if(kSrcPolicy)
351 cout << " (RELATIVE TO INITIAL SOURCE POSITION)";
352 else
353 cout << " (ABSOLUTE POSITION IN THE CAMERA)";
354 cout << ", at " << thetime.GetSqlDateTime() << endl;
355 cout << "De-rotation flag " << kRotate << endl;
356 cout << "Source celestial coordiantes (rad): RA = " << fRA << ", DEC = " << fDEC << endl;
357 }
358 cout << "***********" << endl << endl;
359
360 return kTRUE;
361}
Note: See TracBrowser for help on using the repository browser.