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

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