source: trunk/MagicSoft/Mars/mtemp/meth/macros/analysis.C@ 5258

Last change on this file since 5258 was 4354, checked in by commichau, 22 years ago
*** empty log message ***
File size: 14.6 KB
Line 
1
2/////////////////////////////////////////////////////////////////////////////
3// //
4// analysis.C //
5// //
6// Author(s): S.C. Commichau, 6/2004 //
7// //
8/////////////////////////////////////////////////////////////////////////////
9
10
11// This macro calculates source independent and dependent image parameter
12// distributions from data and writes them into a file called hillas.root.
13// To look at the distributions use the macro analysis_read.C
14// Up to now no derotation and mispointing correction is taken into account!
15
16const TString inpathcal = "/ihp/MAGIC/Period016/rootdata2/2004_04_22/";
17//const TString inpath = "/ihp/MAGIC/Period016/rootdata2/2004_04_23/";
18const TString inpath = "/ihp/MAGIC/Period016/rootdata2/2004_04_23/";
19
20
21const TString dname = "/ihp/MAGIC/Period016/rootdata2/2004_04_23/20040423_2375*_D_OffMrk421-6_E.root";
22
23//const TString dname = "/ihp/MAGIC/Period016/rootdata2/2004_04_23/20040423_2343*_D_Mrk421_E.root";
24
25
26const Int_t pedrun = 23745;//23427; // Pedestal file
27const Int_t calrun1 = 23205; // Calibration files
28const Int_t calrun2 = 0;
29const Int_t datrun1 = 23429; // Data files
30
31const Bool_t usedisplay = kFALSE;
32
33// Range of slices that are taken into account
34const Byte_t hifirst = 1;
35const Byte_t hilast = 14;
36const Byte_t lofirst = 3;
37const Byte_t lolast = 14;
38
39
40void analysis()
41{
42
43 /*******************************************
44 * Define some parameters for the analysis *
45 *******************************************/
46
47 // Change the following line toggle pedestal display
48 Bool_t DrawPedestals = kFALSE;//kTRUE;
49
50 // Set calibration mode
51 // kDefault, kNone, kDummy, kFfactor, kPinDiode,...
52 MCalibrate::CalibrationMode_t CalMode=MCalibrate::kDefault;
53
54 // Cleaning: Core and Tailcut
55 TArrayF clvl(2);
56 clvl[0] = 3.5;
57 clvl[1] = 3.0;
58
59 // Set bad pixels
60 const Short_t x[16] = {330,395,329,396,389,
61 323,388,322,384,385,
62 386,387,321,320,319,
63 394};
64
65 // Set output filename
66 TString oname = "hillas.root";
67
68
69 /***************************************
70 * Define general tasks and containers *
71 ***************************************/
72
73 MCalibrationQECam qecam;
74 MBadPixelsCam badcam;
75 MGeomCamMagic geomcam;
76 MGeomApply geomapl;
77 MExtractFixedWindow extractor;
78 extractor.SetRange(hifirst, hilast, lofirst, lolast);
79 //MExtractFixedWindowPeakSearch extractor;
80 //extractor.SetWindows(8,8);
81 //MExtractSlidingWindow extractor;
82
83
84 // Iterate over runs
85 MRunIter pruns;
86 MRunIter cruns;
87 MRunIter druns;
88
89 MReadMarsFile dread("Events", dname);
90 dread.DisableAutoScheme();
91
92 // Pedestal runs
93 pruns.AddRun(pedrun,inpath);
94
95 // Calibration runs
96 if (calrun2==0)
97 cruns.AddRun(calrun1,inpathcal);
98 else
99 cruns.AddRuns(calrun1,calrun2,inpath);
100
101 // Data
102/* if (1)//datrun2==0)
103 druns.AddRun(datrun1,inpath);
104 else
105 druns.AddRuns(datrun1,datrun2,inpath);//,datrun3,datrun4,datrun5,datrun6,inpath);
106*/
107
108 // Look for the pulser colors
109 MCalibrationCam::PulserColor_t color;
110
111 if (calrun1 < 20000)
112 color = MCalibrationCam::kCT1;
113 else
114 color = FindColor((MDirIter*)&cruns);
115
116 if (color==MCalibrationCam::kNONE)
117 {
118 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
119
120 while (1)
121 {
122 timer.TurnOn();
123 TString input = Getline("Couldn't find the pulser colour, type 'q' to exit, "
124 "'green', 'blue', 'uv' or 'ct1' to go on: ");
125 timer.TurnOff();
126
127 if (input=="q\n")
128 return;
129
130 if (input=="green\n")
131 {
132 color = MCalibrationCam::kGREEN;
133 break;
134 }
135 if (input=="blue\n")
136 {
137 color = MCalibrationCam::kBLUE;
138 break;
139 }
140 if (input=="uv\n" || input=="UV\n")
141 {
142 color = MCalibrationCam::kUV;
143 break;
144 }
145 if (input=="ct1\n" || input=="CT1\n")
146 {
147 color = MCalibrationCam::kCT1;
148 break;
149 }
150 }
151 }
152
153
154 /**********************************
155 * 1st Loop: Pedestal Computation *
156 **********************************/
157
158 // If pixels have to be excluded from the beginning, read
159 // an ASCII-file with the corresponding pixel numbers (MBadPixelsCam)
160 badcam.AsciiRead("badpixels.txt");
161
162 MJPedestal pedloop;
163 pedloop.SetInput(&pruns);
164 pedloop.SetExtractor(&extractor);
165 pedloop.SetBadPixels(badcam);
166
167 if (!pedloop.Process())
168 return;
169
170 // Print the content of MPedestalCam
171 MPedestalCam *pedcam = (MPedestalCam*)pedloop.GetPedestalCam();
172 pedcam->Print();
173
174 // Draw Pedestals and RMS into a camera - ADC counts
175 if (DrawPedestals)
176 {
177 MHCamera *disp0 = new MHCamera(geomcam, "Pedestals", "");
178 MHCamera *disp1 = new MHCamera(geomcam, "PedestalRMS", "");
179
180 disp0->SetMinimum(0);
181 disp1->SetMinimum(0);
182
183 TCanvas *c1 = new TCanvas("Pedestals", "Pedestals", 50, 50, 400, 400);
184 TCanvas *c2 = new TCanvas("RMS", "RMS", 460, 50, 400, 400);
185
186 c1->SetBorderMode(0);
187 c1->SetFillColor(0);
188
189 c2->SetBorderMode(0);
190 c2->SetFillColor(0);
191
192 gStyle->SetOptStat(1101);
193 gStyle->SetStatColor(0);
194 gStyle->SetStatBorderSize(1);
195
196 gPad->SetBorderMode(0);
197
198 disp0.SetPrettyPalette();
199 disp1.SetPrettyPalette();
200
201 c1->cd();
202 disp0.Draw();
203
204 gPad->cd(1);
205
206 disp0.SetCamContent(*pedcam,0);
207 disp0.SetCamError(*pedcam,1);
208
209 c2->cd();
210 disp1.Draw();
211
212 disp1.SetCamContent(*pedcam,2);
213 disp1.SetCamError(*pedcam,3);
214 }
215
216
217 /*************************
218 * 2nd Loop: Calibration *
219 *************************/
220
221 MJCalibration calloop;
222 //calloop.SetColor(color);
223 calloop.SetInput(&cruns);
224 calloop.SetExtractor(&extractor);
225 //calloop.SetQECam(qecam);
226 calloop.SetBadPixels(pedloop.GetBadPixels());
227
228 // Choose the arrival time Extractor:
229 //MExtractTimeHighestIntegral timeext;
230 //MExtractTimeSpline timeext;
231 // Set Ranges or Windows
232 //timeext.SetRange(0,14,6,14);
233 //calloop.SetTimeExtractor(&timeext);
234
235 // To calibrate times enable the following line
236 //Bool_t UseTimes = kFALSE;
237 // Bool_t calibrateTimes = kFALSE;
238 //calloop.SetRelTimeCalibration(UseTimes);
239
240 if (!calloop.Process(pedloop.GetPedestalCam()))
241 return;
242
243 /**********************************
244 * 3rd Loop: Pedestal Calibration *
245 **********************************/
246
247 // Create a Parameter and a Task List
248 MParList plist3;
249 MTaskList tlist3;
250 plist3.AddToList(&tlist3);
251
252 // Parameter Containers
253 MCerPhotEvt evtphot;
254 MPedPhotCam pedphot;
255 MExtractedSignalCam sigcam;
256
257 // Add also some containers from previous loops
258 plist3.AddToList(&geomcam);
259 plist3.AddToList(&pedloop.GetPedestalCam());
260
261 plist3.AddToList(&calloop.GetCalibrationCam());
262 plist3.AddToList(&calloop.GetQECam()); // I have to check this!
263 //plist3.AddToList(&calloop.GetRelTimeCam());
264 //plist3.AddToList(&calloop.GetBadPixels());
265 plist3.AddToList(&badcam);
266 plist3.AddToList(&sigcam);
267 plist3.AddToList(&evtphot);
268 plist3.AddToList(&pedphot);
269
270
271 // Tasks
272 MReadMarsFile read3("Events");
273 static_cast<MRead&>(read3).AddFiles(pruns);
274 read3.DisableAutoScheme();
275
276 //MExtractSignal extsig;
277 //extsig.SetRange(hifirst, hilast, lofirst, lolast);
278 MCalibrate calphot;
279 calphot.SetCalibrationMode(CalMode);
280 MPedPhotCalc pedphotcalc;
281
282 tlist3.AddToList(&read3);
283 tlist3.AddToList(&geomapl);
284 tlist3.AddToList(&extractor);
285 //tlist3.AddToList(&extsig);
286 tlist3.AddToList(&calphot);
287 tlist3.AddToList(&pedphotcalc);
288
289
290 // Create and setup the eventloop
291 MProgressBar bar3;
292 bar3.SetWindowName("Calibrating Pedestals...");
293
294 MEvtLoop evtloop3;
295 evtloop3.SetProgressBar(&bar3);
296 evtloop3.SetParList(&plist3);
297
298 if (!evtloop3.Eventloop())
299 return;
300
301 tlist3.PrintStatistics();
302
303
304 /****************************
305 * 4th Loop: Photonize Data *
306 ****************************/
307
308 // Create a Parameter and a Task List
309 MParList plist4;
310 MTaskList tlist4;
311 plist4.AddToList(&tlist4);
312
313 // Parameter Containers
314 MHillas hillas;
315 //MDCA dca;
316 MHillasSrc hillassrc;
317 MSrcPosCam source;
318 source.SetX(0.0*189/0.6);
319 source.SetY(0.0*189/0.6);
320 MRawRunHeader runhead;
321 MArrivalTimeCam timecam;
322
323 // False Source stuff
324 MTime time;
325 MObservatory obsv;
326 MPointingPos poin;
327
328 hillassrc.SetSrcPos(&source);
329
330
331 // Add also some containers from previous loops
332 //plist4.AddToList(&timecam);
333 plist4.AddToList(&geomcam);
334 plist4.AddToList(&pedloop.GetPedestalCam());
335 plist4.AddToList(&calloop.GetCalibrationCam());
336 //plist4.AddToList(&badcam);
337 plist4.AddToList(&calloop.GetQECam());
338 //plist4.AddToList(&calloop.GetRelTimeCam());
339 //plist4.AddToList(&calloop.GetBadPixels());
340 plist4.AddToList(&source);
341 plist4.AddToList(&hillas);
342 plist4.AddToList(&hillassrc);
343 //plist4.AddToList(&dca);
344 plist4.AddToList(&evtphot);
345 plist4.AddToList(&pedphot);
346
347 plist4.AddToList(&time);
348 plist4.AddToList(&obsv);
349 plist4.AddToList(&poin);
350
351
352 //MArrivalTime times;
353 //plist4.AddToList(&times);
354
355 // Tasks
356/* MReadMarsFile read4("Events");
357 static_cast<MRead&>(read4).AddFiles(druns);
358 read4.DisableAutoScheme();
359*/
360 // Set bad pixels
361 MBlindPixelCalc blind;
362
363 const TArrayS bp(16,(Short_t*)x);
364 blind.SetPixelIndices(bp);
365
366 MImgCleanStd clean(clvl[0],clvl[1]);
367 MHillasCalc hilcalc;
368 MHillasSrcCalc hilsrc;
369 //MDCACalc dcacalc;
370
371
372
373
374 MWriteRootFile write(oname,"RECREATE"); //Alternative options: NEW, UPDATE
375
376 //write.AddContainer("MDCA" , "Parameters");
377 write.AddContainer("MHillas" , "Parameters");
378 write.AddContainer("MHillasSrc" , "Parameters");
379 write.AddContainer("MHillasExt" , "Parameters");
380 write.AddContainer("MNewImagePar" , "Parameters");
381 write.AddContainer("MRawEvtHeader" , "Parameters");
382 write.AddContainer("MRawRunHeader" , "Parameters");
383 write.AddContainer("MConcentration" , "Parameters");
384 write.AddContainer("MSrcPosCam" , "Parameters");
385
386
387 // False Source stuff
388 MHFalseSource fshist;
389 fshist.SetAlphaCut(12.5);
390 //fshist.SetBgMean(35);
391 MFillH hfill(&fshist, "MHillas");
392 //MFillH hfill("MHFalseSource", "MHillas");
393
394 //MFalseSrcCalc falsecalc;
395 //falsecalc.SetAlphaCut(12.5);
396 //falsecalc.SetDistCuts(.25,1.0);
397 //falsecalc.SetSizeCut(1000);
398 //falsecalc.SetWidthCuts(.04,.14);
399 //falsecalc.SetLengthCuts(.14,.26);
400 //falsecalc.SetOutput("ON.root");//outname);
401 //falsecalc.SetRefPoint(0,0);
402
403 //tlist4.AddToList(&read4);
404 tlist4.AddToList(&dread);
405 tlist4.AddToList(&geomapl);
406 tlist4.AddToList(&extractor);
407 //tlist4.AddToList(&blind);
408 //tlist4.AddToList(&timecalc);
409 tlist4.AddToList(&calphot);
410 //tlist4.AddToList(&timecal);
411 tlist4.AddToList(&clean);
412 tlist4.AddToList(&hilcalc);
413 //tlist4.AddToList(&dcacalc);
414 tlist4.AddToList(&hilsrc);
415// tlist4.AddToList(&falsecalc);
416 //tlist4.AddToList(&hfill);
417 tlist4.AddToList(&write);
418
419 // Create and setup the eventloop
420 MProgressBar bar4;
421 bar4.SetWindowName("Analyzing Data...");
422
423 MEvtLoop evtloop4;
424 evtloop4.SetProgressBar(&bar4);
425 evtloop4.SetParList(&plist4);
426
427 if (!evtloop4.Eventloop())
428 return;
429
430 tlist4.PrintStatistics();
431
432 //fshist.Draw();
433
434
435}
436
437
438Bool_t HandleInput()
439{
440 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
441
442 while (1)
443 {
444 // While reading the input process GUI events asynchronously
445 timer.TurnOn();
446 TString input = Getline("Type 'q' to exit, <return> to go on: ");
447 timer.TurnOff();
448
449 if (input=="q\n")
450 return kFALSE;
451
452 if (input=="\n")
453 return kTRUE;
454 };
455
456 return kFALSE;
457}
458
459
460// Find the correct pulser color, done by M. Gaug
461MCalibrationCam::PulserColor_t FindColor(MDirIter* run)
462{
463
464 MCalibrationCam::PulserColor_t col = MCalibrationCam::kNONE;
465
466 TString filenames;
467
468 while (!(filenames=run->Next()).IsNull())
469 {
470
471 filenames.ToLower();
472
473 if (filenames.Contains("green"))
474 if (col == MCalibrationCam::kNONE)
475 {
476 cout << "Found colour: Green in " << filenames << endl;
477 col = MCalibrationCam::kGREEN;
478 }
479 else if (col != MCalibrationCam::kGREEN)
480 {
481 cout << "Different colour found in " << filenames << "... abort" << endl;
482 return MCalibrationCam::kNONE;
483 }
484
485 if (filenames.Contains("blue"))
486 if (col == MCalibrationCam::kNONE)
487 {
488 cout << "Found colour: Blue in " << filenames << endl;
489 col = MCalibrationCam::kBLUE;
490 }
491 else if (col != MCalibrationCam::kBLUE)
492 {
493 cout << "Different colour found in " << filenames << "... abort" << endl;
494 return MCalibrationCam::kNONE;
495 }
496
497 if (filenames.Contains("uv"))
498 if (col == MCalibrationCam::kNONE)
499 {
500 cout << "Found colour: Uv in " << filenames << endl;
501 col = MCalibrationCam::kUV;
502 }
503 else if (col != MCalibrationCam::kUV)
504 {
505 cout << "Different colour found in " << filenames << "... abort" << endl;
506 return MCalibrationCam::kNONE;
507 }
508
509 if (filenames.Contains("ct1"))
510 if (col == MCalibrationCam::kNONE)
511 {
512 cout << "Found colour: Ct1 in " << filenames << endl;
513 col = MCalibrationCam::kCT1;
514 }
515 else if (col != MCalibrationCam::kCT1)
516 {
517 cout << "Different colour found in " << filenames << "... abort" << endl;
518 return MCalibrationCam::kNONE;
519 }
520
521 }
522
523
524 if (col == MCalibrationCam::kNONE)
525 cout << "No colour found in filenames of runs: " << ((MRunIter*)run)->GetRunsAsString()
526 << "... abort" << endl;
527
528 return col;
529}
Note: See TracBrowser for help on using the repository browser.