Index: /trunk/MagicSoft/Mars/mtemp/MApplyPadding.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MApplyPadding.cc (revision 1677) +++ /trunk/MagicSoft/Mars/mtemp/MApplyPadding.cc (revision 1677) @@ -0,0 +1,287 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 10/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MApplyPadding // +// // +// This task applies padding to a given Sigmabar target value. // +// The task checks whether the data stream it is applied to has to be // +// padded or not and if so, Gaussian noise with the difference in Sigma // +// is produced and added to the particular event. Number of photons and // +// error on this number is altered. // +// // +// There are three ways to define the sigmabar value to which the events // +// are padded: // +// // +// 1) Set a fixed level with SetTargetLevel // +// // +// 2) Give a TH1D* which defines the Sigmabar as function of Theta // +// with SetDefiningHistogram // +// The given histogram should have the same binning in Theta as // +// is used in the analysis // +// // +// 3) Do nothing, then PreProcess will try to read in a (workaround!) // +// propriety format ASCII database for the CT1 test. // +// the name of this file is set by SetDatabaseFile // +// Better know what you are doing or use methods 1 or 2. // +// // +// This implementation is still PRELIMINARY and requires some workarounds // +// put in SPECIFICALLY FOR THE CT1 TESTS, since a database to access is // +// missing. It is not the FINAL MAGIC VERSION. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MApplyPadding.h" + +#include + +#include "TH1.h" +#include "TH2.h" +#include "TRandom.h" + +#include "MBinning.h" +#include "MSigmabar.h" +#include "MMcEvt.hxx" +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MGeomCam.h" +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" +#include "MPedestalPix.h" + +ClassImp(MApplyPadding); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MApplyPadding::MApplyPadding(const char *name, const char *title) : fRunType(0), fGroup(0), fUseHistogram(kTRUE), fFixedSigmabar(0.0) +{ + fName = name ? name : "MApplyPadding"; + fTitle = title ? title : "Task to apply padding"; + Print(); +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MApplyPadding::~MApplyPadding() +{ + //nothing yet +} + +// -------------------------------------------------------------------------- +// +// You can provide a TH1D* histogram containing the target Sigmabar in +// bins of theta. Be sure to use the same binning as for the analysis +// +Bool_t MApplyPadding::SetDefiningHistogram(TH1D *histo) +{ + fHSigmabarMax = histo; + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// check if MEvtHeader exists in the Parameter list already. +// if not create one and add them to the list +// +Bool_t MApplyPadding::PreProcess(MParList *pList) +{ + fRnd = new TRandom3(0); + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << dbginf << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar"); + if (!fSigmabar) + { + *fLog << dbginf << "MSigmabar not found... aborting." << endl; + return kFALSE; + } + + // Get Theta Binning + MBinning* binstheta = (MBinning*)pList->FindObject("BinningTheta"); + if (!binstheta) + { + *fLog << err << dbginf << "BinningTheta not found... aborting." << endl; + return kFALSE; + } + + // Create fSigmabarMax histogram + // (only if no fixed Sigmabar target value or a histogram have already been + // provided) + if ((!fUseHistogram) && (fHSigmabarMax==NULL)) { + + fHSigmabarMax = new TH1D(); + fHSigmabarMax->SetNameTitle("fHSigmabarMax","Sigmabarmax for this analysis"); + TAxis &x = *fHSigmabarMax->GetXaxis(); +#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03) + TString xtitle = x.GetTitle(); +#endif + fHSigmabarMax->SetBins(binstheta->GetNumBins(), 0, 1); + // Set the binning of the current histogram to the binning + // in one of the two given histograms + x.Set(binstheta->GetNumBins(), binstheta->GetEdges()); +#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03) + x.SetTitle(xtitle); +#endif + + // ------------------------------------------------- + // read in SigmabarParams + // workaround--proprietary file format--CT1test only BEGIN + // ------------------------------------------------- + + FILE *f; + if( !(f =fopen(fDatabaseFilename, "r")) ) { + *fLog << err << dbginf << "Database file " << fDatabaseFilename << "was not found... (specify with MApplyPadding::SetDatabaseFile) aborting." << endl; + return kFALSE; + } + char line[80]; + Float_t sigmabarMin, sigmabarMax, thetaMin, thetaMax, ra, dec2; + Int_t type, group, mjd, nr; + while ( fgets(line, sizeof(line), f) != NULL) { + if ((line[0]!='#')) { + sscanf(line,"%d %d %f %f %d %d %f %f %f %f",&type, &group, &ra, &dec2, &mjd, &nr, &sigmabarMin,&sigmabarMax,&thetaMin,&thetaMax); + if ((group==fGroup)||(type==1)) //selected ON group or OFF + { + // find out which bin(s) we have to look at + for (Int_t i=fHSigmabarMax->GetXaxis()->FindBin(thetaMin); + iGetXaxis()->FindBin(thetaMax)+1; i++) + if (sigmabarMax > fHSigmabarMax->GetBinContent(i)) + fHSigmabarMax->SetBinContent(i, sigmabarMax); + } + } + }//while + + // workaround--proprietary file format--CT1test only END + + // fHSigmabarMax->DrawClone(); + // fTest = new TH2D("fTest", "Test if padding works", 201, -0.05, 2.05, 201, -0.05, 2.05); + + } //!fUseHistogram + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculate Sigmabar for current event +// Then apply padding +// +// 1) have current event +// 2) get sigmabar(theta) +// 3) pad event +// +Bool_t MApplyPadding::Process() +{ + // Calculate sigmabar of event + fSigmabar->Calc(*fCam, *fPed); + Double_t mySig = fSigmabar->GetSigmabar(); + + // Get sigmabar which we have to pad to + Double_t otherSig; + if (fUseHistogram) { + Int_t binNumber = fHSigmabarMax->GetXaxis()->FindBin(fMcEvt->GetTheta()*kRad2Deg); + otherSig = fHSigmabarMax->GetBinContent(binNumber); + } else { + otherSig = fFixedSigmabar; + } + + // Determine quadratic difference other-mine + Double_t quadraticDiff = otherSig*otherSig - mySig*mySig; + + if (quadraticDiff < 0) { + *fLog << err << dbginf << "Event has higher Sigmabar="< 0 + if (quadraticDiff > 0) { + + MPedestalCam newPed; + newPed.InitSize(fPed->GetSize()); + + const UInt_t npix = fPed->GetSize(); // Total number of pixels + for (UInt_t i=0; ioperator[](i); + if (!pix.IsPixelUsed()) + continue; + pix.SetNumPhotons(pix.GetNumPhotons() + + sqrt(quadraticDiff)* + fRnd->Gaus(0.0, 1.0)/ + fCam->GetPixRatio(pix.GetPixId()) + ); + // error: add sigma of padded noise quadratically + Double_t error = pix.GetErrorPhot(); + pix.SetErrorPhot(sqrt(error*error + quadraticDiff)); + + MPedestalPix ppix = fPed->operator[](i); + MPedestalPix npix; + npix.SetSigma(sqrt(ppix.GetSigma()*ppix.GetSigma() + quadraticDiff)); + newPed[i]=npix; + } //for + // Calculate Sigmabar again and crosscheck + fSigmabar->Calc(*fCam, newPed); + //mySig = fSigmabar->GetSigmabar(); + // fTest->Fill(otherSig,mySig); + return kTRUE; + } //if + return kFALSE; +} + +Bool_t MApplyPadding::PostProcess() +{ + // fTest->DrawClone(); + return kTRUE; +} Index: unk/MagicSoft/Mars/mtemp/MSigmabar.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MSigmabar.cc (revision 1676) +++ (revision ) @@ -1,150 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Robert Wagner 10/2002 -! -! Copyright: MAGIC Software Development, 2002 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MSigmabar // -// // -// This is the storage container to hold information about the mean sigma // -// (aka Sigmabar) of all pedestals // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MSigmabar.h" - -#include - -#include "MLog.h" -#include "MLogManip.h" - -#include "MGeomCam.h" -#include "MPedestalCam.h" -#include "MGeomPix.h" -#include "MPedestalPix.h" - -ClassImp(MSigmabar); - -MSigmabar::MSigmabar(const char *name, const char *title) -{ - fName = name ? name : "MSigmabar"; - fTitle = title ? title : "Storage container for Sigmabar"; - - fSigmabar = 0.0; - fSigmabarInner = 0.0; - fSigmabarOuter = 0.0; - fRatioA = 0.0; - - fCalcPixNum=kTRUE; -} - -MSigmabar::~MSigmabar() -{ - // do nothing special. -} - - -// -------------------------------------------------------------------------- -// -// Actual calculation of sigmabar. This is done for each of the six sectors -// separately due to their possibly different HV behavior. Also inner and -// outer pixels are treated separately -// -// Preliminary! Works for CT1 test, for real MAGIC crosschecks still have -// to be done. Also implementation details will be updated, like -// determination of sector to which a respective pixel belongs -// -Bool_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped) -{ - const UInt_t npix = ped.GetSize(); - Float_t innerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - Float_t outerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - Int_t innerPixels[6] = {0,0,0,0,0,0}; - Int_t outerPixels[6] = {0,0,0,0,0,0}; - Int_t currentSector; - Float_t angle; - - for (UInt_t i=0; i