/* ======================================================================== *\ ! ! * ! * 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 11/2003 ! ! Copyright: MAGIC Software Development, 2000-2001 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MCalibrationCam // // Hold the whole Calibration results of the camera: // // 1) MCalibrationCam initializes a TClonesArray whose elements are // pointers to MCalibrationPix Containers // 2) It initializes a pointer to an MCalibrationBlindPix container // 3) It initializes a pointer to an MCalibrationPINDiode container // ///////////////////////////////////////////////////////////////////////////// #include "MCalibrationCam.h" #include #include #include #include "MLog.h" #include "MLogManip.h" #include "MGeomCam.h" #include "MGeomPix.h" #include "MCalibrationPix.h" #include "MCalibrationConfig.h" #include "MCalibrationBlindPix.h" #include "MCalibrationPINDiode.h" #include "MHCalibrationPixel.h" ClassImp(MCalibrationCam); using namespace std; const Int_t MCalibrationCam::gkBlindPixelId = 559; const Int_t MCalibrationCam::gkPINDiodeId = 9999; // -------------------------------------------------------------------------- // // Default constructor. // // Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry // Later, a call to MCalibrationCam::InitSize(Int_t size) has to be performed // // Creates an MCalibrationBlindPix container // Creates an MCalibrationPINDiode container // MCalibrationCam::MCalibrationCam(const char *name, const char *title) : fOffsets(NULL), fSlopes(NULL), fOffvsSlope(NULL) { fName = name ? name : "MCalibrationCam"; fTitle = title ? title : "Storage container for the Calibration Information in the camera"; fPixels = new TClonesArray("MCalibrationPix",1); fBlindPixel = new MCalibrationBlindPix(); fPINDiode = new MCalibrationPINDiode(); Clear(); } // -------------------------------------------------------------------------- // // Delete the TClonesArray of MCalibrationPix containers // Delete the MCalibrationPINDiode and the MCalibrationBlindPix // // Delete the histograms if they exist // MCalibrationCam::~MCalibrationCam() { // // delete fPixels should delete all Objects stored inside // delete fPixels; delete fBlindPixel; delete fPINDiode; if (fOffsets) delete fOffsets; if (fSlopes) delete fSlopes; if (fOffvsSlope) delete fOffvsSlope; } // ------------------------------------------------------------------- // // This function simply allocates memory via the ROOT command: // (TObject**) TStorage::ReAlloc(fCont, newSize * sizeof(TObject*), // fSize * sizeof(TObject*)); // newSize corresponds to size in our case // fSize is the old size (in most cases: 1) // void MCalibrationCam::InitSize(const UInt_t i) { // // check if we have already initialized to size // if (CheckBounds(i)) return; fPixels->ExpandCreate(i); } // -------------------------------------------------------------------------- // // This function returns the current size of the TClonesArray // independently if the MCalibrationPix is filled with values or not. // // It is the size of the array fPixels. // Int_t MCalibrationCam::GetSize() const { return fPixels->GetEntriesFast(); } // -------------------------------------------------------------------------- // // Check if position i is inside the current bounds of the TClonesArray // Bool_t MCalibrationCam::CheckBounds(Int_t i) const { return i < GetSize(); } // -------------------------------------------------------------------------- // // Get i-th pixel (pixel number) // MCalibrationPix &MCalibrationCam::operator[](Int_t i) { return *static_cast(fPixels->UncheckedAt(i)); } // -------------------------------------------------------------------------- // // Get i-th pixel (pixel number) // MCalibrationPix &MCalibrationCam::operator[](Int_t i) const { return *static_cast(fPixels->UncheckedAt(i)); } // -------------------------------------- // void MCalibrationCam::Clear(Option_t *o) { fPixels->ForEach(TObject, Clear)(); fBlindPixel->Clear(); fPINDiode->Clear(); fMeanPhotInsidePlexiglass = -1.; fMeanPhotErrInsidePlexiglass = -1.; fMeanPhotOutsidePlexiglass = -1.; fMeanPhotErrOutsidePlexiglass = -1.; fNumExcludedPixels = 0; CLRBIT(fFlags,kBlindPixelMethodValid); CLRBIT(fFlags,kPINDiodeMethodValid); CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable); CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); return; } void MCalibrationCam::SetBlindPixelMethodValid(const Bool_t b) { if (b) SETBIT(fFlags, kBlindPixelMethodValid); else CLRBIT(fFlags, kBlindPixelMethodValid); } void MCalibrationCam::SetPINDiodeMethodValid(const Bool_t b) { if (b) SETBIT(fFlags, kPINDiodeMethodValid); else CLRBIT(fFlags, kPINDiodeMethodValid); } Bool_t MCalibrationCam::IsBlindPixelMethodValid() const { return TESTBIT(fFlags,kBlindPixelMethodValid); } Bool_t MCalibrationCam::IsPINDiodeMethodValid() const { return TESTBIT(fFlags,kPINDiodeMethodValid); } Bool_t MCalibrationCam::IsNumPhotInsidePlexiglassAvailable() const { return TESTBIT(fFlags,kNumPhotInsidePlexiglassAvailable); } Bool_t MCalibrationCam::IsNumPhotOutsidePlexiglassAvailable() const { return TESTBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); } // -------------------------------------------------------------------------- // // Print first the well fitted pixels // and then the ones which are not FitValid // void MCalibrationCam::Print(Option_t *o) const { *fLog << all << GetDescriptor() << ":" << endl; int id = 0; *fLog << all << "Succesfully calibrated pixels:" << endl; *fLog << all << endl; TIter Next(fPixels); MCalibrationPix *pix; while ((pix=(MCalibrationPix*)Next())) { if (pix->IsChargeValid() && !pix->IsExcluded()) { *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() << " Nr Phe's: " << pix->GetPheFFactorMethod() << endl; id++; } } *fLog << all << id << " succesful pixels :-))" << endl; id = 0; *fLog << all << endl; *fLog << all << "Pixels with errors:" << endl; *fLog << all << endl; TIter Next2(fPixels); while ((pix=(MCalibrationPix*)Next2())) { if (!pix->IsChargeValid() && !pix->IsExcluded()) { *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() << endl; id++; } } *fLog << all << id << " pixels with errors :-((" << endl; *fLog << all << endl; *fLog << all << "Excluded pixels:" << endl; *fLog << all << endl; TIter Next3(fPixels); while ((pix=(MCalibrationPix*)Next3())) if (pix->IsExcluded()) *fLog << all << pix->GetPixId() << endl; *fLog << all << fNumExcludedPixels << " excluded pixels " << endl; } // -------------------------------------------------------------------------- // // Return true if pixel is inside bounds of the TClonesArray fPixels // Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const { if (!CheckBounds(idx)) return kFALSE; return kTRUE; } // -------------------------------------------------------------------------- // // Return true if pixel has already been fitted once (independent of the result) // Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const { if (!CheckBounds(idx)) return kFALSE; return (*this)[idx].IsFitted(); } // -------------------------------------------------------------------------- // // Sets the user ranges of all histograms such that // empty bins at the edges are not used. Additionally, it rebins the // histograms such that in total, 50 bins are used. // void MCalibrationCam::CutEdges() { fBlindPixel->GetHist()->CutAllEdges(); fPINDiode->GetHist()->CutAllEdges(); TIter Next(fPixels); MCalibrationPix *pix; while ((pix=(MCalibrationPix*)Next())) { pix->GetHist()->CutAllEdges(); } return; } // The types are as follows: // // 0: Fitted Charge // 1: Error of fitted Charge // 2: Sigma of fitted Charge // 3: Error of Sigma of fitted Charge // 4: Returned probability of Gauss fit to Charge distribution // 5: Mean arrival time // 6: Sigma of the arrival time // 7: Chi-square of the Gauss fit to the arrival times // 8: Pedestal // 9: Pedestal RMS // 10: Reduced Sigma Square // 11: Number of Photo-electrons after the F-Factor method // 12: Error on the Number of Photo-electrons after the F-Factor method // 13: Mean conversion factor after the F-Factor method // 14: Error on the conversion factor after the F-Factor method // 15: Number of Photons after the Blind Pixel method // 16: Mean conversion factor after the Blind Pixel method // Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const { if (idx > GetSize()) return kFALSE; Float_t arearatio = cam.GetPixRatio(idx); if (arearatio == 0) return kFALSE; switch (type) { case 0: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetCharge(); break; case 1: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetErrCharge(); break; case 2: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetSigmaCharge(); break; case 3: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetErrSigmaCharge(); break; case 4: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetChargeProb(); break; case 5: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetRSigmaCharge(); break; case 6: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetErrRSigmaCharge(); break; case 7: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge(); break; case 8: if ((*this)[idx].IsExcluded()) return kFALSE; // relative error RsigmaCharge square val = (*this)[idx].GetErrRSigmaCharge()* (*this)[idx].GetErrRSigmaCharge() / ((*this)[idx].GetRSigmaCharge() * (*this)[idx].GetRSigmaCharge() ); // relative error Charge square val += (*this)[idx].GetErrCharge() * (*this)[idx].GetErrCharge() / ((*this)[idx].GetCharge() * (*this)[idx].GetCharge() ); // calculate relative error out of squares val = TMath::Sqrt(val) ; // multiply with value to get absolute error val *= (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge(); break; case 9: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetPheFFactorMethod(); break; case 10: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetPheFFactorMethodError(); break; case 11: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetMeanConversionFFactorMethod(); break; case 12: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetErrorConversionFFactorMethod(); break; case 13: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetTotalFFactorFFactorMethod(); break; case 14: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetTotalFFactorErrorFFactorMethod(); break; case 15: if ((*this)[idx].IsExcluded()) return kFALSE; if (arearatio == 1.) val = GetMeanPhotInsidePlexiglass(); else val = GetMeanPhotInsidePlexiglass()*gkCalibrationOutervsInnerPixelArea; break; case 16: if ((*this)[idx].IsExcluded()) return kFALSE; if (arearatio == 1.) val = (double)fMeanPhotInsidePlexiglass; else val = (double)fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; break; case 17: if ((*this)[idx].IsExcluded()) return kFALSE; if (arearatio == 1.) val = (*this)[idx].GetMeanConversionBlindPixelMethod(); else val = (*this)[idx].GetMeanConversionBlindPixelMethod()*gkCalibrationOutervsInnerPixelArea; break; case 18: if ((*this)[idx].IsExcluded()) return kFALSE; if (arearatio == 1.) val = (*this)[idx].GetErrorConversionBlindPixelMethod(); else { val = (*this)[idx].GetErrorConversionBlindPixelMethod()*(*this)[idx].GetErrorConversionBlindPixelMethod() * gkCalibrationOutervsInnerPixelArea * gkCalibrationOutervsInnerPixelArea; val += gkCalibrationOutervsInnerPixelAreaError * gkCalibrationOutervsInnerPixelAreaError * (*this)[idx].GetMeanConversionBlindPixelMethod() *(*this)[idx].GetMeanConversionBlindPixelMethod(); val = TMath::Sqrt(val); } break; case 19: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetTotalFFactorBlindPixelMethod(); break; case 20: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod(); break; case 21: if ((*this)[idx].IsExcluded()) return kFALSE; if (arearatio == 1.) val = GetMeanPhotOutsidePlexiglass(); else val = GetMeanPhotOutsidePlexiglass()*gkCalibrationOutervsInnerPixelArea; break; case 22: if ((*this)[idx].IsExcluded()) return kFALSE; if (arearatio == 1.) val = (double)fMeanPhotOutsidePlexiglass; else val = (double)fMeanPhotOutsidePlexiglass*gkCalibrationOutervsInnerPixelArea; break; case 23: if ((*this)[idx].IsExcluded()) return kFALSE; if (arearatio == 1.) val = (*this)[idx].GetMeanConversionPINDiodeMethod(); else val = (*this)[idx].GetMeanConversionPINDiodeMethod()*gkCalibrationOutervsInnerPixelArea; break; case 24: if ((*this)[idx].IsExcluded()) return kFALSE; if (arearatio == 1.) val = (*this)[idx].GetErrorConversionPINDiodeMethod(); else { val = (*this)[idx].GetErrorConversionPINDiodeMethod()*(*this)[idx].GetErrorConversionPINDiodeMethod() * gkCalibrationOutervsInnerPixelArea * gkCalibrationOutervsInnerPixelArea; val += gkCalibrationOutervsInnerPixelAreaError * gkCalibrationOutervsInnerPixelAreaError * (*this)[idx].GetMeanConversionPINDiodeMethod() *(*this)[idx].GetMeanConversionPINDiodeMethod(); val = TMath::Sqrt(val); } break; case 25: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetTotalFFactorBlindPixelMethod(); break; case 26: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod(); break; case 27: if ((*this)[idx].IsExcluded()) val = 1.; else return kFALSE; break; case 28: if ((*this)[idx].IsExcluded()) return kFALSE; if (!(*this)[idx].IsFitted()) val = 1; else return kFALSE; break; case 29: if ((*this)[idx].IsExcluded()) return kFALSE; if (!(*this)[idx].IsFitted()) return kFALSE; if (!(*this)[idx].IsChargeValid()) val = 1; else return kFALSE; break; case 30: if ((*this)[idx].IsExcluded()) return kFALSE; if ((*this)[idx].IsOscillating()) val = 1; else return kFALSE; break; case 31: if ((*this)[idx].IsExcluded()) return kFALSE; if ((*this)[idx].IsHiGainSaturation()) val = 1; else return kFALSE; break; case 32: if ((*this)[idx].IsExcluded()) return kFALSE; if ((*this)[idx].IsFFactorMethodValid()) val = 1; else return kFALSE; break; case 33: if ((*this)[idx].IsExcluded()) return kFALSE; if ((*this)[idx].IsBlindPixelMethodValid()) val = 1; else return kFALSE; break; case 34: if ((*this)[idx].IsExcluded()) return kFALSE; if ((*this)[idx].IsPINDiodeMethodValid()) val = 1; else return kFALSE; break; case 35: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetPed(); break; case 36: if ((*this)[idx].IsExcluded()) return kFALSE; val = 1.; // val = (*this)[idx].GetPedError(); break; case 37: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetPedRms(); break; case 38: if ((*this)[idx].IsExcluded()) return kFALSE; val = 1.; // val = (*this)[idx].GetPedRmsError(); break; case 39: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetAbsTimeMean(); break; case 40: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetAbsTimeMeanErr(); break; case 41: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetAbsTimeRms(); break; case 42: if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetAbsTimeMeanErr()/TMath::Sqrt(2.); break; default: return kFALSE; } return val!=-1.; } // -------------------------------------------------------------------------- // // What MHCamera needs in order to draw an individual pixel in the camera // void MCalibrationCam::DrawPixelContent(Int_t idx) const { (*this)[idx].Draw(); } // -------------------------------------------------------------------------- // // // Bool_t MCalibrationCam::CalcNumPhotInsidePlexiglass() { if (!fBlindPixel->IsFitOK()) return kFALSE; const Float_t mean = fBlindPixel->GetLambda(); const Float_t merr = fBlindPixel->GetErrLambda(); // // Start calculation of number of photons // // The blind pixel has exactly one cm^2 area (with negligible error), // thus we omi division by 1. // fMeanPhotInsidePlexiglass = mean * gkCalibrationInnerPixelArea; // Start calculation of number of photons relative Variance (!!) fMeanPhotErrInsidePlexiglass = merr*merr/mean/mean; fMeanPhotErrInsidePlexiglass += gkCalibrationInnerPixelAreaError*gkCalibrationInnerPixelAreaError / gkCalibrationInnerPixelArea/gkCalibrationInnerPixelArea; switch (fColor) { case kECGreen: fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEGreen; fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEGreenError*gkCalibrationBlindPixelQEGreenError / gkCalibrationBlindPixelQEGreen / gkCalibrationBlindPixelQEGreen; fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttGreen); // correct for absorption // attenuation has negligible error break; case kECBlue: fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEBlue; fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEBlueError*gkCalibrationBlindPixelQEBlueError / gkCalibrationBlindPixelQEBlue / gkCalibrationBlindPixelQEBlue; fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttBlue); // correct for absorption // attenuation has negligible error break; case kECUV: fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEUV; fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEUVError*gkCalibrationBlindPixelQEUVError / gkCalibrationBlindPixelQEUV / gkCalibrationBlindPixelQEUV; fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttUV); // correct for absorption // attenuation has negligible error break; case kECCT1: default: fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQECT1; fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQECT1Error*gkCalibrationBlindPixelQECT1Error / gkCalibrationBlindPixelQECT1 / gkCalibrationBlindPixelQECT1; fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttCT1); // correct for absorption // attenuation has negligible error break; } *fLog << inf << endl; *fLog << inf << " Mean number of Photons for an Inner Pixel (inside Plexiglass): " << fMeanPhotInsidePlexiglass << endl; if (fMeanPhotInsidePlexiglass > 0.) SETBIT(fFlags,kNumPhotInsidePlexiglassAvailable); else { CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable); return kFALSE; } if (fMeanPhotErrInsidePlexiglass < 0.) { *fLog << warn << " Relative Variance on number of Photons for an Inner Pixel (inside Plexiglass): " << fMeanPhotErrInsidePlexiglass << endl; CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable); return kFALSE; } // Finish calculation of errors -> convert from relative variance to absolute error fMeanPhotErrInsidePlexiglass = TMath::Sqrt(fMeanPhotErrInsidePlexiglass); fMeanPhotErrInsidePlexiglass *= fMeanPhotInsidePlexiglass; *fLog << inf << " Error on number of Photons for an Inner Pixel (inside Plexiglass): " << fMeanPhotErrInsidePlexiglass << endl; *fLog << inf << endl; TIter Next(fPixels); MCalibrationPix *pix; while ((pix=(MCalibrationPix*)Next())) { if(pix->IsChargeValid()) { const Float_t charge = pix->GetCharge(); const Float_t ratio = fGeomCam->GetPixRatio(pix->GetPixId()); const Float_t chargeerr = pix->GetErrCharge(); Float_t nphot; Float_t nphoterr; if (ratio == 1.) { nphot = fMeanPhotInsidePlexiglass; nphoterr = fMeanPhotErrInsidePlexiglass; } else { nphot = fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; nphoterr = fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea *fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; nphoterr += fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError *fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError; nphoterr = TMath::Sqrt(nphoterr); } const Float_t conversion = nphot/charge; Float_t conversionerr; conversionerr = nphoterr/charge * nphoterr/charge ; conversionerr += chargeerr/charge * chargeerr/charge * conversion*conversion; conversionerr = TMath::Sqrt(conversionerr); const Float_t conversionsigma = 0.; pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma); if (conversionerr/conversion < 0.1) pix->SetBlindPixelMethodValid(); } } return kTRUE; } Bool_t MCalibrationCam::CalcNumPhotOutsidePlexiglass() { if (!fPINDiode->IsChargeFitValid()) return kFALSE; const Float_t mean = fPINDiode->GetCharge(); const Float_t merr = fPINDiode->GetErrCharge(); // Start calculation of number of photons fMeanPhotOutsidePlexiglass = mean * gkCalibrationInnerPixelvsPINDiodeArea; // Start calculation of number of photons relative Variance (!!) fMeanPhotErrOutsidePlexiglass = merr*merr/mean/mean; fMeanPhotErrOutsidePlexiglass += gkCalibrationInnerPixelvsPINDiodeAreaError*gkCalibrationInnerPixelvsPINDiodeAreaError / gkCalibrationInnerPixelvsPINDiodeArea/gkCalibrationInnerPixelvsPINDiodeArea; switch (fColor) { case kECGreen: fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEGreen; fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEGreenError*gkCalibrationPINDiodeQEGreenError / gkCalibrationPINDiodeQEGreen/gkCalibrationPINDiodeQEGreen; break; case kECBlue: fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEBlue; fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEBlueError*gkCalibrationPINDiodeQEBlueError / gkCalibrationPINDiodeQEBlue/gkCalibrationPINDiodeQEBlue; break; case kECUV: fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEUV; fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEUVError*gkCalibrationPINDiodeQEUVError / gkCalibrationPINDiodeQEUV/gkCalibrationPINDiodeQEUV; break; case kECCT1: default: fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQECT1; fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQECT1Error*gkCalibrationPINDiodeQECT1Error / gkCalibrationPINDiodeQECT1/gkCalibrationPINDiodeQECT1; break; } *fLog << inf << endl; *fLog << inf << " Mean number of Photons for an Inner Pixel (outside Plexiglass): " << fMeanPhotOutsidePlexiglass << endl; if (fMeanPhotOutsidePlexiglass > 0.) SETBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); else { CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); return kFALSE; } if (fMeanPhotErrOutsidePlexiglass < 0.) { *fLog << warn << "Relative Variance on number of Photons for an Inner Pixel (outside Plexiglass): " << fMeanPhotErrOutsidePlexiglass << endl; CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); return kFALSE; } // Finish calculation of errors -> convert from relative variance to absolute error fMeanPhotErrOutsidePlexiglass = TMath::Sqrt(fMeanPhotErrOutsidePlexiglass); fMeanPhotErrOutsidePlexiglass *= fMeanPhotOutsidePlexiglass; *fLog << inf << " Error on number of Photons for an Inner Pixel (outside Plexiglass): " << fMeanPhotErrOutsidePlexiglass << endl; *fLog << inf << endl; TIter Next(fPixels); MCalibrationPix *pix; while ((pix=(MCalibrationPix*)Next())) { if (pix->IsChargeValid()) { const Float_t charge = pix->GetCharge(); const Float_t ratio = fGeomCam->GetPixRatio(pix->GetPixId()); const Float_t chargeerr = pix->GetErrCharge(); Float_t nphot; Float_t nphoterr; if (ratio == 1.) { nphot = fMeanPhotInsidePlexiglass; nphoterr = fMeanPhotErrInsidePlexiglass; } else { nphot = fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; nphoterr = fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea *fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; nphoterr += fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError *fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError; nphoterr = TMath::Sqrt(nphoterr); } const Float_t conversion = nphot/charge; Float_t conversionerr; conversionerr = nphoterr/charge * nphoterr/charge ; conversionerr += chargeerr/charge * chargeerr/charge * conversion*conversion; if (conversionerr > 0.) conversionerr = TMath::Sqrt(conversionerr); const Float_t conversionsigma = 0.; pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma); if (conversionerr/conversion < 0.1) pix->SetBlindPixelMethodValid(); } } return kTRUE; } Bool_t MCalibrationCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { if (ipx < 0 || !IsPixelFitted(ipx)) return kFALSE; if (!IsNumPhotInsidePlexiglassAvailable()) if (!CalcNumPhotInsidePlexiglass()) return kFALSE; mean = (*this)[ipx].GetMeanConversionBlindPixelMethod(); err = (*this)[ipx].GetErrorConversionBlindPixelMethod(); sigma = (*this)[ipx].GetSigmaConversionBlindPixelMethod(); return kTRUE; } Bool_t MCalibrationCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { if (ipx < 0 || !IsPixelFitted(ipx)) return kFALSE; Float_t conv = (*this)[ipx].GetMeanConversionFFactorMethod(); if (conv < 0.) return kFALSE; mean = conv; err = (*this)[ipx].GetErrorConversionFFactorMethod(); sigma = (*this)[ipx].GetSigmaConversionFFactorMethod(); return kTRUE; } //----------------------------------------------------------------------------------- // // Calculates the conversion factor between the integral of FADCs slices // (as defined in the signal extractor MExtractSignal.cc) // and the number of photons reaching the plexiglass for one Inner Pixel // // FIXME: The PINDiode is still not working and so is the code // Bool_t MCalibrationCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { if (ipx < 0 || !IsPixelFitted(ipx)) return kFALSE; if (!IsNumPhotOutsidePlexiglassAvailable()) if (!CalcNumPhotOutsidePlexiglass()) return kFALSE; mean = (*this)[ipx].GetMeanConversionPINDiodeMethod(); err = (*this)[ipx].GetErrorConversionPINDiodeMethod(); sigma = (*this)[ipx].GetSigmaConversionPINDiodeMethod(); return kFALSE; } //----------------------------------------------------------------------------------- // // Calculates the best combination of the three used methods possible // between the integral of FADCs slices // (as defined in the signal extractor MExtractSignal.cc) // and the number of photons reaching one Inner Pixel. // The procedure is not yet defined. // // FIXME: The PINDiode is still not working and so is the code // Bool_t MCalibrationCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { if (ipx < 0 || !IsPixelFitted(ipx)) return kFALSE; return kFALSE; } void MCalibrationCam::DrawHiLoFits() { if (!fOffsets) fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-600.,400.); if (!fSlopes) fSlopes = new TH1D("mm","Slopes of the HiGain LoGain Fit",100,-2.,2.); if (!fOffvsSlope) fOffvsSlope = new TH2D("aa","Slopes vs Offsets of the HiGain LoGain Fit",100,-600.,400.,100,-2.,2.); TIter Next(fPixels); MCalibrationPix *pix; MHCalibrationPixel *hist; while ((pix=(MCalibrationPix*)Next())) { hist = pix->GetHist(); hist->FitHiGainvsLoGain(); fOffsets->Fill(hist->GetOffset(),1.); fSlopes->Fill(hist->GetSlope(),1.); fOffvsSlope->Fill(hist->GetOffset(),hist->GetSlope(),1.); } TCanvas *c1 = new TCanvas(); c1->Divide(1,3); c1->cd(1); fOffsets->Draw(); gPad->Modified(); gPad->Update(); c1->cd(2); fSlopes->Draw(); gPad->Modified(); gPad->Update(); c1->cd(3); fOffvsSlope->Draw("col1"); gPad->Modified(); gPad->Update(); }