source: branches/MarsISDCBranchBasedOn17887/mtemp/PaddingExample.C@ 18815

Last change on this file since 18815 was 5509, checked in by wittek, 21 years ago
*** empty log message ***
File size: 19.2 KB
Line 
1
2void InitBinnings(MParList *plist)
3{
4 gLog << "InitBinnings" << endl;
5
6 //--------------------------------------------
7 MBinning *binse = new MBinning("BinningE");
8 //binse->SetEdgesLog(30, 1.0e2, 1.0e5);
9
10 //This is Daniel's binning in energy:
11 binse->SetEdgesLog(14, 296.296, 86497.6);
12 plist->AddToList(binse);
13
14 //--------------------------------------------
15
16 MBinning *binssize = new MBinning("BinningSize");
17 binssize->SetEdgesLog(50, 10, 1.0e5);
18 plist->AddToList(binssize);
19
20 MBinning *binsdistc = new MBinning("BinningDist");
21 binsdistc->SetEdges(50, 0, 1.4);
22 plist->AddToList(binsdistc);
23
24 MBinning *binswidth = new MBinning("BinningWidth");
25 binswidth->SetEdges(50, 0, 1.0);
26 plist->AddToList(binswidth);
27
28 MBinning *binslength = new MBinning("BinningLength");
29 binslength->SetEdges(50, 0, 1.0);
30 plist->AddToList(binslength);
31
32 MBinning *binsalpha = new MBinning("BinningAlpha");
33 binsalpha->SetEdges(100, -100, 100);
34 plist->AddToList(binsalpha);
35
36 MBinning *binsasym = new MBinning("BinningAsym");
37 binsasym->SetEdges(50, -1.5, 1.5);
38 plist->AddToList(binsasym);
39
40 MBinning *binsm3l = new MBinning("BinningM3Long");
41 binsm3l->SetEdges(50, -1.5, 1.5);
42 plist->AddToList(binsm3l);
43
44 MBinning *binsm3t = new MBinning("BinningM3Trans");
45 binsm3t->SetEdges(50, -1.5, 1.5);
46 plist->AddToList(binsm3t);
47
48
49 //.....
50 MBinning *binsb = new MBinning("BinningSigmabar");
51 binsb->SetEdges( 100, 0.0, 50.0);
52 plist->AddToList(binsb);
53
54 MBinning *binssig = new MBinning("BinningSigma");
55 binssig->SetEdges( 100, 0.0, 120.0);
56 plist->AddToList(binssig);
57
58 MBinning *binth = new MBinning("BinningTheta");
59 // this is Daniel's binning in theta
60 //Double_t yedge[8] =
61 // {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99};
62 // this is our binning
63 //Double_t yedge[9] =
64 // {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
65 //TArrayD yed;
66 //yed.Set(9,yedge);
67 //binth->SetEdges(yed);
68
69 //binth->SetEdges( 1, 0.0, 90.0);
70 binth->SetEdgesCos( 10, 0.0, 90.0);
71 plist->AddToList(binth);
72
73 //MBinning *bincosth = new MBinning("BinningCosTheta");
74 //Double_t yedge[9] =
75 // {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
76 //Double_t zedge[9];
77 //for (Int_t i=0; i<9; i++)
78 //{
79 //zedge[8-i] = cos(yedge[i]/kRad2Deg);
80 //}
81 //TArrayD zed;
82 //zed.Set(9,zedge);
83 //bincosth->SetEdges(zed);
84 //plist->AddToList(bincosth);
85
86 MBinning *binsdiff = new MBinning("BinningDiffsigma2");
87 binsdiff->SetEdges(100, -500.0, 1500.0);
88 plist->AddToList(binsdiff);
89
90 // robert ----------------------------------------------
91 MBinning *binsalphaf = new MBinning("BinningAlphaFlux");
92 binsalphaf->SetEdges(100, -100, 100);
93 plist->AddToList(binsalphaf);
94
95 MBinning *binsdifftime = new MBinning("BinningTimeDiff");
96 binsdifftime->SetEdges(50, 0, 10);
97 plist->AddToList(binsdifftime);
98
99 MBinning *binstime = new MBinning("BinningTime");
100 binstime->SetEdges(50, 44500, 61000);
101 plist->AddToList(binstime);
102}
103
104
105void DeleteBinnings(MParList *plist)
106{
107 gLog << "DeleteBinnings" << endl;
108
109 TObject *bin;
110
111 //--------------------------------------------
112 bin = plist->FindObject("BinningE");
113 if (bin) delete bin;
114
115 //--------------------------------------------
116
117 bin = plist->FindObject("BinningSize");
118 if (bin) delete bin;
119
120 bin = plist->FindObject("BinningDist");
121 if (bin) delete bin;
122
123 bin = plist->FindObject("BinningWidth");
124 if (bin) delete bin;
125
126 bin = plist->FindObject("BinningLength");
127 if (bin) delete bin;
128
129 bin = plist->FindObject("BinningAlpha");
130 if (bin) delete bin;
131
132 bin = plist->FindObject("BinningAsym");
133 if (bin) delete bin;
134
135 bin = plist->FindObject("BinningM3Long");
136 if (bin) delete bin;
137
138 bin = plist->FindObject("BinningM3Trans");
139 if (bin) delete bin;
140
141 //.....
142 bin = plist->FindObject("BinningSigmabar");
143 if (bin) delete bin;
144
145 bin = plist->FindObject("BinningSigma");
146 if (bin) delete bin;
147
148 bin = plist->FindObject("BinningTheta");
149 if (bin) delete bin;
150
151 bin = plist->FindObject("BinningCosTheta");
152 if (bin) delete bin;
153
154 bin = plist->FindObject("BinningDiffsigma2");
155 if (bin) delete bin;
156
157
158 // robert ----------------------------------------------
159 bin = plist->FindObject("BinningAlphaFlux");
160 if (bin) delete bin;
161
162 bin = plist->FindObject("BinningTimeDiff");
163 if (bin) delete bin;
164
165 bin = plist->FindObject("BinningTime");
166 if (bin) delete bin;
167}
168
169
170
171//************************************************************************
172void PaddingExample()
173{
174 gLog.SetNoColors();
175
176 if (gRandom)
177 delete gRandom;
178 gRandom = new TRandom3(0);
179
180 //-----------------------------------------------
181 TString tag = "080000";
182
183 const char *offfile = "*OffCrab11*";
184
185 const char *onfile = "*CrabNebula*";
186
187 const char *mcfile = "calibrated_MCdata2";
188
189 //-----------------------------------------------
190
191 // path for input for Mars
192
193 if (tag == "080000")
194 {
195 TString inPathON = "/.magic/magicserv01/MAGIC/calibrateddata/2004_09_21/";
196 TString inPathOFF = "/.magic/magicserv01/MAGIC/calibrateddata/2004_09_21/";
197 TString inPathMC = "/.magic/data21a/mase/Mars/Mars041103/DataCalibUV/";
198
199 }
200
201 // path for output from Mars
202 TString outPath = "/.magic/data21a/wittek/PaddedData/";
203// outPath += tag;
204// outPath += "/";
205
206 //-----------------------------------------------
207
208
209 //************************************************************************
210
211 // Job A :
212 // - produce MHSigmaTheta plots for ON, OFF and MC data
213 // - write out (or read in) these MHSigmaTheta plots
214 // - read ON (or OFF or MC) data
215 // - pad the events;
216 // - write root file for ON (or OFF or MC) data (ON1.root, ...);
217
218 Bool_t JobA = kTRUE;
219 Bool_t GPadON = kTRUE; // \ generate Pad histograms
220 Bool_t GPadOFF = kFALSE; // | and write them onto disk
221 Bool_t GPadMC = kFALSE; // /
222 Bool_t Merge = kFALSE; // read the Pad histograms, merge them
223 // and write them onto disk
224 Bool_t Wout = kFALSE; // \ read in merged padding histograms and
225 // | write out root file of padded data
226 // / (ON1.root, OFF1.root or MC1.root)
227 TString typeInput("ON");
228 //TString typeInput("OFF");
229 //TString typeInput("MC");
230
231
232
233
234 //************************************************************************
235
236
237 //---------------------------------------------------------------------
238 // Job A
239 //=========
240
241 // - produce the histograms "sigmabar versus Theta", etc.
242 // for ON, OFF and MC data (to be used for the padding)
243 //
244 // - write root file of padded ON (OFF, MC) events (ON1.root, ...)
245 // (after the standard cuts, before the g/h separation)
246
247
248 if (JobA)
249 {
250 gLog << "=====================================================" << endl;
251 gLog << "Macro ONOFFAnalysis : Start of Job A" << endl;
252 gLog << "" << endl;
253 gLog << "Macro ONOFFAnalysis : JobA, GPadON, GPadOFF, GPadMC, Merge, Wout = "
254 << (JobA ? "kTRUE" : "kFALSE") << ", "
255 << (GPadON ? "kTRUE" : "kFALSE") << ", "
256 << (GPadOFF ? "kTRUE" : "kFALSE") << ", "
257 << (GPadMC ? "kTRUE" : "kFALSE") << ", "
258 << (Merge ? "kTRUE" : "kFALSE") << ", "
259 << (Wout ? "kTRUE" : "kFALSE") << endl;
260
261 //--------------------------------------------------
262
263
264 // name of MPedPhotCam container
265 TString fNamePedPhotCam("MPedPhotCam");
266
267
268 //************************************************************
269 // generate histograms to be used in the padding
270 //
271 // read ON, OFF and MC data files
272 // generate (or read in) the padding histograms for ON, OFF and MC data
273 //
274
275 MPad pad;
276 pad.SetName("MPad");
277
278 //--------------------------------------------------
279 // names of ON and OFF files to be read
280 // for generating the histograms to be used in the padding
281 TString fileON = inPathON;
282 fileON += onfile;
283 fileON += ".root";
284
285 TString fileOFF = inPathOFF;
286 fileOFF += offfile;
287 fileOFF += ".root";
288
289 TString fileMC = inPathMC;
290 fileMC += mcfile;
291 fileMC += ".root";
292
293 //--------------------------------------------------
294 // name of files to contain the paddding histograms of ON, OFF and MC data
295 TString NamePadON(outPath);
296 NamePadON += "PadON";
297 NamePadON += ".root";
298
299 TString NamePadOFF(outPath);
300 NamePadOFF += "PadOFF";
301 NamePadOFF += ".root";
302
303 TString NamePadMC(outPath);
304 NamePadMC += "PadMC";
305 NamePadMC += ".root";
306
307 // name of file to conatin the merged histograms for the padding
308 TString outNameSigTh = outPath;
309 outNameSigTh += "SigmaTheta";
310 outNameSigTh += ".root";
311
312 //--------------------------------------------------
313
314 if (GPadON || GPadOFF || GPadMC)
315 {
316 // generate the padding histograms
317 gLog << "=====================================================" << endl;
318 gLog << "Start generating the padding histograms" << endl;
319
320
321 gLog << "fileON, fileOFF, fileMC = " << fileON << ", "
322 << fileOFF << ", " << fileMC << endl;
323
324
325
326 //--------------------------------------------------
327 MMakePadHistograms makepad;
328 makepad.SetMaxEvents(-1);
329 makepad.SetNamePedPhotCam(fNamePedPhotCam);
330 makepad.SetPedestalLevel(2.0);
331 makepad.SetUseInterpolation(kTRUE);
332 makepad.SetProcessPedestal(kTRUE);
333 makepad.SetProcessTime(kFALSE);
334
335 //-----------------------------------------
336 // ON events
337
338 if (GPadON)
339 {
340 makepad.SetDataType("ON");
341 makepad.SetNameInputFile(fileON);
342 makepad.SetNameOutputFile(NamePadON);
343 makepad.MakeHistograms();
344 }
345
346 //-----------------------------------------
347 // OFF events
348
349 if (GPadOFF)
350 {
351 makepad.SetDataType("OFF");
352 makepad.SetNameInputFile(fileOFF);
353 makepad.SetNameOutputFile(NamePadOFF);
354 makepad.MakeHistograms();
355 }
356
357 //-----------------------------------------
358 // MC events
359
360 if (GPadMC)
361 {
362 makepad.SetDataType("MC");
363 makepad.SetNameInputFile(fileMC);
364 makepad.SetNameOutputFile(NamePadMC);
365 makepad.MakeHistograms();
366 }
367
368 //-----------------------------------------
369
370
371 gLog << "" << endl;
372 gLog << "End of generating the padding histograms" << endl;
373 gLog << "=====================================================" << endl;
374 }
375
376 //************************************************************
377
378 if (Merge)
379 {
380 gLog << "=====================================================" << endl;
381 gLog << "Start of merging the padding histograms" << endl;
382 gLog << "" << endl;
383
384 //pad.MergeONOFFMC(NamePadON, NamePadOFF, NamePadMC, outNameSigTh);
385 //pad.MergeONOFFMC(NamePadON, "", NamePadMC, outNameSigTh);
386 pad.MergeONOFFMC(NamePadON, NamePadOFF, "", outNameSigTh);
387 //pad.MergeONOFFMC("", NamePadOFF, NamePadMC, outNameSigTh);
388
389 gLog << "" << endl;
390 gLog << "End of merging the padding histograms" << endl;
391 gLog << "=====================================================" << endl;
392 }
393 // end of Merge
394
395
396
397 //************************************************************
398
399 if (Wout)
400 {
401 // read the target padding histograms ---------------------------
402 pad.ReadPaddingDist(outNameSigTh);
403
404
405 gLog << "=====================================================" << endl;
406 gLog << "Start the padding" << endl;
407
408 //--------------------------------------------------
409 // type of data to be padded
410 gLog << "typeInput = " << typeInput << endl;
411
412 //-------------------------------------------
413 // name of input root file
414 if (typeInput == "ON")
415 TString filenamein(fileON);
416 else if (typeInput == "OFF")
417 TString filenamein(fileOFF);
418 else if (typeInput == "MC")
419 TString filenamein(fileMC);
420 gLog << "data to be padded : " << filenamein << endl;
421
422 // name of output root file
423 TString outNameImage = outPath;
424 outNameImage += "Hillas";
425 outNameImage += typeInput;
426 outNameImage += "1.root";
427 gLog << "padded data to be written onto : " << outNameImage << endl;
428
429 //-----------------------------------------------------------
430 pad.SetDataType(typeInput);
431 pad.SetNamePedPhotCam(fNamePedPhotCam);
432
433 MTaskList tliston;
434 MParList pliston;
435
436 // geometry is needed in MHHillas... classes
437 MGeomCam *fGeom =
438 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
439
440 //-------------------------------------------
441 // create the tasks which should be executed
442 //
443
444 MReadReports read;
445 read.AddTree("Drive");
446 read.AddTree("Events","MTime.",kTRUE);
447 static_cast<MRead&>(read).AddFile(filenamein);
448
449 read.AddToBranchList("MReportDrive.*");
450 read.AddToBranchList("MRawEvtHeader.*");
451
452
453// MReadMarsFile read("Events", filenamein);
454// read.DisableAutoScheme();
455
456 MGeomApply apply;
457
458 // a way to find out whether one is dealing with MC :
459 //MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5); // MC
460 //f1.SetName("Select MC");
461 //MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5); // data
462 //f2.SetName("Select Data");
463
464
465
466 MBadPixelsCalc badcalc;
467 badcalc.SetNamePedPhotContainer(fNamePedPhotCam);
468 badcalc.SetPedestalLevel(2.0);
469 badcalc.SetName("BadCalc");
470
471 MBadPixelsTreat badtreat;
472 badtreat.SetNamePedPhotCam(fNamePedPhotCam);
473 badtreat.SetUseInterpolation(kTRUE);
474 badtreat.SetProcessPedestal(kTRUE);
475 badtreat.SetProcessTimes(kFALSE);
476 badtreat.SetName("BadTreat");
477
478 MPointingPosCalc mpointcalc;
479
480 MFSelBasic selbasic;
481 selbasic.SetCuts(20.0, 0.0, 25.8419, 90.0, 126.0);
482
483 MContinue contbasic(&selbasic);
484 contbasic.SetName("SelBasic");
485
486 MHBadPixels bad("BadPixels");
487 bad.SetNamePedPhotCam(fNamePedPhotCam);
488 MFillH fillbad("BadPixels[MHBadPixels]", "MBadPixelsCam");
489 fillbad.SetName("FillBad");
490
491 MHSigmaTheta sigth("SigmaTheta");
492 sigth.SetNamePedPhotCam(fNamePedPhotCam);
493 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "");
494 fillsigtheta.SetName("FillSigTheta");
495
496
497 MImgCleanStd clean(3.0, 2.5);
498 clean.SetMethod(MImgCleanStd::kScaled);
499 clean.SetCleanRings(1);
500 clean.SetName("Clean");
501 clean.SetNamePedPhotCam(fNamePedPhotCam);
502
503 // calculation of image parameters ---------------------
504 TString fHilName = "MHillas";
505 TString fHilNameExt = "MHillasExt";
506 TString fHilNameSrc = "MHillasSrc";
507 TString fImgParName = "MNewImagePar";
508
509 MHillasCalc hcalc;
510 hcalc.SetNameHillas(fHilName);
511 hcalc.SetNameHillasExt(fHilNameExt);
512 hcalc.SetNameNewImagePar(fImgParName);
513
514 MHillasSrcCalc hsrccalc("MSrcPosCam", fHilNameSrc);
515 hsrccalc.SetInput(fHilName);
516
517 MFillH hfill1("MHHillas", fHilName);
518 hfill1.SetName("HHillas");
519
520 MFillH hfill2("MHStarMap", fHilName);
521 hfill2.SetName("HStarMap");
522
523 MFillH hfill3("MHHillasExt", fHilNameSrc);
524 hfill3.SetName("HHillasExt");
525
526 MFillH hfill4("MHHillasSrc", fHilNameSrc);
527 hfill4.SetName("HHillasSrc");
528
529 MFillH hfill5("MHNewImagePar", fImgParName);
530 hfill5.SetName("HNewImagePar");
531 // --------------------------------------------------
532
533 MFSelStandard selstandard(fHilNameSrc);
534 selstandard.SetHillasName(fHilName);
535 selstandard.SetImgParName(fImgParName);
536 selstandard.SetCuts(500, 2, 20, 0.0, 5.0, 0.001, 0.001);
537 MContinue contstandard(&selstandard);
538 contstandard.SetName("SelStandard");
539
540
541 /*
542 MWriteRootFile write(outNameImage);
543
544 write.AddContainer("MRawRunHeader", "RunHeaders");
545 //write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
546 //write.AddContainer("MTime", "Events");
547 write.AddContainer("MPointingPos", "Events");
548 write.AddContainer("MSrcPosCam", "Events");
549 write.AddContainer("MHillas", "Events");
550 write.AddContainer("MHillasExt", "Events");
551 write.AddContainer("MHillasSrc", "Events");
552 write.AddContainer("MNewImagePar", "Events");
553 write.AddContainer("MTime", "Events");
554 write.AddContainer("MImagePar", "Events");
555
556 write.AddContainer("MReportDrive", "Drive", kFALSE);
557 write.AddContainer("MTimeDrive", "Drive", kFALSE); //new
558 */
559
560 //*****************************
561 // entries in MParList
562
563 pliston.AddToList(&tliston);
564 InitBinnings(&pliston);
565 pliston.AddToList(&bad);
566 pliston.AddToList(&sigth);
567
568 //*****************************
569 // entries in MTaskList
570
571 tliston.AddToList(&read);
572 //tliston.AddToList(&f1);
573 //tliston.AddToList(&f2);
574 tliston.AddToList(&apply);
575 tliston.AddToList(&mpointcalc, "Events");
576
577// tliston.AddToList(&badcalc); done in callisto
578// tliston.AddToList(&badtreat); done in callisto
579
580 tliston.AddToList(&contbasic, "Events");
581 //tliston.AddToList(&pad, "Events");
582
583
584 tliston.AddToList(&fillbad, "Events");
585 tliston.AddToList(&fillsigtheta, "Events");
586 tliston.AddToList(&clean, "Events");
587
588 tliston.AddToList(&hcalc, "Events");
589 tliston.AddToList(&hsrccalc, "Events");
590
591 tliston.AddToList(&hfill1, "Events");
592 tliston.AddToList(&hfill2, "Events");
593 tliston.AddToList(&hfill3, "Events");
594 tliston.AddToList(&hfill4, "Events");
595 tliston.AddToList(&hfill5, "Events");
596
597 tliston.AddToList(&contstandard, "Events");
598 //tliston.AddToList(&write);
599
600
601 //*****************************
602
603 //-------------------------------------------
604 // Execute event loop
605 //
606 MProgressBar bar;
607 MEvtLoop evtloop;
608 evtloop.SetParList(&pliston);
609 //evtloop.ReadEnv(env, "", printEnv);
610 evtloop.SetProgressBar(&bar);
611 // evtloop.Write();
612
613 Int_t maxevents = -1;
614 //Int_t maxevents = 2000;
615 if ( !evtloop.Eventloop(maxevents) )
616 return;
617
618 tliston.PrintStatistics(0, kTRUE);
619
620
621 //-------------------------------------------
622 // Display the histograms
623
624 pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone();
625 pliston.FindObject("BadPixels", "MHBadPixels")->DrawClone();
626
627 pliston.FindObject("MHHillas")->DrawClone();
628 pliston.FindObject("MHHillasExt")->DrawClone();
629 pliston.FindObject("MHHillasSrc")->DrawClone();
630 pliston.FindObject("MHNewImagePar")->DrawClone();
631 pliston.FindObject("MHStarMap")->DrawClone();
632
633 //DeleteBinnings(&pliston);
634
635 gLog << "End of padding" << endl;
636 gLog << "=====================================================" << endl;
637 }
638
639
640 gLog << "Macro ONOFFAnalysis : End of Job A" << endl;
641 gLog << "===================================================" << endl;
642 }
643
644}
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
Note: See TracBrowser for help on using the repository browser.