/* ======================================================================== *\ ! $Name: not supported by cvs2svn $:$Id: MHCalibrationChargeCam.cc,v 1.62 2008-06-30 09:36:36 tbretz Exp $ ! -------------------------------------------------------------------------- ! ! * ! * 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 ! Author(s): Thomas Bretz ! ! Copyright: MAGIC Software Development, 2000-2007 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MHCalibrationChargeCam // // Fills the extracted signals of MExtractedSignalCam into the MHCalibrationPix-classes // MHCalibrationChargeHiGainPix and MHCalibrationChargeLoGainPix for every: // // - Pixel, stored in the TOrdCollection's MHCalibrationCam::fHiGainArray and // MHCalibrationCam::fLoGainArray // // - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), // stored in the TOrdCollection's MHCalibrationCam::fAverageHiGainAreas and // MHCalibrationCam::fAverageLoGainAreas // // - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), // stored in the TOrdCollection's MHCalibrationCam::fAverageHiGainSectors and // MHCalibrationCam::fAverageLoGainSectors // // Every signal is taken from MExtractedSignalCam and filled into a histogram and // an array, in order to perform a Fourier analysis (see MHGausEvents). // The signals are moreover averaged on an event-by-event basis and written into // the corresponding average pixels. // // Additionally, the (FADC slice) position of the maximum is stored in an Absolute // Arrival Time histogram. This histogram serves for a rough cross-check if the // signal does not lie at or outside the edges of the extraction window. // // The Charge histograms are fitted to a Gaussian, mean and sigma with its errors // and the fit probability are extracted. If none of these values are NaN's and // if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%), // the fit is declared valid. // Otherwise, the fit is repeated within ranges of the previous mean // +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit()) // In case this does not make the fit valid, the histogram means and RMS's are // taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set: // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) or // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) and // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) // // Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas // from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup) // // If more than fNumHiGainSaturationLimit (default: 15%) of the overall FADC // slices show saturation, the following flag is set: // - MCalibrationChargePix::SetHiGainSaturation(); // In that case, the calibration constants are derived from the low-gain results. // // If more than fNumLoGainSaturationLimit (default: 1%) of the overall // low-gain FADC slices saturate, the following flags are set: // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation ) and // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnsuitableRun ) // // The class also fills arrays with the signal vs. event number, creates a fourier // spectrum and investigates if the projected fourier components follow an exponential // distribution. In case that the probability of the exponential fit is less than // MHGausEvents::fProbLimit (default: 0.5%), the following flags are set: // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) or // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) and // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) // // This same procedure is performed for the average pixels. // // The following results are written into MCalibrationChargeCam: // // - MCalibrationPix::SetHiGainSaturation() // - MCalibrationPix::SetHiGainMean() // - MCalibrationPix::SetHiGainMeanErr() // - MCalibrationPix::SetHiGainSigma() // - MCalibrationPix::SetHiGainSigmaErr() // - MCalibrationPix::SetHiGainProb() // - MCalibrationPix::SetHiGainNumPickup() // // - MCalibrationPix::SetLoGainMean() // - MCalibrationPix::SetLoGainMeanErr() // - MCalibrationPix::SetLoGainSigma() // - MCalibrationPix::SetLoGainSigmaErr() // - MCalibrationPix::SetLoGainProb() // - MCalibrationPix::SetLoGainNumPickup() // // - MCalibrationChargePix::SetAbsTimeMean() // - MCalibrationChargePix::SetAbsTimeRms() // // For all averaged areas, the fitted sigma is multiplied with the square root of // the number involved pixels in order to be able to compare it to the average of // sigmas in the camera. // // Class Version 2: // + Float_t fNumLoGainBlackoutLimit; // Rel. amount blackout logain events until pixel is declared unsuitable // // Class Version 3: // - Byte_t fFirstHiGain; // First used slice High Gain // - Byte_t fLastHiGain; // Last used slice High Gain // - Byte_t fFirstLoGain; // First used slice Low Gain // - Byte_t fLastLoGain; // Last used slice Low Gain // ///////////////////////////////////////////////////////////////////////////// #include "MHCalibrationChargeCam.h" #include #include #include #include #include #include #include #include #include #include #include #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MHCalibrationChargePix.h" #include "MHCalibrationPix.h" #include "MCalibrationChargeCam.h" #include "MCalibrationChargePix.h" #include "MGeomCam.h" #include "MGeomPix.h" #include "MBadPixelsCam.h" #include "MBadPixelsPix.h" #include "MRawRunHeader.h" #include "MRawEvtPixelIter.h" #include "MPedestalSubtractedEvt.h" #include "MExtractedSignalCam.h" #include "MExtractedSignalPix.h" ClassImp(MHCalibrationChargeCam); using namespace std; const Int_t MHCalibrationChargeCam::fgChargeHiGainNbins = 500; const Axis_t MHCalibrationChargeCam::fgChargeHiGainFirst = -98.; const Axis_t MHCalibrationChargeCam::fgChargeHiGainLast = 1902.; const Int_t MHCalibrationChargeCam::fgChargeLoGainNbins = 500; const Axis_t MHCalibrationChargeCam::fgChargeLoGainFirst = -99.; const Axis_t MHCalibrationChargeCam::fgChargeLoGainLast = 901.; const Float_t MHCalibrationChargeCam::fgProbLimit = 0.00000001; const TString MHCalibrationChargeCam::gsHistName = "Charge"; const TString MHCalibrationChargeCam::gsHistTitle = "Signals"; const TString MHCalibrationChargeCam::gsHistXTitle = "Signal [FADC counts]"; const TString MHCalibrationChargeCam::gsHistYTitle = "Nr. events"; const TString MHCalibrationChargeCam::gsAbsHistName = "AbsTime"; const TString MHCalibrationChargeCam::gsAbsHistTitle = "Abs. Arr. Times"; const TString MHCalibrationChargeCam::gsAbsHistXTitle = "Time [FADC slices]"; const TString MHCalibrationChargeCam::gsAbsHistYTitle = "Nr. events"; const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.15; const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005; const Float_t MHCalibrationChargeCam::fgNumLoGainBlackoutLimit = 0.05; const Float_t MHCalibrationChargeCam::fgLoGainBlackoutLimit = 3.5; const Float_t MHCalibrationChargeCam::fgLoGainPickupLimit = 3.5; const Float_t MHCalibrationChargeCam::fgTimeLowerLimit = 1.; const Float_t MHCalibrationChargeCam::fgTimeUpperLimit = 3.; const TString MHCalibrationChargeCam::fgReferenceFile = "mjobs/calibrationref.rc"; // -------------------------------------------------------------------------- // // Default Constructor. // // Sets: // - all pointers to NULL // // Initializes: // - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit // - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit // - fTimeLowerLimit to fgTimeLowerLimit // - fTimeUpperLimit to fgTimeUpperLimit // - fNumLoGainBlackoutLimit to fgNumLoGainBlackoutLimit // // - fNbins to fgChargeHiGainNbins // - fFirst to fgChargeHiGainFirst // - fLast to fgChargeHiGainLast // // - fLoGainNbins to fgChargeLoGainNbins // - fLoGainFirst to fgChargeLoGainFirst // - fLoGainLast to fgChargeLoGainLast // // - fHistName to gsHistName // - fHistTitle to gsHistTitle // - fHistXTitle to gsHistXTitle // - fHistYTitle to gsHistYTitle // // - fAbsHistName to gsAbsHistName // - fAbsHistTitle to gsAbsHistTitle // - fAbsHistXTitle to gsAbsHistXTitle // - fAbsHistYTitle to gsAbsHistYTitle // MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title) : fRawEvt(NULL), fSignal(NULL) { fName = name ? name : "MHCalibrationChargeCam"; fTitle = title ? title : "Class to fill the calibration histograms "; SetNumHiGainSaturationLimit(fgNumHiGainSaturationLimit); SetNumLoGainSaturationLimit(fgNumLoGainSaturationLimit); SetNumLoGainBlackoutLimit (fgNumLoGainBlackoutLimit); SetTimeLowerLimit(); SetTimeUpperLimit(); SetBinning(fgChargeHiGainNbins, fgChargeHiGainFirst, fgChargeHiGainLast ); SetBinningLoGain(fgChargeLoGainNbins, fgChargeLoGainFirst, fgChargeLoGainLast ); SetProbLimit(fgProbLimit); SetHistName (gsHistName .Data()); SetHistTitle (gsHistTitle .Data()); SetHistXTitle(gsHistXTitle.Data()); SetHistYTitle(gsHistYTitle.Data()); SetAbsHistName (gsAbsHistName .Data()); SetAbsHistTitle (gsAbsHistTitle .Data()); SetAbsHistXTitle(gsAbsHistXTitle.Data()); SetAbsHistYTitle(gsAbsHistYTitle.Data()); SetReferenceFile(); fInnerRefCharge = 278.; fOuterRefCharge = 282.; } // -------------------------------------------------------------------------- // // Creates new MHCalibrationChargeCam only with the averaged areas: // the rest has to be retrieved directly, e.g. via: // MHCalibrationChargeCam *cam = MParList::FindObject("MHCalibrationChargeCam"); // - cam->GetAverageSector(5).DrawClone(); // - (*cam)[100].DrawClone() // TObject *MHCalibrationChargeCam::Clone(const char *) const { MHCalibrationChargeCam *cam = new MHCalibrationChargeCam(); // // Copy the data members // cam->fColor = fColor; cam->fRunNumbers = fRunNumbers; cam->fPulserFrequency = fPulserFrequency; cam->fFlags = fFlags; cam->fNbins = fNbins; cam->fFirst = fFirst; cam->fLast = fLast; cam->fLoGainNbins = fLoGainNbins; cam->fLoGainFirst = fLoGainFirst; cam->fLoGainLast = fLoGainLast; cam->fReferenceFile = fReferenceFile; cam->fInnerRefCharge = fInnerRefCharge; cam->fOuterRefCharge = fOuterRefCharge; // // Copy the MArrays // cam->fAverageAreaRelSigma = fAverageAreaRelSigma; cam->fAverageAreaRelSigmaVar = fAverageAreaRelSigmaVar; cam->fAverageAreaSat = fAverageAreaSat; cam->fAverageAreaSigma = fAverageAreaSigma; cam->fAverageAreaSigmaVar = fAverageAreaSigmaVar; cam->fAverageAreaNum = fAverageAreaNum; cam->fAverageSectorNum = fAverageSectorNum; if (!IsAverageing()) return cam; const Int_t navhi = fAverageHiGainAreas->GetSize(); for (int i=0; ifAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i); if (IsLoGain()) { const Int_t navlo = fAverageLoGainAreas->GetSize(); for (int i=0; ifAverageLoGainAreas->AddAt(GetAverageLoGainArea(i).Clone(),i); } return cam; } // -------------------------------------------------------------------------- // // Gets the pointers to: // - MRawEvtData // Bool_t MHCalibrationChargeCam::SetupHists(const MParList *pList) { fRawEvt = (MPedestalSubtractedEvt*)pList->FindObject("MPedestalSubtractedEvt"); if (!fRawEvt) { *fLog << err << dbginf << "MPedestalSubtractedEvt not found... aborting." << endl; return kFALSE; } return kTRUE; } // -------------------------------------------------------------------------- // // Gets or creates the pointers to: // - MExtractedSignalCam // - MCalibrationChargeCam // - MBadPixelsCam // // Initializes the number of used FADC slices from MExtractedSignalCam // into MCalibrationChargeCam and test for changes in that variable // // Calls: // - InitHiGainArrays() // - InitLoGainArrays() // // Sets: // - fSumhiarea to nareas // - fSumloarea to nareas // - fTimehiarea to nareas // - fTimeloarea to nareas // - fSumhisector to nsectors // - fSumlosector to nsectors // - fTimehisector to nsectors // - fTimelosector to nsectors // - fSathiarea to nareas // - fSatloarea to nareas // - fSathisector to nsectors // - fSatlosector to nsectors // Bool_t MHCalibrationChargeCam::ReInitHists(MParList *pList) { fSignal = (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam")); if (!fSignal) { *fLog << err << "MExtractedSignalCam not found... abort." << endl; return kFALSE; } if (!InitCams(pList,"Charge")) return kFALSE; const Int_t npixels = fGeom->GetNumPixels(); const Int_t nsectors = fGeom->GetNumSectors(); const Int_t nareas = fGeom->GetNumAreas(); // // In case of the intense blue, double the range // if (fGeom->InheritsFrom("MGeomCamMagic")) if ( fColor == MCalibrationCam::kBLUE) SetBinning(fLoGainNbins, fLoGainFirst, 2.*fLoGainLast - fLoGainFirst); InitHiGainArrays(npixels,nareas,nsectors); InitLoGainArrays(npixels,nareas,nsectors); fSumhiarea .Set(nareas); fSumloarea .Set(nareas); fTimehiarea .Set(nareas); fTimeloarea .Set(nareas); fSumhisector.Set(nsectors); fSumlosector.Set(nsectors); fTimehisector.Set(nsectors); fTimelosector.Set(nsectors); fSathiarea .Set(nareas); fSatloarea .Set(nareas); fSathisector.Set(nsectors); fSatlosector.Set(nsectors); return kTRUE; } // -------------------------------------------------------------------------- // // Retrieve: // - fRunHeader->GetNumSamplesHiGain(); // // Initializes the High Gain Arrays: // // - For every entry in the expanded arrays: // * Initialize an MHCalibrationChargePix // * Set Binning from fNbins, fFirst and fLast // * Set Binning of Abs Times histogram from fAbsNbins, fAbsFirst and fAbsLast // * Set Histgram names and titles from fHistName and fHistTitle // * Set Abs Times Histgram names and titles from fAbsHistName and fAbsHistTitle // * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle // * Set X-axis and Y-axis titles of Abs Times Histogram from fAbsHistXTitle and fAbsHistYTitle // * Call InitHists // // void MHCalibrationChargeCam::InitHiGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors) { TH1F *h; const Int_t higainsamples = fRunHeader->GetNumSamplesHiGain(); if (fHiGainArray->GetSize()==0) { for (Int_t i=0; iAddAt(new MHCalibrationChargePix(Form("%sHiGainPix%04d",fHistName.Data(),i), Form("%s High Gain Pixel%04d",fHistTitle.Data(),i)),i); MHCalibrationChargePix &pix = (MHCalibrationChargePix&)(*this)[i]; pix.SetBinning(fNbins, fFirst, fLast); pix.SetProbLimit(fProbLimit); pix.SetBinningAbsTime(higainsamples, -0.5, higainsamples-0.5); InitHists(pix,(*fBadPixels)[i], i); h = pix.GetHAbsTime(); h->SetName (Form("H%sHiGainPix%04d",fAbsHistName.Data(),i)); h->SetTitle(Form("%s High Gain Pixel %04d",fAbsHistTitle.Data(),i)); h->SetXTitle(fAbsHistXTitle.Data()); h->SetYTitle(fAbsHistYTitle.Data()); } } if (fAverageHiGainAreas->GetSize()==0) { for (Int_t j=0; jAddAt(new MHCalibrationChargePix(Form("%sHiGainArea%d",fHistName.Data(),j), Form("%s High Gain Area Idx %d",fHistTitle.Data(),j)),j); MHCalibrationChargePix &pix = (MHCalibrationChargePix&)GetAverageHiGainArea(j); pix.SetBinning(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas), fFirst, fLast); pix.SetBinningAbsTime(higainsamples, -0.5, higainsamples-0.5); InitHists(pix, fCam->GetAverageBadArea(j),j); h = pix.GetHAbsTime(); h->SetName (Form("H%sHiGainArea%d",fAbsHistName.Data(),j)); h->SetTitle(Form("%s%s%d",fAbsHistTitle.Data(), " averaged on event-by-event basis High Gain Area Idx ",j)); h->SetXTitle(fAbsHistXTitle.Data()); h->SetYTitle(fAbsHistYTitle.Data()); } } if (fAverageHiGainSectors->GetSize()==0) { for (Int_t j=0; jAddAt(new MHCalibrationChargePix(Form("%sHiGainSector%02d",fHistName.Data(),j), Form("%s High Gain Sector %02d",fHistTitle.Data(),j)),j); MHCalibrationChargePix &pix = (MHCalibrationChargePix&)GetAverageHiGainSector(j); pix.SetBinning(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas), fFirst, fLast); pix.SetBinningAbsTime(higainsamples, -0.5, higainsamples-0.5); InitHists(pix, fCam->GetAverageBadSector(j),j); h = pix.GetHAbsTime(); h->SetName (Form("H%sHiGainSector%02d",fAbsHistName.Data(),j)); h->SetTitle(Form("%s%s%02d",fAbsHistTitle.Data(), " averaged on event-by-event basis High Gain Area Sector ",j)); h->SetXTitle(fAbsHistXTitle.Data()); h->SetYTitle(fAbsHistYTitle.Data()); } } } //-------------------------------------------------------------------------------------- // // Return, if IsLoGain() is kFALSE // // Retrieve: // - fRunHeader->GetNumSamplesHiGain(); // // Initializes the Low Gain Arrays: // // - For every entry in the expanded arrays: // * Initialize an MHCalibrationChargePix // * Set Binning from fNbins, fFirst and fLast // * Set Binning of Abs Times histogram from fAbsNbins, fAbsFirst and fAbsLast // * Set Histgram names and titles from fHistName and fHistTitle // * Set Abs Times Histgram names and titles from fAbsHistName and fAbsHistTitle // * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle // * Set X-axis and Y-axis titles of Abs Times Histogram from fAbsHistXTitle and fAbsHistYTitle // * Call InitHists // void MHCalibrationChargeCam::InitLoGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors) { if (!IsLoGain()) return; const Int_t logainsamples = fRunHeader->GetNumSamplesLoGain(); TH1F *h; if (fLoGainArray->GetSize()==0 ) { for (Int_t i=0; iAddAt(new MHCalibrationChargePix(Form("%sLoGainPix%04d",fHistName.Data(),i), Form("%s Low Gain Pixel %04d",fHistTitle.Data(),i)),i); MHCalibrationChargePix &pix = (MHCalibrationChargePix&)(*this)(i); pix.SetBinning(fLoGainNbins, fLoGainFirst, fLoGainLast); pix.SetProbLimit(fProbLimit); pix.SetBinningAbsTime(logainsamples, -0.5, logainsamples-0.5); pix.SetPickupLimit(fgLoGainPickupLimit); pix.SetBlackoutLimit(fgLoGainBlackoutLimit); InitHists(pix,(*fBadPixels)[i], i); h = pix.GetHAbsTime(); h->SetName (Form("H%sLoGainPix%04d",fAbsHistName.Data(),i)); h->SetTitle(Form("%s Low Gain Pixel %04d",fAbsHistTitle.Data(),i)); h->SetXTitle(fAbsHistXTitle.Data()); h->SetYTitle(fAbsHistYTitle.Data()); } } if (fAverageLoGainAreas->GetSize()==0) { for (Int_t j=0; jAddAt(new MHCalibrationChargePix(Form("%sLoGainArea%d",fHistName.Data(),j), Form("%s Low Gain Area Idx %d",fHistTitle.Data(),j)),j); MHCalibrationChargePix &pix = (MHCalibrationChargePix&)GetAverageLoGainArea(j); pix.SetBinning(fLoGainNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas), fLoGainFirst, fLoGainLast); pix.SetBinningAbsTime(logainsamples, -0.5, logainsamples-0.5); InitHists(pix, fCam->GetAverageBadArea(j),j); h = pix.GetHAbsTime(); h->SetName (Form("H%sLoGainArea%02d",fAbsHistName.Data(),j)); h->SetTitle(Form("%s%s%02d",fAbsHistTitle.Data(), " averaged on event-by-event basis Low Gain Area Idx ",j)); h->SetXTitle(fAbsHistXTitle.Data()); h->SetYTitle(fAbsHistYTitle.Data()); } } if (fAverageLoGainSectors->GetSize()==0 && IsLoGain()) { for (Int_t j=0; jAddAt(new MHCalibrationChargePix(Form("%sLoGainSector%02d",fHistName.Data(),j), Form("%s Low Gain Sector %02d",fHistTitle.Data(),j)),j); MHCalibrationChargePix &pix = (MHCalibrationChargePix&)GetAverageLoGainSector(j); pix.SetBinning(fLoGainNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas), fLoGainFirst, fLoGainLast); pix.SetBinningAbsTime(logainsamples, -0.5, logainsamples-0.5); InitHists(pix, fCam->GetAverageBadSector(j),j); h = pix.GetHAbsTime(); h->SetName (Form("H%sLoGainSector%02d",fAbsHistName.Data(),j)); h->SetTitle(Form("%s%s%02d",fAbsHistTitle.Data(), " averaged on event-by-event basis Low Gain Area Sector ",j)); h->SetXTitle(fAbsHistXTitle.Data()); h->SetYTitle(fAbsHistYTitle.Data()); } } } // -------------------------------------------------------------------------- // // Retrieves from MExtractedSignalCam: // - first used LoGain FADC slice // // Retrieves from MGeomCam: // - number of pixels // - number of pixel areas // - number of sectors // // For all TOrdCollection's (including the averaged ones), the following steps are performed: // // 1) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with: // - MExtractedSignalPix::GetExtractedSignalHiGain(); // - MExtractedSignalPix::GetExtractedSignalLoGain(); // // 2) Set number of saturated slices (MHCalibrationChargePix::AddSaturated()) with: // - MExtractedSignalPix::IsHiGainSaturated(); // - MExtractedSignalPix::GetNumLoGainSaturated(); // // 3) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with: // - MRawEvtPixelIter::GetIdxMaxHiGainSample(); // - MRawEvtPixelIter::GetIdxMaxLoGainSample(first slice); // Bool_t MHCalibrationChargeCam::FillHists(const MParContainer *par, const Stat_t w) { MExtractedSignalCam *signal = (MExtractedSignalCam*)par; if (!signal) { *fLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl; return kFALSE; } const UInt_t npixels = fGeom->GetNumPixels(); const UInt_t nareas = fGeom->GetNumAreas(); const UInt_t nsectors = fGeom->GetNumSectors(); const Int_t hifirst = fSignal->GetFirstUsedSliceHiGain(); const Int_t hilast = fSignal->GetLastUsedSliceHiGain(); const Int_t lofirst = fSignal->GetFirstUsedSliceLoGain(); const Int_t lolast = fSignal->GetLastUsedSliceLoGain(); fSumhiarea .Reset(); fSumloarea .Reset(); fTimehiarea .Reset(); fTimeloarea .Reset(); fSumhisector.Reset(); fSumlosector.Reset(); fTimehisector.Reset(); fTimelosector.Reset(); fSathiarea .Reset(); fSatloarea .Reset(); fSathisector.Reset(); fSatlosector.Reset(); for (UInt_t i=0; iGetNumPixels(); const Int_t nhi = fRunHeader->GetNumSamplesHiGain(); for (int pixid=0; pixidGetMaxPos(pixid, hifirst, hilast)+hifirst; histhi.FillAbsTime(timehi); const Int_t aidx = (*fGeom)[pixid].GetAidx(); const Int_t sector = (*fGeom)[pixid].GetSector(); fTimehiarea [aidx] += timehi; fTimehisector[sector] += timehi; if (IsLoGain()) { MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(pixid); const Float_t timelo = fRawEvt->GetMaxPos(pixid, nhi+lofirst, nhi+lolast)+lofirst; histlo.FillAbsTime(timelo); fTimeloarea[aidx] += timelo; fTimelosector[sector] += timelo; } } for (UInt_t j=0; j 0.5 ? 1 : 0); hipix.FillAbsTime (fTimehiarea[j]/npix); if (IsLoGain()) { MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainArea(j); if (IsOscillations()) lopix.FillHistAndArray(fSumloarea [j]/npix); else lopix.FillHist(fSumloarea [j]/npix); lopix.AddSaturated ((Float_t)fSatloarea [j]/npix > 0.5 ? 1 : 0); lopix.FillAbsTime (fTimeloarea[j]/npix); } } for (UInt_t j=0; j 0.5 ? 1 : 0); hipix.FillAbsTime (fTimehisector[j]/npix); if (IsLoGain()) { MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainSector(j); if (IsOscillations()) lopix.FillHistAndArray(fSumlosector [j]/npix); else lopix.FillHist(fSumlosector [j]/npix); lopix.AddSaturated ((Float_t)fSatlosector[j]/npix > 0.5 ? 1 : 0); lopix.FillAbsTime (fTimelosector[j]/npix); } } return kTRUE; } // -------------------------------------------------------------------------- // // For all TOrdCollection's (including the averaged ones), the following steps are performed: // // 1) Returns if the pixel is excluded. // 2) Tests saturation. In case yes, set the flag: MCalibrationPix::SetHiGainSaturation() // or the flag: MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturated ) // 3) Store the absolute arrival times in the MCalibrationChargePix's. If flag // MCalibrationPix::IsHiGainSaturation() is set, the Low-Gain arrival times are stored, // otherwise the Hi-Gain ones. // 4) Calls to MHCalibrationCam::FitHiGainArrays() and MCalibrationCam::FitLoGainArrays() // with the flags: // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) // Bool_t MHCalibrationChargeCam::FinalizeHists() { *fLog << endl; TH1F *h = NULL; const Int_t hifirst = fSignal->GetFirstUsedSliceHiGain(); const Int_t hilast = fSignal->GetLastUsedSliceHiGain(); const Int_t lofirst = fSignal->GetFirstUsedSliceLoGain(); const Int_t lolast = fSignal->GetLastUsedSliceLoGain(); for (Int_t i=0; iGetSize(); i++) { MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i]; if (histhi.IsExcluded()) continue; MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i]; const Int_t numsat = histhi.GetSaturated(); pix.SetNumSaturated(numsat); if (numsat > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) { pix.SetHiGainSaturation(); if (IsOscillations()) histhi.CreateFourierSpectrum(); continue; } MBadPixelsPix &bad = (*fBadPixels)[i]; h = histhi.GetHGausHist(); Stat_t overflow = h->GetBinContent(h->GetNbinsX()+1); if (overflow > fOverflowLimit*histhi.GetHGausHist()->GetEntries()) { *fLog << warn << "HiGain Hist-overflow " << overflow << " times in: " << histhi.GetName() << " (w/o saturation!) " << endl; bad.SetUncalibrated( MBadPixelsPix::kHiGainOverFlow ); } overflow = h->GetBinContent(0); if (overflow > fOverflowLimit*histhi.GetHGausHist()->GetEntries()) { *fLog << warn << "HiGain Hist-underflow " << overflow << " times in pix: " << histhi.GetName() << " (w/o saturation!) " << endl; bad.SetUncalibrated( MBadPixelsPix::kHiGainOverFlow ); } FinalizeAbsTimes(histhi, pix, bad, hifirst, hilast); } if (IsLoGain()) for (Int_t i=0; iGetSize(); i++) { MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i); if (histlo.IsExcluded()) continue; MBadPixelsPix &bad = (*fBadPixels)[i]; if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) { *fLog << warn << "Pixel " << setw(4) << i << ": More than " << Form("%.1f%%", fNumLoGainSaturationLimit*100) << " lo-gains saturated." << endl; bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); if (IsOscillations()) histlo.CreateFourierSpectrum(); continue; } MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i] ; if (!pix.IsHiGainSaturation()) continue; h = histlo.GetHGausHist(); Stat_t overflow = h->GetBinContent(h->GetNbinsX()+1); if (overflow > fOverflowLimit*histlo.GetHGausHist()->GetEntries()) { *fLog << warn << "LoGain Hist-overflow " << overflow << " times in: " << histlo.GetName() << " (w/o saturation!) " << endl; bad.SetUncalibrated( MBadPixelsPix::kLoGainOverFlow ); } overflow = h->GetBinContent(0); if (overflow > fOverflowLimit*histlo.GetHGausHist()->GetEntries()) { *fLog << warn << "LoGain Hist-underflow " << overflow << " times in: " << histlo.GetName() << " (w/o saturation!) " << endl; bad.SetUncalibrated( MBadPixelsPix::kLoGainOverFlow ); } FinalizeAbsTimes(histlo, pix, bad, lofirst, lolast); } for (Int_t j=0; jGetSize(); j++) { MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainArea(j); MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(j); if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) { pix.SetHiGainSaturation(); if (IsOscillations()) histhi.CreateFourierSpectrum(); continue; } MBadPixelsPix &bad = fCam->GetAverageBadArea(j); FinalizeAbsTimes(histhi, pix, bad, hifirst, hilast); } if (IsLoGain()) for (Int_t j=0; jGetSize(); j++) { MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainArea(j); if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) { *fLog << warn << "Area " << setw(4) << j << ": More than " << Form("%.1f%%", fNumLoGainSaturationLimit*100) << " lo-gains saturated." << endl; if (IsOscillations()) histlo.CreateFourierSpectrum(); continue; } MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(j); if (pix.IsHiGainSaturation()) { MBadPixelsPix &bad = fCam->GetAverageBadArea(j); FinalizeAbsTimes(histlo, pix, bad, lofirst, lolast); } } for (Int_t j=0; jGetSize(); j++) { MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainSector(j); MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(j); if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) { pix.SetHiGainSaturation(); if (IsOscillations()) histhi.CreateFourierSpectrum(); continue; } MBadPixelsPix &bad = fCam->GetAverageBadSector(j); FinalizeAbsTimes(histhi, pix, bad, hifirst, hilast); } if (IsLoGain()) for (Int_t j=0; jGetSize(); j++) { MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j); MBadPixelsPix &bad = fCam->GetAverageBadSector(j); if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) { *fLog << warn << "Sector " << setw(4) << j << ": More than " << Form("%.1f%%", fNumLoGainSaturationLimit*100) << " lo-gains saturated." << endl; bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); if (IsOscillations()) histlo.CreateFourierSpectrum(); continue; } MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(j); if (pix.IsHiGainSaturation()) FinalizeAbsTimes(histlo, pix, bad, lofirst, lolast); } // // Perform the fitting for the High Gain (done in MHCalibrationCam) // FitHiGainArrays(*fCam, *fBadPixels, MBadPixelsPix::kHiGainNotFitted, MBadPixelsPix::kHiGainOscillating); // // Perform the fitting for the Low Gain (done in MHCalibrationCam) // if (IsLoGain()) FitLoGainArrays(*fCam, *fBadPixels, MBadPixelsPix::kLoGainNotFitted, MBadPixelsPix::kLoGainOscillating); // // Check for pixels which have the high-gain saturated, but the low-gain // switch not applied in sufficient cases. Have to exclude these pixels, // although they not abnormal. They simply cannot be calibrated... // for (Int_t i=0; iGetSize(); i++) { MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i); if (histlo.IsExcluded()) continue; MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i]; if (!pix.IsHiGainSaturation()) continue; // // Now,treat only low-gain calibrated pixels: // const Double_t lim = fNumLoGainBlackoutLimit*histlo.GetHGausHist()->GetEntries(); if (histlo.GetBlackout() <= lim) continue; MBadPixelsPix &bad = (*fBadPixels)[i]; bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun); bad.SetUncalibrated(MBadPixelsPix::kLoGainBlackout); } return kTRUE; } // -------------------------------------------------------------------------------- // // Fill the absolute time results into MCalibrationChargePix // // Check absolute time validity: // - Mean arrival time is at least fTimeLowerLimit slices from the lower edge // - Mean arrival time is at least fUpperLimit slices from the upper edge // void MHCalibrationChargeCam::FinalizeAbsTimes(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad, Int_t first, Int_t last) { const Float_t mean = hist.GetAbsTimeMean(); const Float_t rms = hist.GetAbsTimeRms(); pix.SetAbsTimeMean(mean); pix.SetAbsTimeRms(rms); const Float_t lowerlimit = (Float_t)first+1;// + fTimeLowerLimit; const Float_t upperlimit = (Float_t)last -1;// - fTimeUpperLimit; // FIXME: instead of checking whether the maximum is in the first or // last extracted slice we should check whether the extractor // was able to properly extract the signal!!! if (meanupperlimit) { *fLog << warn << hist.GetName() << ": Mean Arr.Time: " << Form("%4.1f > %4.1f (last used fadc sl-1)", mean, upperlimit) << endl; bad.SetUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins); } } // -------------------------------------------------------------------------- // // Sets all pixels to MBadPixelsPix::kUnsuitableRun, if following flags are set: // - MBadPixelsPix::kLoGainSaturation // // Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set: // - if MBadPixelsPix::kHiGainNotFitted and !MCalibrationPix::IsHiGainSaturation() // - if MBadPixelsPix::kHiGainOscillating and !MCalibrationPix::IsHiGainSaturation() // - if MBadPixelsPix::kLoGainNotFitted and MCalibrationPix::IsLoGainSaturation() // - if MBadPixelsPix::kLoGainOscillating and MCalibrationPix::IsLoGainSaturation() // void MHCalibrationChargeCam::FinalizeBadPixels() { for (Int_t i=0; iGetSize(); i++) { MBadPixelsPix &bad = (*fBadPixels)[i]; MCalibrationPix &pix = (*fCam)[i]; if (bad.IsUncalibrated(MBadPixelsPix::kHiGainNotFitted)) if (!pix.IsHiGainSaturation()) bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun); if (bad.IsUncalibrated(MBadPixelsPix::kLoGainNotFitted)) if (pix.IsHiGainSaturation()) bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun); if (bad.IsUncalibrated(MBadPixelsPix::kLoGainSaturation)) bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun); if (IsOscillations()) { if (bad.IsUncalibrated(MBadPixelsPix::kHiGainOscillating)) if (!pix.IsHiGainSaturation()) bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun); if (bad.IsUncalibrated(MBadPixelsPix::kLoGainOscillating)) if (pix.IsHiGainSaturation()) bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun); } } } // -------------------------------------------------------------------------- // // Calls MHCalibrationPix::DrawClone() for pixel idx // void MHCalibrationChargeCam::DrawPixelContent(Int_t idx) const { (*this)[idx].DrawClone(); } // ----------------------------------------------------------------------------- // // Default draw: // // Displays the averaged areas, both High Gain and Low Gain // // Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options // void MHCalibrationChargeCam::Draw(const Option_t *opt) { const Int_t nareas = fAverageHiGainAreas->GetSize(); if (nareas == 0) return; TString option(opt); option.ToLower(); if (!option.Contains("datacheck")) { MHCalibrationCam::Draw(opt); return; } // // From here on , the datacheck - Draw // TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this); pad->SetBorderMode(0); pad->Divide(1,nareas); // // Loop over inner and outer pixels // for (Int_t i=0; icd(i+1); MHCalibrationChargePix &hipix = (MHCalibrationChargePix&)GetAverageHiGainArea(i); // // Ask for Hi-Gain saturation // if (hipix.GetSaturated() > fNumHiGainSaturationLimit*hipix.GetHGausHist()->GetEntries() && IsLoGain()) { if (i>=fAverageLoGainAreas->GetSize()) continue; MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainArea(i); DrawDataCheckPixel(lopix,i ? fOuterRefCharge : fInnerRefCharge); } else DrawDataCheckPixel(hipix,i ? fOuterRefCharge : fInnerRefCharge); } } // ----------------------------------------------------------------------------- // // Draw the average pixel for the datacheck: // // Displays the averaged areas, both High Gain and Low Gain // // Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options // void MHCalibrationChargeCam::DrawDataCheckPixel(MHCalibrationChargePix &pix, const Float_t refline) { if (pix.IsEmpty()) return; TVirtualPad *pad = gPad; pad->Divide(1,2, 1e-10, 1e-10); pad->cd(1); gPad->SetBorderMode(0); gPad->SetTicks(); if (!pix.IsEmpty() && !pix.IsOnlyOverflow() && !pix.IsOnlyUnderflow()) gPad->SetLogy(); TH1F *hist = pix.GetHGausHist(); TAxis *xaxe = hist->GetXaxis(); TAxis *yaxe = hist->GetYaxis(); xaxe->CenterTitle(); yaxe->CenterTitle(); xaxe->SetTitleSize(0.07); yaxe->SetTitleSize(0.07); xaxe->SetTitleOffset(0.7); yaxe->SetTitleOffset(0.55); xaxe->SetLabelSize(0.06); yaxe->SetLabelSize(0.06); xaxe->SetTitle(hist->GetXaxis()->GetTitle()); yaxe->SetTitle(hist->GetYaxis()->GetTitle()); xaxe->SetRange(hist->GetMaximumBin()-50, hist->GetMaximumBin()+50); gStyle->SetOptFit(); hist->Draw(); TF1 *fit = pix.GetFGausFit(); if (fit) { switch (fColor) { case MCalibrationCam::kGREEN: fit->SetLineColor(kGreen); break; case MCalibrationCam::kBLUE: fit->SetLineColor(kBlue); break; case MCalibrationCam::kUV: fit->SetLineColor(51); break; case MCalibrationCam::kCT1: fit->SetLineColor(6); break; default: fit->SetLineColor(kRed); } fit->Draw("same"); } DisplayRefLines(hist,refline); pad->cd(2); gPad->SetBorderMode(0); gPad->SetTicks(); pix.CreateGraphEvents(); TGraph *gr = pix.GetGraphEvents(); if (gr) { TH1F *null2 = gr->GetHistogram(); null2->SetMinimum(pix.GetMean()-10.*pix.GetSigma()); null2->SetMaximum(pix.GetMean()+10.*pix.GetSigma()); null2->SetStats(kFALSE); // // set the labels bigger // TAxis *xaxe2 = null2->GetXaxis(); TAxis *yaxe2 = null2->GetYaxis(); xaxe2->CenterTitle(); yaxe2->CenterTitle(); xaxe2->SetTitleSize(0.07); yaxe2->SetTitleSize(0.07); xaxe2->SetTitleOffset(0.7); yaxe2->SetTitleOffset(0.55); xaxe2->SetLabelSize(0.06); yaxe2->SetLabelSize(0.06); xaxe2->SetRangeUser(0, pix.GetEvents()->GetSize()/pix.GetEventFrequency()); } pix.DrawEvents(); // newpad->cd(3); // pix.DrawPowerSpectrum(*newpad,4); } void MHCalibrationChargeCam::DisplayRefLines(const TH1F *hist, const Float_t refline) const { TLine *line = new TLine(refline, 0, refline, hist->GetMaximum()); line->SetLineColor(51); line->SetLineStyle(2); line->SetLineWidth(3); line->SetBit(kCanDelete); line->Draw(); TLegend *leg = new TLegend(0.8,0.01,0.99,0.15); leg->AddEntry(line, "10 Leds UV", "l"); leg->SetBit(kCanDelete); leg->Draw(); } Int_t MHCalibrationChargeCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { Int_t rc = MHCalibrationCam::ReadEnv(env,prefix,print); if (rc==kERROR) return kERROR; if (IsEnvDefined(env, prefix, "HiGainNbins", print)) { fNbins = GetEnvValue(env, prefix, "HiGainNbins", fNbins); rc = kTRUE; } if (IsEnvDefined(env, prefix, "HiGainFirst", print)) { fFirst = GetEnvValue(env, prefix, "HiGainFirst", fFirst); rc = kTRUE; } if (IsEnvDefined(env, prefix, "HiGainLast", print)) { fLast = GetEnvValue(env, prefix, "HiGainLast", fLast); rc = kTRUE; } if (IsEnvDefined(env, prefix, "LoGainNbins", print)) { fLoGainNbins = GetEnvValue(env, prefix, "LoGainNbins", fLoGainNbins); rc = kTRUE; } if (IsEnvDefined(env, prefix, "LoGainFirst", print)) { fLoGainFirst = GetEnvValue(env, prefix, "LoGainFirst", fLoGainFirst); rc = kTRUE; } if (IsEnvDefined(env, prefix, "LoGainLast", print)) { fLoGainLast = GetEnvValue(env, prefix, "LoGainLast", fLoGainLast); rc = kTRUE; } if (IsEnvDefined(env, prefix, "TimeLowerLimit", print)) { SetTimeLowerLimit(GetEnvValue(env, prefix, "TimeLowerLimit", fTimeLowerLimit)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "TimeUpperLimit", print)) { SetTimeUpperLimit(GetEnvValue(env, prefix, "TimeUpperLimit", fTimeUpperLimit)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "ReferenceFile", print)) { SetReferenceFile(GetEnvValue(env,prefix,"ReferenceFile",fReferenceFile.Data())); rc = kTRUE; } if (IsEnvDefined(env, prefix, "NumHiGainSaturationLimit", print)) { SetNumHiGainSaturationLimit(GetEnvValue(env, prefix, "NumHiGainSaturationLimit", fNumHiGainSaturationLimit)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "NumLoGainSaturationLimit", print)) { SetNumLoGainSaturationLimit(GetEnvValue(env, prefix, "NumLoGainSaturationLimit", fNumLoGainSaturationLimit)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "NumLoGainBlackoutLimit", print)) { SetNumLoGainBlackoutLimit(GetEnvValue(env, prefix, "NumLoGainBlackoutLimit", fNumLoGainBlackoutLimit)); rc = kTRUE; } TEnv refenv(fReferenceFile); fInnerRefCharge = refenv.GetValue("InnerRefCharge",fInnerRefCharge); fOuterRefCharge = refenv.GetValue("OuterRefCharge",fOuterRefCharge); return rc; }