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

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