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

Last change on this file since 3947 was 3947, checked in by rico, 21 years ago
*** empty log message ***
File size: 26.5 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 "MReadMarsFile.h"
13#include "MGeomApply.h"
14#include "MPedCalcPedRun.h"
15#include "MEvtLoop.h"
16#include "MGeomCamMagic.h"
17#include "MExtractedSignalCam.h"
18#include "MCalibrationChargeCam.h"
19#include "MHCalibrationChargeCam.h"
20#include "MHCalibrationRelTimeCam.h"
21#include "MExtractSignal.h"
22#include "MCalibrationChargeCalc.h"
23#include "MFCosmics.h"
24#include "MContinue.h"
25#include "MFillH.h"
26#include "MLog.h"
27#include "MCerPhotEvt.h"
28#include "MPedPhotCam.h"
29#include "MCalibrate.h"
30#include "MPedPhotCalc.h"
31#include "MHillas.h"
32#include "MRawRunHeader.h"
33#include "MSrcPosCam.h"
34#include "MBlindPixelCalc.h"
35#include "MImgCleanStd.h"
36#include "MHillasSrcCalc.h"
37#include "MHillasCalc.h"
38#include "MWriteRootFile.h"
39#include "MProgressBar.h"
40#include "MArgs.h"
41#include "MRunIter.h"
42#include "MJPedestal.h"
43#include "MJCalibration.h"
44
45#include <iostream>
46#include <fstream>
47#include <stdlib.h>
48
49using namespace std;
50
51Bool_t readDatacards(TString& filename);
52void makeHillas();
53
54// initial and final time slices to be used in signal extraction
55const Byte_t hifirst = 1;
56const Byte_t hilast = 14;
57const Byte_t lofirst = 4;
58const Byte_t lolast = 13;
59
60// declaration of variables read from datacards
61TString outname;
62TString idirname;
63MRunIter caliter;
64MRunIter pediter;
65MRunIter datiter;
66ULong_t nmaxevents=999999999;
67Short_t calflag=1;
68Float_t lcore = 3.0;
69Float_t ltail = 1.5;
70Int_t nfiles = 0;
71
72const TString defaultcard="input.datacard";
73static void Usage()
74{
75 gLog <<endl;
76 gLog << "Usage is:" << endl;
77 gLog << " makeHillas [-h] [-?] <datacards>" << endl << endl;
78 gLog << " <datacards>: datacards file name (dafault input.datacards)" << endl;
79 gLog << " -?/-h: This help" << endl << endl;
80}
81
82int main(int argc, char **argv)
83{
84 // evaluate arguments
85 MArgs arg(argc, argv);
86 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
87 {
88 Usage();
89 return -1;
90 }
91
92 TString datacard = arg.GetArgumentStr(0);
93 if(!datacard.Length())
94 datacard = defaultcard;
95
96 if(!readDatacards(datacard))
97 {
98 cout << "Error reading datacards. Stoping" << endl;
99 return -1;
100 }
101 makeHillas();
102}
103
104/*************************************************************/
105void makeHillas()
106{
107#if 0
108 MStatusDisplay *display = new MStatusDisplay;
109 display->SetUpdateTime(3000);
110 display->Resize(850,700);
111
112 gStyle->SetOptStat(1111);
113 gStyle->SetOptFit();
114#endif
115
116 /************************************/
117 /* FIRST LOOP: PEDESTAL COMPUTATION */
118 /************************************/
119
120 MParList plist1;
121 MTaskList tlist1;
122 plist1.AddToList(&tlist1);
123
124 // containers
125 MPedestalCam pedcam;
126 MBadPixelsCam badcam;
127 // badcam.AsciiRead("badpixels.dat");
128
129 plist1.AddToList(&pedcam);
130 plist1.AddToList(&badcam);
131
132 //tasks
133 MReadMarsFile read("Events");
134 static_cast<MRead&>(read).AddFiles(pediter);
135 read.DisableAutoScheme();
136
137 MGeomApply geomapl;
138 MPedCalcPedRun pedcalc;
139
140 tlist1.AddToList(&read);
141 tlist1.AddToList(&geomapl);
142 tlist1.AddToList(&pedcalc);
143
144 // Create and setup the eventloop
145 MEvtLoop pedloop;
146 pedloop.SetParList(&plist1);
147 //pedloop.SetDisplay(display);
148
149 cout << "*************************" << endl;
150 cout << "** COMPUTING PEDESTALS **" << endl;
151 cout << "*************************" << endl;
152
153 if (!pedloop.Eventloop())
154 return;
155
156 tlist1.PrintStatistics();
157
158 /*****************************/
159 /* SECOND LOOP: CALIBRATION */
160 /*****************************/
161
162 MParList plist2;
163 MTaskList tlist2;
164 plist2.AddToList(&tlist2);
165 plist2.AddToList(&pedcam);
166 plist2.AddToList(&badcam);
167
168 // new containers
169 MGeomCamMagic geomcam;
170 MExtractedSignalCam sigcam;
171 MCalibrationChargeCam calcam;
172 MHCalibrationChargeCam histcharge;
173 MHCalibrationRelTimeCam histtime;
174
175
176 plist2.AddToList(&geomcam);
177 plist2.AddToList(&sigcam);
178 plist2.AddToList(&calcam);
179 plist2.AddToList(&histcharge);
180
181 //tasks
182 MReadMarsFile read2("Events");
183 static_cast<MRead&>(read2).AddFiles(caliter);
184 read2.DisableAutoScheme();
185
186 MExtractSignal sigcalc;
187 sigcalc.SetRange(hifirst,hilast,lofirst,lolast);
188 MCalibrationChargeCalc calcalc;
189 MFCosmics cosmics;
190 MContinue cont(&cosmics);
191
192
193 calcalc.SetPulserColor(MCalibrationCam::kCT1);
194 MFillH fillcam ("MHCalibrationChargeCam" , "MExtractedSignalCam");
195
196 tlist2.AddToList(&read2);
197 tlist2.AddToList(&geomapl);
198 tlist2.AddToList(&sigcalc);
199 tlist2.AddToList(&cont);
200 tlist2.AddToList(&fillcam);
201 tlist2.AddToList(&calcalc);
202
203 // Create and setup the eventloop
204 MEvtLoop calloop;
205 calloop.SetParList(&plist2);
206 //calloop.SetDisplay(display);
207
208 cout << "***************************" << endl;
209 cout << "** COMPUTING CALIBRATION **" << endl;
210 cout << "***************************" << endl;
211
212 if (!calloop.Eventloop())
213 return;
214
215 tlist2.PrintStatistics();
216
217 MLog gauglog;
218 gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
219 calcam.SetLogStream(&gauglog);
220 calcam.Print();
221 calcam.SetLogStream(&gLog);
222#if 0
223 // Create histograms to display
224 MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges");
225 MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges");
226 MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit");
227 MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas");
228 MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge");
229 MHCamera disp6 (geomcam, "Cal;FFactorPh", "Nr. of Photo-electrons (F-Factor Method)");
230 MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor to photons (F-Factor Method)");
231 MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
232 MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)");
233 MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor to photons (Blind Pixel Method)");
234 MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
235 MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)");
236 MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor tp photons (PIN Diode Method)");
237 MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
238 MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded");
239 MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted");
240 MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results");
241 MHCamera disp18 (geomcam, "Cal;HiGainOscillating", "Oscillating Pixels HI Gain");
242 MHCamera disp19 (geomcam, "Cal;LoGainOscillating", "Oscillating Pixels LO Gain");
243 MHCamera disp20 (geomcam, "Cal;HiGainPickup", "Number Pickup events Hi Gain");
244 MHCamera disp21 (geomcam, "Cal;LoGainPickup", "Number Pickup events Lo Gain");
245 MHCamera disp22 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain");
246 MHCamera disp23 (geomcam, "Cal;FFactorValid", "Pixels with valid F-Factor calibration");
247 MHCamera disp24 (geomcam, "Cal;BlindPixelValid", "Pixels with valid BlindPixel calibration");
248 MHCamera disp25 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration");
249
250 MHCamera disp26 (geomcam, "Cal;Ped", "Pedestals");
251 MHCamera disp27 (geomcam, "Cal;PedRms", "Pedestal RMS");
252
253 MHCamera disp28 (geomcam, "time;Time", "Rel. Arrival Times");
254 MHCamera disp29 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times");
255 MHCamera disp30 (geomcam, "time;TimeProb", "Probability of Time Fit");
256 MHCamera disp31 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
257 MHCamera disp32 (geomcam, "time;Oscillating", "Oscillating Pixels");
258
259 MHCamera disp33 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
260 MHCamera disp34 (geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times");
261
262 // Fitted charge means and sigmas
263 disp1.SetCamContent(calcam, 0);
264 disp1.SetCamError( calcam, 1);
265 disp2.SetCamContent(calcam, 2);
266 disp2.SetCamError( calcam, 3);
267
268 // Fit probabilities
269 disp3.SetCamContent(calcam, 4);
270
271 // Reduced Sigmas and reduced sigmas per charge
272 disp4.SetCamContent(calcam, 5);
273 disp4.SetCamError( calcam, 6);
274 disp5.SetCamContent(calcam, 7);
275 disp5.SetCamError( calcam, 8);
276
277 // F-Factor Method
278 disp6.SetCamContent(calcam, 9);
279 disp6.SetCamError( calcam, 10);
280 disp7.SetCamContent(calcam, 11);
281 disp7.SetCamError( calcam, 12);
282 disp8.SetCamContent(calcam, 13);
283 disp8.SetCamError( calcam, 14);
284
285 // Blind Pixel Method
286 disp9.SetCamContent(calcam, 15);
287 disp9.SetCamError( calcam, 16);
288 disp10.SetCamContent(calcam,17);
289 disp10.SetCamError( calcam,18);
290 disp11.SetCamContent(calcam,19);
291 disp11.SetCamError( calcam,20);
292
293 // PIN Diode Method
294 disp12.SetCamContent(calcam,21);
295 disp12.SetCamError( calcam,22);
296 disp13.SetCamContent(calcam,23);
297 disp13.SetCamError( calcam,24);
298 disp14.SetCamContent(calcam,25);
299 disp14.SetCamError( calcam,26);
300
301 // Pixels with defects
302 disp15.SetCamContent(calcam,27);
303 disp16.SetCamContent(calcam,28);
304 disp17.SetCamContent(badcam,9);
305 disp18.SetCamContent(badcam,16);
306 disp19.SetCamContent(badcam,15);
307 disp20.SetCamContent(calcam,29);
308 disp21.SetCamContent(calcam,30);
309
310 // Lo Gain calibration
311 disp22.SetCamContent(calcam,31);
312
313 // Valid flags
314 disp23.SetCamContent(calcam,32);
315 disp24.SetCamContent(calcam,33);
316 disp25.SetCamContent(calcam,34);
317
318 // Pedestals
319 disp26.SetCamContent(calcam,35);
320 disp26.SetCamError( calcam,36);
321 disp27.SetCamContent(calcam,37);
322 disp27.SetCamError( calcam,38);
323
324 // Relative Times
325 disp28.SetCamContent(histtime,0);
326 disp28.SetCamError( histtime,1);
327 disp29.SetCamContent(histtime,2);
328 disp29.SetCamError( histtime,3);
329 disp30.SetCamContent(histtime,4);
330 disp31.SetCamContent(histtime,5);
331 disp32.SetCamContent(histtime,6);
332
333 // Absolute Times
334 disp33.SetCamContent(calcam,39);
335 disp33.SetCamError( calcam,40);
336 disp34.SetCamContent(calcam,41);
337
338 disp1.SetYTitle("Charge [FADC units]");
339 disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
340 disp3.SetYTitle("P_{Charge} [1]");
341
342 disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
343 disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
344
345 disp6.SetYTitle("Nr. Photo-electrons [1]");
346 disp7.SetYTitle("Conversion Factor [Ph/FADC Count]");
347 disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] ");
348
349 disp9.SetYTitle("Photon flux [ph/mm^2]");
350 disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
351 disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
352
353 disp12.SetYTitle("Photon flux [ph/mm^2]");
354 disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
355 disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
356
357 disp15.SetYTitle("[1]");
358 disp16.SetYTitle("[1]");
359 disp17.SetYTitle("[1]");
360 disp18.SetYTitle("[1]");
361 disp19.SetYTitle("[1]");
362 disp20.SetYTitle("[1]");
363 disp21.SetYTitle("[1]");
364 disp22.SetYTitle("[1]");
365 disp23.SetYTitle("[1]");
366 disp24.SetYTitle("[1]");
367 disp25.SetYTitle("[1]");
368
369 disp26.SetYTitle("Ped [FADC Counts ]");
370 disp27.SetYTitle("RMS_{Ped} [FADC Counts ]");
371
372 disp28.SetYTitle("Time Offset [ns]");
373 disp29.SetYTitle("Timing resolution [ns]");
374 disp30.SetYTitle("P_{Time} [1]");
375
376 disp31.SetYTitle("[1]");
377 disp32.SetYTitle("[1]");
378
379 disp33.SetYTitle("Mean Abs. Time [FADC slice]");
380 disp34.SetYTitle("RMS Abs. Time [FADC slices]");
381
382 gStyle->SetOptStat(1111);
383 gStyle->SetOptFit();
384
385 // Charges
386 TCanvas &c1 = display->AddTab("Fit.Charge");
387 c1.Divide(2, 3);
388
389 CamDraw(c1, disp1,calcam,1, 2 , 2);
390 CamDraw(c1, disp2,calcam,2, 2 , 2);
391
392 // Fit Probability
393 TCanvas &c2 = display->AddTab("Fit.Prob");
394 c2.Divide(1,3);
395
396 CamDraw(c2, disp3,calcam,1, 1 , 4);
397
398 // Reduced Sigmas
399 TCanvas &c3 = display->AddTab("Red.Sigma");
400 c3.Divide(2,3);
401
402 CamDraw(c3, disp4,calcam,1, 2 , 2);
403 CamDraw(c3, disp5,calcam,2, 2 , 2);
404
405 // F-Factor Method
406 TCanvas &c4 = display->AddTab("F-Factor");
407 c4.Divide(3,3);
408
409 CamDraw(c4, disp6,calcam,1, 3 , 2);
410 CamDraw(c4, disp7,calcam,2, 3 , 2);
411 CamDraw(c4, disp8,calcam,3, 3 , 2);
412
413 // Blind Pixel Method
414 TCanvas &c5 = display->AddTab("BlindPix");
415 c5.Divide(3, 3);
416
417 CamDraw(c5, disp9,calcam,1, 3 , 9);
418 CamDraw(c5, disp10,calcam,2, 3 , 2);
419 CamDraw(c5, disp11,calcam,3, 3 , 2);
420
421 // PIN Diode Method
422 TCanvas &c6 = display->AddTab("PINDiode");
423 c6.Divide(3,3);
424
425 CamDraw(c6, disp12,calcam,1, 3 , 9);
426 CamDraw(c6, disp13,calcam,2, 3 , 2);
427 CamDraw(c6, disp14,calcam,3, 3 , 2);
428
429 // Defects
430 TCanvas &c7 = display->AddTab("Defects");
431 c7.Divide(4,2);
432
433 CamDraw(c7, disp15,calcam,1,4, 0);
434 CamDraw(c7, disp16,calcam,2,4, 0);
435 CamDraw(c7, disp20,calcam,3,4, 0);
436 CamDraw(c7, disp21,calcam,4,4, 0);
437
438 // BadCam
439 TCanvas &c8 = display->AddTab("Defects");
440 c8.Divide(3,2);
441
442 CamDraw(c8, disp17,badcam,1,3, 0);
443 CamDraw(c8, disp18,badcam,2,3, 0);
444 CamDraw(c8, disp19,badcam,3,3, 0);
445
446 // Valid flags
447 TCanvas &c9 = display->AddTab("Validity");
448 c9.Divide(4,2);
449
450 CamDraw(c9, disp22,calcam,1,4,0);
451 CamDraw(c9, disp23,calcam,2,4,0);
452 CamDraw(c9, disp24,calcam,3,4,0);
453 CamDraw(c9, disp25,calcam,4,4,0);
454
455 // Pedestals
456 TCanvas &c10 = display->AddTab("Pedestals");
457 c10.Divide(2,3);
458
459 CamDraw(c10,disp26,calcam,1,2,1);
460 CamDraw(c10,disp27,calcam,2,2,2);
461
462 // Rel. Times
463 TCanvas &c11 = display->AddTab("Fitted Rel. Times");
464 c11.Divide(3,3);
465
466 CamDraw(c11,disp28,calcam,1,3,2);
467 CamDraw(c11,disp29,calcam,2,3,2);
468 CamDraw(c11,disp30,calcam,3,3,4);
469
470 // Time Defects
471 TCanvas &c12 = display->AddTab("Time Def.");
472 c12.Divide(2,2);
473
474 CamDraw(c12, disp31,calcam,1,2, 0);
475 CamDraw(c12, disp32,calcam,2,2, 0);
476
477 // Abs. Times
478 TCanvas &c13 = display->AddTab("Abs. Times");
479 c13.Divide(2,3);
480
481 CamDraw(c13,disp33,calcam,1,2,2);
482 CamDraw(c13,disp34,calcam,2,2,2);
483#endif
484 /************************************************************************/
485 /* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */
486 /************************************************************************/
487 MParList plist3;
488 MTaskList tlist3;
489 plist3.AddToList(&tlist3);
490
491 // containers
492 MCerPhotEvt nphot;
493 MPedPhotCam nphotrms;
494
495 plist3.AddToList(&geomcam);
496 plist3.AddToList(&pedcam);
497 plist3.AddToList(&calcam);
498 // plist3.AddToList(&badcam);
499 plist3.AddToList(&sigcam);
500 plist3.AddToList(&nphot);
501 plist3.AddToList(&nphotrms);
502
503
504 MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault;
505 if(calflag==0)
506 calMode=MCalibrate::kNone;
507
508
509 //tasks
510 MReadMarsFile read3("Events");
511 static_cast<MRead&>(read3).AddFiles(pediter);
512 read3.DisableAutoScheme();
513
514 MExtractSignal extsig;
515 extsig.SetRange(hifirst,hilast,lofirst,lolast);
516 MCalibrate photcalc(calMode);
517 MPedPhotCalc photrmscalc;
518
519 tlist3.AddToList(&read3);
520 tlist3.AddToList(&geomapl);
521 tlist3.AddToList(&extsig);
522 tlist3.AddToList(&photcalc);
523 tlist3.AddToList(&photrmscalc);
524
525 // Create and setup the eventloop
526 MEvtLoop evtloop3;
527 evtloop3.SetParList(&plist3);
528 if (!evtloop3.Eventloop())
529 return;
530
531 tlist3.PrintStatistics();
532
533 /************************************************************************/
534 /* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */
535 /************************************************************************/
536
537 MParList plist4;
538 MTaskList tlist4;
539 plist4.AddToList(&tlist4);
540
541 // containers
542 MHillas hillas;
543 MSrcPosCam source;
544 MRawRunHeader runhead;
545
546 plist4.AddToList(&geomcam);
547 plist4.AddToList(&pedcam);
548 plist4.AddToList(&calcam);
549 // plist4.AddToList(&badcam);
550 plist4.AddToList(&nphot);
551 plist4.AddToList(&nphotrms);
552 plist4.AddToList(&source);
553 plist4.AddToList(&hillas);
554 plist4.AddToList(&runhead);
555
556 //tasks
557 MReadMarsFile read4("Events");
558 static_cast<MRead&>(read4).AddFiles(datiter);
559 read4.DisableAutoScheme();
560
561 // set bad pixels
562 MBlindPixelCalc blind;
563 MBlindPixelCalc blind2;
564 const Short_t x[16] = {330,395,329,396,389,
565 323,388,322,384,385,
566 386,387,321,320,319,
567 394};
568 const TArrayS bp(16,(Short_t*)x);
569 blind.SetPixelIndices(bp);
570 blind2.SetPixelIndices(bp);
571
572 MImgCleanStd clean(lcore,ltail);
573 MHillasCalc hcalc;
574 MHillasSrcCalc csrc1;
575
576 MWriteRootFile write(outname,"RECREATE");
577
578 // write.AddContainer("MGeomCam" , "RunHeaders");
579 // write.AddContainer("MRawRunHeader" , "RunHeaders");
580 // write.AddContainer("MSrcPosCam" , "RunHeaders");
581 // write.AddContainer("MCalibrationChargeCam" , "RunHeaders");
582 // write.AddContainer("MPedPhotCam" , "RunHeaders"); // Attention, was in Events - Tree!!
583 // write.AddContainer("MPedestalCam" , "RunHeaders");
584 // write.AddContainer("MHCalibrationRelTimeCam","RunHeaders");
585
586 // write.AddContainer("MCerPhotEvt" , "Events");
587 // write.AddContainer("MRawEvtHeader" , "Events");
588 // write.AddContainer("MBadPixelsCam" , "Events");
589 // write.AddContainer("MPedPhotCam" , "Events");
590
591 write.AddContainer("MHillas" , "Parameters");
592 write.AddContainer("MHillasSrc" , "Parameters");
593 write.AddContainer("MHillasExt" , "Parameters");
594 write.AddContainer("MNewImagePar" , "Parameters");
595 write.AddContainer("MRawEvtHeader" , "Parameters");
596 write.AddContainer("MRawRunHeader" , "Parameters");
597 write.AddContainer("MConcentration" , "Parameters");
598
599 tlist4.AddToList(&read4);
600 tlist4.AddToList(&geomapl);
601 tlist4.AddToList(&extsig);
602 tlist4.AddToList(&photcalc);
603 //tlist4.AddToList(&blind);
604 tlist4.AddToList(&clean);
605 //tlist4.AddToList(&blind2);
606 tlist4.AddToList(&hcalc);
607 // tlist4.AddToList(&srcposcalc);
608 tlist4.AddToList(&csrc1);
609 tlist4.AddToList(&write);
610
611 // Create and setup the eventloop
612 MEvtLoop datloop;
613 datloop.SetParList(&plist4);
614 // MProgressBar bar;
615 // datloop.SetProgressBar(&bar);
616
617 cout << "*************************************************************" << endl;
618 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
619 cout << "*************************************************************" << endl;
620
621 if (!datloop.Eventloop(nmaxevents))
622 return;
623
624 tlist4.PrintStatistics();
625 nphotrms.Print();
626
627}
628
629#if 0
630void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
631{
632
633 c.cd(i);
634 gPad->SetBorderMode(0);
635 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
636 // obj1->AddNotify(evt);
637
638 c.cd(i+j);
639 gPad->SetBorderMode(0);
640 obj1->Draw();
641 ((MHCamera*)obj1)->SetPrettyPalette();
642
643 if (fit != 0)
644 {
645 c.cd(i+2*j);
646 gPad->SetBorderMode(0);
647 TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName());
648
649// obj2->Sumw2();
650 obj2->Draw();
651 obj2->SetBit(kCanDelete);
652
653 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
654 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
655 const Double_t integ = obj2->Integral("width")/2.5066283;
656 const Double_t mean = obj2->GetMean();
657 const Double_t rms = obj2->GetRMS();
658 const Double_t width = max-min;
659
660 if (rms == 0. || width == 0. )
661 return;
662
663 switch (fit)
664 {
665 case 1:
666 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
667 sgaus->SetBit(kCanDelete);
668 sgaus->SetParNames("Area","#mu","#sigma");
669 sgaus->SetParameters(integ/rms,mean,rms);
670 sgaus->SetParLimits(0,0.,integ);
671 sgaus->SetParLimits(1,min,max);
672 sgaus->SetParLimits(2,0,width/1.5);
673 obj2->Fit("sgaus","QLR");
674 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
675 break;
676
677 case 2:
678 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
679 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
680 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
681 dgaus->SetBit(kCanDelete);
682 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
683 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
684 integ/width/2.,(max+mean)/2.,width/4.);
685 // The left-sided Gauss
686 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
687 dgaus->SetParLimits(1,min+(width/10.),mean);
688 dgaus->SetParLimits(2,0,width/2.);
689 // The right-sided Gauss
690 dgaus->SetParLimits(3,0,integ);
691 dgaus->SetParLimits(4,mean,max-(width/10.));
692 dgaus->SetParLimits(5,0,width/2.);
693 obj2->Fit("dgaus","QLRM");
694 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
695 break;
696
697 case 3:
698 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
699 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
700 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
701 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
702 tgaus->SetBit(kCanDelete);
703 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
704 "A_{2}","#mu_{2}","#sigma_{2}",
705 "A_{3}","#mu_{3}","#sigma_{3}");
706 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
707 integ/width/3.,(max+mean)/2.,width/4.,
708 integ/width/3.,mean,width/2.);
709 // The left-sided Gauss
710 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
711 tgaus->SetParLimits(1,min+(width/10.),mean);
712 tgaus->SetParLimits(2,width/15.,width/2.);
713 // The right-sided Gauss
714 tgaus->SetParLimits(3,0.,integ);
715 tgaus->SetParLimits(4,mean,max-(width/10.));
716 tgaus->SetParLimits(5,width/15.,width/2.);
717 // The Gauss describing the outliers
718 tgaus->SetParLimits(6,0.,integ);
719 tgaus->SetParLimits(7,min,max);
720 tgaus->SetParLimits(8,width/4.,width/1.5);
721 obj2->Fit("tgaus","QLRM");
722 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
723 break;
724 case 4:
725 obj2->Fit("pol0","Q");
726 obj2->GetFunction("pol0")->SetLineColor(kYellow);
727 break;
728 case 9:
729 break;
730 default:
731 obj2->Fit("gaus","Q");
732 obj2->GetFunction("gaus")->SetLineColor(kYellow);
733 break;
734 }
735
736 TArrayI s0(3);
737 s0[0] = 6;
738 s0[1] = 1;
739 s0[2] = 2;
740
741 TArrayI s1(3);
742 s1[0] = 3;
743 s1[1] = 4;
744 s1[2] = 5;
745
746 TArrayI inner(1);
747 inner[0] = 0;
748
749 TArrayI outer(1);
750 outer[0] = 1;
751
752 // Just to get the right (maximum) binning
753 TH1D *half[4];
754 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
755 half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
756 half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
757 half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
758
759 for (int i=0; i<4; i++)
760 {
761 half[i]->SetLineColor(kRed+i);
762 half[i]->SetDirectory(0);
763 half[i]->SetBit(kCanDelete);
764 half[i]->Draw("same");
765 }
766
767 gPad->Modified();
768 gPad->Update();
769
770 }
771}
772#endif
773
774Bool_t readDatacards(TString& filename)
775{
776 ifstream ifun(filename.Data());
777 if(!ifun)
778 {
779 cout << "File " << filename << " not found" << endl;
780 return kFALSE;
781 }
782
783 TString word;
784
785 while(ifun >> word)
786 {
787 // skip comments
788 if(word[0]=='/' && word[1]=='/')
789 {
790 while(ifun.get()!='\n'); // skip line
791 continue;
792 }
793
794 // number of events
795 if(strcmp(word.Data(),"NEVENTS")==0)
796 ifun >> nmaxevents;
797
798
799 // input file directory
800 if(strcmp(word.Data(),"IDIR")==0)
801 {
802 if(idirname.Length())
803 cout << "readDataCards Warning: overriding input directory file name" << endl;
804 ifun >> idirname;
805 }
806
807 // pedestal runs
808 if(strcmp(word.Data(),"PRUNS")==0)
809 {
810 if(pediter.GetNumRuns())
811 cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl;
812 ifun >> word;
813 pediter.AddRuns(word.Data(),idirname.Data());
814 }
815
816 // calibration runs
817 if(strcmp(word.Data(),"CRUNS")==0)
818 {
819 if(caliter.GetNumRuns())
820 cout << "readDataCards Warning: adding calibration runs to the existing list" << endl;
821 ifun >> word;
822 caliter.AddRuns(word.Data(),idirname.Data());
823 }
824
825 // data runs
826 if(strcmp(word.Data(),"DRUNS")==0)
827 {
828 if(datiter.GetNumRuns())
829 cout << "readDataCards Warning: adding data runs to the existing list" << endl;
830 ifun >> word;
831 datiter.AddRuns(word.Data(),idirname.Data());
832 }
833
834 // output file name
835 if(strcmp(word.Data(),"OUTFILE")==0)
836 {
837 if(outname.Length())
838 cout << "readDataCards Warning: overriding output file name" << endl;
839 ifun >> outname;
840 }
841
842 // calibration flag
843 if(strcmp(word.Data(),"CALFLAG")==0)
844 ifun >> calflag;
845
846 // cleaning level
847 if(strcmp(word.Data(),"CLEANLEVEL")==0)
848 {
849 ifun >> lcore;
850 ifun >> ltail;
851 }
852 }
853
854 pediter.Reset();
855 caliter.Reset();
856 datiter.Reset();
857 TString pfile;
858
859 // Dump read values
860 cout << "************************************************" << endl;
861 cout << "* Datacards read from file " << filename << endl;
862 cout << "************************************************" << endl;
863 cout << "Pedestal file (s): " << endl;
864 while(!(pfile=pediter.Next()).IsNull())
865 cout << pfile << endl;
866 cout << "Calibration file (s): " << endl;
867 while(!(pfile=caliter.Next()).IsNull())
868 cout << pfile << endl;
869 cout << "Data file (s): " << endl;
870 while(!(pfile=datiter.Next()).IsNull())
871 cout << pfile << endl;
872 cout << "Maximum number of events: " << nmaxevents << endl;
873 cout << "Output file name: " << outname << endl;
874 cout << "Calibration flag: " << calflag << endl;
875 cout << "Cleaning level: ("<<lcore<<","<<ltail<<")" << endl;
876 cout << "***********" << endl << endl;
877
878 if(!pediter.GetNumEntries())
879 {
880 cout << "No pedestal file name specified" << endl;
881 return kFALSE;
882 }
883 if(!caliter.GetNumEntries())
884 {
885 cout << "No calibration file name specified" << endl;
886 return kFALSE;
887 }
888 if(!datiter.GetNumEntries())
889 {
890 cout << "No data file name specified" << endl;
891 return kFALSE;
892 }
893 if(!outname.Length())
894 {
895 cout << "No output file name specified" << endl;
896 return kFALSE;
897 }
898
899
900 return kTRUE;
901}
Note: See TracBrowser for help on using the repository browser.