/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Markus Gaug 02/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // MCalibrationChargeCalc // // Task to calculate the calibration conversion factors and quantum efficiencies // from the fit results to the summed FADC slice distributions delivered by // MCalibrationChargeCam, MCalibrationChargePix, MCalibrationChargeBlindPix and // MCalibrationChargePINDiode, calculated and filled by MHCalibrationChargeCam, // MHCalibrationChargePix, MHCalibrationChargeBlindPix and MHCalibrationChargePINDiode. // // PreProcess(): Initialize pointers to MCalibrationChargeCam, MCalibrationChargeBlindPix // MCalibrationChargePINDiode and MCalibrationQECam // // Initialize pulser light wavelength // // ReInit(): MCalibrationCam::InitSize(NumPixels) is called from MGeomApply (which allocates // memory in a TClonesArray of type MCalibrationChargePix) // Initializes pointer to MBadPixelsCam // // Process(): Nothing to be done, histograms getting filled by MHCalibrationChargeCam // // PostProcess(): - FinalizePedestals() // - FinalizeCharges() // - FinalizeFFactorMethod() // - FinalizeBadPixels() // - FinalizeBlindPixel() // - FinalizePINDiode() // - FinalizeFFactorQECam() // - FinalizeBlindPixelQECam() // - FinalizePINDiodeQECam() // // Input Containers: // MCalibrationChargeCam // MCalibrationChargeBlindPix // MCalibrationChargePINDiode // MCalibrationQECam // MPedestalCam // MBadPixelsCam // MGeomCam // MTime // // Output Containers: // MCalibrationChargeCam // MCalibrationChargeBlindPix // MCalibrationChargePINDiode // MCalibrationQECam // MBadPixelsCam // // // Preliminary description of the calibration in photons (email from 12/02/04) // // Why calibrating in photons: // =========================== // // At the Barcelona meeting in 2002, we decided to calibrate the camera in // photons. This for the following reasons: // // * The physical quantity arriving at the camera are photons. This is // the direct physical information from the air shower. The photons // have a flux and a spectrum. // // * The photon fluxes depend mostly on the shower energy (with // corrections deriving from the observation conditions), while the photon // spectra depend mostly on the observation conditions: zenith angle, // quality of the air, also the impact parameter of the shower. // // * The photomultiplier, in turn, has different response properties // (quantum efficiencies) for photons of different colour. (Moreover, // different pixels have slightly different quantum efficiencies). // The resulting number of photo-electrons is then amplified (linearly) // with respect to the photo-electron flux. // // * In the ideal case, one would like to disentagle the effects // of the observation conditions from the primary particle energy (which // one likes to measure). To do so, one needs: // // 1) A reliable calibration relating the FADC counts to the photo-electron // flux -> This is accomplished with the F-Factor method. // // 2) A reliable calibration of the wavelength-dependent quantum efficiency // -> This is accomplished with the combination of the three methods, // together with QE-measurements performed by David in order to do // the interpolation. // // 3) A reliable calibration of the observation conditions. This means: // - Tracing the atmospheric conditions -> LIDAR // - Tracing the observation zenith angle -> Drive System // // 4) Some knowlegde about the impact parameter: // - This is the only part which cannot be accomplished well with a // single telescope. We would thus need to convolute the spectrum // over the distribution of impact parameters. // // // How an ideal calibration would look like: // ========================================= // // We know from the combined PIN-Diode and Blind-Pixel Method the response of // each pixel to well-measured light fluxes in three representative // wavelengths (green, blue, UV). We also know the response to these light // fluxes in photo-electrons. Thus, we can derive: // // - conversion factors to photo-electrons // - conversion factors to photons in three wavelengths. // // Together with David's measurements and some MC-simulation, we should be // able to derive tables for typical Cherenkov-photon spectra - convoluted // with the impact parameters and depending on the athmospheric conditions // and the zenith angle (the "outer parameters"). // // From these tables we can create "calibration tables" containing some // effective quantum efficiency depending on these outer parameters and which // are different for each pixel. // // In an ideal MCalibrate, one would thus have to convert first the FADC // slices to Photo-electrons and then, depending on the outer parameters, // look up the effective quantum efficiency and get the mean number of // photons which is then used for the further analysis. // // How the (first) MAGIC calibration should look like: // =================================================== // // For the moment, we have only one reliable calibration method, although // with very large systematic errors. This is the F-Factor method. Knowing // that the light is uniform over the whole camera (which I would not at all // guarantee in the case of the CT1 pulser), one could in principle already // perform a relative calibration of the quantum efficiencies in the UV. // However, the spread in QE at UV is about 10-15% (according to the plot // that Abelardo sent around last time. The spread in photo-electrons is 15% // for the inner pixels, but much larger (40%) for the outer ones. // // I'm not sure if we can already say that we have measured the relative // difference in quantum efficiency for the inner pixels and produce a first // QE-table for each pixel. To so, I would rather check in other wavelengths // (which we can do in about one-two weeks when the optical transmission of // the calibration trigger is installed). // // Thus, for the moment being, I would join Thomas proposal to calibrate in // photo-electrons and apply one stupid average quantum efficiency for all // pixels. This keeping in mind that we will have much preciser information // in about one to two weeks. // // // What MCalibrate should calculate and what should be stored: // =========================================================== // // It is clear that in the end, MCerPhotEvt will store photons. // MCalibrationCam stores the conversionfactors to photo-electrons and also // some tables of how to apply the conversion to photons, given the outer // parameters. This is not yet implemented and not even discussed. // // To start, I would suggest that we define the "average quantum efficiency" // (maybe something like 25+-3%) and apply them equally to all // photo-electrons. Later, this average factor can be easily replaced by a // pixel-dependent factor and later by a (pixel-dependent) table. // // // ////////////////////////////////////////////////////////////////////////////// #include "MCalibrationChargeCalc.h" #include #include #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MRawRunHeader.h" #include "MRawEvtPixelIter.h" #include "MGeomCam.h" #include "MGeomPix.h" #include "MPedestalCam.h" #include "MPedestalPix.h" #include "MCalibrationChargeCam.h" #include "MCalibrationChargePix.h" #include "MCalibrationChargePINDiode.h" #include "MCalibrationChargeBlindPix.h" #include "MExtractedSignalCam.h" #include "MExtractedSignalPix.h" #include "MExtractedSignalBlindPixel.h" #include "MExtractedSignalPINDiode.h" #include "MBadPixelsCam.h" #include "MBadPixelsPix.h" #include "MCalibrationQECam.h" #include "MCalibrationQEPix.h" #include "MCalibrationCam.h" ClassImp(MCalibrationChargeCalc); using namespace std; const Float_t MCalibrationChargeCalc::fgChargeLimit = 3.; const Float_t MCalibrationChargeCalc::fgChargeErrLimit = 0.; const Float_t MCalibrationChargeCalc::fgChargeRelErrLimit = 1.; const Float_t MCalibrationChargeCalc::fgLambdaErrLimit = 0.2; const Float_t MCalibrationChargeCalc::fgLambdaCheckLimit = 0.2; const Float_t MCalibrationChargeCalc::fgPheErrLimit = 5.; // -------------------------------------------------------------------------- // // Default constructor. // // Sets all pointers to NULL // // Calls AddToBranchList for: // - MRawEvtData.fHiGainPixId // - MRawEvtData.fLoGainPixId // - MRawEvtData.fHiGainFadcSamples // - MRawEvtData.fLoGainFadcSamples // // Initializes: // - fChargeLimit to fgChargeLimit // - fChargeErrLimit to fgChargeErrLimit // - fChargeRelErrLimit to fgChargeRelErrLimit // - fLambdaCheckLimit to fgLambdaCheckLimit // - fLambdaErrLimit to fgLambdaErrLimit // - fPheErrLimit to fgPheErrLimit // - fPulserColor to MCalibrationCam::kCT1 // // Calls: // - Clear() // MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title) : fBadPixels(NULL), fCam(NULL), fBlindPixel(NULL), fPINDiode(NULL), fQECam(NULL), fGeom(NULL), fPedestals(NULL), fEvtTime(NULL) { fName = name ? name : "MCalibrationChargeCalc"; fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam "; AddToBranchList("MRawEvtData.fHiGainPixId"); AddToBranchList("MRawEvtData.fLoGainPixId"); AddToBranchList("MRawEvtData.fHiGainFadcSamples"); AddToBranchList("MRawEvtData.fLoGainFadcSamples"); SetChargeLimit(); SetChargeErrLimit(); SetChargeRelErrLimit(); SetLambdaCheckLimit(); SetLambdaErrLimit(); SetPheErrLimit(); SetPulserColor(MCalibrationCam::kNONE); Clear(); } // -------------------------------------------------------------------------- // // Sets: // - all variables to 0., // - all flags to kFALSE // void MCalibrationChargeCalc::Clear(const Option_t *o) { fNumHiGainSamples = 0.; fNumLoGainSamples = 0.; fSqrtHiGainSamples = 0.; fSqrtLoGainSamples = 0.; SkipHiLoGainCalibration( kFALSE ); } // ----------------------------------------------------------------------------------- // // The following container are searched for and execution aborted if not in MParList: // - MPedestalCam // // The following containers are searched and created if they were not found: // // - MCalibrationQECam // - MBadPixelsCam // // The following output containers are only searched, but not created. If they // cannot be found, the corresponding calibration part is only skipped. // // - MTime // Int_t MCalibrationChargeCalc::PreProcess(MParList *pList) { // // Containers that have to be there. // fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); if (!fPedestals) { *fLog << err << "MPedestalCam not found... aborting" << endl; return kFALSE; } // // Containers that are created in case that they are not there. // fQECam = (MCalibrationQECam*)pList->FindCreateObj("MCalibrationQECam"); if (!fQECam) { *fLog << err << "Cannot find nor create MCalibrationQECam... aborting" << endl; return kFALSE; } fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam"); if (!fBadPixels) { *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl; return kFALSE; } fEvtTime = (MTime*)pList->FindObject("MTime"); // // Check the pulser colour --> FIXME: this solution is only valid until the arrival of the DM's // if (fPulserColor == MCalibrationCam::kNONE) { *fLog << endl; *fLog << err << GetDescriptor() << ": No Pulser colour has been chosen. Since the installation of the IFAE pulser box," << " you HAVE to provide the LEDs colour, otherwise there is no calibration. " << endl; *fLog << "See e.g. the macro calibration.C " << endl; return kFALSE; } return kTRUE; } // -------------------------------------------------------------------------- // // Search for the following input containers and abort if not existing: // - MGeomCam // - MCalibrationChargeCam // // Search for the following input containers and give a warning if not existing: // - MCalibrationChargeBlindPix // - MCalibrationChargePINDiode // // It retrieves the following variables from MCalibrationChargeCam: // // - fNumHiGainSamples // - fNumLoGainSamples // // It defines the PixId of every pixel in: // // - MCalibrationChargeCam // - MCalibrationQECam // // It sets all pixels in excluded which have the flag fBadBixelsPix::IsBad() set in: // // - MCalibrationChargePix // - MCalibrationQEPix // // Sets the pulser colour and tests if it has not changed w.r.t. fPulserColor in: // // - MCalibrationChargeCam // - MCalibrationChargeBlindPix (if existing) // - MCalibrationChargePINDiode (if existing) // Bool_t MCalibrationChargeCalc::ReInit(MParList *pList ) { fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); if (!fGeom) { *fLog << err << "No MGeomCam found... aborting." << endl; return kFALSE; } fCam = (MCalibrationChargeCam*)pList->FindObject("MCalibrationChargeCam"); if (!fCam) { *fLog << err << "Cannot find MCalibrationChargeCam... aborting" << endl; *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationChargeCam before..." << endl; return kFALSE; } // // Optional Containers // fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindObject("MCalibrationChargeBlindPix"); if (!fBlindPixel) *fLog << warn << GetDescriptor() << ": MCalibrationChargeBlindPix not found... no blind pixel method! " << endl; fPINDiode = (MCalibrationChargePINDiode*)pList->FindObject("MCalibrationChargePINDiode"); if (!fPINDiode) *fLog << warn << GetDescriptor() << "MCalibrationChargePINDiode not found... no PIN Diode method! " << endl; // // Initialize the pulser colours // if (fCam->GetPulserColor() == MCalibrationCam::kNONE) { fCam->SetPulserColor( fPulserColor ); if (fBlindPixel) fBlindPixel->SetColor( fPulserColor ); if (fPINDiode) fPINDiode->SetColor( fPulserColor ); } if (fPulserColor != fCam->GetPulserColor()) { *fLog << err << GetDescriptor() << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeCam" << endl; *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl; return kFALSE; } if (fBlindPixel) if (fPulserColor != fBlindPixel->GetColor()) { *fLog << err << GetDescriptor() << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeBlindPix." << endl; *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl; return kFALSE; } if (fPINDiode) if (fPulserColor != fPINDiode->GetColor()) { *fLog << err << GetDescriptor() << ": Pulser colour has changed w.r.t. last file in MCalibrationChargePINDiode." << endl; *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl; return kFALSE; } fNumHiGainSamples = fCam->GetNumHiGainFADCSlices(); fNumLoGainSamples = fCam->GetNumLoGainFADCSlices(); fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); UInt_t npixels = fGeom->GetNumPixels(); for (UInt_t i=0; iIsValid()) { *fLog << warn << GetDescriptor() << ": MCalibrationChargePINDiode is declared not valid... no PIN Diode method! " << endl; fPINDiode = NULL; } if (fBlindPixel) if (!fBlindPixel->IsValid()) { *fLog << warn << GetDescriptor() << ": MCalibrationChargeBlindPix is declared not valid... no Blind Pixel method! " << endl; fBlindPixel = NULL; } // // First loop over pixels, call FinalizePedestals and FinalizeCharges // Int_t nvalid = 0; for (Int_t pixid=0; pixidGetSize(); pixid++) { MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[pixid]; // // Check if the pixel has been excluded from the fits // if (pix.IsExcluded()) continue; MPedestalPix &ped = (*fPedestals)[pixid]; MBadPixelsPix &bad = (*fBadPixels)[pixid]; FinalizePedestals(ped,pix); if (FinalizeCharges(pix,bad)) nvalid++; } // // The Michele check ... // if (nvalid == 0) { *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. " << "Did you forget to fill the histograms " << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl; *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run " << "instead of a calibration run " << endl; return kFALSE; } for (UInt_t aidx=0; aidxGetNumAreas(); aidx++) { const MPedestalPix &ped = fPedestals->GetAverageArea(aidx); MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx); FinalizePedestals(ped,pix); FinalizeCharges(pix, fCam->GetAverageBadArea(aidx)); } for (UInt_t sector=0; sectorGetNumSectors(); sector++) { const MPedestalPix &ped = fPedestals->GetAverageSector(sector); MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(sector); FinalizePedestals(ped,pix); FinalizeCharges(pix, fCam->GetAverageBadSector(sector)); } // // Finalize Bad Pixels // FinalizeBadPixels(); // // Finalize F-Factor method // if (!FinalizeFFactorMethod()) { *fLog << warn << "Could not calculate the photons flux from the F-Factor method " << endl; fCam->SetFFactorMethodValid(kFALSE); return kFALSE; } else fCam->SetFFactorMethodValid(kTRUE); // // Finalize Blind Pixel // if (FinalizeBlindPixel()) fQECam->SetBlindPixelMethodValid(kTRUE); else fQECam->SetBlindPixelMethodValid(kFALSE); // // Finalize PIN Diode // if (FinalizePINDiode()) fQECam->SetPINDiodeMethodValid(kTRUE); else fQECam->SetPINDiodeMethodValid(kFALSE); // // Finalize QE Cam // FinalizeFFactorQECam(); FinalizeBlindPixelQECam(); FinalizePINDiodeQECam(); fCam ->SetReadyToSave(); fQECam ->SetReadyToSave(); fBadPixels->SetReadyToSave(); if (fBlindPixel) fBlindPixel->SetReadyToSave(); if (fPINDiode) fPINDiode->SetReadyToSave(); *fLog << inf << endl; *fLog << GetDescriptor() << ": Calibration statistics:" << endl; *fLog << dec << setfill(' '); UInt_t countinner = 0; UInt_t countouter = 0; for (Int_t i=0; iGetSize(); i++) { MBadPixelsPix &bad = (*fBadPixels)[i]; if (bad.IsOK()) { if (fGeom->GetPixRatio(i) == 1.) countinner++; else countouter++; } } *fLog << " " << setw(7) << "Successfully calibrated Pixels: " << "Inner: " << countinner << " Outer: " << countouter << endl; PrintUnsuitable(MBadPixelsPix::kUnsuitableRun, "Bad Pixels: "); PrintUnsuitable(MBadPixelsPix::kUnreliableRun, "Unreliable Pixels: "); *fLog << inf << endl; *fLog << GetDescriptor() << ": Errors statistics:" << endl; PrintUncalibrated(MBadPixelsPix::kChargeIsPedestal, Form("%s%2.1f%s","Signal less than ",fChargeLimit," Pedestal RMS: ")); PrintUncalibrated(MBadPixelsPix::kChargeErrNotValid, Form("%s%2.1f%s","Signal Error smaller than ",fChargeErrLimit,": ")); PrintUncalibrated(MBadPixelsPix::kChargeRelErrNotValid, Form("%s%2.1f%s","Signal Error bigger than ",fChargeRelErrLimit," times Mean Signal: ")); PrintUncalibrated(MBadPixelsPix::kChargeSigmaNotValid, "Signal Sigma smaller than Pedestal RMS: "); PrintUncalibrated(MBadPixelsPix::kLoGainSaturation, "Pixels with Low Gain Saturation: "); PrintUncalibrated(MBadPixelsPix::kMeanTimeInFirstBin, Form("%s%2.1f%s","Mean Abs. Arr. Time in First ",1.," Bin(s): ")); PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins, Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",2.," Bin(s): ")); PrintUncalibrated(MBadPixelsPix::kHiGainOscillating, "Pixels with changing Hi Gain signal over time: "); PrintUncalibrated(MBadPixelsPix::kLoGainOscillating, "Pixels with changing Lo Gain signal over time: "); PrintUncalibrated(MBadPixelsPix::kDeviatingNumPhes, "Pixels with deviating number of phes: "); PrintUncalibrated(MBadPixelsPix::kHiGainNotFitted, "Pixels with unsuccesful Gauss fit to the Hi Gain: "); PrintUncalibrated(MBadPixelsPix::kLoGainNotFitted, "Pixels with unsuccesful Gauss fit to the Lo Gain: "); return kTRUE; } // ---------------------------------------------------------------------------------- // // Retrieves pedestal and pedestal RMS from MPedestalPix // Retrieves total entries from MPedestalCam // Sets pedestal*fNumHiGainSamples and pedestal*fNumLoGainSamples in MCalibrationChargePix // Sets pedRMS *fSqrtHiGainSamples and pedRMS *fSqrtLoGainSamples in MCalibrationChargePix // // If the flag MCalibrationPix::IsHiGainSaturation() is set, call also: // - MCalibrationChargePix::CalcLoGainPedestal() // void MCalibrationChargeCalc::FinalizePedestals(const MPedestalPix &ped, MCalibrationChargePix &cal) { // // get the pedestals // const Float_t pedes = ped.GetPedestal(); const Float_t prms = ped.GetPedestalRms(); const Float_t num = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries()); // // set them in the calibration camera // if (cal.IsHiGainSaturation()) { cal.SetPedestal(pedes* fNumLoGainSamples, prms * fSqrtLoGainSamples, prms * fNumLoGainSamples / num); cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples); } else { cal.SetPedestal(pedes* fNumHiGainSamples, prms * fSqrtHiGainSamples, prms * fNumHiGainSamples / num); } } // ---------------------------------------------------------------------------------------------------- // // Check fit results validity. Bad Pixels flags are set if: // // 1) Pixel has a mean smaller than fChargeLimit*PedRMS ( Flag: MBadPixelsPix::kChargeIsPedestal) // 2) Pixel has a mean error smaller than fChargeErrLimit ( Flag: MBadPixelsPix::kChargeErrNotValid) // 3) Pixel has mean smaller than fChargeRelVarLimit times its mean error // ( Flag: MBadPixelsPix::kChargeRelErrNotValid) // 4) Pixel has a sigma bigger than its Pedestal RMS ( Flag: MBadPixelsPix::kChargeSigmaNotValid ) // // Further returns if flags: MBadPixelsPix::kUnsuitableRun is set // // Calls MCalibrationChargePix::CalcReducedSigma() and sets flag: MBadPixelsPix::kChargeIsPedestal // if not succesful. // // Calls MCalibrationChargePix::CalcFFactorMethod() and sets flag: MBadPixelsPix::kDeviatingNumPhes) // if not succesful. // Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad) { if (cal.GetMean() < fChargeLimit*cal.GetPedRms()) { *fLog << warn << GetDescriptor() << ": Fitted Charge: " << cal.GetMean() << " is smaller than " << fChargeLimit << " Pedestal RMS: " << cal.GetPedRms() << " in Pixel " << cal.GetPixId() << endl; bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal); } if (cal.GetMeanErr() < fChargeErrLimit) { *fLog << warn << GetDescriptor() << ": Error of Fitted Charge: " << cal.GetMeanErr() << " is smaller than " << fChargeErrLimit << " in Pixel " << cal.GetPixId() << endl; bad.SetUncalibrated( MBadPixelsPix::kChargeErrNotValid ); } if (cal.GetMean() < fChargeRelErrLimit*cal.GetMeanErr()) { *fLog << warn << GetDescriptor() << ": Fitted Charge: " << cal.GetMean() << " is smaller than " << fChargeRelErrLimit << "* its error: " << cal.GetMeanErr() << " in Pixel " << cal.GetPixId() << endl; bad.SetUncalibrated( MBadPixelsPix::kChargeRelErrNotValid ); } if (cal.GetSigma() < cal.GetPedRms()) { *fLog << warn << GetDescriptor() << ": Sigma of Fitted Charge: " << cal.GetSigma() << " smaller than Pedestal RMS: " << cal.GetPedRms() << " in Pixel " << cal.GetPixId() << endl; bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid ); } if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun)) return kFALSE; if (!cal.CalcReducedSigma()) { *fLog << warn << GetDescriptor() << ": Could not calculate reduced sigmas of pixel: " << cal.GetPixId() << endl; bad.SetUncalibrated(MBadPixelsPix::kChargeIsPedestal); return kFALSE; } if (!cal.CalcFFactorMethod()) { *fLog << warn << GetDescriptor() << ": Could not calculate F-Factor of pixel: " << cal.GetPixId() << endl; bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes); return kFALSE; } return kTRUE; } // ----------------------------------------------------------------------------------- // // Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set: // - MBadPixelsPix::kChargeIsPedestal // - MBadPixelsPix::kChargeErrNotValid // - MBadPixelsPix::kChargeRelErrNotValid // - MBadPixelsPix::kChargeSigmaNotValid // - MBadPixelsPix::kMeanTimeInFirstBin // - MBadPixelsPix::kMeanTimeInLast2Bins // // Sets pixel to MBadPixelsPix::kUnreliableRun, if one of the following flags is set: // - MBadPixelsPix::kDeviatingNumPhes // void MCalibrationChargeCalc::FinalizeBadPixels() { for (Int_t i=0; iGetSize(); i++) { MBadPixelsPix &bad = (*fBadPixels)[i]; if (bad.IsUncalibrated( MBadPixelsPix::kChargeIsPedestal)) bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); if (bad.IsUncalibrated( MBadPixelsPix::kChargeErrNotValid )) bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); if (bad.IsUncalibrated( MBadPixelsPix::kChargeRelErrNotValid )) bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); if (bad.IsUncalibrated( MBadPixelsPix::kChargeSigmaNotValid )) bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin )) bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins )) bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingNumPhes )) bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); } } // ------------------------------------------------------------------------ // // // First loop: Calculate a mean and mean RMS of photo-electrons per area index // Include only pixels which are not MBadPixelsPix::kUnsuitableRun or // MBadPixelsPix::kUnreliableRun (see FinalizeBadPixels()) and set // MCalibrationChargePix::SetFFactorMethodValid(kFALSE) in that case. // // Second loop: Get weighted mean number of photo-electrons and its RMS including // only pixels with flag MCalibrationChargePix::IsFFactorMethodValid() // and further exclude those deviating by more than fPheErrLimit mean // sigmas from the mean (obtained in first loop). Set // MBadPixelsPix::kDeviatingNumPhes if excluded. // // Set weighted mean and variance of photo-electrons per area index in: // average area pixels of MCalibrationChargeCam (obtained from: // MCalibrationChargeCam::GetAverageArea() ) // // Set weighted mean and variance of photo-electrons per sector in: // average sector pixels of MCalibrationChargeCam (obtained from: // MCalibrationChargeCam::GetAverageSector() ) // Bool_t MCalibrationChargeCalc::FinalizeFFactorMethod() { const UInt_t npixels = fGeom->GetNumPixels(); const UInt_t nareas = fGeom->GetNumAreas(); const UInt_t nsectors = fGeom->GetNumSectors(); Float_t lowlim [nareas]; Float_t upplim [nareas]; Float_t areavars [nareas]; Float_t areaweights [nareas], sectorweights [nsectors]; Float_t areaphes [nareas], sectorphes [nsectors]; Int_t numareavalid[nareas], numsectorvalid[nsectors]; memset(lowlim ,0, nareas * sizeof(Float_t)); memset(upplim ,0, nareas * sizeof(Float_t)); memset(areaphes ,0, nareas * sizeof(Float_t)); memset(areavars ,0, nareas * sizeof(Float_t)); memset(areaweights ,0, nareas * sizeof(Float_t)); memset(numareavalid ,0, nareas * sizeof(Int_t )); memset(sectorweights ,0, nsectors * sizeof(Float_t)); memset(sectorphes ,0, nsectors * sizeof(Float_t)); memset(numsectorvalid,0, nsectors * sizeof(Int_t )); // // First loop: Get mean number of photo-electrons and the RMS // The loop is only to recognize later pixels with very deviating numbers // for (UInt_t i=0; i 0.) { areaphes [aidx] += nphe; areavars [aidx] += nvar; numareavalid[aidx] ++; } } for (UInt_t i=0; i upplim[aidx] ) { *fLog << warn << GetDescriptor() << ": Deviating number of photo-electrons: " << Form("%4.2f",nphe) << " out of accepted limits: [" << Form("%4.2f%s%4.2f",lowlim[aidx],",",upplim[aidx]) << "] in pixel " << i << endl; bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes ); bad.SetUnsuitable ( MBadPixelsPix::kUnreliableRun ); continue; } const Float_t weight = 1./nvar; areaweights [aidx] += weight; areaphes [aidx] += weight*nphe; numareavalid [aidx] ++; sectorweights [sector] += weight; sectorphes [sector] += weight*nphe; numsectorvalid[sector] ++; } for (UInt_t aidx=0; aidxGetAverageArea(aidx); if (areaweights[aidx] <= 0. || areaphes[aidx] <= 0.) { *fLog << warn << " Mean number of phe's from area index " << aidx << " cannot be calculated: " << " Sum of weights: " << areaweights[aidx] << " Sum of weighted phes: " << areaphes[aidx] << endl; apix.SetFFactorMethodValid(kFALSE); continue; } *fLog << inf << "Replacing number photo-electrons of average area idx " << aidx << ": " << Form("%5.3f%s%5.3f",apix.GetPheFFactorMethod()," +- ",apix.GetPheFFactorMethodErr()) << endl; *fLog << inf << " by average number of photo-electrons from area idx " << aidx << ": " << Form("%5.3f%s%5.3f",areaphes[aidx] / areaweights[aidx]," +- ", TMath::Sqrt(1./areaweights[aidx])) << endl; apix.SetPheFFactorMethod ( areaphes[aidx]/ areaweights[aidx] ); apix.SetPheFFactorMethodVar( 1. / areaweights[aidx] ); apix.SetFFactorMethodValid ( kTRUE ); } for (UInt_t sector=0; sectorGetAverageSector(sector); if (sectorweights[sector] <= 0. || sectorphes[sector] <= 0.) { *fLog << warn << " Mean number of phe's from sector " << sector << " cannot be calculated: " << " Sum of weights: " << sectorweights[sector] << " Sum of weighted phes: " << sectorphes[sector] << endl; spix.SetFFactorMethodValid(kFALSE); continue; } *fLog << inf << "Replacing number photo-electrons of average sector " << sector << ": " << Form("%5.3f%s%5.3f",spix.GetPheFFactorMethod()," +- ",spix.GetPheFFactorMethodErr()) << endl; *fLog << inf << " by average number photo-electrons from sector " << sector << ": " << Form("%5.3f%s%5.3f",sectorphes[sector]/ sectorweights[sector]," +- ", TMath::Sqrt(1./sectorweights[sector])) << endl; spix.SetPheFFactorMethod ( sectorphes[sector]/ sectorweights[sector] ); spix.SetPheFFactorMethodVar( 1. / sectorweights[sector] ); spix.SetFFactorMethodValid ( kTRUE ); } return kTRUE; } // ------------------------------------------------------------------------ // // Returns kFALSE if pointer to MCalibrationChargeBlindPix is NULL // // The check returns kFALSE if: // // 1) fLambda and fLambdaCheck are separated relatively to each other by more than fLambdaCheckLimit // 2) BlindPixel has an fLambdaErr greater than fLambdaErrLimit // // Calls: // - MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass() // Bool_t MCalibrationChargeCalc::FinalizeBlindPixel() { if (!fBlindPixel) return kFALSE; const Float_t lambda = fBlindPixel->GetLambda(); const Float_t lambdaerr = fBlindPixel->GetLambdaErr(); const Float_t lambdacheck = fBlindPixel->GetLambdaCheck(); if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) < fLambdaCheckLimit) { *fLog << warn << GetDescriptor() << ": Lambda and Lambda-Check differ by more than " << fLambdaCheckLimit << " in the Blind Pixel " << endl; return kFALSE; } if (lambdaerr < fLambdaErrLimit) { *fLog << warn << GetDescriptor() << ": Error of Fitted Lambda is greater than " << fLambdaErrLimit << " in Blind Pixel " << endl; return kFALSE; } if (!fBlindPixel->CalcFluxInsidePlexiglass()) { *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, " << "will skip Blind Pixel Calibration " << endl; return kFALSE; } return kTRUE; } // ------------------------------------------------------------------------ // // Returns kFALSE if pointer to MCalibrationChargePINDiode is NULL // // The check returns kFALSE if: // // 1) PINDiode has a fitted charge smaller than fChargeLimit*PedRMS // 2) PINDiode has a fit error smaller than fChargeErrLimit // 3) PINDiode has a fitted charge smaller its fChargeRelErrLimit times its charge error // 4) PINDiode has a charge sigma smaller than its Pedestal RMS // // Calls: // - MCalibrationChargePINDiode::CalcFluxOutsidePlexiglass() // Bool_t MCalibrationChargeCalc::FinalizePINDiode() { if (!fPINDiode) return kFALSE; if (fPINDiode->GetMean() < fChargeLimit*fPINDiode->GetPedRms()) { *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than " << fChargeLimit << " Pedestal RMS in PINDiode " << endl; return kFALSE; } if (fPINDiode->GetMeanErr() < fChargeErrLimit) { *fLog << warn << GetDescriptor() << ": Error of Fitted Charge is smaller than " << fChargeErrLimit << " in PINDiode " << endl; return kFALSE; } if (fPINDiode->GetMean() < fChargeRelErrLimit*fPINDiode->GetMeanErr()) { *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than " << fChargeRelErrLimit << "* its error in PINDiode " << endl; return kFALSE; } if (fPINDiode->GetSigma() < fPINDiode->GetPedRms()) { *fLog << warn << GetDescriptor() << ": Sigma of Fitted Charge smaller than Pedestal RMS in PINDiode " << endl; return kFALSE; } if (!fPINDiode->CalcFluxOutsidePlexiglass()) { *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, " << "will skip PIN Diode Calibration " << endl; return kFALSE; } return kTRUE; } // ------------------------------------------------------------------------ // // Calculate the average number of photons outside the plexiglass with the // formula: // // av.Num.photons(area index) = av.Num.Phes(area index) // / MCalibrationQEPix::GetDefaultQE(fPulserColor) // / MCalibrationQECam::GetPlexiglassQE() // // Calculate the variance on the average number of photons. // // Loop over pixels: // // - Continue, if not MCalibrationChargePix::IsFFactorMethodValid() and set: // MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor) // // - Call MCalibrationChargePix::CalcMeanFFactor(av.Num.photons) and set: // MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor) if not succesful // // - Calculate the quantum efficiency with the formula: // // QE = ( Num.Phes / av.Num.photons ) * MGeomCam::GetPixRatio() // // - Set QE in MCalibrationQEPix::SetQEFFactor ( QE, fPulserColor ); // - Set Variance of QE in MCalibrationQEPix::SetQEFFactorVar ( Variance, fPulserColor ); // - Set bit MCalibrationQEPix::SetFFactorMethodValid(kTRUE,fPulserColor) // // - Call MCalibrationQEPix::UpdateFFactorMethod() // void MCalibrationChargeCalc::FinalizeFFactorQECam() { MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCam->GetAverageArea(0); MCalibrationQEPix &qepix = (MCalibrationQEPix&) fQECam->GetAverageArea(0); const Float_t avphotons = avpix.GetPheFFactorMethod() / qepix.GetQEFFactor(fPulserColor) / fQECam->GetPlexiglassQE(); const Float_t avphotrelvar = avpix.GetPheFFactorMethodRelVar() + qepix.GetQEFFactorRelVar(fPulserColor) + fQECam->GetPlexiglassQERelVar(); const UInt_t npixels = fGeom->GetNumPixels(); for (UInt_t i=0; iGetPixRatio(i); const Float_t qe = pix.GetPheFFactorMethod() / photons ; if (!pix.CalcMeanFFactor( photons , avphotrelvar )) { pix.SetFFactorMethodValid(kFALSE); qepix.SetFFactorMethodValid(kFALSE, fPulserColor); (*fBadPixels)[i].SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes ); } const Float_t qerelvar = avphotrelvar + pix.GetPheFFactorMethodRelVar(); qepix.SetQEFFactor ( qe , fPulserColor ); qepix.SetQEFFactorVar ( qerelvar*qe*qe, fPulserColor ); qepix.SetFFactorMethodValid( kTRUE , fPulserColor ); if (!qepix.UpdateFFactorMethod()) *fLog << warn << GetDescriptor() << ": Cannot update Quantum efficiencies with the F-Factor Method" << endl; } } // ------------------------------------------------------------------------ // // Loop over pixels: // // - Continue, if not MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable() and set: // MCalibrationQEPix::SetBlindPixelMethodValid(kFALSE,fPulserColor) // // - Calculate the quantum efficiency with the formula: // // QE = Num.Phes / MCalibrationChargeBlindPix::GetFluxInsidePlexiglass() // / MGeomPix::GetA() * MCalibrationQECam::GetPlexiglassQE() // // - Set QE in MCalibrationQEPix::SetQEBlindPixel ( QE, fPulserColor ); // - Set Variance of QE in MCalibrationQEPix::SetQEBlindPixelVar ( Variance, fPulserColor ); // - Set bit MCalibrationQEPix::SetBlindPixelMethodValid(kTRUE,fPulserColor) // // - Call MCalibrationQEPix::UpdateBlindPixelMethod() // void MCalibrationChargeCalc::FinalizeBlindPixelQECam() { const UInt_t npixels = fGeom->GetNumPixels(); // // With the knowledge of the overall photon flux, calculate the // quantum efficiencies after the Blind Pixel and PIN Diode method // for (UInt_t i=0; iIsFluxInsidePlexiglassAvailable()) { qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor); continue; } MBadPixelsPix &bad = (*fBadPixels)[i]; if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun)) { qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor); continue; } MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i]; MGeomPix &geo = (*fGeom)[i]; const Float_t qe = pix.GetPheFFactorMethod() / fBlindPixel->GetFluxInsidePlexiglass() / geo.GetA() * fQECam->GetPlexiglassQE(); const Float_t qerelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar() + fQECam->GetPlexiglassQERelVar() + pix.GetPheFFactorMethodRelVar(); qepix.SetQEBlindPixel ( qe , fPulserColor ); qepix.SetQEBlindPixelVar ( qerelvar*qe*qe, fPulserColor ); qepix.UpdateBlindPixelMethod(); } } // ------------------------------------------------------------------------ // // Loop over pixels: // // - Continue, if not MCalibrationChargePINDiode::IsFluxOutsidePlexiglassAvailable() and set: // MCalibrationQEPix::SetPINDiodeMethodValid(kFALSE,fPulserColor) // // - Calculate the quantum efficiency with the formula: // // QE = Num.Phes / MCalibrationChargePINDiode::GetFluxOutsidePlexiglass() / MGeomPix::GetA() // // - Set QE in MCalibrationQEPix::SetQEPINDiode ( QE, fPulserColor ); // - Set Variance of QE in MCalibrationQEPix::SetQEPINDiodeVar ( Variance, fPulserColor ); // - Set bit MCalibrationQEPix::SetPINDiodeMethodValid(kTRUE,fPulserColor) // // - Call MCalibrationQEPix::UpdatePINDiodeMethod() // void MCalibrationChargeCalc::FinalizePINDiodeQECam() { const UInt_t npixels = fGeom->GetNumPixels(); // // With the knowledge of the overall photon flux, calculate the // quantum efficiencies after the PIN Diode method // for (UInt_t i=0; iIsFluxOutsidePlexiglassAvailable()) { qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor); continue; } MBadPixelsPix &bad = (*fBadPixels)[i]; if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun)) { qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor); continue; } MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i]; MGeomPix &geo = (*fGeom)[i]; const Float_t qe = pix.GetPheFFactorMethod() / fPINDiode->GetFluxOutsidePlexiglass() / geo.GetA(); const Float_t qerelvar = fPINDiode->GetFluxOutsidePlexiglassRelVar() + pix.GetPheFFactorMethodRelVar(); qepix.SetQEPINDiode ( qe , fPulserColor ); qepix.SetQEPINDiodeVar ( qerelvar*qe*qe, fPulserColor ); qepix.UpdateBlindPixelMethod(); } } // ----------------------------------------------------------------------------------------------- // // Print out statistics about BadPixels of type UnsuitableType_t // void MCalibrationChargeCalc::PrintUnsuitable(MBadPixelsPix::UnsuitableType_t typ, const char *text) const { UInt_t countinner = 0; UInt_t countouter = 0; for (Int_t i=0; iGetSize(); i++) { MBadPixelsPix &bad = (*fBadPixels)[i]; if (bad.IsUnsuitable(typ)) { if (fGeom->GetPixRatio(i) == 1.) countinner++; else countouter++; } } *fLog << " " << setw(7) << text << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl; } // ----------------------------------------------------------------------------------------------- // // Print out statistics about BadPixels of type UncalibratedType_t // void MCalibrationChargeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const { UInt_t countinner = 0; UInt_t countouter = 0; for (Int_t i=0; iGetSize(); i++) { MBadPixelsPix &bad = (*fBadPixels)[i]; if (bad.IsUncalibrated(typ)) { if (fGeom->GetPixRatio(i) == 1.) countinner++; else countouter++; } } *fLog << " " << setw(7) << text << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl; }