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

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