Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3083) +++ /trunk/MagicSoft/Mars/Changelog (revision 3084) @@ -29,4 +29,20 @@ spectrum + * mfilter/MCosmics.[h,cc] + * mfilter/Makefile + * mfilter/FilterLinkDef.h + - new filter removing cosmics, the same as in MCalibrationCalc + where it was removed now. + Call: MCosmics cosmics; + MContinue cont(&cosmics); + tlist.AddToList(&cont); + + * mcalib/MCalibrationCalc.[h,cc] + - removed cosmics rejection from there + + * macros/calibration.C + * mjobs/MJCalibration.cc + - changed cosmics rejection to the filter algorithm + 2004/02/09: Markus Gaug Index: /trunk/MagicSoft/Mars/macros/calibration.C =================================================================== --- /trunk/MagicSoft/Mars/macros/calibration.C (revision 3083) +++ /trunk/MagicSoft/Mars/macros/calibration.C (revision 3084) @@ -219,4 +219,11 @@ MArrivalTimeCalc timecalc; MCalibrationCalc calcalc; + + // + // Apply a filter against cosmics + // (was directly in MCalibrationCalc in earlier versions) + // + MFCosmics cosmics; + MContinue cont(&cosmics); // @@ -241,5 +248,5 @@ // (This is a preliminary feature) // - //calcalc.ExcludePixelsFromAsciiFile("badpixels_all.dat"); + calcalc.ExcludePixelsFromAsciiFile("badpixels.dat"); // @@ -248,10 +255,4 @@ // // calcalc.SkipBlindPixelFit(); - - // - // In case, you want to skip the cosmics rejection - // (NOT RECOMMENDED!!!) - // - // calcalc.SkipCosmicsRejection(); // @@ -276,4 +277,5 @@ // // tlist2.AddToList(&timecalc); + tlist2.AddToList(&cont); tlist2.AddToList(&calcalc); @@ -301,5 +303,5 @@ // just one example how to get the plots of individual pixels // - calcam[356].DrawClone(); + calcam[543].DrawClone(); // Create histograms to display @@ -312,8 +314,8 @@ MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)"); MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)"); - MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Nr. of Photons inside plexiglass (Blind Pixel Method)"); + MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)"); MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)"); MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)"); - MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Nr. of Photons outside plexiglass (PIN Diode Method)"); + MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)"); MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor (PIN Diode Method)"); MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)"); @@ -416,9 +418,9 @@ disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); - disp9.SetYTitle("Nr. Photons [1]"); + disp9.SetYTitle("Photon flux [ph/mm^2]"); disp10.SetYTitle("Conversion Factor [Phot/FADC Count]"); disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); - disp12.SetYTitle("Nr. Photons [1]"); + disp12.SetYTitle("Photon flux [ph/mm^2]"); disp13.SetYTitle("Conversion Factor [Phot/FADC Count]"); disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); @@ -465,9 +467,9 @@ // F-Factor Method TCanvas &c4 = display->AddTab("F-Factor"); - c4.Divide(3,3); - - CamDraw(c4, disp6,calcam,1, 3 , 2); - CamDraw(c4, disp7,calcam,2, 3 , 2); - CamDraw(c4, disp8,calcam,3, 3 , 2); + c4.Divide(2,3); + + CamDraw(c4, disp6,calcam,1, 2 , 2); + CamDraw(c4, disp7,calcam,2, 2 , 2); + // CamDraw(c4, disp8,calcam,3, 3 , 2); // Blind Pixel Method Index: /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc (revision 3083) +++ /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc (revision 3084) @@ -43,7 +43,5 @@ // Optionally exclude pixels from calibration // -// Process: Optionally, a cut on cosmics can be performed -// -// Every MCalibrationPix holds a histogram class, +// Process: Every MCalibrationPix holds a histogram class, // MHCalibrationPixel which itself hold histograms of type: // HCharge(npix) (distribution of summed FADC time slice @@ -144,5 +142,4 @@ SETBIT(fFlags, kUseBlindPixelFit); - SETBIT(fFlags, kUseCosmicsRejection); SETBIT(fFlags, kUseQualityChecks); SETBIT(fFlags, kHiLoGainCalibration); @@ -156,5 +153,4 @@ fEvents = 0; - fCosmics = 0; fNumHiGainSamples = 0; @@ -391,60 +387,6 @@ MRawEvtPixelIter pixel(fRawEvt); - // - // Perform cosmics cut - // - if (TESTBIT(fFlags,kUseCosmicsRejection)) - { - - Int_t cosmicpix = 0; - - // - // 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 UInt_t pixid = pixel.GetPixelId(); - - MExtractedSignalPix &sig = (*fSignals)[pixid]; - MPedestalPix &ped = (*fPedestals)[pixid]; - Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples; - Float_t sumhi = sig.GetExtractedSignalHiGain(); - - // - // We consider a pixel as presumably due to cosmics - // if its sum of FADC slices is lower than 3 pedestal RMS - // - if (sumhi < 3.*pedrms ) - cosmicpix++; - } - - // - // If the camera contains more than 230 - // (this is the number of outer pixels plus about 50 inner ones) - // presumed pixels due to cosmics, then the event is discarted. - // This procedure is more or less equivalent to keeping only events - // with at least 350 pixels with high signals. - // - if (cosmicpix > 230.) - { - fCosmics++; - return kCONTINUE; - } - - pixel.Reset(); - } - fEvents++; - // Float_t referencetime = 0.; - // // Create a (second) loop to do fill the calibration histograms @@ -624,14 +566,8 @@ { *fLog << err << GetDescriptor() - << ": This run contains only cosmics or pedestals, " - << "cannot find events with more than 350 illuminated pixels. " << endl; + << ": This run contains no calibration data! " << endl; return kFALSE; } - if (fEvents < fCosmics) - *fLog << warn << GetDescriptor() - << ": WARNING: Run contains more cosmics or pedestals than calibration events " << endl; - - *fLog << inf << GetDescriptor() << ": Cut Histogram Edges" << endl; @@ -784,12 +720,7 @@ fCalibrations->SetReadyToSave(); - if (GetNumExecutions()==0) - return kTRUE; - - *fLog << inf << endl; - *fLog << dec << setfill(' ') << fCosmics << " Events presumably cosmics" << endl; - return kTRUE; } + Bool_t MCalibrationCalc::CalcSignalBlindPixel(Byte_t *ptr, Float_t &signal, const Float_t ped) const Index: /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.h (revision 3083) +++ /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.h (revision 3084) @@ -54,5 +54,4 @@ Int_t fEvents; // Number of events - Int_t fCosmics; // Number of events due to supposed cosmics Byte_t fNumHiGainSamples; @@ -67,5 +66,5 @@ enum { kUseBlindPixelFit, kUsePinDiodeFit, - kUseCosmicsRejection, kUseQualityChecks, + kUseQualityChecks, kHiLoGainCalibration, kHiGainOverFlow, kLoGainOverFlow }; @@ -96,6 +95,4 @@ void SkipPinDiodeFit(Bool_t b=kTRUE) {b ? CLRBIT(fFlags, kUsePinDiodeFit) : SETBIT(fFlags, kUsePinDiodeFit);} - void SkipCosmicsRejection(Bool_t b=kTRUE) - {b ? CLRBIT(fFlags, kUseCosmicsRejection) : SETBIT(fFlags, kUseCosmicsRejection);} void SkipQualityChecks(Bool_t b=kTRUE) {b ? CLRBIT(fFlags, kUseQualityChecks) : SETBIT(fFlags, kUseQualityChecks);} Index: /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 3083) +++ /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 3084) @@ -26,4 +26,5 @@ #pragma link C++ class MFSelFinal+; #pragma link C++ class MFSoftwareTrigger+; +#pragma link C++ class MFCosmics+; #pragma link C++ class MFEnergySlope+; Index: /trunk/MagicSoft/Mars/mfilter/MFCosmics.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFCosmics.cc (revision 3084) +++ /trunk/MagicSoft/Mars/mfilter/MFCosmics.cc (revision 3084) @@ -0,0 +1,219 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MFCosmics +// +// Filter to reject cosmics by the criterion that at least +// fMaxEmptyPixels pixels have values of lower than 3 Pedestal RMS. +// fMaxEmptyPixels is set to 230 by default which is slightly higher +// than the number of outer pixels in MAGIC (for the case that +// the outer pixels have some defect). +// +// ProProcess: Search for MPedestalCam, MExtractedSignalCam +// +// ReInit: Initialize number of used FADC slices +// +// Process: if fMaxEmptyPixels pixels lower than 3 pedRMS, +// the event is supposed to +// be a cosmic and kContinue is returned +// +// +// Input Containers: +// MRawEvtData +// MPedestalCam +// MExtractedSignalCam +// +// Output Containers: +// +////////////////////////////////////////////////////////////////////////////// +#include "MFCosmics.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MRawEvtPixelIter.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" + +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" + +ClassImp(MFCosmics); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MFCosmics::MFCosmics(const char *name, const char *title) + : fPedestals(NULL), fSignals(NULL), + fRawEvt(NULL), fMaxEmptyPixels(230) +{ + + fName = name ? name : "MFCosmics"; + fTitle = title ? title : "Filter to reject cosmics"; +} + +// -------------------------------------------------------------------------- +// +// The PreProcess searches for the following input containers: +// - MRawEvtData +// - MPedestalCam +// - MExtractedSignalCam +// +Int_t MFCosmics::PreProcess(MParList *pList) +{ + + fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); + if (!fRawEvt) + { + *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl; + return kFALSE; + } + + fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPedestals) + { + *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl; + return kFALSE; + } + + fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); + if (!fSignals) + { + *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl; + return kFALSE; + } + + memset(fCut, 0, sizeof(fCut)); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// The ReInit searches the following input containers for information: +// - MExtractedSignalCam +// +Bool_t MFCosmics::ReInit(MParList *pList ) +{ + + fSqrtHiGainSamples = TMath::Sqrt((Float_t) fSignals->GetNumUsedHiGainFADCSlices()); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Retrieve the integral of the FADC time slices and compare them to the +// pedestal values. +// +Int_t MFCosmics::Process() +{ + + fResult = CosmicsRejection(); + + fCut[fResult ? 0 : 1]++; + return kTRUE; +} + +// --------------------------------------------------------- +// +// Cosmics rejection: +// +// Requiring less than fMaxEmptyPixels pixels to have values +// lower than 3 Pedestal RMS. +// +// fMaxEmptyPixels is set to 230 by default which is slightly higher +// than the number of outer pixels in MAGIC (for the case that +// the outer pixels have some defect). +// +Bool_t MFCosmics::CosmicsRejection() +{ + + MRawEvtPixelIter pixel(fRawEvt); + + Int_t cosmicpix = 0; + + // + // Create a first loop to sort out the cosmics ... + // + while (pixel.Next()) + { + + const UInt_t idx = pixel.GetPixelId(); + + MExtractedSignalPix &sig = (*fSignals)[idx]; + MPedestalPix &ped = (*fPedestals)[idx]; + const Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples; + const Float_t sumhi = sig.GetExtractedSignalHiGain(); + + // + // We consider a pixel as presumably due to cosmics + // if its sum of FADC slices is lower than 3 pedestal RMS + // + if (sumhi < 3.*pedrms ) + cosmicpix++; + } + + // + // If the camera contains more than fMaxEmptyPixels + // presumed pixels due to cosmics, then the event is discarted. + // + if (cosmicpix > fMaxEmptyPixels) + return kTRUE; + + return kFALSE; +} + +Int_t MFCosmics::PostProcess() +{ + + if (GetNumExecutions()==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + + *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ; + *fLog << (int)(fCut[1]*100/GetNumExecutions()) ; + *fLog << "%) Evts skipped due to: Cosmics Rejection applied " ; + *fLog << " (with fMaxEmptyPixels = " << fMaxEmptyPixels << ")" << endl; + + *fLog << " " << fCut[0] << " (" << (int)(fCut[0]*100/GetNumExecutions()) ; + *fLog << "%) Evts survived the cosmics rejection!" << endl; + *fLog << endl; + + return kTRUE; +} + Index: /trunk/MagicSoft/Mars/mfilter/MFCosmics.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFCosmics.h (revision 3084) +++ /trunk/MagicSoft/Mars/mfilter/MFCosmics.h (revision 3084) @@ -0,0 +1,54 @@ +#ifndef MARS_MFCosmics +#define MARS_MFCosmics + +///////////////////////////////////////////////////////////////////////////// +// // +// MFCosmics // +// // +// Filter against cosmics (used especially by the calibration // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MFilter +#include "MFilter.h" +#endif + +class MRawEvtData; + +class MPedestalCam; +class MExtractedSignalCam; + +class MFCosmics : public MFilter +{ +private: + + MPedestalCam *fPedestals; // Pedestals of all pixels in the camera + MExtractedSignalCam *fSignals; // Calibration events of all pixels in the camera + + MRawEvtData *fRawEvt; // raw event data (time slices) + + Int_t fCut[2]; + Bool_t fResult; + Int_t fMaxEmptyPixels; // Maximum number of empty pixels before declaring event as cosmic + Float_t fSqrtHiGainSamples; // Square root of the number of used Hi-Gain Samples + + Bool_t ReInit(MParList *pList); + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + Bool_t CosmicsRejection(); + + Bool_t IsExpressionTrue() const { return fResult; } + +public: + + MFCosmics(const char *name=NULL, const char *title=NULL); + + void SetMaxEmptyPixels(const Int_t n) { fMaxEmptyPixels = n; } + Int_t GetMaxEmptyPixels() const { return fMaxEmptyPixels; } + + ClassDef(MFCosmics, 0) // Filter to perform a cosmics rejection +}; + +#endif Index: /trunk/MagicSoft/Mars/mfilter/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mfilter/Makefile (revision 3083) +++ /trunk/MagicSoft/Mars/mfilter/Makefile (revision 3084) @@ -52,4 +52,5 @@ MFSelFinal.cc \ MFSoftwareTrigger.cc \ + MFCosmics.cc \ MFEnergySlope.cc Index: /trunk/MagicSoft/Mars/mtools/MFFT.cc =================================================================== --- /trunk/MagicSoft/Mars/mtools/MFFT.cc (revision 3083) +++ /trunk/MagicSoft/Mars/mtools/MFFT.cc (revision 3084) @@ -68,5 +68,5 @@ #include "MFFT.h" -#include "TMath.h" +#include #include "MLog.h" @@ -930,3 +930,2 @@ } - Index: /trunk/MagicSoft/Mars/mtools/MFFT.h =================================================================== --- /trunk/MagicSoft/Mars/mtools/MFFT.h (revision 3083) +++ /trunk/MagicSoft/Mars/mtools/MFFT.h (revision 3084) @@ -76,5 +76,5 @@ TArrayF* PowerSpectrumDensity(const TArrayF *array); TArrayD* PowerSpectrumDensity(const TArrayD *array); - + TArrayF* RealFunctionSpectrum(const TArrayF *data);