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

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