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

Last change on this file since 4409 was 4392, checked in by blanch, 21 years ago
*** empty log message ***
File size: 15.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 "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;
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.SetExtractor(&extractor);
172 calloop.SetInput(&caliter);
173 calloop.SetBadPixels(pedcalloop.GetBadPixels());
174 if(calflag==2)
175 calloop.SetUseBlindPixel();
176 if(calflag>0)
177 if (!calloop.Process(pedcalloop.GetPedestalCam()))
178 return;
179
180
181 /************************************************************************/
182 /* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */
183 /************************************************************************/
184 // First Compute the pedestals
185 MJPedestal pedloop;
186 pedloop.SetInput(&pediter);
187
188 if (!pedloop.Process())
189 return;
190
191 MParList plist3;
192 MTaskList tlist3;
193 plist3.AddToList(&tlist3);
194
195 // containers
196 MCerPhotEvt nphot;
197 MPedPhotCam nphotrms;
198 MExtractedSignalCam sigcam;
199
200 plist3.AddToList(&geomcam);
201 plist3.AddToList(&pedloop.GetPedestalCam());
202 plist3.AddToList(&calloop.GetCalibrationCam());
203 plist3.AddToList(&calloop.GetQECam());
204 plist3.AddToList(&calloop.GetRelTimeCam());
205 plist3.AddToList(&calloop.GetBadPixels());
206 plist3.AddToList(&sigcam);
207 plist3.AddToList(&nphot);
208 plist3.AddToList(&nphotrms);
209
210
211 MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault;
212 if(calflag==0)
213 calMode=MCalibrate::kNone;
214 if(calflag==-1)
215 calMode=MCalibrate::kDummy;
216
217 //tasks
218 MReadMarsFile read3("Events");
219 static_cast<MRead&>(read3).AddFiles(pediter);
220 read3.DisableAutoScheme();
221
222 MCalibrate photcalc(calMode);
223 MPedPhotCalc photrmscalc;
224
225 tlist3.AddToList(&read3);
226 tlist3.AddToList(&geomapl);
227 tlist3.AddToList(&extractor);
228 tlist3.AddToList(&photcalc);
229 tlist3.AddToList(&photrmscalc);
230
231 // Create and setup the eventloop
232 MEvtLoop evtloop3;
233 evtloop3.SetParList(&plist3);
234 if (!evtloop3.Eventloop())
235 return;
236
237 tlist3.PrintStatistics();
238
239 /************************************************************************/
240 /* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */
241 /************************************************************************/
242
243 MParList plist4;
244 MTaskList tlist4;
245 plist4.AddToList(&tlist4);
246
247 // containers
248 MHillas hillas;
249 MNewImagePar newimagepar;
250 MSrcPosCam source;
251 MRawRunHeader runhead;
252
253 MArrivalTimeCam timecam;
254
255
256 // islands
257 MIslands isl;
258 MIslands isl2;
259 MIslands isl3;
260
261 isl.SetName("MIslands");
262 isl2.SetName("MIslands2");
263 isl3.SetName("MIslands3");
264
265 plist4.AddToList(&timecam);
266 plist4.AddToList(&isl);
267
268 if (islflag == 2)
269 plist4.AddToList(&isl2);
270 if (islflag == 3)
271 plist4.AddToList(&isl3);
272
273 plist4.AddToList(&geomcam);
274 plist4.AddToList(&pedloop.GetPedestalCam());
275 plist4.AddToList(&calloop.GetCalibrationCam());
276 plist4.AddToList(&calloop.GetQECam());
277 plist4.AddToList(&calloop.GetRelTimeCam());
278 plist4.AddToList(&calloop.GetBadPixels());
279 plist4.AddToList(&nphot);
280 plist4.AddToList(&nphotrms);
281 plist4.AddToList(&source);
282 plist4.AddToList(&hillas);
283 plist4.AddToList(&newimagepar);
284 plist4.AddToList(&runhead);
285
286 // cuts
287 MF cut(filter);
288
289 //tasks
290 MReadMarsFile read4("Events");
291 static_cast<MRead&>(read4).AddFiles(datiter);
292 read4.DisableAutoScheme();
293
294 MImgCleanStd clean(lcore,ltail);
295 clean.SetCleanRings(lrings);
296 MImgCleanStd::CleaningMethod_t cleanMeth= MImgCleanStd::kStandard;
297 if(cleanflag==2)
298 cleanMeth=MImgCleanStd::kDemocratic;
299 clean.SetMethod(cleanMeth);
300
301 MArrivalTimeCalc2 timecalc;
302 MIslandCalc island;
303 island.SetOutputName("MIslands");
304 island.SetAlgorithm(kalgorithm);
305
306 MBadPixelsTreat interpolatebadpixels;
307 interpolatebadpixels.SetUseInterpolation();
308
309 MIslandClean islclean(lnew);
310 islclean.SetInputName("MIslands");
311 islclean.SetMethod(kmethod);
312
313 MIslandCalc island2;
314 island2.SetOutputName("MIslands2");
315 island2.SetAlgorithm(kalgorithm);
316
317 MIslandCalc island3;
318 island3.SetOutputName("MIslands3");
319
320
321 MHillasCalc hcalc;
322 MHillasSrcCalc csrc1;
323
324 MContinue applycut(&cut);
325 applycut.SetInverted(kTRUE);
326 MWriteRootFile write(outname,"RECREATE");
327
328 MHillasDisplay* disphillas=NULL;
329
330 write.AddContainer("MHillas" , "Parameters");
331 write.AddContainer("MHillasSrc" , "Parameters");
332 write.AddContainer("MHillasExt" , "Parameters");
333 write.AddContainer("MNewImagePar" , "Parameters");
334 write.AddContainer("MRawEvtHeader" , "Parameters");
335 write.AddContainer("MRawRunHeader" , "Parameters");
336 write.AddContainer("MConcentration" , "Parameters");
337 write.AddContainer("MSrcPosCam" , "Parameters");
338 write.AddContainer("MIslands" , "Parameters");
339
340 if (islflag == 2)
341 write.AddContainer("MIslands2" , "Parameters");
342 if (islflag == 3)
343 write.AddContainer("MIslands3" , "Parameters");
344
345
346 if(display)
347 {
348 disphillas = new MHillasDisplay(&nphot,&geomcam);
349 disphillas->SetIslandsName("MIslands");
350 if (islflag == 2)
351 disphillas->SetIslandsName("MIslands2");
352 if (islflag == 3)
353 disphillas->SetIslandsName("MIslands3");
354 }
355
356 tlist4.AddToList(&read4);
357 tlist4.AddToList(&geomapl);
358 tlist4.AddToList(&extractor);
359 tlist4.AddToList(&photcalc);
360 if(calflag==11 || calflag==21)
361 tlist4.AddToList(&interpolatebadpixels);
362 tlist4.AddToList(&clean);
363 tlist4.AddToList(&timecalc);
364 tlist4.AddToList(&island);
365
366 if (islflag == 2)
367 {
368 tlist4.AddToList(&islclean);
369 tlist4.AddToList(&island2);
370 }
371
372 if (islflag == 3)
373 {
374 tlist4.AddToList(&islclean);
375 tlist4.AddToList(&island3);
376 }
377
378 tlist4.AddToList(&hcalc);
379 tlist4.AddToList(&csrc1);
380 if(filter.Length())
381 tlist4.AddToList(&applycut);
382 tlist4.AddToList(&write);
383 if(display)
384 {
385 disphillas->SetPSFile();
386 disphillas->SetPSFileName(psfilename);
387 if(display==2)
388 disphillas->SetPause(kFALSE);
389 tlist4.AddToList(disphillas);
390 }
391
392 // Create and setup the eventloop
393 MEvtLoop datloop;
394 datloop.SetParList(&plist4);
395
396 cout << "*************************************************************" << endl;
397 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
398 cout << "*************************************************************" << endl;
399
400 if (!datloop.Eventloop(nmaxevents))
401 return;
402
403 tlist4.PrintStatistics();
404
405}
406//-------------------------------------------------------------------------------
407
408Bool_t readDatacards(TString& filename)
409{
410 ifstream ifun(filename.Data());
411 if(!ifun)
412 {
413 cout << "File " << filename << " not found" << endl;
414 return kFALSE;
415 }
416
417 TString word;
418
419 while(ifun >> word)
420 {
421 // skip comments
422 if(word[0]=='/' && word[1]=='/')
423 {
424 while(ifun.get()!='\n'); // skip line
425 continue;
426 }
427
428 // number of events
429 if(strcmp(word.Data(),"NEVENTS")==0)
430 ifun >> nmaxevents;
431
432
433 // input file directory
434 if(strcmp(word.Data(),"IDIR")==0)
435 {
436 if(idirname.Length())
437 cout << "readDataCards Warning: overriding input directory file name" << endl;
438 ifun >> idirname;
439 }
440
441 // pedestal runs
442 if(strcmp(word.Data(),"PRUNS")==0)
443 {
444 if(pediter.GetNumRuns())
445 cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl;
446 ifun >> word;
447 pediter.AddRuns(word.Data(),idirname.Data());
448 }
449
450 // pedestal runs for calibration
451 if(strcmp(word.Data(),"PCRUNS")==0)
452 {
453 if(pedcaliter.GetNumRuns())
454 cout << "readDataCards Warning: adding pedestal runs for calibration to the existing list" << endl;
455 ifun >> word;
456 pedcaliter.AddRuns(word.Data(),idirname.Data());
457 }
458
459 // calibration runs
460 if(strcmp(word.Data(),"CRUNS")==0)
461 {
462 if(caliter.GetNumRuns())
463 cout << "readDataCards Warning: adding calibration runs to the existing list" << endl;
464 ifun >> word;
465 caliter.AddRuns(word.Data(),idirname.Data());
466 }
467
468 // data runs
469 if(strcmp(word.Data(),"DRUNS")==0)
470 {
471 if(datiter.GetNumRuns())
472 cout << "readDataCards Warning: adding data runs to the existing list" << endl;
473 ifun >> word;
474 datiter.AddRuns(word.Data(),idirname.Data());
475 }
476
477 // output file name
478 if(strcmp(word.Data(),"OUTFILE")==0)
479 {
480 if(outname.Length())
481 cout << "readDataCards Warning: overriding output file name" << endl;
482 ifun >> outname;
483 }
484
485 // exclusion cut
486 if(strcmp(word.Data(),"FILTER")==0)
487 {
488 if(filter.Length())
489 cout << "readDataCards Warning: overriding existing cut" << endl;
490
491 char ch;
492 while((ch=ifun.get())!='\n')
493 filter.Append(ch);
494 }
495
496 // display flag
497 if(strcmp(word.Data(),"DISPLAY")==0)
498 ifun >> display;
499
500 // ps file name
501 if(strcmp(word.Data(),"PSFILENAME")==0)
502 ifun >> psfilename;
503
504 // calibration flag
505 if(strcmp(word.Data(),"CALFLAG")==0)
506 ifun >> calflag;
507
508
509 // cleaning level
510 if(strcmp(word.Data(),"CLEANLEVEL")==0)
511 {
512 ifun >> lcore;
513 ifun >> ltail;
514 ifun >> lrings;
515 ifun >> cleanflag;
516 }
517
518 if(strcmp(word.Data(),"ISLFLAG")==0)
519 {
520 ifun >> islflag;
521
522 // if (islflag == 1 || islflag == 2)
523 ifun >> kalgorithm;
524 }
525
526 // island cleaning
527 if (islflag == 2){
528 if(strcmp(word.Data(),"ISLANDCLEAN")==0)
529 {
530 ifun >> kmethod;
531 ifun >> lnew;
532 }
533 }
534 }
535
536 pediter.Reset();
537 caliter.Reset();
538 datiter.Reset();
539 TString pfile;
540
541 // Dump read values
542 cout << "************************************************" << endl;
543 cout << "* Datacards read from file " << filename << endl;
544 cout << "************************************************" << endl;
545 cout << "Pedestal file (s): " << endl;
546 while(!(pfile=pediter.Next()).IsNull())
547 cout << pfile << endl;
548 cout << "Calibration file (s): " << endl;
549 while(!(pfile=caliter.Next()).IsNull())
550 cout << pfile << endl;
551 cout << "Data file (s): " << endl;
552 while(!(pfile=datiter.Next()).IsNull())
553 cout << pfile << endl;
554 cout << "Maximum number of events: " << nmaxevents << endl;
555 if(filter.Length())
556 cout << "Applying rejection cut: " << filter << endl;
557 cout << "Output file name: " << outname << endl;
558 if(display)
559 cout << "Generating PS file: " << psfilename << endl;
560 cout << "Calibration: ";
561 if(calflag==0)
562 cout << "Pixel area proportional intercalibration" << endl;
563 else if(calflag==-1)
564 cout << "No calibration whatsoever" << endl;
565 else if(calflag==1)
566 cout << "Default calibration" << endl;
567 else if(calflag==11)
568 cout << "Default calibration + bad pixels interpolation" << endl;
569 cout << "Cleaning level: ("<<lcore<<","<<ltail<<") - " << lrings << "ring" << endl;
570 cout << "Cleaning methode: "<< cleanflag << endl;
571 if (islflag == 1 || islflag == 2)
572 cout << "Island calcultation..." << "using algorithm #" << kalgorithm <<endl;
573 if (islflag == 2)
574 {
575 cout << "Island Cleaning: "<< kmethod <<" method "<< lnew << " new threshold" << endl;
576 }
577 cout << "***********" << endl << endl;
578
579 if(!pediter.GetNumEntries())
580 {
581 cout << "No pedestal file name specified" << endl;
582 return kFALSE;
583 }
584 if(!caliter.GetNumEntries() && calflag>0)
585 {
586 cout << "No calibration file name specified" << endl;
587 return kFALSE;
588 }
589 if(!datiter.GetNumEntries())
590 {
591 cout << "No data file name specified" << endl;
592 return kFALSE;
593 }
594 if(!outname.Length())
595 {
596 cout << "No output file name specified" << endl;
597 return kFALSE;
598 }
599
600
601 return kTRUE;
602}
603
Note: See TracBrowser for help on using the repository browser.