| 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): Thomas Bretz 12/2000 <mailto:tbretz@atsro.uni-wuerzburh.de>
|
|---|
| 19 | ! Author(s): Markus Gaug 02/2004 <mailto:markus@ifae.es>
|
|---|
| 20 | !
|
|---|
| 21 | ! Copyright: MAGIC Software Development, 2000-2004
|
|---|
| 22 | !
|
|---|
| 23 | !
|
|---|
| 24 | \* ======================================================================== */
|
|---|
| 25 |
|
|---|
| 26 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 27 | //
|
|---|
| 28 | // MHPedestalCam
|
|---|
| 29 | //
|
|---|
| 30 | // Fills the extracted pedestals of MExtractedSignalCam into
|
|---|
| 31 | // the MHCalibrationPix-classes MHPedestalPix for every:
|
|---|
| 32 | //
|
|---|
| 33 | // - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray
|
|---|
| 34 | // or MHCalibrationCam::fHiGainArray, respectively
|
|---|
| 35 | //
|
|---|
| 36 | // - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
|
|---|
| 37 | // stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
|
|---|
| 38 | // MHCalibrationCam::fAverageHiGainAreas
|
|---|
| 39 | //
|
|---|
| 40 | // - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
|
|---|
| 41 | // stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
|
|---|
| 42 | // and MHCalibrationCam::fAverageHiGainSectors
|
|---|
| 43 | //
|
|---|
| 44 | // Every pedestal is directly taken from MExtractedSignalCam, filled by the
|
|---|
| 45 | // appropriate extractor.
|
|---|
| 46 | //
|
|---|
| 47 | // The pedestals are filled into a histogram and an array, in order to perform
|
|---|
| 48 | // a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an
|
|---|
| 49 | // event-by-event basis and written into the corresponding average pixels.
|
|---|
| 50 | //
|
|---|
| 51 | // The histograms are fitted to a Gaussian, mean and sigma with its errors
|
|---|
| 52 | // and the fit probability are extracted. If none of these values are NaN's and
|
|---|
| 53 | // if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%),
|
|---|
| 54 | // the fit is declared valid.
|
|---|
| 55 | // Otherwise, the fit is repeated within ranges of the previous mean
|
|---|
| 56 | // +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
|
|---|
| 57 | // In case this does not make the fit valid, the histogram means and RMS's are
|
|---|
| 58 | // taken directly (see MHCalibrationPix::BypassFit()).
|
|---|
| 59 | //
|
|---|
| 60 | // Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas
|
|---|
| 61 | // from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup)
|
|---|
| 62 | //
|
|---|
| 63 | // The number of summed FADC slices is taken to re-normalize
|
|---|
| 64 | // the result to a pedestal per pixel with the formulas (see MHPedestalPix::Renorm()):
|
|---|
| 65 | // - Mean Pedestal / slice = Mean Pedestal / Number slices
|
|---|
| 66 | // - Mean Pedestal Error / slice = Mean Pedestal Error / Number slices
|
|---|
| 67 | // - Sigma Pedestal / slice = Sigma Pedestal / Sqrt (Number slices)
|
|---|
| 68 | // - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices)
|
|---|
| 69 | //
|
|---|
| 70 | // The class also fills arrays with the signal vs. event number, creates a fourier
|
|---|
| 71 | // spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the
|
|---|
| 72 | // projected fourier components follow an exponential distribution.
|
|---|
| 73 | //
|
|---|
| 74 | // This same procedure is performed for the average pixels.
|
|---|
| 75 | //
|
|---|
| 76 | // The following results are written into MPedestalCam:
|
|---|
| 77 | //
|
|---|
| 78 | // - MCalibrationPix::SetMean()
|
|---|
| 79 | // - MCalibrationPix::SetMeanErr()
|
|---|
| 80 | // - MCalibrationPix::SetSigma()
|
|---|
| 81 | // - MCalibrationPix::SetSigmaErr()
|
|---|
| 82 | // - MCalibrationPix::SetProb()
|
|---|
| 83 | // - MCalibrationPix::SetNumPickup()
|
|---|
| 84 | //
|
|---|
| 85 | // For all averaged areas, the fitted sigma is multiplied with the square root of
|
|---|
| 86 | // the number involved pixels in order to be able to compare it to the average of
|
|---|
| 87 | // sigmas in the camera.
|
|---|
| 88 | //
|
|---|
| 89 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 90 | #include "MHPedestalCam.h"
|
|---|
| 91 | #include "MHPedestalPix.h"
|
|---|
| 92 |
|
|---|
| 93 | #include "MBadPixelsCam.h"
|
|---|
| 94 | #include "MBadPixelsPix.h"
|
|---|
| 95 |
|
|---|
| 96 | #include "MLog.h"
|
|---|
| 97 | #include "MLogManip.h"
|
|---|
| 98 |
|
|---|
| 99 | #include "MParList.h"
|
|---|
| 100 |
|
|---|
| 101 | #include "MExtractedSignalCam.h"
|
|---|
| 102 | #include "MExtractedSignalPix.h"
|
|---|
| 103 |
|
|---|
| 104 | #include "MPedestalCam.h"
|
|---|
| 105 | #include "MPedestalPix.h"
|
|---|
| 106 | #include "MCalibrationPix.h"
|
|---|
| 107 |
|
|---|
| 108 | #include "MGeomCam.h"
|
|---|
| 109 | #include "MGeomPix.h"
|
|---|
| 110 |
|
|---|
| 111 | #include "MCalibrationPedCam.h"
|
|---|
| 112 |
|
|---|
| 113 | #include "TH1.h"
|
|---|
| 114 | #include "TH1F.h"
|
|---|
| 115 |
|
|---|
| 116 | #include <TOrdCollection.h>
|
|---|
| 117 |
|
|---|
| 118 | ClassImp(MHPedestalCam);
|
|---|
| 119 |
|
|---|
| 120 | using namespace std;
|
|---|
| 121 |
|
|---|
| 122 | const Int_t MHPedestalCam::fgNbins = 100;
|
|---|
| 123 | const Axis_t MHPedestalCam::fgFirst = -50.;
|
|---|
| 124 | const Axis_t MHPedestalCam::fgLast = 50.;
|
|---|
| 125 | const TString MHPedestalCam::gsHistName = "Pedestal";
|
|---|
| 126 | const TString MHPedestalCam::gsHistTitle = "Pedestal";
|
|---|
| 127 | const TString MHPedestalCam::gsHistXTitle = "Charge [FADC slices]";
|
|---|
| 128 | const TString MHPedestalCam::gsHistYTitle = "Nr. events";
|
|---|
| 129 | // --------------------------------------------------------------------------
|
|---|
| 130 | //
|
|---|
| 131 | // Default constructor.
|
|---|
| 132 | //
|
|---|
| 133 | // Initializes:
|
|---|
| 134 | // - fExtractHiGainSlices to 0.
|
|---|
| 135 | // - the event frequency to 1200 Hz.
|
|---|
| 136 | // - the fRenorm flag to kFALSE
|
|---|
| 137 | //
|
|---|
| 138 | // - fNbins to fgNbins
|
|---|
| 139 | // - fFirst to fgFirst
|
|---|
| 140 | // - fLast to fgLast
|
|---|
| 141 | //
|
|---|
| 142 | // - fHistName to gsHistName
|
|---|
| 143 | // - fHistTitle to gsHistTitle
|
|---|
| 144 | // - fHistXTitle to gsHistXTitle
|
|---|
| 145 | // - fHistYTitle to gsHistYTitle
|
|---|
| 146 | //
|
|---|
| 147 | MHPedestalCam::MHPedestalCam(const char *name, const char *title)
|
|---|
| 148 | : fNumEvents(0), fExtractHiGainSlices(0.)
|
|---|
| 149 | {
|
|---|
| 150 |
|
|---|
| 151 | fName = name ? name : "MHPedestalCam";
|
|---|
| 152 | fTitle = title ? title : "Class to fill the pedestal histograms";
|
|---|
| 153 |
|
|---|
| 154 | SetPulserFrequency(1200);
|
|---|
| 155 | SetRenorm(kTRUE);
|
|---|
| 156 | SetLoGain(kFALSE);
|
|---|
| 157 |
|
|---|
| 158 | SetNbins(fgNbins);
|
|---|
| 159 | SetFirst(fgFirst);
|
|---|
| 160 | SetLast (fgLast );
|
|---|
| 161 |
|
|---|
| 162 | SetHistName (gsHistName .Data());
|
|---|
| 163 | SetHistTitle (gsHistTitle .Data());
|
|---|
| 164 | SetHistXTitle(gsHistXTitle.Data());
|
|---|
| 165 | SetHistYTitle(gsHistYTitle.Data());
|
|---|
| 166 | }
|
|---|
| 167 |
|
|---|
| 168 | // --------------------------------------------------------------------------
|
|---|
| 169 | //
|
|---|
| 170 | // Calls MHCalibrationCam::ResetHists()
|
|---|
| 171 | //
|
|---|
| 172 | // Resets:
|
|---|
| 173 | // - fSum
|
|---|
| 174 | // - fSum2
|
|---|
| 175 | // - fAreaSum
|
|---|
| 176 | // - fAreaSum2
|
|---|
| 177 | // - fAreaVal
|
|---|
| 178 | // - fSectorSum
|
|---|
| 179 | // - fSectorSum2
|
|---|
| 180 | // - fSectorVal
|
|---|
| 181 | //
|
|---|
| 182 | void MHPedestalCam::ResetHists()
|
|---|
| 183 | {
|
|---|
| 184 |
|
|---|
| 185 | MHCalibrationCam::ResetHists();
|
|---|
| 186 |
|
|---|
| 187 | fNumEvents = 0;
|
|---|
| 188 |
|
|---|
| 189 | // If the size is yet set, set the size
|
|---|
| 190 | if (fSum.GetSize()>0)
|
|---|
| 191 | {
|
|---|
| 192 | // Reset contents of arrays.
|
|---|
| 193 | fSum.Reset();
|
|---|
| 194 | fSum2.Reset();
|
|---|
| 195 |
|
|---|
| 196 | fAreaSum. Reset();
|
|---|
| 197 | fAreaSum2.Reset();
|
|---|
| 198 | fAreaVal.Reset();
|
|---|
| 199 |
|
|---|
| 200 | fSectorSum. Reset();
|
|---|
| 201 | fSectorSum2.Reset();
|
|---|
| 202 | fSectorVal.Reset();
|
|---|
| 203 | }
|
|---|
| 204 | }
|
|---|
| 205 |
|
|---|
| 206 | // --------------------------------------------------------------------------
|
|---|
| 207 | //
|
|---|
| 208 | // Creates new MHCalibrationChargeCam only with the averaged areas:
|
|---|
| 209 | // the rest has to be retrieved directly, e.g. via:
|
|---|
| 210 | // MHPedestalCam *cam = MParList::FindObject("MHPedestalCam");
|
|---|
| 211 | // - cam->GetAverageSector(5).DrawClone();
|
|---|
| 212 | // - (*cam)[100].DrawClone()
|
|---|
| 213 | //
|
|---|
| 214 | TObject *MHPedestalCam::Clone(const char *) const
|
|---|
| 215 | {
|
|---|
| 216 |
|
|---|
| 217 | MHPedestalCam *cam = new MHPedestalCam();
|
|---|
| 218 |
|
|---|
| 219 | //
|
|---|
| 220 | // Copy the data members
|
|---|
| 221 | //
|
|---|
| 222 | cam->fRunNumbers = fRunNumbers;
|
|---|
| 223 | cam->fPulserFrequency = fPulserFrequency;
|
|---|
| 224 | cam->fFlags = fFlags;
|
|---|
| 225 | cam->fNbins = fNbins;
|
|---|
| 226 | cam->fFirst = fFirst;
|
|---|
| 227 | cam->fLast = fLast;
|
|---|
| 228 |
|
|---|
| 229 | if (!IsAverageing())
|
|---|
| 230 | return cam;
|
|---|
| 231 |
|
|---|
| 232 | const Int_t navhi = fAverageHiGainAreas->GetSize();
|
|---|
| 233 |
|
|---|
| 234 | for (int i=0; i<navhi; i++)
|
|---|
| 235 | cam->fAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i);
|
|---|
| 236 |
|
|---|
| 237 | return cam;
|
|---|
| 238 | }
|
|---|
| 239 |
|
|---|
| 240 | // --------------------------------------------------------------------------
|
|---|
| 241 | //
|
|---|
| 242 | // Searches pointer to:
|
|---|
| 243 | // - MPedestalCam
|
|---|
| 244 | // - MExtractedSignalCam
|
|---|
| 245 | //
|
|---|
| 246 | // Searches or creates:
|
|---|
| 247 | // - MCalibrationPedCam
|
|---|
| 248 | //
|
|---|
| 249 | // Retrieves from MExtractedSignalCam:
|
|---|
| 250 | // - number of used High Gain FADC slices (MExtractedSignalCam::GetNumUsedHiGainFADCSlices())
|
|---|
| 251 | //
|
|---|
| 252 | // Initializes, if empty to MGeomCam::GetNumPixels():
|
|---|
| 253 | // - MHCalibrationCam::fHiGainArray
|
|---|
| 254 | //
|
|---|
| 255 | // Initializes, if empty to MGeomCam::GetNumAreas() for:
|
|---|
| 256 | // - MHCalibrationCam::fAverageHiGainAreas
|
|---|
| 257 | //
|
|---|
| 258 | // Initializes, if empty to MGeomCam::GetNumSectors() for:
|
|---|
| 259 | // - MHCalibrationCam::fAverageHiGainSectors
|
|---|
| 260 | //
|
|---|
| 261 | // Calls MHCalibrationCam::InitPedHists() for every entry in:
|
|---|
| 262 | // - MHCalibrationCam::fHiGainArray
|
|---|
| 263 | // - MHCalibrationCam::fAverageHiGainAreas
|
|---|
| 264 | // - MHCalibrationCam::fAverageHiGainSectors
|
|---|
| 265 | //
|
|---|
| 266 | // Sets Titles and Names for the Histograms
|
|---|
| 267 | // - MHCalibrationCam::fAverageHiGainAreas
|
|---|
| 268 | //
|
|---|
| 269 | Bool_t MHPedestalCam::ReInitHists(MParList *pList)
|
|---|
| 270 | {
|
|---|
| 271 |
|
|---|
| 272 | MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
|
|---|
| 273 | if (!signal)
|
|---|
| 274 | {
|
|---|
| 275 | *fLog << err << "Cannot find MExtractedSignalCam... abort." << endl;
|
|---|
| 276 | return kFALSE;
|
|---|
| 277 | }
|
|---|
| 278 |
|
|---|
| 279 |
|
|---|
| 280 | fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
|
|---|
| 281 |
|
|---|
| 282 | if (!fPedestals)
|
|---|
| 283 | *fLog << warn << "Cannot find MPedestalCam ..." << endl;
|
|---|
| 284 |
|
|---|
| 285 | const Int_t npixels = fGeom->GetNumPixels();
|
|---|
| 286 | const Int_t nsectors = fGeom->GetNumSectors();
|
|---|
| 287 | const Int_t nareas = fGeom->GetNumAreas();
|
|---|
| 288 |
|
|---|
| 289 | fCam = (MCalibrationCam*)pList->FindObject("MCalibrationPedCam");
|
|---|
| 290 | if (!fCam)
|
|---|
| 291 | {
|
|---|
| 292 | fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationPedCam"));
|
|---|
| 293 | if (!fCam)
|
|---|
| 294 | return kFALSE;
|
|---|
| 295 |
|
|---|
| 296 | fCam->Init(*fGeom);
|
|---|
| 297 | }
|
|---|
| 298 |
|
|---|
| 299 | Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices();
|
|---|
| 300 |
|
|---|
| 301 | if (sliceshi == 0.)
|
|---|
| 302 | {
|
|---|
| 303 | *fLog << err << "Number of used signal slices in MExtractedSignalCam is zero ... abort."
|
|---|
| 304 | << endl;
|
|---|
| 305 | return kFALSE;
|
|---|
| 306 | }
|
|---|
| 307 |
|
|---|
| 308 | if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices )
|
|---|
| 309 | {
|
|---|
| 310 | *fLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam ... abort."
|
|---|
| 311 | << endl;
|
|---|
| 312 | return kFALSE;
|
|---|
| 313 | }
|
|---|
| 314 |
|
|---|
| 315 | fExtractHiGainSlices = sliceshi;
|
|---|
| 316 |
|
|---|
| 317 | *fLog << endl;
|
|---|
| 318 | *fLog << inf << GetDescriptor()
|
|---|
| 319 | << ": Number of found High Gain signal slices: " << fExtractHiGainSlices << endl;
|
|---|
| 320 |
|
|---|
| 321 | InitHiGainArrays(npixels,nareas,nsectors);
|
|---|
| 322 |
|
|---|
| 323 | if (fSum.GetSize()==0)
|
|---|
| 324 | {
|
|---|
| 325 | fSum. Set(npixels);
|
|---|
| 326 | fSum2. Set(npixels);
|
|---|
| 327 | fAreaSum. Set(nareas);
|
|---|
| 328 | fAreaSum2. Set(nareas);
|
|---|
| 329 | fAreaVal. Set(nareas);
|
|---|
| 330 | fSectorSum. Set(nsectors);
|
|---|
| 331 | fSectorSum2.Set(nsectors);
|
|---|
| 332 | fSectorVal.Set(nsectors);
|
|---|
| 333 | }
|
|---|
| 334 |
|
|---|
| 335 | fNumEvents = 0;
|
|---|
| 336 |
|
|---|
| 337 | return kTRUE;
|
|---|
| 338 | }
|
|---|
| 339 |
|
|---|
| 340 |
|
|---|
| 341 | // -------------------------------------------------------------------------------
|
|---|
| 342 | //
|
|---|
| 343 | // Retrieves pointer to MExtractedSignalCam:
|
|---|
| 344 | //
|
|---|
| 345 | // Retrieves from MGeomCam:
|
|---|
| 346 | // - number of pixels
|
|---|
| 347 | // - number of pixel areas
|
|---|
| 348 | // - number of sectors
|
|---|
| 349 | //
|
|---|
| 350 | // Fills HiGain histograms (MHGausEvents::FillHistAndArray()), respectively
|
|---|
| 351 | // with the signals MExtractedSignalCam::GetExtractedSignalHiGain().
|
|---|
| 352 | //
|
|---|
| 353 | Bool_t MHPedestalCam::FillHists(const MParContainer *par, const Stat_t w)
|
|---|
| 354 | {
|
|---|
| 355 |
|
|---|
| 356 | MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
|
|---|
| 357 | if (!signal)
|
|---|
| 358 | {
|
|---|
| 359 | *fLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
|
|---|
| 360 | return kFALSE;
|
|---|
| 361 | }
|
|---|
| 362 |
|
|---|
| 363 | const UInt_t npixels = fGeom->GetNumPixels();
|
|---|
| 364 | const UInt_t nareas = fGeom->GetNumAreas();
|
|---|
| 365 | const UInt_t nsectors = fGeom->GetNumSectors();
|
|---|
| 366 |
|
|---|
| 367 | MArrayF sumareahi(nareas);
|
|---|
| 368 | MArrayF sumsectorhi(nsectors);
|
|---|
| 369 | MArrayI numareahi(nareas);
|
|---|
| 370 | MArrayI numsectorhi(nsectors);
|
|---|
| 371 |
|
|---|
| 372 | for (UInt_t i=0; i<npixels; i++)
|
|---|
| 373 | {
|
|---|
| 374 |
|
|---|
| 375 | MHCalibrationPix &histhi = (*this)[i];
|
|---|
| 376 |
|
|---|
| 377 | if (histhi.IsExcluded())
|
|---|
| 378 | continue;
|
|---|
| 379 |
|
|---|
| 380 | const MExtractedSignalPix &pix = (*signal)[i];
|
|---|
| 381 |
|
|---|
| 382 | const Float_t pedhi = pix.GetExtractedSignalHiGain();
|
|---|
| 383 |
|
|---|
| 384 | const Int_t aidx = (*fGeom)[i].GetAidx();
|
|---|
| 385 | const Int_t sector = (*fGeom)[i].GetSector();
|
|---|
| 386 |
|
|---|
| 387 | histhi.FillHistAndArray(pedhi) ;
|
|---|
| 388 |
|
|---|
| 389 | fSum[i] += pedhi;
|
|---|
| 390 | fAreaSum[aidx] += pedhi;
|
|---|
| 391 | fSectorSum[sector] += pedhi;
|
|---|
| 392 |
|
|---|
| 393 | const Float_t sqrsum = pedhi*pedhi;
|
|---|
| 394 | fSum2[i] += sqrsum;
|
|---|
| 395 | fAreaSum2[aidx] += sqrsum;
|
|---|
| 396 | fSectorSum2[sector] += sqrsum;
|
|---|
| 397 |
|
|---|
| 398 | sumareahi [aidx] += pedhi;
|
|---|
| 399 | numareahi [aidx] ++;
|
|---|
| 400 | sumsectorhi[sector] += pedhi;
|
|---|
| 401 | numsectorhi[sector] ++;
|
|---|
| 402 |
|
|---|
| 403 | }
|
|---|
| 404 |
|
|---|
| 405 | for (UInt_t j=0; j<nareas; j++)
|
|---|
| 406 | {
|
|---|
| 407 | MHCalibrationPix &histhi = GetAverageHiGainArea(j);
|
|---|
| 408 | histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
|
|---|
| 409 | }
|
|---|
| 410 |
|
|---|
| 411 | for (UInt_t j=0; j<nsectors; j++)
|
|---|
| 412 | {
|
|---|
| 413 | MHCalibrationPix &histhi = GetAverageHiGainSector(j);
|
|---|
| 414 | histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
|
|---|
| 415 |
|
|---|
| 416 | }
|
|---|
| 417 |
|
|---|
| 418 | fNumEvents++;
|
|---|
| 419 |
|
|---|
| 420 | return kTRUE;
|
|---|
| 421 | }
|
|---|
| 422 |
|
|---|
| 423 | // --------------------------------------------------------------------------
|
|---|
| 424 | //
|
|---|
| 425 | // Calls:
|
|---|
| 426 | // - MHCalibrationCam::FitHiGainArrays() with Bad Pixels flags 0
|
|---|
| 427 | //
|
|---|
| 428 | Bool_t MHPedestalCam::FinalizeHists()
|
|---|
| 429 | {
|
|---|
| 430 |
|
|---|
| 431 | if (fNumEvents <= 1)
|
|---|
| 432 | {
|
|---|
| 433 | *fLog << err << GetDescriptor()
|
|---|
| 434 | << ": Number of processed events smaller or equal to 1" << endl;
|
|---|
| 435 | return kFALSE;
|
|---|
| 436 | }
|
|---|
| 437 |
|
|---|
| 438 | FitHists();
|
|---|
| 439 |
|
|---|
| 440 | if (fRenorm)
|
|---|
| 441 | RenormResults();
|
|---|
| 442 |
|
|---|
| 443 | return kTRUE;
|
|---|
| 444 |
|
|---|
| 445 | }
|
|---|
| 446 |
|
|---|
| 447 | void MHPedestalCam::FitHists()
|
|---|
| 448 | {
|
|---|
| 449 |
|
|---|
| 450 | for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
|
|---|
| 451 | {
|
|---|
| 452 |
|
|---|
| 453 | MHCalibrationPix &hist = (*this)[i];
|
|---|
| 454 | MCalibrationPix &pix = (*fCam)[i];
|
|---|
| 455 |
|
|---|
| 456 | //
|
|---|
| 457 | // 1) Store calculated means and variances in the low-gain slices
|
|---|
| 458 | //
|
|---|
| 459 | pix.SetLoGainMean ( fSum[i] / fNumEvents );
|
|---|
| 460 | const Double_t diff = fSum2[i] - fSum[i]*fSum[i]/fNumEvents;
|
|---|
| 461 | pix.SetLoGainSigma ( diff > 0. ? TMath::Sqrt( diff / (fNumEvents-1) ) : 0.);
|
|---|
| 462 | pix.SetLoGainMeanVar ( pix.GetLoGainSigma() * pix.GetLoGainSigma() / fNumEvents );
|
|---|
| 463 | pix.SetLoGainSigmaVar( pix.GetLoGainMeanVar() / 4. );
|
|---|
| 464 |
|
|---|
| 465 | if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
|
|---|
| 466 | continue;
|
|---|
| 467 |
|
|---|
| 468 | //
|
|---|
| 469 | // 2) Fit the Hi Gain histograms with a Gaussian
|
|---|
| 470 | //
|
|---|
| 471 | hist.FitGaus();
|
|---|
| 472 | //
|
|---|
| 473 | // 3) Take histogram means and RMS
|
|---|
| 474 | //
|
|---|
| 475 | hist.BypassFit();
|
|---|
| 476 | //
|
|---|
| 477 | // 4) Check for oscillations
|
|---|
| 478 | //
|
|---|
| 479 | if (IsOscillations())
|
|---|
| 480 | hist.CreateFourierSpectrum();
|
|---|
| 481 | //
|
|---|
| 482 | // 5) Retrieve the results and store them in this class
|
|---|
| 483 | //
|
|---|
| 484 | pix.SetHiGainMean ( hist.GetMean() );
|
|---|
| 485 | pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
|
|---|
| 486 | pix.SetHiGainRms ( hist.GetHistRms() );
|
|---|
| 487 | pix.SetHiGainSigma ( hist.GetSigma() );
|
|---|
| 488 | pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
|
|---|
| 489 | pix.SetHiGainProb ( hist.GetProb() );
|
|---|
| 490 | pix.SetHiGainNumBlackout( hist.GetBlackout() );
|
|---|
| 491 | pix.SetHiGainNumPickup ( hist.GetPickup() );
|
|---|
| 492 | }
|
|---|
| 493 |
|
|---|
| 494 | if (!IsAverageing())
|
|---|
| 495 | return;
|
|---|
| 496 |
|
|---|
| 497 | for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
|
|---|
| 498 | {
|
|---|
| 499 |
|
|---|
| 500 | MHCalibrationPix &hist = GetAverageHiGainArea(j);
|
|---|
| 501 | MCalibrationPix &pix = fCam->GetAverageArea(j);
|
|---|
| 502 |
|
|---|
| 503 | if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
|
|---|
| 504 | continue;
|
|---|
| 505 |
|
|---|
| 506 | //
|
|---|
| 507 | // 2) Fit the Hi Gain histograms with a Gaussian
|
|---|
| 508 | //
|
|---|
| 509 | hist.FitGaus();
|
|---|
| 510 | //
|
|---|
| 511 | // 3) Take histogram means and RMS
|
|---|
| 512 | //
|
|---|
| 513 | hist.BypassFit();
|
|---|
| 514 | //
|
|---|
| 515 | // 4) Check for oscillations
|
|---|
| 516 | //
|
|---|
| 517 | if (IsOscillations())
|
|---|
| 518 | hist.CreateFourierSpectrum();
|
|---|
| 519 | //
|
|---|
| 520 | // 5) Retrieve the results and store them in this class
|
|---|
| 521 | //
|
|---|
| 522 | pix.SetHiGainMean ( hist.GetMean() );
|
|---|
| 523 | pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
|
|---|
| 524 | pix.SetHiGainRms ( hist.GetHistRms() );
|
|---|
| 525 | pix.SetHiGainSigma ( hist.GetSigma() );
|
|---|
| 526 | pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
|
|---|
| 527 | pix.SetHiGainProb ( hist.GetProb() );
|
|---|
| 528 | pix.SetHiGainNumBlackout( hist.GetBlackout() );
|
|---|
| 529 | pix.SetHiGainNumPickup ( hist.GetPickup() );
|
|---|
| 530 | //
|
|---|
| 531 | // 6) Store calculated means and variances in the low-gain slices
|
|---|
| 532 | //
|
|---|
| 533 | const ULong_t aevts = fNumEvents * fAreaVal[j];
|
|---|
| 534 | if (aevts <= 1)
|
|---|
| 535 | continue;
|
|---|
| 536 |
|
|---|
| 537 | pix.SetLoGainMean ( fAreaSum[j] / aevts );
|
|---|
| 538 | const Double_t diff = fAreaSum2[j] - fAreaSum[j]*fAreaSum[j]/aevts ;
|
|---|
| 539 | pix.SetLoGainSigma( diff > 0. ? TMath::Sqrt( diff / (aevts-1) ) : 0.);
|
|---|
| 540 | }
|
|---|
| 541 |
|
|---|
| 542 | for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
|
|---|
| 543 | {
|
|---|
| 544 | MHCalibrationPix &hist = GetAverageHiGainSector(j);
|
|---|
| 545 | MCalibrationPix &pix = fCam->GetAverageSector(j);
|
|---|
| 546 |
|
|---|
| 547 | if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
|
|---|
| 548 | continue;
|
|---|
| 549 |
|
|---|
| 550 | //
|
|---|
| 551 | // 2) Fit the Hi Gain histograms with a Gaussian
|
|---|
| 552 | //
|
|---|
| 553 | hist.FitGaus();
|
|---|
| 554 | //
|
|---|
| 555 | // 3) Take histogram means and RMS
|
|---|
| 556 | //
|
|---|
| 557 | hist.BypassFit();
|
|---|
| 558 | //
|
|---|
| 559 | // 4) Check for oscillations
|
|---|
| 560 | //
|
|---|
| 561 | if (IsOscillations())
|
|---|
| 562 | hist.CreateFourierSpectrum();
|
|---|
| 563 | //
|
|---|
| 564 | // 5) Retrieve the results and store them in this class
|
|---|
| 565 | //
|
|---|
| 566 | pix.SetHiGainMean ( hist.GetMean() );
|
|---|
| 567 | pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
|
|---|
| 568 | pix.SetHiGainRms ( hist.GetHistRms() );
|
|---|
| 569 | pix.SetHiGainSigma ( hist.GetSigma() );
|
|---|
| 570 | pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
|
|---|
| 571 | pix.SetHiGainProb ( hist.GetProb() );
|
|---|
| 572 | pix.SetHiGainNumBlackout( hist.GetBlackout() );
|
|---|
| 573 | pix.SetHiGainNumPickup ( hist.GetPickup() );
|
|---|
| 574 | //
|
|---|
| 575 | // 6) Store calculated means and variances in the low-gain slices
|
|---|
| 576 | //
|
|---|
| 577 | const ULong_t sevts = fNumEvents * fSectorVal[j];
|
|---|
| 578 | if (sevts <= 1)
|
|---|
| 579 | continue;
|
|---|
| 580 |
|
|---|
| 581 | pix.SetLoGainMean ( fSectorSum[j] / sevts );
|
|---|
| 582 | const Double_t diff = fSectorSum2[j] - fSectorSum[j]*fSectorSum[j]/sevts ;
|
|---|
| 583 | pix.SetLoGainSigma( diff > 0. ? TMath::Sqrt( diff / (sevts-1) ) : 0.);
|
|---|
| 584 | }
|
|---|
| 585 | }
|
|---|
| 586 |
|
|---|
| 587 | // --------------------------------------------------------------------------
|
|---|
| 588 | //
|
|---|
| 589 | // Renormalizes the pedestal fit results by the following formulae:
|
|---|
| 590 | //
|
|---|
| 591 | // - Mean Pedestal / slice -> Mean Pedestal / Number slices
|
|---|
| 592 | // - Mean Pedestal Error / slice -> Mean Pedestal Error / Number slices
|
|---|
| 593 | // - Sigma Pedestal / slice -> Sigma Pedestal / Sqrt (Number slices)
|
|---|
| 594 | // - Sigma Pedestal Error / slice -> Sigma Pedestal Error / Sqrt (Number slices)
|
|---|
| 595 | //
|
|---|
| 596 | void MHPedestalCam::RenormResults()
|
|---|
| 597 | {
|
|---|
| 598 |
|
|---|
| 599 | //
|
|---|
| 600 | // One never knows...
|
|---|
| 601 | //
|
|---|
| 602 | if (fExtractHiGainSlices <= 0)
|
|---|
| 603 | return;
|
|---|
| 604 |
|
|---|
| 605 | const Float_t sqslices = TMath::Sqrt(fExtractHiGainSlices);
|
|---|
| 606 |
|
|---|
| 607 | for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
|
|---|
| 608 | {
|
|---|
| 609 |
|
|---|
| 610 | MCalibrationPix &pix = (*fCam)[i];
|
|---|
| 611 | pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
|
|---|
| 612 | pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
|
|---|
| 613 | //
|
|---|
| 614 | // Mean error goes with PedestalRMS/Sqrt(entries) -> scale with sqrt(slices)
|
|---|
| 615 | //
|
|---|
| 616 | pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
|
|---|
| 617 | pix.SetLoGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
|
|---|
| 618 | //
|
|---|
| 619 | // Sigma goes like PedestalRMS -> scale with sqrt(slices)
|
|---|
| 620 | //
|
|---|
| 621 | pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
|
|---|
| 622 | pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
|
|---|
| 623 | //
|
|---|
| 624 | // Sigma error goes like PedestalRMS/2.(entries) -> scale with sqrt(slices)
|
|---|
| 625 | //
|
|---|
| 626 | pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
|
|---|
| 627 | pix.SetLoGainSigmaVar ( pix.GetLoGainSigmaVar() / fExtractHiGainSlices );
|
|---|
| 628 | }
|
|---|
| 629 |
|
|---|
| 630 | if (!IsAverageing())
|
|---|
| 631 | return;
|
|---|
| 632 |
|
|---|
| 633 | for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
|
|---|
| 634 | {
|
|---|
| 635 |
|
|---|
| 636 | MCalibrationPix &pix = fCam->GetAverageArea(j);
|
|---|
| 637 | pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
|
|---|
| 638 | pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
|
|---|
| 639 | pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
|
|---|
| 640 | pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
|
|---|
| 641 | pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
|
|---|
| 642 | pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
|
|---|
| 643 | }
|
|---|
| 644 |
|
|---|
| 645 | for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
|
|---|
| 646 | {
|
|---|
| 647 | MCalibrationPix &pix = fCam->GetAverageSector(j);
|
|---|
| 648 | pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
|
|---|
| 649 | pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
|
|---|
| 650 | pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
|
|---|
| 651 | pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
|
|---|
| 652 | pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
|
|---|
| 653 | pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
|
|---|
| 654 | }
|
|---|
| 655 |
|
|---|
| 656 | }
|
|---|
| 657 |
|
|---|
| 658 |
|
|---|
| 659 | // ------------------------------------------------------------------
|
|---|
| 660 | //
|
|---|
| 661 | // The types are as follows:
|
|---|
| 662 | //
|
|---|
| 663 | // Fitted values:
|
|---|
| 664 | // ==============
|
|---|
| 665 | //
|
|---|
| 666 | // 0: Fitted Charge
|
|---|
| 667 | // 1: Error of fitted Charge
|
|---|
| 668 | // 2: Sigma of fitted Charge
|
|---|
| 669 | // 3: Error of Sigma of fitted Charge
|
|---|
| 670 | //
|
|---|
| 671 | //
|
|---|
| 672 | // Useful variables derived from the fit results:
|
|---|
| 673 | // =============================================
|
|---|
| 674 | //
|
|---|
| 675 | // 4: Returned probability of Gauss fit to Charge distribution
|
|---|
| 676 | // 5: Relative differenc of calculated pedestal (per slice) and fitted (per slice)
|
|---|
| 677 | // 6: Error of the Relative differenc of calculated pedestal (per slice) and fitted (per slice)
|
|---|
| 678 | // 7: Relative difference of the error of the mean pedestal (per slice) - calculated and fitted
|
|---|
| 679 | // 8: Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
|
|---|
| 680 | // 9: Error of Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
|
|---|
| 681 | // 10: Relative difference of the error of the pedestal RMS (per slice) - calculated and fitted
|
|---|
| 682 | //
|
|---|
| 683 | // Localized defects:
|
|---|
| 684 | // ==================
|
|---|
| 685 | //
|
|---|
| 686 | // 11: Gaus fit not OK
|
|---|
| 687 | // 12: Fourier spectrum not OK
|
|---|
| 688 | //
|
|---|
| 689 | Bool_t MHPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
|
|---|
| 690 | {
|
|---|
| 691 |
|
|---|
| 692 | if (fHiGainArray->GetSize() <= idx)
|
|---|
| 693 | return kFALSE;
|
|---|
| 694 |
|
|---|
| 695 | if ((*this)[idx].IsExcluded())
|
|---|
| 696 | return kFALSE;
|
|---|
| 697 |
|
|---|
| 698 | const Float_t ped = (*fPedestals)[idx].GetPedestal();
|
|---|
| 699 | const Float_t rms = (*fPedestals)[idx].GetPedestalRms();
|
|---|
| 700 |
|
|---|
| 701 | const Float_t entsqr = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
|
|---|
| 702 |
|
|---|
| 703 | const Float_t pederr = rms/entsqr;
|
|---|
| 704 | const Float_t rmserr = rms/entsqr/2.;
|
|---|
| 705 |
|
|---|
| 706 | const Float_t mean = (*this)[idx].GetMean();
|
|---|
| 707 | const Float_t meanerr = (*this)[idx].GetMeanErr();
|
|---|
| 708 | const Float_t sigma = (*this)[idx].GetSigma() ;
|
|---|
| 709 | const Float_t sigmaerr = (*this)[idx].GetSigmaErr();
|
|---|
| 710 |
|
|---|
| 711 | switch (type)
|
|---|
| 712 | {
|
|---|
| 713 | case 0:
|
|---|
| 714 | val = mean;
|
|---|
| 715 | break;
|
|---|
| 716 | case 1:
|
|---|
| 717 | val = meanerr;
|
|---|
| 718 | break;
|
|---|
| 719 | case 2:
|
|---|
| 720 | val = sigma;
|
|---|
| 721 | break;
|
|---|
| 722 | case 3:
|
|---|
| 723 | val = sigmaerr;
|
|---|
| 724 | break;
|
|---|
| 725 | case 4:
|
|---|
| 726 | val = (*this)[idx].GetProb();
|
|---|
| 727 | break;
|
|---|
| 728 | case 5:
|
|---|
| 729 | val = 2.*(mean-ped)/(ped+mean);
|
|---|
| 730 | break;
|
|---|
| 731 | case 6:
|
|---|
| 732 | val = TMath::Sqrt((pederr*pederr + meanerr*meanerr) * (ped*ped + mean*mean))
|
|---|
| 733 | *2./(ped+mean)/(ped+mean);
|
|---|
| 734 | break;
|
|---|
| 735 | case 7:
|
|---|
| 736 | val = 2.*(meanerr-pederr)/(pederr + meanerr);
|
|---|
| 737 | break;
|
|---|
| 738 | case 8:
|
|---|
| 739 | val = 2.*(sigma-rms)/(sigma+rms);
|
|---|
| 740 | break;
|
|---|
| 741 | case 9:
|
|---|
| 742 | val = TMath::Sqrt((rmserr*rmserr + sigmaerr*sigmaerr) * (rms*rms + sigma*sigma))
|
|---|
| 743 | *2./(rms+sigma)/(rms+sigma);
|
|---|
| 744 | break;
|
|---|
| 745 | case 10:
|
|---|
| 746 | val = 2.*(sigmaerr - rmserr)/(sigmaerr + rmserr);
|
|---|
| 747 | break;
|
|---|
| 748 | case 11:
|
|---|
| 749 | if (!(*this)[idx].IsGausFitOK())
|
|---|
| 750 | val = 1.;
|
|---|
| 751 | break;
|
|---|
| 752 | case 12:
|
|---|
| 753 | if (!(*this)[idx].IsFourierSpectrumOK())
|
|---|
| 754 | val = 1.;
|
|---|
| 755 | break;
|
|---|
| 756 | default:
|
|---|
| 757 | return kFALSE;
|
|---|
| 758 | }
|
|---|
| 759 | return kTRUE;
|
|---|
| 760 | }
|
|---|
| 761 |
|
|---|
| 762 | // --------------------------------------------------------------------------
|
|---|
| 763 | //
|
|---|
| 764 | // Calls MHGausEvents::DrawClone() for pixel idx
|
|---|
| 765 | //
|
|---|
| 766 | void MHPedestalCam::DrawPixelContent(Int_t idx) const
|
|---|
| 767 | {
|
|---|
| 768 | (*this)[idx].DrawClone();
|
|---|
| 769 | }
|
|---|