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

Last change on this file since 5138 was 4561, checked in by rico, 20 years ago
*** empty log message ***
File size: 15.8 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 "MBadPixelsTreat.h"
13#include "MReadMarsFile.h"
14#include "MGeomApply.h"
15#include "MPedCalcPedRun.h"
16#include "MEvtLoop.h"
17#include "MGeomCamMagic.h"
18#include "MExtractedSignalCam.h"
19#include "MCalibrationChargeCam.h"
20#include "MHCalibrationChargeCam.h"
21#include "MHCalibrationRelTimeCam.h"
22#include "MExtractor.h"
23#include "MExtractFixedWindow.h"
24#include "MExtractSlidingWindow.h"
25#include "MExtractSignal.h"
26#include "MCalibrationChargeCalc.h"
27#include "MFCosmics.h"
28#include "MContinue.h"
29#include "MLog.h"
30#include "MCerPhotEvt.h"
31#include "MPedPhotCam.h"
32#include "MCalibrate.h"
33#include "MPedPhotCalc.h"
34#include "MHillas.h"
35#include "MNewImagePar.h"
36#include "MRawRunHeader.h"
37#include "MSrcPosCam.h"
38#include "MImgCleanStd.h"
39#include "MHillasSrcCalc.h"
40#include "MHillasCalc.h"
41#include "MArrivalTimeCam.h"
42#include "MArrivalTimeCalc2.h"
43#include "MIslands.h"
44#include "MIslandCalc.h"
45#include "MIslandClean.h"
46#include "MWriteRootFile.h"
47#include "MArgs.h"
48#include "MRunIter.h"
49#include "MJPedestal.h"
50#include "MJCalibration.h"
51#include "MHillasDisplay.h"
52#include "MF.h"
53#include "MContinue.h"
54
55#include "TApplication.h"
56
57#include <iostream>
58#include <fstream>
59#include <stdlib.h>
60
61using namespace std;
62
63Bool_t readDatacards(TString& filename);
64void makeHillas();
65
66// initial and final time slices to be used in signal extraction
67const Byte_t hifirst = 1;
68const Byte_t hilast = 14;
69const Byte_t lofirst = 3;
70const Byte_t lolast = 14;
71
72// declaration of variables read from datacards
73TString outname;
74const TString outpath = "./";
75TString idirname;
76TString filter;
77TString psfilename("makehillas.ps");
78MRunIter pedcaliter;
79MRunIter caliter;
80MRunIter pediter;
81MRunIter datiter;
82UInt_t display = 0;
83ULong_t nmaxevents= 999999999;
84Short_t calflag = 1;
85Bool_t caltimeflag = kFALSE;
86Short_t cleanflag = 1;
87UShort_t lrings = 1;
88Float_t lcore = 3.0;
89Float_t ltail = 1.5;
90Int_t islflag = 0;
91Float_t lnew = 40;
92Int_t kmethod = 1;
93Int_t kalgorithm = 1;
94Int_t nfiles = 0;
95
96const TString defaultcard="makehillas.datacard";
97/*************************************************************/
98static void Usage()
99{
100 gLog <<endl;
101 gLog << "Usage is:" << endl;
102 gLog << " makeHillas [-h] [-?] <datacards>" << endl << endl;
103 gLog << " <datacards>: datacards file name (dafault input.datacards)" << endl;
104 gLog << " -?/-h: This help" << endl << endl;
105}
106/*************************************************************/
107int main(int argc, char **argv)
108{
109 // create a TApplication to be able to
110 TApplication app("Application",0,0);
111
112 // evaluate arguments
113 MArgs arg(argc, argv);
114 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
115 {
116 Usage();
117 return -1;
118 }
119
120 TString datacard = arg.GetArgumentStr(0);
121 if(!datacard.Length())
122 datacard = defaultcard;
123
124 if(!readDatacards(datacard))
125 {
126 cout << "Error reading datacards. Stoping" << endl;
127 return -1;
128 }
129
130 makeHillas();
131}
132
133/*************************************************************/
134void makeHillas()
135{
136 // Set the general tasks/containers
137 MExtractFixedWindow extractor;
138 extractor.SetRange(hifirst,hilast,lofirst,lolast);
139
140 //MExtractSlidingWindow extractor;
141 //extractor.SetRange(hifirst,hilast,lofirst,lolast);
142 //extractor.SetWindowSize(4,4);
143
144 MGeomCamMagic geomcam;
145 MGeomApply geomapl;
146
147 /************************************/
148 /* FIRST LOOP: PEDESTAL COMPUTATION */
149 /************************************/
150
151 // If you want to exclude pixels from the beginning, read
152 // an ascii-file with the corr. pixel numbers (see MBadPixelsCam)
153 //badcam.AsciiRead("badpixels.dat");
154
155 MJPedestal pedcalloop;
156 pedcalloop.SetInput(&pediter);
157 pedcalloop.SetExtractor(&extractor);
158 // pedloop.SetBadPixels(badcam);
159 //pedcalloop.SetOutputPath(outpath.Data());
160
161
162 if (!pedcalloop.Process())
163 return;
164
165 /*****************************/
166 /* SECOND LOOP: CALIBRATION */
167 /*****************************/
168
169 MJCalibration calloop;
170
171 calloop.SetRelTimeCalibration(caltimeflag);
172 calloop.SetExtractor(&extractor);
173 calloop.SetInput(&caliter);
174 calloop.SetBadPixels(pedcalloop.GetBadPixels());
175 if(calflag==2)
176 calloop.SetUseBlindPixel();
177 if(calflag>0)
178 if (!calloop.Process(pedcalloop.GetPedestalCam()))
179 return;
180
181
182 /************************************************************************/
183 /* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */
184 /************************************************************************/
185 // First Compute the pedestals
186 MJPedestal pedloop;
187 pedloop.SetInput(&pediter);
188
189 if (!pedloop.Process())
190 return;
191
192 MParList plist3;
193 MTaskList tlist3;
194 plist3.AddToList(&tlist3);
195
196 // containers
197 MCerPhotEvt nphot;
198 MPedPhotCam nphotrms;
199 MExtractedSignalCam sigcam;
200
201 plist3.AddToList(&geomcam);
202 plist3.AddToList(&pedloop.GetPedestalCam());
203 plist3.AddToList(&calloop.GetCalibrationCam());
204 plist3.AddToList(&calloop.GetQECam());
205 plist3.AddToList(&calloop.GetRelTimeCam());
206 plist3.AddToList(&calloop.GetBadPixels());
207 plist3.AddToList(&sigcam);
208 plist3.AddToList(&nphot);
209 plist3.AddToList(&nphotrms);
210
211
212 MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault;
213 if(calflag==0)
214 calMode=MCalibrate::kNone;
215 if(calflag==-1)
216 calMode=MCalibrate::kDummy;
217
218 //tasks
219 MReadMarsFile read3("Events");
220 static_cast<MRead&>(read3).AddFiles(pediter);
221 read3.DisableAutoScheme();
222
223 MCalibrate photcalc(calMode);
224 MPedPhotCalc photrmscalc;
225
226 tlist3.AddToList(&read3);
227 tlist3.AddToList(&geomapl);
228 tlist3.AddToList(&extractor);
229 tlist3.AddToList(&photcalc);
230 tlist3.AddToList(&photrmscalc);
231
232 // Create and setup the eventloop
233 MEvtLoop evtloop3;
234 evtloop3.SetParList(&plist3);
235 if (!evtloop3.Eventloop())
236 return;
237
238 tlist3.PrintStatistics();
239
240 /************************************************************************/
241 /* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */
242 /************************************************************************/
243
244 MParList plist4;
245 MTaskList tlist4;
246 plist4.AddToList(&tlist4);
247
248 // containers
249 MHillas hillas;
250 MNewImagePar newimagepar;
251 MSrcPosCam source;
252 MRawRunHeader runhead;
253
254 MArrivalTimeCam timecam;
255
256
257 // islands
258 MIslands isl;
259 MIslands isl2;
260 MIslands isl3;
261
262 isl.SetName("MIslands");
263 isl2.SetName("MIslands2");
264 isl3.SetName("MIslands3");
265
266 plist4.AddToList(&timecam);
267 plist4.AddToList(&isl);
268
269 if (islflag == 2)
270 plist4.AddToList(&isl2);
271 if (islflag == 3)
272 plist4.AddToList(&isl3);
273
274 plist4.AddToList(&geomcam);
275 plist4.AddToList(&pedloop.GetPedestalCam());
276 plist4.AddToList(&calloop.GetCalibrationCam());
277 plist4.AddToList(&calloop.GetQECam());
278 plist4.AddToList(&calloop.GetRelTimeCam());
279 plist4.AddToList(&calloop.GetBadPixels());
280 plist4.AddToList(&nphot);
281 plist4.AddToList(&nphotrms);
282 plist4.AddToList(&source);
283 plist4.AddToList(&hillas);
284 plist4.AddToList(&newimagepar);
285 plist4.AddToList(&runhead);
286
287 // cuts
288 MF cut(filter);
289
290 //tasks
291 MReadMarsFile read4("Events");
292 static_cast<MRead&>(read4).AddFiles(datiter);
293 read4.DisableAutoScheme();
294
295 MImgCleanStd clean(lcore,ltail);
296 clean.SetCleanRings(lrings);
297 MImgCleanStd::CleaningMethod_t cleanMeth= MImgCleanStd::kStandard;
298 if(cleanflag==2)
299 cleanMeth=MImgCleanStd::kDemocratic;
300 clean.SetMethod(cleanMeth);
301
302 MArrivalTimeCalc2 timecalc;
303 MIslandCalc island;
304 island.SetOutputName("MIslands");
305 island.SetAlgorithm(kalgorithm);
306
307 MBadPixelsTreat interpolatebadpixels;
308 interpolatebadpixels.SetUseInterpolation();
309 interpolatebadpixels.SetProcessRMS();
310 // interpolatebadpixels.SetSloppyTreatment();
311
312 MIslandClean islclean(lnew);
313 islclean.SetInputName("MIslands");
314 islclean.SetMethod(kmethod);
315
316 MIslandCalc island2;
317 island2.SetOutputName("MIslands2");
318 island2.SetAlgorithm(kalgorithm);
319
320 MIslandCalc island3;
321 island3.SetOutputName("MIslands3");
322
323
324 MHillasCalc hcalc;
325 MHillasSrcCalc csrc1;
326
327 MContinue applycut(&cut);
328 applycut.SetInverted(kTRUE);
329 MWriteRootFile write(outname,"RECREATE");
330
331 MHillasDisplay* disphillas=NULL;
332
333 write.AddContainer("MHillas" , "Parameters");
334 write.AddContainer("MHillasSrc" , "Parameters");
335 write.AddContainer("MHillasExt" , "Parameters");
336 write.AddContainer("MNewImagePar" , "Parameters");
337 write.AddContainer("MRawEvtHeader" , "Parameters");
338 write.AddContainer("MRawRunHeader" , "Parameters");
339 write.AddContainer("MTime" , "Parameters");
340 write.AddContainer("MConcentration" , "Parameters");
341 write.AddContainer("MSrcPosCam" , "Parameters");
342 write.AddContainer("MIslands" , "Parameters");
343
344 if (islflag == 2)
345 write.AddContainer("MIslands2" , "Parameters");
346 if (islflag == 3)
347 write.AddContainer("MIslands3" , "Parameters");
348
349
350 if(display)
351 {
352 disphillas = new MHillasDisplay(&nphot,&geomcam);
353 disphillas->SetIslandsName("MIslands");
354 if (islflag == 2)
355 disphillas->SetIslandsName("MIslands2");
356 if (islflag == 3)
357 disphillas->SetIslandsName("MIslands3");
358 }
359
360 tlist4.AddToList(&read4);
361 tlist4.AddToList(&geomapl);
362 tlist4.AddToList(&extractor);
363 tlist4.AddToList(&photcalc);
364 if(calflag==11 || calflag==21)
365 tlist4.AddToList(&interpolatebadpixels);
366 tlist4.AddToList(&clean);
367 tlist4.AddToList(&timecalc);
368 tlist4.AddToList(&island);
369
370 if (islflag == 2)
371 {
372 tlist4.AddToList(&islclean);
373 tlist4.AddToList(&island2);
374 }
375
376 if (islflag == 3)
377 {
378 tlist4.AddToList(&islclean);
379 tlist4.AddToList(&island3);
380 }
381
382 tlist4.AddToList(&hcalc);
383 tlist4.AddToList(&csrc1);
384 if(filter.Length())
385 tlist4.AddToList(&applycut);
386 tlist4.AddToList(&write);
387 if(display)
388 {
389 disphillas->SetPSFile();
390 disphillas->SetPSFileName(psfilename);
391 if(display==2)
392 disphillas->SetPause(kFALSE);
393 tlist4.AddToList(disphillas);
394 }
395
396 // Create and setup the eventloop
397 MEvtLoop datloop;
398 datloop.SetParList(&plist4);
399
400 cout << "*************************************************************" << endl;
401 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
402 cout << "*************************************************************" << endl;
403
404 if (!datloop.Eventloop(nmaxevents))
405 return;
406
407 tlist4.PrintStatistics();
408
409}
410//-------------------------------------------------------------------------------
411
412Bool_t readDatacards(TString& filename)
413{
414 ifstream ifun(filename.Data());
415 if(!ifun)
416 {
417 cout << "File " << filename << " not found" << endl;
418 return kFALSE;
419 }
420
421 TString word;
422
423 while(ifun >> word)
424 {
425 // skip comments
426 if(word[0]=='/' && word[1]=='/')
427 {
428 while(ifun.get()!='\n'); // skip line
429 continue;
430 }
431
432 // number of events
433 if(strcmp(word.Data(),"NEVENTS")==0)
434 ifun >> nmaxevents;
435
436
437 // input file directory
438 if(strcmp(word.Data(),"IDIR")==0)
439 {
440 if(idirname.Length())
441 cout << "readDataCards Warning: overriding input directory file name" << endl;
442 ifun >> idirname;
443 }
444
445 // pedestal runs
446 if(strcmp(word.Data(),"PRUNS")==0)
447 {
448 if(pediter.GetNumRuns())
449 cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl;
450 ifun >> word;
451 pediter.AddRuns(word.Data(),idirname.Data());
452 }
453
454 // pedestal runs for calibration
455 if(strcmp(word.Data(),"PCRUNS")==0)
456 {
457 if(pedcaliter.GetNumRuns())
458 cout << "readDataCards Warning: adding pedestal runs for calibration to the existing list" << endl;
459 ifun >> word;
460 pedcaliter.AddRuns(word.Data(),idirname.Data());
461 }
462
463 // calibration runs
464 if(strcmp(word.Data(),"CRUNS")==0)
465 {
466 if(caliter.GetNumRuns())
467 cout << "readDataCards Warning: adding calibration runs to the existing list" << endl;
468 ifun >> word;
469 caliter.AddRuns(word.Data(),idirname.Data());
470 }
471
472 // data runs
473 if(strcmp(word.Data(),"DRUNS")==0)
474 {
475 if(datiter.GetNumRuns())
476 cout << "readDataCards Warning: adding data runs to the existing list" << endl;
477 ifun >> word;
478 datiter.AddRuns(word.Data(),idirname.Data());
479 }
480
481 // output file name
482 if(strcmp(word.Data(),"OUTFILE")==0)
483 {
484 if(outname.Length())
485 cout << "readDataCards Warning: overriding output file name" << endl;
486 ifun >> outname;
487 }
488
489 // exclusion cut
490 if(strcmp(word.Data(),"FILTER")==0)
491 {
492 if(filter.Length())
493 cout << "readDataCards Warning: overriding existing cut" << endl;
494
495 char ch;
496 while((ch=ifun.get())!='\n')
497 filter.Append(ch);
498 }
499
500 // display flag
501 if(strcmp(word.Data(),"DISPLAY")==0)
502 ifun >> display;
503
504 // ps file name
505 if(strcmp(word.Data(),"PSFILENAME")==0)
506 ifun >> psfilename;
507
508 // calibration flag
509 if(strcmp(word.Data(),"CALFLAG")==0)
510 ifun >> calflag;
511
512 // calibration flag
513 if(strcmp(word.Data(),"CALTIME")==0)
514 ifun >> caltimeflag;
515
516 // cleaning level
517 if(strcmp(word.Data(),"CLEANLEVEL")==0)
518 {
519 ifun >> lcore;
520 ifun >> ltail;
521 if(ifun.get()!='\n'){
522 ifun.unget();
523 ifun >> lrings;
524 if(ifun.get()!='\n'){
525 ifun.unget();
526 ifun >> cleanflag;
527 }
528 }
529 }
530
531 if(strcmp(word.Data(),"ISLFLAG")==0)
532 {
533 ifun >> islflag;
534
535 // if (islflag == 1 || islflag == 2)
536 ifun >> kalgorithm;
537 }
538
539 // island cleaning
540 if (islflag == 2){
541 if(strcmp(word.Data(),"ISLANDCLEAN")==0)
542 {
543 ifun >> kmethod;
544 ifun >> lnew;
545 }
546 }
547 }
548
549 pediter.Reset();
550 caliter.Reset();
551 datiter.Reset();
552 TString pfile;
553
554 // Dump read values
555 cout << "************************************************" << endl;
556 cout << "* Datacards read from file " << filename << endl;
557 cout << "************************************************" << endl;
558 cout << "Pedestal file (s): " << endl;
559 while(!(pfile=pediter.Next()).IsNull())
560 cout << pfile << endl;
561 cout << "Calibration file (s): " << endl;
562 while(!(pfile=caliter.Next()).IsNull())
563 cout << pfile << endl;
564 cout << "Data file (s): " << endl;
565 while(!(pfile=datiter.Next()).IsNull())
566 cout << pfile << endl;
567 cout << "Maximum number of events: " << nmaxevents << endl;
568 if(filter.Length())
569 cout << "Applying rejection cut: " << filter << endl;
570 cout << "Output file name: " << outname << endl;
571 if(display)
572 cout << "Generating PS file: " << psfilename << endl;
573 cout << "Calibration: ";
574 if(calflag==0)
575 cout << "Pixel area proportional intercalibration" << endl;
576 else if(calflag==-1)
577 cout << "No calibration whatsoever" << endl;
578 else if(calflag==1)
579 cout << "Default calibration" << endl;
580 else if(calflag==11)
581 cout << "Default calibration + bad pixels interpolation" << endl;
582 cout << "Cleaning level: ("<<lcore<<","<<ltail<<") - " << lrings << "ring" << endl;
583 cout << "Cleaning methode: "<< cleanflag << endl;
584 if (islflag == 1 || islflag == 2)
585 cout << "Island calcultation..." << "using algorithm #" << kalgorithm <<endl;
586 if (islflag == 2)
587 {
588 cout << "Island Cleaning: "<< kmethod <<" method "<< lnew << " new threshold" << endl;
589 }
590 cout << "***********" << endl << endl;
591
592 if(!pediter.GetNumEntries())
593 {
594 cout << "No pedestal file name specified" << endl;
595 return kFALSE;
596 }
597 if(!caliter.GetNumEntries() && calflag>0)
598 {
599 cout << "No calibration file name specified" << endl;
600 return kFALSE;
601 }
602 if(!datiter.GetNumEntries())
603 {
604 cout << "No data file name specified" << endl;
605 return kFALSE;
606 }
607 if(!outname.Length())
608 {
609 cout << "No output file name specified" << endl;
610 return kFALSE;
611 }
612
613
614 return kTRUE;
615}
616
Note: See TracBrowser for help on using the repository browser.