source: trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillas.C@ 5160

Last change on this file since 5160 was 4429, checked in by jlopez, 21 years ago
*** empty log message ***
File size: 14.1 KB
Line 
1/*********************************/
2/* Compute the hillas parameters */
3/*********************************/
4
5Bool_t readDatacards(TString& filename);
6
7// initial and final time slices to be used in signal extraction
8const Byte_t hifirst = 0;
9const Byte_t hilast = 13;
10const Byte_t lofirst = 3;
11const Byte_t lolast = 12;
12
13// declaration of variables read from datacards
14TString outname;
15TString idirname;
16MRunIter pedcaliter;
17MRunIter caliter;
18MRunIter pediter;
19MRunIter datiter;
20ULong_t nmaxevents=999999999;
21Short_t calflag=1;
22Float_t lcore = 3.0;
23Float_t ltail = 1.5;
24Int_t islflag = 0;
25Float_t lnew = 40;
26Int_t kmethod = 1;
27Int_t nfiles = 0;
28
29const TString defaultcard="input.datacard";
30
31//Slow control varialbles
32Short_t slowflag=0;
33Bool_t isSlowControlAvailable = kFALSE;
34TString hvConfigFile="/mnt/users/jlopez/Mars/Files4Mars/Config/HVSettings_FF35q.conf";
35TString continuosLightFile="/local_disk/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
36
37
38/*************************************************************/
39int makeHillas(TString datacard = defaultcard)
40{
41
42 if(!datacard.Length())
43 datacard = defaultcard;
44
45 if(!readDatacards(datacard))
46 {
47 cout << "Error reading datacards. Stoping" << endl;
48 return -1;
49 }
50
51 // Set the general tasks/containers
52 MExtractFixedWindow extractor;
53 extractor.SetRange(hifirst,hilast,lofirst,lolast);
54
55 MCalibrationQECam qecam;
56 MBadPixelsCam badcam;
57 MGeomCamMagic geomcam;
58 MGeomApply geomapl;
59
60 cout << endl;
61 cout << "/****************************************************/" << endl;
62 cout << "/* FIRST LOOP: PEDESTAL FOR CALIBRATION COMPUTATION */" << endl;
63 cout << "/****************************************************/" << endl;
64 cout << endl;
65
66 // If you want to exclude pixels from the beginning, read
67 // an ascii-file with the corr. pixel numbers (see MBadPixelsCam)
68 ifstream fin("badpixels.dat");
69 badcam.AsciiRead(fin);
70 fin.close();
71
72 MJPedestal pedloop0;
73 pedloop0.SetOutputPath("./");
74 pedloop0.SetInput(&pedcaliter);
75 pedloop0.SetExtractor(&extractor);
76 pedloop0.SetBadPixels(badcam);
77
78 if (!pedloop0.Process())
79 return;
80
81 cout << endl;
82 cout << "/*****************************/" << endl;
83 cout << "/* SECOND LOOP: CALIBRATION */" << endl;
84 cout << "/*****************************/" << endl;
85 cout << endl;
86
87 MJCalibration calloop;
88
89 calloop.SetOutputPath("./");
90 calloop.SetExtractor(&extractor);
91 calloop.SetInput(&caliter);
92 calloop.SetQECam(qecam);
93 calloop.SetBadPixels(pedloop0.GetBadPixels());
94 if (!calloop.Process(pedloop0.GetPedestalCam()))
95 return;
96
97 cout << endl;
98 cout << "/*************************************************/" << endl;
99 cout << "/* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */" << endl;
100 cout << "/*************************************************/" << endl;
101 cout << endl;
102
103 // First Compute the pedestals
104 MJPedestal pedloop;
105 pedloop.SetOutputPath("./");
106 pedloop.SetInput(&pediter);
107
108 if (!pedloop.Process())
109 return;
110
111 MPedestalCam pedcam = pedloop.GetPedestalCam();
112
113 // Now convert this pedestals in photons
114 MParList plist3;
115 MTaskList tlist3;
116 plist3.AddToList(&tlist3);
117
118 // containers
119 MCerPhotEvt nphot;
120 MPedPhotCam nphotrms;
121 MExtractedSignalCam sigcam;
122
123 plist3.AddToList(&geomcam);
124 plist3.AddToList(&pedcam);
125 plist3.AddToList(&calloop.GetCalibrationCam());
126 plist3.AddToList(&qecam);
127 plist3.AddToList(&badcam);
128 plist3.AddToList(&sigcam);
129 plist3.AddToList(&nphot);
130 plist3.AddToList(&nphotrms);
131
132
133 MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault;
134 if(calflag==0)
135 calMode=MCalibrate::kNone;
136 if(calflag==-1)
137 calMode=MCalibrate::kDummy;
138
139 //tasks
140 MReadMarsFile read3("Events");
141 static_cast<MRead&>(read3).AddFiles(pediter);
142 read3.DisableAutoScheme();
143
144// MReadReports read3;
145// read3.AddTree("Events","MTime.",kFALSE);
146// static_cast<MRead&>(read3).AddFiles(pediter);
147
148 // MExtractSignal extsig;
149 // extsig.SetRange(hifirst,hilast,lofirst,lolast);
150 MCalibrate photcalc(calMode);
151 MPedPhotCalc photrmscalc;
152
153 tlist3.AddToList(&read3);
154 tlist3.AddToList(&geomapl);
155 tlist3.AddToList(&extractor);
156 tlist3.AddToList(&photcalc);
157 tlist3.AddToList(&photrmscalc);
158
159 // Create and setup the eventloop
160 MEvtLoop evtloop3;
161 evtloop3.SetParList(&plist3);
162 if (!evtloop3.Eventloop())
163 return;
164
165 tlist3.PrintStatistics();
166
167 cout << "/**********************************************/" << endl;
168 cout << "/* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */" << endl;
169 cout << "/**********************************************/" << endl;
170
171 MParList plist4;
172 MTaskList tlist4;
173 plist4.AddToList(&tlist4);
174
175 // containers
176 MHillas hillas;
177 MSrcPosCam source;
178 MRawRunHeader runhead;
179
180 MArrivalTimeCam timecam;
181
182 MIslands isl;
183 isl.SetName("MIslands1");
184
185 MIslands isl2;
186 isl2.SetName("MIslands2");
187
188 MReportDrive drive;
189 MStarLocalCam starcam;
190
191 if (islflag == 1 || islflag == 2)
192 {
193 plist4.AddToList(&timecam);
194 plist4.AddToList(&isl);
195 }
196
197 if (islflag == 2)
198 {
199 plist4.AddToList(&isl2);
200 }
201
202 if (isSlowControlAvailable)
203 {
204 plist4.AddToList(&starcam);
205 }
206
207
208 plist4.AddToList(&geomcam);
209 plist4.AddToList(&pedcam);
210 plist4.AddToList(&calloop.GetCalibrationCam());
211 plist4.AddToList(&qecam);
212 plist4.AddToList(&badcam);
213 plist4.AddToList(&nphot);
214 plist4.AddToList(&nphotrms);
215 plist4.AddToList(&source);
216 plist4.AddToList(&hillas);
217 plist4.AddToList(&runhead);
218
219 //tasks
220// MReadMarsFile read4("Events");
221// static_cast<MRead&>(read4).AddFiles(datiter);
222// read4.DisableAutoScheme();
223
224 MReadReports read4;
225 read4.AddTree("Events","MTime.",kTRUE);
226 read4.AddTree("Currents");
227 read4.AddTree("Camera");
228 read4.AddTree("Drive");
229 static_cast<MRead&>(read4).AddFiles(datiter);
230 read4.AddToBranchList("MReportCurrents.*");
231
232 // Task that use Slow control information
233 MFHVNotNominal fhvnotnominal;
234 fhvnotnominal.SetHVNominalValues(hvConfigFile);
235 fhvnotnominal.SetMaxNumPixelsDeviated(40);
236
237 MContinue hvnotnominal(&fhvnotnominal);
238
239 MCalibrateDC dccalib;
240 dccalib.SetFileName(continuosLightFile);
241
242 const Int_t numblind = 5;
243 const Short_t w[numblind] = { 47, 124, 470, 475, 571};
244 const TArrayS blindpixels(numblind,(Short_t*)w);
245 Float_t ringinterest = 100; //[mm]
246 Float_t tailcut = 3.5;
247 UInt_t integratedevents = 4;
248
249 MFindStars findstars;
250 findstars.SetBlindPixels(blindpixels);
251 findstars.SetRingInterest(ringinterest);
252 findstars.SetDCTailCut(tailcut);
253 findstars.SetNumIntegratedEvents(integratedevents);
254 findstars.SetMinuitPrintOutLevel(-1);
255
256 MSrcPosFromStars srcposfromstar;
257
258 MBadPixelsTreat interpolatebadpixels;
259 interpolatebadpixels.SetUseInterpolation();
260
261 MImgCleanStd clean(lcore,ltail);
262
263 MArrivalTimeCalc2 timecalc;
264
265 MIslandCalc island;
266 island.SetOutputName("MIslands1");
267
268 MIslandClean islclean(lnew);
269 islclean.SetInputName("MIslands1");
270 islclean.SetMethod(kmethod);
271
272 MIslandCalc island2;
273 island2.SetOutputName("MIslands2");
274
275
276 MHillasCalc hcalc;
277 MHillasSrcCalc csrc1;
278
279 MWriteRootFile write(outname,"RECREATE");
280
281 write.AddContainer("MTime" , "Parameters");
282 write.AddContainer("MHillas" , "Parameters");
283 write.AddContainer("MHillasSrc" , "Parameters");
284 write.AddContainer("MHillasExt" , "Parameters");
285 write.AddContainer("MNewImagePar" , "Parameters");
286 write.AddContainer("MRawEvtHeader" , "Parameters");
287 write.AddContainer("MRawRunHeader" , "Parameters");
288 write.AddContainer("MConcentration" , "Parameters");
289 write.AddContainer("MSrcPosCam" , "Parameters");
290
291 if (islflag == 1 || islflag == 2)
292 write.AddContainer("MIslands1" , "Parameters");
293 if (islflag == 2)
294 write.AddContainer("MIslands2" , "Parameters");
295 if (isSlowControlAvailable)
296 {
297 write.AddContainer("MStarLocalCam" , "Parameters");
298 write.AddContainer("MReportDrive" , "Parameters");
299 }
300
301 tlist4.AddToList(&read4);
302 tlist4.AddToList(&geomapl);
303 if (isSlowControlAvailable)
304 {
305 tlist4.AddToList(&hvnotnominal,"Events");
306 tlist4.AddToList(&dccalib,"Currents");
307 tlist4.AddToList(&findstars,"Currents");
308 tlist4.AddToList(&srcposfromstar,"Events");
309 }
310
311 tlist4.AddToList(&extractor,"Events");
312 tlist4.AddToList(&photcalc,"Events");
313 if(calflag==11)
314 tlist4.AddToList(&interpolatebadpixels);
315 tlist4.AddToList(&clean,"Events");
316
317 if (islflag == 1 || islflag == 2)
318 {
319 tlist4.AddToList(&timecalc,"Events");
320 tlist4.AddToList(&island,"Events");
321 }
322
323 if (islflag == 2)
324 {
325 tlist4.AddToList(&islclean,"Events");
326 tlist4.AddToList(&island2,"Events");
327 }
328
329 //tlist4.AddToList(&blind2,"Events");
330 tlist4.AddToList(&hcalc,"Events");
331 tlist4.AddToList(&csrc1,"Events");
332 tlist4.AddToList(&write,"Events");
333
334 // Create and setup the eventloop
335 MEvtLoop datloop;
336 datloop.SetParList(&plist4);
337
338 cout << endl;
339 cout << "******************************************************" << endl;
340 cout << "* COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) *" << endl;
341 cout << "******************************************************" << endl;
342 cout << endl;
343
344 if (!datloop.Eventloop(nmaxevents))
345 return;
346
347 tlist4.PrintStatistics();
348
349}
350//-------------------------------------------------------------------------------
351
352Bool_t readDatacards(TString& filename)
353{
354 ifstream ifun(filename.Data());
355 if(!ifun)
356 {
357 cout << "File " << filename << " not found" << endl;
358 return kFALSE;
359 }
360
361 TString word;
362
363 while(ifun >> word)
364 {
365
366 // skip comments
367 if(word[0]=='/' && word[1]=='/')
368 {
369 while(ifun.get()!='\n'); // skip line
370 continue;
371 }
372
373 // number of events
374 if(strcmp(word.Data(),"NEVENTS")==0)
375 ifun >> nmaxevents;
376
377
378 // input file directory
379 if(strcmp(word.Data(),"IDIR")==0)
380 {
381 if(idirname.Length())
382 cout << "readDataCards Warning: overriding input directory file name" << endl;
383 ifun >> idirname;
384 }
385
386 // pedestal runs for calibration
387 if(strcmp(word.Data(),"PCRUNS")==0)
388 {
389 if(pedcaliter.GetNumRuns())
390 cout << "readDataCards Warning: adding pedestal runs for calibration to the existing list" << endl;
391 ifun >> word;
392 pedcaliter.AddRuns(word.Data(),idirname.Data());
393 }
394
395 // calibration runs
396 if(strcmp(word.Data(),"CRUNS")==0)
397 {
398 if(caliter.GetNumRuns())
399 cout << "readDataCards Warning: adding calibration runs to the existing list" << endl;
400 ifun >> word;
401 caliter.AddRuns(word.Data(),idirname.Data());
402 }
403
404 // pedestal runs for data
405 if(strcmp(word.Data(),"PRUNS")==0)
406 {
407 if(pediter.GetNumRuns())
408 cout << "readDataCards Warning: adding pedestal runs for data to the existing list" << endl;
409 ifun >> word;
410 pediter.AddRuns(word.Data(),idirname.Data());
411 }
412
413 // data runs
414 if(strcmp(word.Data(),"DRUNS")==0)
415 {
416 if(datiter.GetNumRuns())
417 cout << "readDataCards Warning: adding data runs to the existing list" << endl;
418 ifun >> word;
419 datiter.AddRuns(word.Data(),idirname.Data());
420 }
421
422 // output file name
423 if(strcmp(word.Data(),"OUTFILE")==0)
424 {
425 if(outname.Length())
426 cout << "readDataCards Warning: overriding output file name" << endl;
427 ifun >> outname;
428 }
429
430 // calibration flag
431 if(strcmp(word.Data(),"CALFLAG")==0)
432 ifun >> calflag;
433
434 // cleaning level
435 if(strcmp(word.Data(),"CLEANLEVEL")==0)
436 {
437 ifun >> lcore;
438 ifun >> ltail;
439 }
440
441 if(strcmp(word.Data(),"ISLFLAG")==0)
442 {
443 ifun >> islflag;
444 }
445
446 // island cleaning
447 if (islflag == 2){
448 if(strcmp(word.Data(),"ISLANDCLEAN")==0)
449 {
450 ifun >> kmethod;
451 ifun >> lnew;
452 }
453 }
454
455 // slow control data available
456 if(strcmp(word.Data(),"SLOWCRT")==0)
457 {
458 ifun >> slowflag;
459 if (slowflag == 1) isSlowControlAvailable = kTRUE;
460 }
461
462 // hv configuration file
463 if(strcmp(word.Data(),"HVCONFFILE")==0)
464 {
465 ifun >> hvConfigFile;
466 }
467
468 // countinous light file to calibrate the dc data
469 if(strcmp(word.Data(),"CLFILE")==0)
470 {
471 ifun >> continuosLightFile;
472 }
473 }
474
475 pedcaliter.Reset();
476 pediter.Reset();
477 caliter.Reset();
478 datiter.Reset();
479 TString pfile;
480
481 // Dump read values
482 cout << "************************************************" << endl;
483 cout << "* Datacards read from file " << filename << endl;
484 cout << "************************************************" << endl;
485 cout << "Pedestal file (s) for calibration: " << endl;
486 while(kTRUE)
487 {
488 pfile=pedcaliter.Next();
489 if (pfile.IsNull())
490 break;
491 cout << pfile << endl;
492 }
493 cout << "Calibration file (s): " << endl;
494 while(kTRUE)
495 {
496 pfile=caliter.Next();
497 if (pfile.IsNull())
498 break;
499 cout << pfile << endl;
500 }
501 cout << "Pedestal file (s) for data: " << endl;
502 while(kTRUE)
503 {
504 pfile=pediter.Next();
505 if (pfile.IsNull())
506 break;
507 cout << pfile << endl;
508 }
509 cout << "Data file (s): " << endl;
510 while(kTRUE)
511 {
512 pfile=datiter.Next();
513 if (pfile.IsNull())
514 break;
515 cout << pfile << endl;
516 }
517 cout << "Maximum number of events: " << nmaxevents << endl;
518 cout << "Output file name: " << outname << endl;
519 cout << "Calibration flag: " << calflag << endl;
520 cout << "Cleaning level: ("<<lcore<<","<<ltail<<")" << endl;
521 if (islflag == 1 || islflag == 2)
522 cout << "Island calcultation..." << endl;
523 if (islflag == 2)
524 {
525 cout << "Island Cleaning: "<< kmethod <<" method "<< lnew << " new threshold" << endl;
526 }
527 cout << "***********" << endl << endl;
528
529 if(!pediter.GetNumEntries())
530 {
531 cout << "No pedestal file name specified" << endl;
532 return kFALSE;
533 }
534 if(!caliter.GetNumEntries())
535 {
536 cout << "No calibration file name specified" << endl;
537 return kFALSE;
538 }
539 if(!datiter.GetNumEntries())
540 {
541 cout << "No data file name specified" << endl;
542 return kFALSE;
543 }
544 if(!outname.Length())
545 {
546 cout << "No output file name specified" << endl;
547 return kFALSE;
548 }
549
550 if (isSlowControlAvailable)
551 {
552 cout << "HV Configuratin file: " << hvConfigFile << endl;
553 cout << "Continous Light file: " << continuosLightFile << endl;
554 }
555
556
557 return kTRUE;
558}
Note: See TracBrowser for help on using the repository browser.