source: trunk/MagicSoft/Mars/manalysis/MPadding.cc@ 1972

Last change on this file since 1972 was 1951, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 17.0 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Robert Wagner, 10/2002 <mailto:magicsoft@rwagner.de>
19! Author(s): Wolfgang Wittek, 01/2003 <mailto:wittek@mppmu.mpg.de>
20! Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2003
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28// //
29// MPadding //
30// (developped from MApplyPadding) //
31// //
32// This task applies padding to a given Sigmabar target value. //
33// The task checks whether the data stream it is applied to has to be //
34// padded or not and if so, Gaussian noise with the difference in Sigma //
35// is produced and added to the particular event. The number of photons, //
36// its error and the pedestal sigmas are altered. //
37// //
38// The padding has to be done before the image cleaning because the //
39// image cleaning depends on the pedestal sigmas. //
40// //
41// There are several ways of defining the sigmabar value to which the //
42// events are padded: //
43// //
44// 1) Set a fixed level (fFixedSigmabar) by calling 'SetTargetLevel'. //
45// //
46// 2) By calling 'SetDefiningHistogram', give a TH1D histogram //
47// (fHSigmabarMax) which defines the Sigmabar as a function of Theta. //
48// //
49// 3) By calling 'SetSigmaThetaHist', give a TH2D histogram //
50// (fHSigmaTheta) which contains the Sigmabar distribution for the //
51// different bins in Theta. For a given event, the sigmabar value to //
52// be used for the padding is thrown from this distribution. //
53// //
54// Workaround : //
55// If none of these options is specified then PreProcess will try to read //
56// in a propriety format ASCII database for the CT1 test. The name of //
57// this file is set by 'SetDatabaseFile'. From the data in this file a //
58// TH1D histogram (fHSigmabarMax) is generated. //
59// //
60// This implementation is still PRELIMINARY and requires some workarounds //
61// put in SPECIFICALLY FOR THE CT1 TESTS, since a database to access is //
62// missing. It is not the FINAL MAGIC VERSION. //
63// //
64// For random numbers gRandom is used. //
65// //
66/////////////////////////////////////////////////////////////////////////////
67#include "MPadding.h"
68
69#include <stdio.h>
70
71#include <TH1.h>
72#include <TH2.h>
73#include <TH3.h>
74#include <TRandom.h>
75#include <TCanvas.h>
76#include <TProfile.h>
77
78#include "MH.h"
79#include "MBinning.h"
80
81#include "MSigmabar.h"
82
83#include "MMcEvt.hxx"
84
85#include "MParList.h"
86
87#include "MLog.h"
88#include "MLogManip.h"
89
90#include "MGeomCam.h"
91
92#include "MCerPhotEvt.h"
93#include "MCerPhotPix.h"
94
95#include "MPedestalCam.h"
96#include "MPedestalPix.h"
97
98ClassImp(MPadding);
99
100// --------------------------------------------------------------------------
101//
102// Default constructor.
103//
104MPadding::MPadding(const char *name, const char *title)
105 : fRunType(0), fGroup(0), fFixedSigmabar(0), fHSigMaxAllocated(kFALSE), fHSigmabarMax(NULL), fHSigmaTheta(NULL)
106{
107 fName = name ? name : "MPadding";
108 fTitle = title ? title : "Task for the padding";
109
110 //--------------------------------------------------------------------
111 // plot pedestal sigmas for testing purposes
112 fHSigmaPedestal = new TH2D("SigPed", "Padded vs orig. sigma",
113 100, 0.0, 5.0, 100, 0.0, 5.0);
114 fHSigmaPedestal->SetXTitle("Orig. Pedestal sigma");
115 fHSigmaPedestal->SetYTitle("Padded Pedestal sigma");
116
117 // plot no.of photons (before vs. after padding) for testing purposes
118 fHPhotons = new TH2D("Photons", "Photons after vs.before padding",
119 100, -10.0, 90.0, 100, -10, 90);
120 fHPhotons->SetXTitle("No.of photons before padding");
121 fHPhotons->SetYTitle("No.of photons after padding");
122
123 // plot of added NSB
124 fHNSB = new TH1D("NSB", "Distribution of added NSB", 100, -10.0, 10.0);
125 fHNSB->SetXTitle("No.of added NSB photons");
126 fHNSB->SetYTitle("No.of pixels");
127
128 fHSigmaOld = new TH2D;
129 fHSigmaOld->SetNameTitle("fHSigmaOld", "Sigma before padding");
130 fHSigmaOld->SetXTitle("Theta");
131 fHSigmaOld->SetYTitle("Sigma");
132
133 fHSigmaNew = new TH2D;
134 fHSigmaNew->SetNameTitle("fHSigmaNew", "Sigma after padding");
135 fHSigmaNew->SetXTitle("Theta");
136 fHSigmaNew->SetYTitle("Sigma");
137}
138
139// --------------------------------------------------------------------------
140//
141// Destructor.
142//
143MPadding::~MPadding()
144{
145 delete fHSigmaPedestal;
146 delete fHPhotons;
147 delete fHNSB;
148 delete fHSigmaOld;
149 delete fHSigmaNew;
150 if (fHSigMaxAllocated && fHSigmabarMax)
151 delete fHSigmabarMax;
152}
153
154// --------------------------------------------------------------------------
155//
156// You can provide a TH1D* histogram containing
157// - the target Sigmabar in bins of theta.
158// Be sure to use the same binning as for the analysis
159//
160Bool_t MPadding::SetDefiningHistogram(TH1D *histo)
161{
162 if (fHSigmabarMax)
163 {
164 *fLog << warn << "MPadding - SigmabarMax already set.";
165 return kFALSE;
166 }
167
168 fHSigmabarMax = histo;
169
170 fFixedSigmabar = 0;
171 fHSigmaTheta = NULL;
172
173 *fLog << inf << "MPadding - Use Defining Histogram.";
174
175 return kTRUE;
176}
177
178// --------------------------------------------------------------------------
179//
180// You can provide a TH2D* histogram containing
181// - the Sigmabar distribution in bins of theta.
182//
183Bool_t MPadding::SetSigmaThetaHist(TH2D *histo)
184{
185 if (fHSigmaTheta)
186 {
187 *fLog << warn << "MPadding - SigmaTheta already set.";
188 return kFALSE;
189 }
190
191 fHSigmaTheta = histo;
192
193 fFixedSigmabar = 0;
194 if (fHSigMaxAllocated)
195 {
196 fHSigMaxAllocated = kFALSE;
197 delete fHSigmabarMax;
198 }
199 fHSigmabarMax = NULL;
200
201 *fLog << inf << "MPadding - Use Sigma Theta Histogram.";
202
203 return kTRUE;
204}
205
206// --------------------------------------------------------------------------
207//
208//
209void MPadding::SetTargetLevel(Double_t sigmabar)
210{
211 fFixedSigmabar = sigmabar;
212
213 fHSigmaTheta = NULL;
214 if (fHSigMaxAllocated)
215 {
216 fHSigMaxAllocated = kFALSE;
217 delete fHSigmabarMax;
218 }
219 fHSigmabarMax = NULL;
220
221 *fLog << inf << "MPadding - Use fixed sigmabar: fFixedSigmabar = ";
222 *fLog << fFixedSigmabar << endl;
223}
224
225// --------------------------------------------------------------------------
226//
227// check if MEvtHeader exists in the Parameter list already.
228// if not create one and add them to the list
229//
230Bool_t MPadding::PreProcess(MParList *pList)
231{
232 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
233 if (!fMcEvt)
234 {
235 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
236 return kFALSE;
237 }
238
239 fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
240 if (!fPed)
241 {
242 *fLog << err << dbginf << "MPedestalCam not found... aborting." << endl;
243 return kFALSE;
244 }
245
246 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
247 if (!fCam)
248 {
249 *fLog << err << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
250 return kFALSE;
251 }
252
253 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
254 if (!fEvt)
255 {
256 *fLog << err << dbginf << "MCerPhotEvt not found... aborting." << endl;
257 return kFALSE;
258 }
259
260 fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
261 if (!fSigmabar)
262 {
263 *fLog << err << dbginf << "MSigmabar not found... aborting." << endl;
264 return kFALSE;
265 }
266
267 // Get Theta Binning
268 MBinning* binstheta = (MBinning*)pList->FindObject("BinningTheta");
269 if (!binstheta)
270 {
271 *fLog << err << dbginf << "BinningTheta not found... aborting." << endl;
272 return kFALSE;
273 }
274
275 // Get Sigma Binning
276 MBinning* binssigma = (MBinning*)pList->FindObject("BinningSigmabar");
277 if (!binssigma)
278 {
279 *fLog << err << dbginf << "BinningSigmabar not found... aborting." << endl;
280 return kFALSE;
281 }
282
283 MH::SetBinning(fHSigmaOld, binstheta, binssigma);
284 MH::SetBinning(fHSigmaNew, binstheta, binssigma);
285
286 //************************************************************************
287 // Create fSigmabarMax histogram
288 // (only if no fixed Sigmabar target value and no histogram have been
289 // provided)
290 //
291 if (fFixedSigmabar==0 && !fHSigmabarMax && !fHSigmaTheta)
292 {
293 *fLog << inf << "MPadding - Creating fSigmabarMax histogram: ";
294 *fLog << "fFixedSigmabar=" << fFixedSigmabar << ", ";
295 *fLog << "fHSigmabarMax = " << fHSigmabarMax << endl;
296
297 // FIXME: Not deleted
298 fHSigmabarMax = new TH1D;
299 fHSigmabarMax->SetNameTitle("fHSigmabarMax", "Sigmabarmax for this analysis");
300
301 fHSigMaxAllocated = kTRUE;
302
303 MH::SetBinning(fHSigmabarMax, binstheta);
304
305 // -------------------------------------------------
306 // read in SigmabarParams
307 // workaround--proprietary file format--CT1test only BEGIN
308 // -------------------------------------------------
309
310 FILE *f=fopen(fDatabaseFilename, "r");
311 if(!f) {
312 *fLog << err << dbginf << "Database file '" << fDatabaseFilename;
313 *fLog << "' was not found (specified by MPadding::SetDatabaseFile) ...aborting." << endl;
314 return kFALSE;
315 }
316
317 TAxis &axe = *fHSigmabarMax->GetXaxis();
318
319 char line[80];
320 while ( fgets(line, sizeof(line), f) != NULL) {
321 if (line[0]=='#')
322 continue;
323
324 Float_t sigmabarMin, sigmabarMax, thetaMin, thetaMax, ra, dec2;
325 Int_t type, group, mjd, nr;
326
327 sscanf(line,"%d %d %f %f %d %d %f %f %f %f",
328 &type, &group, &ra, &dec2, &mjd, &nr,
329 &sigmabarMin, &sigmabarMax, &thetaMin, &thetaMax);
330
331 if (group!=fGroup && type!=1) //selected ON group or OFF
332 continue;
333
334 const Int_t from = axe.FindFixBin(thetaMin);
335 const Int_t to = axe.FindFixBin(thetaMax);
336
337 // find out which bin(s) we have to look at
338 for (Int_t i=from; i<to+1; i++)
339 if (sigmabarMax > fHSigmabarMax->GetBinContent(i))
340 fHSigmabarMax->SetBinContent(i, sigmabarMax);
341 }//while
342
343 } //fFixedSigmabar
344 //************************************************************************
345
346 if (!fHSigmabarMax && !fHSigmaTheta && fFixedSigmabar==0)
347 {
348 *fLog << err << "ERROR: Sigmabar for padding not defined... aborting." << endl;
349 return kFALSE;
350 }
351
352 return kTRUE;
353}
354
355Double_t MPadding::CalcOtherSig(const Double_t mySig, const Double_t theta) const
356{
357 //
358 // Get sigmabar which we have to pad to
359 //
360 const TAxis &axe = *fHSigmabarMax->GetXaxis();
361 const Int_t binnum = axe.FindFixBin(theta);
362 const Bool_t inrange = theta>=axe.GetXmin() && theta<=axe.GetXmax();
363
364 if ((fHSigmabarMax || fHSigmaTheta) && !inrange)
365 {
366 *fLog << err << dbginf;
367 *fLog << "Theta of current event is beyond the limits, Theta = ";
368 *fLog << theta << " ...skipping." <<endl;
369 return -1;
370 }
371
372
373 //
374 // get target sigma for the current Theta from the histogram fHSigmabarMax
375 //
376 if (fHSigmabarMax != NULL)
377 return fHSigmabarMax->GetBinContent(binnum);
378
379 //
380 // for the current Theta,
381 // generate a sigma according to the histogram fHSigmaTheta
382 //
383 if (fHSigmaTheta != NULL)
384 {
385 Double_t otherSig = -1;
386
387 TH1D* fHSigma = fHSigmaTheta->ProjectionY("", binnum, binnum, "");
388
389 if (fHSigma->GetEntries()>0)
390 otherSig = fHSigma->GetRandom();
391
392 delete fHSigma;
393
394 return otherSig;
395 }
396
397 //
398 // use a fixed target sigma
399 //
400 return fFixedSigmabar;
401}
402
403// --------------------------------------------------------------------------
404//
405// Do the padding (mySig ==> otherSig)
406//
407Bool_t MPadding::Padding(const Double_t quadraticDiff, const Double_t theta)
408{
409 const UInt_t npix = fEvt->GetNumPixels();
410
411 // pad only pixels - which are used (before image cleaning)
412 // - and for which the no.of photons is != 0.0
413 //
414 for (UInt_t i=0; i<npix; i++)
415 {
416 MCerPhotPix &pix = (*fEvt)[i];
417 if ( !pix.IsPixelUsed() )
418 continue;
419/*
420 if ( pix.GetNumPhotons() == 0)
421 {
422 *fLog << "MPadding::Process(); no.of photons is 0 for used pixel"
423 << endl;
424 continue;
425 }
426*/
427 const Double_t area = fCam->GetPixRatio(pix.GetPixId());
428
429 // add additional NSB to the number of photons
430 const Double_t NSB = sqrt(quadraticDiff*area) * gRandom->Gaus(0, 1);
431 const Double_t oldphotons = pix.GetNumPhotons();
432 const Double_t newphotons = oldphotons + NSB;
433 pix.SetNumPhotons( newphotons );
434
435 fHNSB->Fill( NSB/sqrt(area) );
436 fHPhotons->Fill( newphotons/sqrt(area), oldphotons/sqrt(area) );
437
438 // error: add sigma of padded noise quadratically
439 const Double_t olderror = pix.GetErrorPhot();
440 const Double_t newerror = sqrt( olderror*olderror + quadraticDiff*area );
441 pix.SetErrorPhot( newerror );
442
443 MPedestalPix &ppix = (*fPed)[i];
444
445 ppix.SetMeanRms(0);
446
447 const Double_t oldsigma = ppix.GetMeanRms();
448 const Double_t newsigma = sqrt( oldsigma*oldsigma + quadraticDiff*area );
449 ppix.SetMeanRms( newsigma );
450
451 fHSigmaPedestal->Fill( oldsigma, newsigma );
452 fHSigmaOld->Fill( theta, oldsigma );
453 fHSigmaNew->Fill( theta, newsigma );
454 } //for
455
456 return kTRUE;
457}
458
459
460// --------------------------------------------------------------------------
461//
462// Calculate Sigmabar for current event
463// Then apply padding
464//
465// 1) have current event
466// 2) get sigmabar(theta)
467// 3) pad event
468//
469Bool_t MPadding::Process()
470{
471 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
472
473 //
474 // Calculate sigmabar of event
475 //
476 Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
477
478 //$$$$$$$$$$$$$$$$$$$$$$$$$$
479 mySig = 0.0; // FIXME?
480 //$$$$$$$$$$$$$$$$$$$$$$$$$$
481
482 const Double_t otherSig = CalcOtherSig(mySig, theta);
483
484 // Skip event if target sigma is zero
485 if (otherSig<=0)
486 return kCONTINUE;
487
488 // Determine quadratic difference other-mine
489 const Double_t quadraticDiff = otherSig*otherSig - mySig*mySig;
490
491 if (quadraticDiff < 0) {
492 *fLog << err << "ERROR - MPadding: Event has higher Sigmabar=" << mySig;
493 *fLog << " than Sigmabarmax=" << otherSig << " @ Theta =" << theta;
494 *fLog << " ...skipping." << endl;
495 return kCONTINUE; //skip
496 }
497
498 if (quadraticDiff == 0)
499 return kTRUE; //no padding necessary.
500
501 //
502 // quadratic difference is > 0, do the padding;
503 //
504 Padding(quadraticDiff, theta);
505
506 // Calculate Sigmabar again and crosscheck
507 //mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
508
509 return kTRUE;
510}
511
512// --------------------------------------------------------------------------
513//
514// Draws some histograms if IsGraphicalOutputEnabled
515//
516Bool_t MPadding::PostProcess()
517{
518 if (!IsGraphicalOutputEnabled())
519 return kTRUE;
520
521 TCanvas &c = *MH::MakeDefCanvas("Padding", "", 600, 900);
522 c.Divide(2,3);
523 gROOT->SetSelectedPad(NULL);
524
525 if (fHSigmabarMax != NULL)
526 {
527 c.cd(1);
528 fHSigmabarMax->DrawClone();
529 }
530 else if (fHSigmaTheta != NULL)
531 {
532 c.cd(1);
533 fHSigmaTheta->DrawClone();
534 }
535
536 c.cd(3);
537 fHSigmaPedestal->DrawClone();
538
539 c.cd(5);
540 fHPhotons->DrawClone();
541
542 c.cd(2);
543 fHNSB->DrawClone();
544
545 c.cd(4);
546 fHSigmaOld->DrawClone();
547
548 c.cd(6);
549 fHSigmaNew->DrawClone();
550
551 return kTRUE;
552}
Note: See TracBrowser for help on using the repository browser.