/* ======================================================================== *\ ! ! * ! * 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 09/2003 ! ! Copyright: MAGIC Software Development, 2000-2001 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // MCalibrationCalc // // // // This is a task which calculates the number of photons from the FADC // // time slices. At the moment it integrates simply the FADC values. // // // // Input Containers: // // MRawEvtData // // // // Output Containers: // // MCalibrationCam // // // // // // The class MCalibrationCam hold one entry of type MCalibrationPix for // // every pixel. It is filled in the following way: // // PreParocess: MalibrationCam::InitSize(577) is called which allocates // // memory in an TClonesArray of type MCalibrationPix and // // all pointers to NULL. // // // // Process: The NULL pointer is tested on every pixel via // // MalibrationCam::IsPixelUsed(npix). // // // // In case, IsPixelUsed returns NULL, // // MalibrationCam::AddPixel(npix) is invoked which creates a // // new MCalibrationPix(npix) in the npix's entry // // of the TClonesArray. // // // // Every MCalibrationPix holds a histogram class, // // MHCalibrationPixel which itself hold histograms of type: // // HCharge(npix) (distribution of summed FADC time slice entries) // // HTime(npix) (distribution of position of maximum) // // HChargevsN(npix) (distribution of charges vs. event number. // // // // PostProcess: All histograms HCharge(npix) are fitted to a Gaussian // // All histograms HTime(npix) are fitted to a Gaussian // // The histogram HBlindPixelCharge (blind pixel) is fitted to a single // // PhE fit // // The histogram HBlindPixelTime (blind pixel) is fitted to a Gaussian // // The histograms of the PIN Diode are fitted to Gaussians // // // // Fits can be excluded via the commands: // // MalibrationCam::SetSkipTimeFits() (skip all time fits) // // MalibrationCam::SetSkipBlindPixelFits() (skip all blind pixel fits) // // MalibrationCam::SetSkipPinDiodeFits() (skip all PIN Diode fits) // // // ////////////////////////////////////////////////////////////////////////////// #include "MCalibrationCalc.h" #include "MCalibrationConfig.h" #include "MCalibrationFits.h" #include "MCalibrationCam.h" #include "MCalibrationPix.h" #include "MCalibrationBlindPix.h" #include "MCalibrationPINDiode.h" #include "MPedestalCam.h" #include "MPedestalPix.h" #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MH.h" #include "MRawRunHeader.h" #include "MRawEvtData.h" // MRawEvtData::GetNumPixels #include "MRawEvtPixelIter.h" #include "MTime.h" #include "TMath.h" ClassImp(MCalibrationCalc); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. b is the number of slices before the maximum slice, // a the number of slices behind the maximum slice which is taken as signal. // MCalibrationCalc::MCalibrationCalc(const char *name, const char *title) : fColor(kEBlue) { fName = name ? name : "MCalibrationCalc"; fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam "; AddToBranchList("MRawEvtData.fHiGainPixId"); AddToBranchList("MRawEvtData.fLoGainPixId"); AddToBranchList("MRawEvtData.fHiGainFadcSamples"); AddToBranchList("MRawEvtData.fLoGainFadcSamples"); SETBIT(fFlags, kUseTimeFits); SETBIT(fFlags, kUseBlindPixelFit); SETBIT(fFlags, kUsePinDiodeFit); } // -------------------------------------------------------------------------- // // The PreProcess searches for the following input containers: // - MRawEvtData // - MPedestalCam // // The following output containers are also searched and created if // they were not found: // // - MHCalibrationBlindPixel // - MCalibrationCam // - MTime // Int_t MCalibrationCalc::PreProcess(MParList *pList) { fHistOverFlow = 0; fEvents = 0; fCosmics = 0; fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); if (!fRawEvt) { *fLog << dbginf << "MRawEvtData not found... aborting." << endl; return kFALSE; } const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); if (!runheader) *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; else if (runheader->GetRunType() == kRTMonteCarlo) { return kTRUE; } fCalibrations = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam"); if (!fCalibrations) { *fLog << err << dbginf << "MCalibrationCam could not be created ... aborting." << endl; return kFALSE; } switch (fColor) { case kEBlue: fCalibrations->SetColor(MCalibrationCam::kECBlue); break; case kEGreen: fCalibrations->SetColor(MCalibrationCam::kECGreen); break; case kEUV: fCalibrations->SetColor(MCalibrationCam::kECUV); } fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); if (!fPedestals) { *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl; return kFALSE; } return kTRUE; } // -------------------------------------------------------------------------- // // The ReInit searches for the following input containers: // - MRawRunHeader // Bool_t MCalibrationCalc::ReInit(MParList *pList ) { fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); if (!fRunHeader) { *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; return kFALSE; } fNumHiGainSamples = fRunHeader->GetNumSamplesHiGain(); fNumLoGainSamples = fRunHeader->GetNumSamplesLoGain(); // // FIXME: The next statement simply does not work: // fRawEvt->GetNumPixels() returns always 0 // fCalibrations->InitSize(fRawEvt->GetNumPixels()); // fCalibrations->InitSize(577); for (Int_t i=0;i<577;i++) { MCalibrationPix &pix = (*fCalibrations)[i]; pix.DefinePixId(i); } return kTRUE; } // -------------------------------------------------------------------------- // // Calculate the integral of the FADC time slices and store them as a new // pixel in the MCerPhotEvt container. // Int_t MCalibrationCalc::Process() { fEvents++; Int_t cosmicpix = 0; MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); MCalibrationPINDiode &pindiode = *(fCalibrations->GetPINDiode()); MRawEvtPixelIter pixel(fRawEvt); // Create a first loop to sort out the cosmics ... // // This is a very primitive check for the number of cosmicpixs // The cut will be applied in the fit, but for the blind pixel, // we need to remove this event // // FIXME: In the future need a much more sophisticated one!!! // while (pixel.Next()) { const Int_t pixid = pixel.GetPixelId(); Int_t sum = pixel.GetSumHiGainSamples(); MPedestalPix &ped = (*fPedestals)[pixid]; Float_t pedes = ped.GetPedestal(); Float_t pedrms = ped.GetPedestalRms(); if ((float)sum < (pedes*fNumHiGainSamples)+(1.5*pedrms) ) cosmicpix++; } if (cosmicpix > 50.) { fCosmics++; return kTRUE; } pixel.Reset(); while (pixel.Next()) { UShort_t sat = 0; UShort_t lowgainoverflow = 0; const Int_t pixid = pixel.GetPixelId(); Byte_t *ptr = pixel.GetHiGainSamples(); Byte_t mid = pixel.GetIdxMaxHiGainSample(); UInt_t max = pixel.GetMaxHiGainSample(); Int_t sum = (max > gkSaturationLimit // overflow ? ? sat++, pixel.GetSumLoGainSamples() // take Low Gain : pixel.GetSumHiGainSamples()); // no overflow if (sat) { ptr = pixel.GetLoGainSamples(); max = pixel.GetMaxLoGainSample(); mid = pixel.GetIdxMaxLoGainSample(); // // FIXME: It seems the conversion HiGain LoGain is already // performed in the data?!? // sum = (max > gkSaturationLimit // overflow of LoGain ??? -> GimmeABreak!!! ? lowgainoverflow++, gkLoGainOverFlow // OUCH (Florian was maybe right) : sum*gkConversionHiLo ); // OUFF (Florian was wrong) !! // : sum ); if (lowgainoverflow) { *fLog << err << dbginf << "Warning: Saturation of LoGain reached in pixel: " << pixid << " " << " sum = " << sum << endl; fHistOverFlow++; } } MPedestalPix &ped = (*fPedestals)[pixid]; MCalibrationPix &pix = (*fCalibrations)[pixid]; Float_t pedes = ped.GetPedestal(); Float_t pedrms = ped.GetPedestalRms(); // // FIXME: This is preliminary, we will change to pedestals per slice!!! // Assume pedestals per time slice ==> multiply with number of slices // pedes *= (sat ? fNumLoGainSamples : fNumHiGainSamples ); pedrms *= (sat ? fNumLoGainSamples : fNumHiGainSamples ); Float_t rsum = (float)sum - pedes; switch(pixid) { case gkCalibrationBlindPixelId: // // FIXME: This works only when the blind pixel ID is much larger than // the rest of the pixels (which is the case right now) // if (!blindpixel.FillCharge(sum)) *fLog << warn << "Overflow or Underflow occurred filling Blind Pixel sum = " << sum << endl; if (!blindpixel.FillTime((int)mid)) *fLog << warn << "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl; if (!blindpixel.FillRChargevsTime(rsum,fEvents)) *fLog << warn << "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl; case gkCalibrationPINDiodeId: if (!pindiode.FillCharge(sum)) *fLog << warn << "Overflow or Underflow occurred filling HCharge: means = " << sum << endl; if (!pindiode.FillTime((int)mid)) *fLog << warn << "Overflow or Underflow occurred filling HTime: time = " << (int)mid << endl; if (!pindiode.FillRChargevsTime(rsum,fEvents)) *fLog << warn << "Overflow or Underflow occurred filling HChargevsN: eventnr = " << fEvents << endl; default: if (!pix.FillCharge(sum)) *fLog << warn << "Could not fill Charge of pixel: " << pixid << " signal = " << sum << endl; // // Fill the reduced charge into the control histo for better visibility // if (!pix.FillRChargevsTime(rsum,fEvents)) *fLog << warn << "Could not fill red. Charge vs. EvtNr of pixel: " << pixid << " signal = " << rsum << " event Nr: " << fEvents << endl; if (!pix.FillTime((int)mid)) *fLog << warn << "Could not fill Time of pixel: " << pixid << " time = " << (int)mid << endl; } /* switch(pixid) */ } /* while (pixel.Next()) */ return kTRUE; } Int_t MCalibrationCalc::PostProcess() { *fLog << inf << endl; *fLog << GetDescriptor() << " Cut Histogram Edges" << endl; // // Cut edges to make fits and viewing of the hists easier // fCalibrations->CutEdges(); // // Get pointer to blind pixel // MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel()); *fLog << GetDescriptor() << " Fitting the Blind Pixel" << endl; // // Fit the blind pixel // if (TESTBIT(fFlags,kUseBlindPixelFit)) { if (blindpixel.FitCharge()) if (!fCalibrations->CalcNrPhotInnerPixel()) *fLog << err << dbginf << "Could not calculate Number of photons from the blind pixel " << endl; else *fLog << err << dbginf << "Could not fit the blind pixel " << endl; if (!blindpixel.FitTime()) *fLog << warn << "Could not the Times of the blind pixel " << endl; blindpixel.Draw(); } *fLog << GetDescriptor() << " Fitting the Normal Pixels" << endl; // // loop over the pedestal events and check if we have calibration // for (Int_t pixid=0; pixidGetSize(); pixid++) { MCalibrationPix &pix = (*fCalibrations)[pixid]; const Float_t ped = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples; const Float_t prms = (*fPedestals)[pixid].GetPedestalRms() * TMath::Sqrt((float)fNumHiGainSamples); pix.SetPedestal(ped,prms); if (TESTBIT(fFlags,kUseTimeFits)) pix.FitTime(); pix.FitCharge(); } fCalibrations->SetReadyToSave(); if (GetNumExecutions()==0) return kTRUE; *fLog << endl; *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl; return kTRUE; }