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

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