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

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