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

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