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

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