source: branches/Mars_use_drstimefiles/manalysis/MPadding.cc@ 18459

Last change on this file since 18459 was 2798, checked in by wittek, 21 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 "MPedPhotCam.h"
96#include "MPedPhotPix.h"
97
98ClassImp(MPadding);
99
100using namespace std;
101
102// --------------------------------------------------------------------------
103//
104// Default constructor.
105//
106MPadding::MPadding(const char *name, const char *title)
107 : fRunType(0), fGroup(0), fFixedSigmabar(0), fHSigMaxAllocated(kFALSE), fHSigmabarMax(NULL), fHSigmaTheta(NULL)
108{
109 fName = name ? name : "MPadding";
110 fTitle = title ? title : "Task for the padding";
111
112 //--------------------------------------------------------------------
113 // plot pedestal sigmas for testing purposes
114 fHSigmaPedestal = new TH2D("SigPed", "Padded vs orig. sigma",
115 100, 0.0, 5.0, 100, 0.0, 5.0);
116 fHSigmaPedestal->SetXTitle("Orig. Pedestal sigma");
117 fHSigmaPedestal->SetYTitle("Padded Pedestal sigma");
118
119 // plot no.of photons (before vs. after padding) for testing purposes
120 fHPhotons = new TH2D("Photons", "Photons after vs.before padding",
121 100, -10.0, 90.0, 100, -10, 90);
122 fHPhotons->SetXTitle("No.of photons before padding");
123 fHPhotons->SetYTitle("No.of photons after padding");
124
125 // plot of added NSB
126 fHNSB = new TH1D("NSB", "Distribution of added NSB", 100, -10.0, 10.0);
127 fHNSB->SetXTitle("No.of added NSB photons");
128 fHNSB->SetYTitle("No.of pixels");
129
130 fHSigmaOld = new TH2D;
131 fHSigmaOld->SetNameTitle("fHSigmaOld", "Sigma before padding");
132 fHSigmaOld->SetXTitle("Theta");
133 fHSigmaOld->SetYTitle("Sigma");
134
135 fHSigmaNew = new TH2D;
136 fHSigmaNew->SetNameTitle("fHSigmaNew", "Sigma after padding");
137 fHSigmaNew->SetXTitle("Theta");
138 fHSigmaNew->SetYTitle("Sigma");
139}
140
141// --------------------------------------------------------------------------
142//
143// Destructor.
144//
145MPadding::~MPadding()
146{
147 delete fHSigmaPedestal;
148 delete fHPhotons;
149 delete fHNSB;
150 delete fHSigmaOld;
151 delete fHSigmaNew;
152 if (fHSigMaxAllocated && fHSigmabarMax)
153 delete fHSigmabarMax;
154}
155
156// --------------------------------------------------------------------------
157//
158// You can provide a TH1D* histogram containing
159// - the target Sigmabar in bins of theta.
160// Be sure to use the same binning as for the analysis
161//
162Bool_t MPadding::SetDefiningHistogram(TH1D *histo)
163{
164 if (fHSigmabarMax)
165 {
166 *fLog << warn << "MPadding - SigmabarMax already set.";
167 return kFALSE;
168 }
169
170 fHSigmabarMax = histo;
171
172 fFixedSigmabar = 0;
173 fHSigmaTheta = NULL;
174
175 *fLog << inf << "MPadding - Use Defining Histogram.";
176
177 return kTRUE;
178}
179
180// --------------------------------------------------------------------------
181//
182// You can provide a TH2D* histogram containing
183// - the Sigmabar distribution in bins of theta.
184//
185Bool_t MPadding::SetSigmaThetaHist(TH2D *histo)
186{
187 if (fHSigmaTheta)
188 {
189 *fLog << warn << "MPadding - SigmaTheta already set.";
190 return kFALSE;
191 }
192
193 fHSigmaTheta = histo;
194
195 fFixedSigmabar = 0;
196 if (fHSigMaxAllocated)
197 {
198 fHSigMaxAllocated = kFALSE;
199 delete fHSigmabarMax;
200 }
201 fHSigmabarMax = NULL;
202
203 *fLog << inf << "MPadding - Use Sigma Theta Histogram.";
204
205 return kTRUE;
206}
207
208// --------------------------------------------------------------------------
209//
210//
211void MPadding::SetTargetLevel(Double_t sigmabar)
212{
213 fFixedSigmabar = sigmabar;
214
215 fHSigmaTheta = NULL;
216 if (fHSigMaxAllocated)
217 {
218 fHSigMaxAllocated = kFALSE;
219 delete fHSigmabarMax;
220 }
221 fHSigmabarMax = NULL;
222
223 *fLog << inf << "MPadding - Use fixed sigmabar: fFixedSigmabar = ";
224 *fLog << fFixedSigmabar << endl;
225}
226
227// --------------------------------------------------------------------------
228//
229// check if MEvtHeader exists in the Parameter list already.
230// if not create one and add them to the list
231//
232Int_t MPadding::PreProcess(MParList *pList)
233{
234 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
235 if (!fMcEvt)
236 {
237 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
238 return kFALSE;
239 }
240
241 fPed = (MPedPhotCam*)pList->FindObject("MPedPhotCam");
242 if (!fPed)
243 {
244 *fLog << err << dbginf << "MPedPhotCam not found... aborting." << endl;
245 return kFALSE;
246 }
247
248 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
249 if (!fCam)
250 {
251 *fLog << err << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
252 return kFALSE;
253 }
254
255 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
256 if (!fEvt)
257 {
258 *fLog << err << dbginf << "MCerPhotEvt not found... aborting." << endl;
259 return kFALSE;
260 }
261
262 fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
263 if (!fSigmabar)
264 {
265 *fLog << err << dbginf << "MSigmabar not found... aborting." << endl;
266 return kFALSE;
267 }
268
269 // Get Theta Binning
270 MBinning* binstheta = (MBinning*)pList->FindObject("BinningTheta");
271 if (!binstheta)
272 {
273 *fLog << err << dbginf << "BinningTheta not found... aborting." << endl;
274 return kFALSE;
275 }
276
277 // Get Sigma Binning
278 MBinning* binssigma = (MBinning*)pList->FindObject("BinningSigmabar");
279 if (!binssigma)
280 {
281 *fLog << err << dbginf << "BinningSigmabar not found... aborting." << endl;
282 return kFALSE;
283 }
284
285 MH::SetBinning(fHSigmaOld, binstheta, binssigma);
286 MH::SetBinning(fHSigmaNew, binstheta, binssigma);
287
288 //************************************************************************
289 // Create fSigmabarMax histogram
290 // (only if no fixed Sigmabar target value and no histogram have been
291 // provided)
292 //
293 if (fFixedSigmabar==0 && !fHSigmabarMax && !fHSigmaTheta)
294 {
295 *fLog << inf << "MPadding - Creating fSigmabarMax histogram: ";
296 *fLog << "fFixedSigmabar=" << fFixedSigmabar << ", ";
297 *fLog << "fHSigmabarMax = " << fHSigmabarMax << endl;
298
299 // FIXME: Not deleted
300 fHSigmabarMax = new TH1D;
301 fHSigmabarMax->SetNameTitle("fHSigmabarMax", "Sigmabarmax for this analysis");
302
303 fHSigMaxAllocated = kTRUE;
304
305 MH::SetBinning(fHSigmabarMax, binstheta);
306
307 // -------------------------------------------------
308 // read in SigmabarParams
309 // workaround--proprietary file format--CT1test only BEGIN
310 // -------------------------------------------------
311
312 FILE *f=fopen(fDatabaseFilename, "r");
313 if(!f) {
314 *fLog << err << dbginf << "Database file '" << fDatabaseFilename;
315 *fLog << "' was not found (specified by MPadding::SetDatabaseFile) ...aborting." << endl;
316 return kFALSE;
317 }
318
319 TAxis &axe = *fHSigmabarMax->GetXaxis();
320
321 char line[80];
322 while ( fgets(line, sizeof(line), f) != NULL) {
323 if (line[0]=='#')
324 continue;
325
326 Float_t sigmabarMin, sigmabarMax, thetaMin, thetaMax, ra, dec2;
327 Int_t type, group, mjd, nr;
328
329 sscanf(line,"%d %d %f %f %d %d %f %f %f %f",
330 &type, &group, &ra, &dec2, &mjd, &nr,
331 &sigmabarMin, &sigmabarMax, &thetaMin, &thetaMax);
332
333 if (group!=fGroup && type!=1) //selected ON group or OFF
334 continue;
335
336 const Int_t from = axe.FindFixBin(thetaMin);
337 const Int_t to = axe.FindFixBin(thetaMax);
338
339 // find out which bin(s) we have to look at
340 for (Int_t i=from; i<to+1; i++)
341 if (sigmabarMax > fHSigmabarMax->GetBinContent(i))
342 fHSigmabarMax->SetBinContent(i, sigmabarMax);
343 }//while
344
345 } //fFixedSigmabar
346 //************************************************************************
347
348 if (!fHSigmabarMax && !fHSigmaTheta && fFixedSigmabar==0)
349 {
350 *fLog << err << "ERROR: Sigmabar for padding not defined... aborting." << endl;
351 return kFALSE;
352 }
353
354 return kTRUE;
355}
356
357Double_t MPadding::CalcOtherSig(const Double_t mySig, const Double_t theta) const
358{
359 //
360 // Get sigmabar which we have to pad to
361 //
362 const TAxis &axe = *fHSigmabarMax->GetXaxis();
363 const Int_t binnum = axe.FindFixBin(theta);
364 const Bool_t inrange = theta>=axe.GetXmin() && theta<=axe.GetXmax();
365
366 if ((fHSigmabarMax || fHSigmaTheta) && !inrange)
367 {
368 *fLog << err << dbginf;
369 *fLog << "Theta of current event is beyond the limits, Theta = ";
370 *fLog << theta << " ...skipping." <<endl;
371 return -1;
372 }
373
374
375 //
376 // get target sigma for the current Theta from the histogram fHSigmabarMax
377 //
378 if (fHSigmabarMax != NULL)
379 return fHSigmabarMax->GetBinContent(binnum);
380
381 //
382 // for the current Theta,
383 // generate a sigma according to the histogram fHSigmaTheta
384 //
385 if (fHSigmaTheta != NULL)
386 {
387 Double_t otherSig = -1;
388
389 TH1D* fHSigma = fHSigmaTheta->ProjectionY("", binnum, binnum, "");
390
391 if (fHSigma->GetEntries()>0)
392 otherSig = fHSigma->GetRandom();
393
394 delete fHSigma;
395
396 return otherSig;
397 }
398
399 //
400 // use a fixed target sigma
401 //
402 return fFixedSigmabar;
403}
404
405// --------------------------------------------------------------------------
406//
407// Do the padding (mySig ==> otherSig)
408//
409Int_t MPadding::Padding(const Double_t quadraticDiff, const Double_t theta)
410{
411 const UInt_t npix = fEvt->GetNumPixels();
412
413 // pad only pixels - which are used (before image cleaning)
414 // - and for which the no.of photons is != 0.0
415 //
416 for (UInt_t i=0; i<npix; i++)
417 {
418 MCerPhotPix &pix = (*fEvt)[i];
419 if ( !pix.IsPixelUsed() )
420 continue;
421/*
422 if ( pix.GetNumPhotons() == 0)
423 {
424 *fLog << "MPadding::Process(); no.of photons is 0 for used pixel"
425 << endl;
426 continue;
427 }
428*/
429 const Double_t area = fCam->GetPixRatio(pix.GetPixId());
430
431 // add additional NSB to the number of photons
432 const Double_t NSB = sqrt(quadraticDiff*area) * gRandom->Gaus(0, 1);
433 const Double_t oldphotons = pix.GetNumPhotons();
434 const Double_t newphotons = oldphotons + NSB;
435 pix.SetNumPhotons( newphotons );
436
437 fHNSB->Fill( NSB/sqrt(area) );
438 fHPhotons->Fill( newphotons/sqrt(area), oldphotons/sqrt(area) );
439
440 // error: add sigma of padded noise quadratically
441 const Double_t olderror = pix.GetErrorPhot();
442 const Double_t newerror = sqrt( olderror*olderror + quadraticDiff*area );
443 pix.SetErrorPhot( newerror );
444
445 MPedPhotPix &ppix = (*fPed)[i];
446
447 ppix.SetRms(0);
448
449 const Double_t oldsigma = ppix.GetRms();
450 const Double_t newsigma = sqrt( oldsigma*oldsigma + quadraticDiff*area );
451 ppix.SetRms( newsigma );
452
453 fHSigmaPedestal->Fill( oldsigma, newsigma );
454 fHSigmaOld->Fill( theta, oldsigma );
455 fHSigmaNew->Fill( theta, newsigma );
456 } //for
457
458 return kTRUE;
459}
460
461
462// --------------------------------------------------------------------------
463//
464// Calculate Sigmabar for current event
465// Then apply padding
466//
467// 1) have current event
468// 2) get sigmabar(theta)
469// 3) pad event
470//
471Int_t MPadding::Process()
472{
473 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
474
475 //
476 // Calculate sigmabar of event
477 //
478 Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
479
480 //$$$$$$$$$$$$$$$$$$$$$$$$$$
481 mySig = 0.0; // FIXME?
482 //$$$$$$$$$$$$$$$$$$$$$$$$$$
483
484 const Double_t otherSig = CalcOtherSig(mySig, theta);
485
486 // Skip event if target sigma is zero
487 if (otherSig<=0)
488 return kCONTINUE;
489
490 // Determine quadratic difference other-mine
491 const Double_t quadraticDiff = otherSig*otherSig - mySig*mySig;
492
493 if (quadraticDiff < 0) {
494 *fLog << err << "ERROR - MPadding: Event has higher Sigmabar=" << mySig;
495 *fLog << " than Sigmabarmax=" << otherSig << " @ Theta =" << theta;
496 *fLog << " ...skipping." << endl;
497 return kCONTINUE; //skip
498 }
499
500 if (quadraticDiff == 0)
501 return kTRUE; //no padding necessary.
502
503 //
504 // quadratic difference is > 0, do the padding;
505 //
506 Padding(quadraticDiff, theta);
507
508 // Calculate Sigmabar again and crosscheck
509 //mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
510
511 return kTRUE;
512}
513
514// --------------------------------------------------------------------------
515//
516// Draws some histograms if IsGraphicalOutputEnabled
517//
518Bool_t MPadding::PostProcess()
519{
520 if (!IsGraphicalOutputEnabled())
521 return kTRUE;
522
523 TCanvas &c = *MH::MakeDefCanvas("Padding", "", 600, 900);
524 c.Divide(2,3);
525 gROOT->SetSelectedPad(NULL);
526
527 if (fHSigmabarMax != NULL)
528 {
529 c.cd(1);
530 fHSigmabarMax->DrawClone();
531 }
532 else if (fHSigmaTheta != NULL)
533 {
534 c.cd(1);
535 fHSigmaTheta->DrawClone();
536 }
537
538 c.cd(3);
539 fHSigmaPedestal->DrawClone();
540
541 c.cd(5);
542 fHPhotons->DrawClone();
543
544 c.cd(2);
545 fHNSB->DrawClone();
546
547 c.cd(4);
548 fHSigmaOld->DrawClone();
549
550 c.cd(6);
551 fHSigmaNew->DrawClone();
552
553 return kTRUE;
554}
Note: See TracBrowser for help on using the repository browser.