source: trunk/MagicSoft/Mars/mtemp/mifae/programs/makeHillas.cc@ 4019

Last change on this file since 4019 was 3965, checked in by rico, 21 years ago
*** empty log message ***
File size: 10.7 KB
Line 
1/*********************************/
2/* Compute the hillas parameters */
3/*********************************/
4
5#include "TString.h"
6#include "TArrayS.h"
7
8#include "MParList.h"
9#include "MTaskList.h"
10#include "MPedestalCam.h"
11#include "MBadPixelsCam.h"
12#include "MReadMarsFile.h"
13#include "MGeomApply.h"
14#include "MPedCalcPedRun.h"
15#include "MEvtLoop.h"
16#include "MGeomCamMagic.h"
17#include "MExtractedSignalCam.h"
18#include "MCalibrationChargeCam.h"
19#include "MHCalibrationChargeCam.h"
20#include "MHCalibrationRelTimeCam.h"
21#include "MExtractor.h"
22#include "MExtractFixedWindow.h"
23#include "MExtractSignal.h"
24#include "MCalibrationChargeCalc.h"
25#include "MFCosmics.h"
26#include "MContinue.h"
27#include "MFillH.h"
28#include "MLog.h"
29#include "MCerPhotEvt.h"
30#include "MPedPhotCam.h"
31#include "MCalibrate.h"
32#include "MPedPhotCalc.h"
33#include "MHillas.h"
34#include "MRawRunHeader.h"
35#include "MSrcPosCam.h"
36#include "MBlindPixelCalc.h"
37#include "MImgCleanStd.h"
38#include "MHillasSrcCalc.h"
39#include "MHillasCalc.h"
40#include "MWriteRootFile.h"
41#include "MProgressBar.h"
42#include "MArgs.h"
43#include "MRunIter.h"
44#include "MJPedestal.h"
45#include "MJCalibration.h"
46
47#include <iostream>
48#include <fstream>
49#include <stdlib.h>
50
51using namespace std;
52
53Bool_t readDatacards(TString& filename);
54void makeHillas();
55
56// initial and final time slices to be used in signal extraction
57const Byte_t hifirst = 0;
58const Byte_t hilast = 13;
59const Byte_t lofirst = 3;
60const Byte_t lolast = 12;
61
62// declaration of variables read from datacards
63TString outname;
64TString idirname;
65MRunIter caliter;
66MRunIter pediter;
67MRunIter datiter;
68ULong_t nmaxevents=999999999;
69Short_t calflag=1;
70Float_t lcore = 3.0;
71Float_t ltail = 1.5;
72Int_t nfiles = 0;
73
74const TString defaultcard="input.datacard";
75/*************************************************************/
76static void Usage()
77{
78 gLog <<endl;
79 gLog << "Usage is:" << endl;
80 gLog << " makeHillas [-h] [-?] <datacards>" << endl << endl;
81 gLog << " <datacards>: datacards file name (dafault input.datacards)" << endl;
82 gLog << " -?/-h: This help" << endl << endl;
83}
84/*************************************************************/
85int main(int argc, char **argv)
86{
87 // evaluate arguments
88 MArgs arg(argc, argv);
89 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
90 {
91 Usage();
92 return -1;
93 }
94
95 TString datacard = arg.GetArgumentStr(0);
96 if(!datacard.Length())
97 datacard = defaultcard;
98
99 if(!readDatacards(datacard))
100 {
101 cout << "Error reading datacards. Stoping" << endl;
102 return -1;
103 }
104 makeHillas();
105}
106
107/*************************************************************/
108void makeHillas()
109{
110 // Set the general tasks/containers
111 MExtractFixedWindow extractor;
112 extractor.SetRange(hifirst,hilast,lofirst,lolast);
113
114 MCalibrationQECam qecam;
115 MBadPixelsCam badcam;
116 MGeomCamMagic geomcam;
117 MGeomApply geomapl;
118
119 /************************************/
120 /* FIRST LOOP: PEDESTAL COMPUTATION */
121 /************************************/
122
123 // If you want to exclude pixels from the beginning, read
124 // an ascii-file with the corr. pixel numbers (see MBadPixelsCam)
125 //badcam.AsciiRead("badpixels.dat");
126
127 MJPedestal pedloop;
128 pedloop.SetInput(&pediter);
129 pedloop.SetExtractor(&extractor);
130 pedloop.SetBadPixels(badcam);
131
132 if (!pedloop.Process())
133 return;
134
135 MPedestalCam pedcam = pedloop.GetPedestalCam();
136
137 /*****************************/
138 /* SECOND LOOP: CALIBRATION */
139 /*****************************/
140
141 MJCalibration calloop;
142
143 calloop.SetExtractor(&extractor);
144 calloop.SetInput(&caliter);
145 calloop.SetQECam(qecam);
146 calloop.SetBadPixels(pedloop.GetBadPixels());
147 if (!calloop.Process(pedcam))
148 return;
149
150 MCalibrationChargeCam calcam = calloop.GetCalibrationCam();
151
152 /************************************************************************/
153 /* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */
154 /************************************************************************/
155 MParList plist3;
156 MTaskList tlist3;
157 plist3.AddToList(&tlist3);
158
159 // containers
160 MCerPhotEvt nphot;
161 MPedPhotCam nphotrms;
162 MExtractedSignalCam sigcam;
163
164 plist3.AddToList(&geomcam);
165 plist3.AddToList(&pedcam);
166 plist3.AddToList(&calcam);
167 plist3.AddToList(&badcam);
168 plist3.AddToList(&sigcam);
169 plist3.AddToList(&nphot);
170 plist3.AddToList(&nphotrms);
171
172
173 MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault;
174 if(calflag==0)
175 calMode=MCalibrate::kNone;
176 if(calflag==-1)
177 calMode=MCalibrate::kDummy;
178
179 //tasks
180 MReadMarsFile read3("Events");
181 static_cast<MRead&>(read3).AddFiles(pediter);
182 read3.DisableAutoScheme();
183
184 // MExtractSignal extsig;
185 // extsig.SetRange(hifirst,hilast,lofirst,lolast);
186 MCalibrate photcalc(calMode);
187 MPedPhotCalc photrmscalc;
188
189 tlist3.AddToList(&read3);
190 tlist3.AddToList(&geomapl);
191 tlist3.AddToList(&extractor);
192 tlist3.AddToList(&photcalc);
193 tlist3.AddToList(&photrmscalc);
194
195 // Create and setup the eventloop
196 MEvtLoop evtloop3;
197 evtloop3.SetParList(&plist3);
198 if (!evtloop3.Eventloop())
199 return;
200
201 tlist3.PrintStatistics();
202
203 /************************************************************************/
204 /* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */
205 /************************************************************************/
206
207 MParList plist4;
208 MTaskList tlist4;
209 plist4.AddToList(&tlist4);
210
211 // containers
212 MHillas hillas;
213 MSrcPosCam source;
214 MRawRunHeader runhead;
215
216 plist4.AddToList(&geomcam);
217 plist4.AddToList(&pedcam);
218 plist4.AddToList(&calcam);
219 // plist4.AddToList(&badcam);
220 plist4.AddToList(&nphot);
221 plist4.AddToList(&nphotrms);
222 plist4.AddToList(&source);
223 plist4.AddToList(&hillas);
224 plist4.AddToList(&runhead);
225
226 //tasks
227 MReadMarsFile read4("Events");
228 static_cast<MRead&>(read4).AddFiles(datiter);
229 read4.DisableAutoScheme();
230
231 // set bad pixels
232 MBlindPixelCalc blind;
233 MBlindPixelCalc blind2;
234 const Short_t x[16] = {330,395,329,396,389,
235 323,388,322,384,385,
236 386,387,321,320,319,
237 394};
238 const TArrayS bp(16,(Short_t*)x);
239 blind.SetPixelIndices(bp);
240 blind2.SetPixelIndices(bp);
241
242 MImgCleanStd clean(lcore,ltail);
243 MHillasCalc hcalc;
244 MHillasSrcCalc csrc1;
245
246 MWriteRootFile write(outname,"RECREATE");
247
248 write.AddContainer("MHillas" , "Parameters");
249 write.AddContainer("MHillasSrc" , "Parameters");
250 write.AddContainer("MHillasExt" , "Parameters");
251 write.AddContainer("MNewImagePar" , "Parameters");
252 write.AddContainer("MRawEvtHeader" , "Parameters");
253 write.AddContainer("MRawRunHeader" , "Parameters");
254 write.AddContainer("MConcentration" , "Parameters");
255 write.AddContainer("MSrcPosCam" , "Parameters");
256
257 tlist4.AddToList(&read4);
258 tlist4.AddToList(&geomapl);
259 tlist4.AddToList(&extractor);
260 tlist4.AddToList(&photcalc);
261 //tlist4.AddToList(&blind);
262 tlist4.AddToList(&clean);
263 //tlist4.AddToList(&blind2);
264 tlist4.AddToList(&hcalc);
265 // tlist4.AddToList(&srcposcalc);
266 tlist4.AddToList(&csrc1);
267 tlist4.AddToList(&write);
268
269 // Create and setup the eventloop
270 MEvtLoop datloop;
271 datloop.SetParList(&plist4);
272
273 cout << "*************************************************************" << endl;
274 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
275 cout << "*************************************************************" << endl;
276
277 if (!datloop.Eventloop(nmaxevents))
278 return;
279
280 tlist4.PrintStatistics();
281}
282//-------------------------------------------------------------------------------
283
284Bool_t readDatacards(TString& filename)
285{
286 ifstream ifun(filename.Data());
287 if(!ifun)
288 {
289 cout << "File " << filename << " not found" << endl;
290 return kFALSE;
291 }
292
293 TString word;
294
295 while(ifun >> word)
296 {
297 // skip comments
298 if(word[0]=='/' && word[1]=='/')
299 {
300 while(ifun.get()!='\n'); // skip line
301 continue;
302 }
303
304 // number of events
305 if(strcmp(word.Data(),"NEVENTS")==0)
306 ifun >> nmaxevents;
307
308
309 // input file directory
310 if(strcmp(word.Data(),"IDIR")==0)
311 {
312 if(idirname.Length())
313 cout << "readDataCards Warning: overriding input directory file name" << endl;
314 ifun >> idirname;
315 }
316
317 // pedestal runs
318 if(strcmp(word.Data(),"PRUNS")==0)
319 {
320 if(pediter.GetNumRuns())
321 cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl;
322 ifun >> word;
323 pediter.AddRuns(word.Data(),idirname.Data());
324 }
325
326 // calibration runs
327 if(strcmp(word.Data(),"CRUNS")==0)
328 {
329 if(caliter.GetNumRuns())
330 cout << "readDataCards Warning: adding calibration runs to the existing list" << endl;
331 ifun >> word;
332 caliter.AddRuns(word.Data(),idirname.Data());
333 }
334
335 // data runs
336 if(strcmp(word.Data(),"DRUNS")==0)
337 {
338 if(datiter.GetNumRuns())
339 cout << "readDataCards Warning: adding data runs to the existing list" << endl;
340 ifun >> word;
341 datiter.AddRuns(word.Data(),idirname.Data());
342 }
343
344 // output file name
345 if(strcmp(word.Data(),"OUTFILE")==0)
346 {
347 if(outname.Length())
348 cout << "readDataCards Warning: overriding output file name" << endl;
349 ifun >> outname;
350 }
351
352 // calibration flag
353 if(strcmp(word.Data(),"CALFLAG")==0)
354 ifun >> calflag;
355
356 // cleaning level
357 if(strcmp(word.Data(),"CLEANLEVEL")==0)
358 {
359 ifun >> lcore;
360 ifun >> ltail;
361 }
362 }
363
364 pediter.Reset();
365 caliter.Reset();
366 datiter.Reset();
367 TString pfile;
368
369 // Dump read values
370 cout << "************************************************" << endl;
371 cout << "* Datacards read from file " << filename << endl;
372 cout << "************************************************" << endl;
373 cout << "Pedestal file (s): " << endl;
374 while(!(pfile=pediter.Next()).IsNull())
375 cout << pfile << endl;
376 cout << "Calibration file (s): " << endl;
377 while(!(pfile=caliter.Next()).IsNull())
378 cout << pfile << endl;
379 cout << "Data file (s): " << endl;
380 while(!(pfile=datiter.Next()).IsNull())
381 cout << pfile << endl;
382 cout << "Maximum number of events: " << nmaxevents << endl;
383 cout << "Output file name: " << outname << endl;
384 cout << "Calibration flag: " << calflag << endl;
385 cout << "Cleaning level: ("<<lcore<<","<<ltail<<")" << endl;
386 cout << "***********" << endl << endl;
387
388 if(!pediter.GetNumEntries())
389 {
390 cout << "No pedestal file name specified" << endl;
391 return kFALSE;
392 }
393 if(!caliter.GetNumEntries())
394 {
395 cout << "No calibration file name specified" << endl;
396 return kFALSE;
397 }
398 if(!datiter.GetNumEntries())
399 {
400 cout << "No data file name specified" << endl;
401 return kFALSE;
402 }
403 if(!outname.Length())
404 {
405 cout << "No output file name specified" << endl;
406 return kFALSE;
407 }
408
409
410 return kTRUE;
411}
Note: See TracBrowser for help on using the repository browser.