/* ======================================================================== *\ ! ! * ! * 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 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MHCalibrationPulseTimeCam // // Fills the extracted signals of MExtractedSignalCam into the MHCalibrationPix-classes // MHCalibrationPulseTimeHiGainPix and MHCalibrationPulseTimeLoGainPix 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 PulseTime 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) // // Unless more than fNumHiGainSaturationLimit (default: 1%) of the overall FADC // slices show saturation, the following flag is set: // - MCalibrationPulseTimePix::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 MCalibrationPulseTimeCam: // // - MCalibrationPix::SetHiGainSaturation() // - MCalibrationPix::SetHiGainMean() // - MCalibrationPix::SetHiGainMeanErr() // - MCalibrationPix::SetHiGainSigma() // - MCalibrationPix::SetHiGainSigmaErr() // - MCalibrationPix::SetHiGainProb() // - MCalibrationPix::SetHiGainNumPickup() // // 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. // ///////////////////////////////////////////////////////////////////////////// #include "MHCalibrationPulseTimeCam.h" #include "MHCalibrationCam.h" #include #include #include #include #include #include #include #include #include #include #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MHCalibrationPix.h" #include "MCalibrationIntensityCam.h" #include "MCalibrationChargeCam.h" #include "MCalibrationChargePix.h" #include "MGeomCam.h" #include "MGeomPix.h" #include "MBadPixelsIntensityCam.h" #include "MBadPixelsCam.h" #include "MBadPixelsPix.h" #include "MRawEvtData.h" #include "MRawRunHeader.h" #include "MRawEvtPixelIter.h" #include "MExtractedSignalCam.h" #include "MExtractedSignalPix.h" #include "MArrayI.h" #include "MArrayD.h" ClassImp(MHCalibrationPulseTimeCam); using namespace std; const Byte_t MHCalibrationPulseTimeCam::fgSaturationLimit = 254; const Byte_t MHCalibrationPulseTimeCam::fgLowerSignalLimit = 100; const Int_t MHCalibrationPulseTimeCam::fgNumPixelsRequired = 2; const Int_t MHCalibrationPulseTimeCam::fgHiGainNbins = 40; const Axis_t MHCalibrationPulseTimeCam::fgHiGainFirst = -0.5; const Axis_t MHCalibrationPulseTimeCam::fgHiGainLast = 19.5; const Float_t MHCalibrationPulseTimeCam::fgProbLimit = 0.001; const TString MHCalibrationPulseTimeCam::gsHistName = "PulseTime"; const TString MHCalibrationPulseTimeCam::gsHistTitle = "Extracted Times"; const TString MHCalibrationPulseTimeCam::gsHistXTitle = "Time [FADC slices]"; const TString MHCalibrationPulseTimeCam::gsHistYTitle = "Nr. events"; const TString MHCalibrationPulseTimeCam::fgReferenceFile = "mjobs/signalref.rc"; // -------------------------------------------------------------------------- // // Default Constructor. // // Sets: // - all pointers to NULL // // - fNbins to fgHiGainNbins // - fFirst to fgHiGainFirst // - fLast to fgHiGainLast // // - fHistName to gsHistName // - fHistTitle to gsHistTitle // - fHistXTitle to gsHistXTitle // - fHistYTitle to gsHistYTitle // // - fSaturationLimit to fgSaturationLimit // - fLowerSignalLimit to fgLowerSignalLimit // - fNumPixelsRequired to fgNumPixelsRequired // MHCalibrationPulseTimeCam::MHCalibrationPulseTimeCam(const char *name, const char *title) : fBadPixels(NULL) { fName = name ? name : "MHCalibrationPulseTimeCam"; fTitle = title ? title : "Class to fill the extracted pulse times for cosmics "; SetNbins(fgHiGainNbins); SetFirst(fgHiGainFirst); SetLast (fgHiGainLast ); SetProbLimit(fgProbLimit); SetHistName (gsHistName .Data()); SetHistTitle (gsHistTitle .Data()); SetHistXTitle(gsHistXTitle.Data()); SetHistYTitle(gsHistYTitle.Data()); SetReferenceFile(); SetLoGain(kFALSE); SetOscillations(kFALSE); SetSaturationLimit(); SetLowerSignalLimit(); SetNumPixelsRequired(); fInnerRefTime = 5.; fOuterRefTime = 5.; } // -------------------------------------------------------------------------- // // Creates new MHCalibrationPulseTimeCam only with the averaged areas: // the rest has to be retrieved directly, e.g. via: // MHCalibrationPulseTimeCam *cam = MParList::FindObject("MHCalibrationPulseTimeCam"); // - cam->GetAverageSector(5).DrawClone(); // - (*cam)[100].DrawClone() // TObject *MHCalibrationPulseTimeCam::Clone(const char *) const { MHCalibrationPulseTimeCam *cam = new MHCalibrationPulseTimeCam(); // // 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->fReferenceFile = fReferenceFile; if (!IsAverageing()) return cam; const Int_t navhi = fAverageHiGainAreas->GetSize(); for (int i=0; ifAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i); return cam; } // -------------------------------------------------------------------------- // // Gets the pointers to: // - MRawEvtData // Bool_t MHCalibrationPulseTimeCam::SetupHists(const MParList *pList) { fBadPixels = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam"); if (!fBadPixels) { *fLog << warn << GetDescriptor() << "MBadPixelsCam not found... " << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // Gets or creates the pointers to: // - MExtractedSignalCam // - MCalibrationPulseTimeCam or MCalibrationIntensityPulseTimeCam // - MBadPixelsCam // // Initializes the number of used FADC slices from MExtractedSignalCam // into MCalibrationPulseTimeCam and test for changes in that variable // // Calls: // - InitHiGainArrays() // // Sets: // - fSumhiarea to nareas // - fSumloarea to nareas // - fSumhisector to nsectors // - fSumlosector to nsectors // Bool_t MHCalibrationPulseTimeCam::ReInitHists(MParList *pList) { MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam")); if (!signal) { *fLog << err << "MExtractedSignalCam not found... abort." << endl; return kFALSE; } if (!InitCams(pList,"PulseTime")) return kFALSE; const Int_t npixels = fGeom->GetNumPixels(); const Int_t nsectors = fGeom->GetNumSectors(); const Int_t nareas = fGeom->GetNumAreas(); InitHiGainArrays(npixels,nareas,nsectors); fSumhiarea .Set(nareas); fSumhisector.Set(nsectors); return kTRUE; } void MHCalibrationPulseTimeCam::InitHiGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors) { if (fHiGainArray->GetSize()==0) { for (Int_t i=0; iAddAt(new MHCalibrationPix(Form("%sHiGainPix%04d",fHistName.Data(),i), Form("%s High Gain Pixel %4d",fHistTitle.Data(),i)),i); MHCalibrationPix &pix = (*this)[i]; pix.SetNbins(fNbins*2); pix.SetFirst(fFirst); pix.SetLast (fLast); MBadPixelsPix &bad = fIntensBad ? (*fIntensBad)[i] : (*fBadPixels)[i]; InitHists(pix,bad,i); if (fCam) (*fCam)[i].SetPixId(i); } } if (!IsAverageing()) return; if (fAverageHiGainAreas->GetSize()==0) { for (Int_t j=0; jAddAt(new MHCalibrationPix(Form("%sHiGainArea%d",fHistName.Data(),j), Form("%s High Gain Area Idx %d",fHistTitle.Data(),j)),j); MHCalibrationPix &pix = GetAverageHiGainArea(j); pix.SetNbins(fNbins*2); pix.SetFirst(fFirst); pix.SetLast (fLast); InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j); } } if (fAverageHiGainSectors->GetSize()==0) { for (Int_t j=0; jAddAt(new MHCalibrationPix(Form("%sHiGainSector%02d",fHistName.Data(),j), Form("%s High Gain Sector %02d",fHistTitle.Data(),j)),j); MHCalibrationPix &pix = GetAverageHiGainSector(j); pix.SetNbins(fNbins); pix.SetFirst(fFirst); pix.SetLast (fLast); InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j); } } } // -------------------------------------------------------------------------- // // 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 PulseTimes histograms (MHGausEvents::FillHistAndArray()) with: // - MExtractedSignalPix::GetExtractedSignalHiGain(); // - MExtractedSignalPix::GetExtractedSignalLoGain(); // Bool_t MHCalibrationPulseTimeCam::FillHists(const MParContainer *par, const Stat_t w) { MRawEvtData *data = (MRawEvtData*)par; if (!data) { *fLog << err << "No argument in MHCalibrationPulseTimeCam::Fill... abort." << endl; return kFALSE; } const UInt_t nareas = fGeom->GetNumAreas(); const UInt_t nsectors = fGeom->GetNumSectors(); fSumhiarea .Reset(); fSumhisector.Reset(); fAverageAreaNum.Reset(); fAverageSectorNum.Reset(); MRawEvtPixelIter pixel(data); while (pixel.Next()) { const Int_t i = pixel.GetPixelId(); if (fBadPixels) { MBadPixelsPix &bad = (*fBadPixels)[i]; if (bad.IsUnsuitable()) continue; } Byte_t *start = pixel.GetHiGainSamples(); Byte_t *end = start + pixel.GetNumHiGainSamples(); Byte_t *p = start; Byte_t max = 0; Int_t maxpos = 0; while (p < end) { if ((*p > max) && (*p < fSaturationLimit)) { max = *p; maxpos = p-start-1; } p++; } start = pixel.GetLoGainSamples(); end = start + pixel.GetNumLoGainSamples(); p = start; while (p < end) { if ((*p > max) && (*p < fSaturationLimit)) { max = *p; maxpos = p-start+pixel.GetNumHiGainSamples() - 1; } p++; } if (max < fLowerSignalLimit) continue; const Float_t time = (Float_t)maxpos; (*this)[i].FillHist(time); const Int_t aidx = (*fGeom)[i].GetAidx(); const Int_t sector = (*fGeom)[i].GetSector(); fSumhiarea[aidx] += time; fSumhisector[sector] += time; fAverageAreaNum[aidx]++; fAverageSectorNum[sector]++; } for (UInt_t j=0; j fNumPixelsRequired) { if (IsOscillations()) GetAverageHiGainArea(j).FillHistAndArray(fSumhiarea[j]/npix); else GetAverageHiGainArea(j).FillHist(fSumhiarea[j]/npix); } } for (UInt_t j=0; j 0) { if (IsOscillations()) GetAverageHiGainSector(j).FillHistAndArray(fSumhisector [j]/npix); else GetAverageHiGainSector(j).FillHist(fSumhisector [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 MCalibrationPulseTimePix'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 MHCalibrationPulseTimeCam::FinalizeHists() { *fLog << endl; MCalibrationCam &calcam = *(fIntensCam ? fIntensCam->GetCam() : fCam); // // Perform the fitting for the High Gain (done in MHCalibrationCam) // for (Int_t i=0; iGetSize(); i++) { MHCalibrationPix &hist = (*this)[i]; if (hist.IsExcluded()) continue; MCalibrationPix &pix = calcam[i]; CalcHists(hist,pix); } if (!IsAverageing()) return kTRUE; for (Int_t j=0; jGetSize(); j++) { MHCalibrationPix &hist = GetAverageHiGainArea(j); MCalibrationPix &pix = calcam.GetAverageArea(j); CalcHists(hist,pix); } for (Int_t j=0; jGetSize(); j++) { MHCalibrationPix &hist = GetAverageHiGainSector(j); MCalibrationPix &pix = calcam.GetAverageSector(j); CalcHists(hist,pix); } return kTRUE; } void MHCalibrationPulseTimeCam::CalcHists(MHCalibrationPix &hist, MCalibrationPix &pix) const { if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow()) { *fLog << warn << GetDescriptor() << ": Only over- or underflow in " << hist.GetName() << endl; return; } hist.BypassFit(); pix.SetHiGainMean ( hist.GetMean() ); pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() ); pix.SetHiGainRms ( hist.GetHistRms() ); pix.SetHiGainSigma ( hist.GetSigma() ); pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() ); if (IsDebug()) { *fLog << dbginf << GetDescriptor() << ": ID " << GetName() << " "<GetSize(); if (nareas == 0) return; TString option(opt); option.ToLower(); 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); MHCalibrationPix &hipix = GetAverageHiGainArea(i); DrawDataCheckPixel(hipix,i ? fOuterRefTime : fInnerRefTime); } } // ----------------------------------------------------------------------------- // // 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 MHCalibrationPulseTimeCam::DrawDataCheckPixel(MHCalibrationPix &pix, const Float_t refline) { gPad->SetTicks(); TH1F *hist = pix.GetHGausHist(); TH1F *null = new TH1F("Null",hist->GetTitle(),100, pix.GetFirst() > 0. ? pix.GetFirst() : 0., pix.GetLast() > pix.GetFirst() ? pix.GetLast() : pix.GetFirst()*2.); null->SetMaximum(1.1*hist->GetMaximum()); null->SetDirectory(NULL); null->SetBit(kCanDelete); null->SetStats(kFALSE); // // set the labels bigger // TAxis *xaxe = null->GetXaxis(); TAxis *yaxe = null->GetYaxis(); xaxe->CenterTitle(); yaxe->CenterTitle(); xaxe->SetTitleSize(0.06); yaxe->SetTitleSize(0.076); xaxe->SetTitleOffset(0.6); yaxe->SetTitleOffset(0.65); xaxe->SetLabelSize(0.06); yaxe->SetLabelSize(0.06); xaxe->SetTitle(hist->GetXaxis()->GetTitle()); yaxe->SetTitle(hist->GetYaxis()->GetTitle()); null->Draw(); hist->Draw("same"); DisplayRefLines(null,refline); return; } void MHCalibrationPulseTimeCam::DisplayRefLines(const TH1F *hist, const Float_t refline) const { TGraph *gr = new TGraph(2); gr->SetPoint(0,refline,0.); gr->SetPoint(1,refline,hist->GetMaximum()); gr->SetBit(kCanDelete); gr->SetLineColor(106); gr->SetLineStyle(2); gr->SetLineWidth(3); gr->Draw("L"); TLegend *leg = new TLegend(0.8,0.35,0.99,0.65); leg->SetBit(kCanDelete); leg->AddEntry(gr,"Reference","l"); leg->Draw(); } Int_t MHCalibrationPulseTimeCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { Bool_t rc = kFALSE; if (MHCalibrationCam::ReadEnv(env,prefix,print)) rc = kTRUE; if (IsEnvDefined(env, prefix, "SaturationLimit", print)) { SetSaturationLimit(GetEnvValue(env, prefix, "SaturationLimit", fSaturationLimit)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "LowerSignalLimit", print)) { SetLowerSignalLimit(GetEnvValue(env,prefix,"LowerSignalLimit",fLowerSignalLimit)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "NumPixelsRequired", print)) { SetNumPixelsRequired(GetEnvValue(env,prefix,"NumPixelsRequired",fNumPixelsRequired)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "ReferenceFile", print)) { SetReferenceFile(GetEnvValue(env,prefix,"ReferenceFile",fReferenceFile.Data())); rc = kTRUE; } TEnv refenv(fReferenceFile); fInnerRefTime = refenv.GetValue("InnerRefTime",fInnerRefTime); fOuterRefTime = refenv.GetValue("OuterRefTime",fOuterRefTime); return rc; }