source: trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc@ 3570

Last change on this file since 3570 was 3374, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 7.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): Abelardo Moralejo, 12/2003 <mailto:moralejo@pd.infn.it>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MMcCalibrationCalc
28//
29// Input Containers:
30// MRawRunHeader
31// MMcFadcHeader
32// MHillas
33// MNewImagePar
34// MMcEvt
35//
36// Output Containers:
37// MCalibrationChargeCam
38//
39/////////////////////////////////////////////////////////////////////////////
40#include "MMcCalibrationCalc.h"
41
42#include <TH1.h>
43
44#include "MLog.h"
45#include "MLogManip.h"
46
47#include "MParList.h"
48
49#include "MCalibrationChargePix.h"
50#include "MCalibrationChargeCam.h"
51
52#include "MGeomCam.h"
53#include "MRawRunHeader.h"
54
55#include "MHillas.h"
56#include "MNewImagePar.h"
57
58#include "MMcEvt.hxx"
59#include "MMcFadcHeader.hxx"
60
61ClassImp(MMcCalibrationCalc);
62
63using namespace std;
64
65MMcCalibrationCalc::MMcCalibrationCalc(const char *name, const char *title)
66{
67 fName = name ? name : "MMcCalibrationCalc";
68 fTitle = title ? title : "Calculate and write conversion factors into MCalibrationCam Container";
69
70 fHistRatio = new TH1F(AddSerialNumber("HistRatio"), "log10(fPassPhotCone/fSize)", 1500, -3., 3.);
71 fHistRatio->SetXTitle("log_{10}(fPassPhotCone / fSize) [phot/ADC count]");
72}
73
74// --------------------------------------------------------------------------
75//
76// Check for the run type. Return kTRUE if it is a MC run or if there
77// is no MC run header (old camera files) kFALSE in case of a different
78// run type
79//
80Bool_t MMcCalibrationCalc::CheckRunType(MParList *pList) const
81{
82 const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
83 if (!run)
84 {
85 *fLog << warn << "Warning - cannot check file type, MRawRunHeader not found." << endl;
86 return kTRUE;
87 }
88
89 return run->IsMonteCarloRun();
90}
91
92// --------------------------------------------------------------------------
93//
94// Make sure, that there is an MCalibrationCam Object in the Parameter List.
95//
96Int_t MMcCalibrationCalc::PreProcess(MParList *pList)
97{
98 fHistRatio->Reset();
99 fADC2Phot = 0;
100
101 fCalCam = (MCalibrationChargeCam*) pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
102 if (!fCalCam)
103 {
104 *fLog << err << AddSerialNumber("MCalibrationChargeCam") << "not found... aborting." << endl;
105 return kFALSE;
106 }
107
108 fHillas = (MHillas*) pList->FindObject(AddSerialNumber("MHillas"));
109 if ( !fHillas)
110 {
111 *fLog << err << AddSerialNumber("MHillas") << "not found... aborting." << endl;
112 return kFALSE;
113 }
114
115 fNew = (MNewImagePar*)pList->FindObject(AddSerialNumber("MNewImagePar"));
116 if (!fNew)
117 {
118 *fLog << err << AddSerialNumber("MNewImagePar") << "not found... aborting." << endl;
119 return kFALSE;
120 }
121
122 fMcEvt = (MMcEvt*) pList->FindObject(AddSerialNumber("MMcEvt"));
123 if (!fMcEvt)
124 {
125 *fLog << err << AddSerialNumber("MMcEvt") << "not found... aborting." << endl;
126 return kFALSE;
127 }
128
129 return kTRUE;
130}
131
132// --------------------------------------------------------------------------
133//
134// Check for the runtype.
135// Search for MGeomCam and MMcFadcHeader.
136//
137Bool_t MMcCalibrationCalc::ReInit(MParList *pList)
138{
139 //
140 // If it is no MC file display error and exit
141 //
142 if (!CheckRunType(pList))
143 {
144 *fLog << err << "MMcCalibrationCalc can only be used with MC files... aborting." << endl;
145 return kFALSE;
146 }
147
148 //
149 // Now check the existence of all necessary containers.
150 //
151 fGeom = (MGeomCam*) pList->FindObject(AddSerialNumber("MGeomCam"));
152 if (!fGeom)
153 {
154 *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl;
155 return kFALSE;
156 }
157
158 fHeaderFadc = (MMcFadcHeader*)pList->FindObject(AddSerialNumber("MMcFadcHeader"));
159 if (!fHeaderFadc)
160 {
161 *fLog << err << AddSerialNumber("MMcFadcHeader") << " not found... aborting." << endl;
162 return kFALSE;
163 }
164
165 for (UInt_t ipix = 0; ipix < fGeom->GetNumPixels(); ipix++)
166 {
167 if (fHeaderFadc->GetPedestalRmsHigh(ipix) > 0 ||
168 fHeaderFadc->GetPedestalRmsLow(ipix) > 0 )
169 {
170 *fLog << err << "Trying to calibrate the data using a Camera file produced with added noise." << endl;
171 *fLog << "Please use a noiseless file for calibration... aborting." << endl << endl;
172 return kFALSE;
173 }
174 }
175
176 return kTRUE;
177}
178
179
180// --------------------------------------------------------------------------
181//
182// Obtain average ratio of photons in camera to image Size.
183//
184Int_t MMcCalibrationCalc::Process()
185{
186 //
187 // Exclude events with some saturated pixel
188 //
189 if (fNew->GetNumSaturatedPixels()>0)
190 return kTRUE;
191
192 //
193 // Exclude events with low Size (larger fluctuations)
194 // FIXME? The present cut (1000 "inner-pixel-counts") is somehow
195 // arbitrary. Might it be optimized?
196 //
197 if (fHillas->GetSize()<1000)
198 return kTRUE;
199
200 fADC2Phot += fMcEvt->GetPassPhotCone()/fHillas->GetSize();
201
202 fHistRatio->Fill(TMath::Log10(fMcEvt->GetPassPhotCone()/fHillas->GetSize()));
203
204 return kTRUE;
205}
206
207// --------------------------------------------------------------------------
208//
209// Fill the MCalibrationCam object
210//
211Int_t MMcCalibrationCalc::PostProcess()
212{
213 const Stat_t n = fHistRatio->GetEntries();
214 if (n<1)
215 {
216 *fLog << err << "No events read... aborting." << endl;
217 return kFALSE;
218 }
219
220 fADC2Phot /= n;
221
222 //
223 // For the calibration we no longer use the mean,
224 // but the peak of the distribution:
225 //
226 const Int_t reach = 2;
227
228 Stat_t summax = 0;
229 Int_t mode = 0;
230
231 // FIXME: Is this necessary? We could use GetMaximumBin instead..
232 for (Int_t ibin = 1+reach; ibin <= fHistRatio->GetNbinsX()-reach; ibin++)
233 {
234 const Stat_t sum = fHistRatio->Integral(ibin-reach, ibin+reach);
235
236 if (sum <= summax)
237 continue;
238
239 summax = sum;
240 mode = ibin;
241 }
242
243 fADC2Phot = TMath::Power(10, fHistRatio->GetBinCenter(mode));
244
245 const Int_t num = fCalCam->GetSize();
246 for (int i=0; i<num; i++)
247 {
248 MCalibrationChargePix &calpix = (*fCalCam)[i];
249
250 const Float_t factor = fADC2Phot*calpix.GetMeanConversionBlindPixelMethod();
251
252 calpix.SetConversionBlindPixelMethod(factor, 0., 0.);
253 }
254
255 return kTRUE;
256}
Note: See TracBrowser for help on using the repository browser.