Index: trunk/MagicSoft/Mars/manalysis/MPadding.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadding.cc (revision 1748) +++ trunk/MagicSoft/Mars/manalysis/MPadding.cc (revision 1748) @@ -0,0 +1,569 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! Wolfgang Wittek 01/2003 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MPadding // +// (developped from 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. The number of photons, // +// its error and the pedestal sigmas are altered. // +// // +// The padding has to be done before the image cleaning because the // +// image cleaning depends on the pedestal sigmas. // +// // +// There are several ways of defining the sigmabar value to which the // +// events are padded: // +// // +// 1) Set a fixed level (fFixedSigmabar) by calling 'SetTargetLevel'. // +// // +// 2) By calling 'SetDefiningHistogram', give a TH1D histogram // +// (fHSigmabarMax) which defines the Sigmabar as a function of Theta. // +// // +// 3) By calling 'SetSigmaThetaHist', give a TH2D histogram // +// (fHSigmaTheta) which contains the Sigmabar distribution for the // +// different bins in Theta. For a given event, the sigmabar value to // +// be used for the padding is thrown from this distribution. // +// // +// Workaround : // +// If none of these options is specified then PreProcess will try to read // +// in a propriety format ASCII database for the CT1 test. The name of // +// this file is set by 'SetDatabaseFile'. From the data in this file a // +// TH1D histogram (fHSigmabarMax) is generated. // +// // +// 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 "MPadding.h" + +#include + +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TProfile.h" +#include "TRandom.h" +#include "TCanvas.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(MPadding); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MPadding::MPadding(const char *name, const char *title) : fRunType(0), fGroup(0), fFixedSigmabar(0.0) +{ + fName = name ? name : "MPadding"; + fTitle = title ? title : "Task for the padding"; + + fFixedSigmabar = 0.0; + fHSigmabarMax = NULL; + fHSigmaTheta = NULL; + + Print(); +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MPadding::~MPadding() +{ + //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 MPadding::SetDefiningHistogram(TH1D *histo) +{ + fHSigmabarMax = histo; + + fFixedSigmabar = 0.0; + fHSigmaTheta = NULL; + + *fLog << "MPadding::SetDefiningHistogram" << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// You can provide a TH2D* histogram containing +// - the Sigmabar distribution in bins of theta. +// +Bool_t MPadding::SetSigmaThetaHist(TH2D *histo) +{ + fHSigmaTheta = histo; + + fFixedSigmabar = 0.0; + fHSigmabarMax = NULL; + + *fLog << "MPadding::SetSigmaThetaHist" << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +void MPadding::SetTargetLevel(Double_t sigmabar) +{ + fFixedSigmabar = sigmabar; + + fHSigmabarMax = NULL; + fHSigmaTheta = NULL; + + *fLog << "MPadding::SetTargetLevel; use fixed sigmabar : fFixedSigmabar = " + << fFixedSigmabar << endl; +} + +// -------------------------------------------------------------------------- +// +// check if MEvtHeader exists in the Parameter list already. +// if not create one and add them to the list +// +Bool_t MPadding::PreProcess(MParList *pList) +{ + fRnd = new TRandom3(0); + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MPadding::PreProcess : MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << dbginf << "MPadding::PreProcess : MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MPadding::PreProcess : MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MPadding::PreProcess : MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar"); + if (!fSigmabar) + { + *fLog << dbginf << "MPadding::PreProcess : MSigmabar not found... aborting." << endl; + return kFALSE; + } + + // Get Theta Binning + MBinning* binstheta = (MBinning*)pList->FindObject("BinningTheta"); + if (!binstheta) + { + *fLog << err << dbginf << "MPadding::PreProcess : BinningTheta not found... aborting." << endl; + return kFALSE; + } + + // Get Sigma Binning + MBinning* binssigma = (MBinning*)pList->FindObject("BinningSigmabar"); + if (!binssigma) + { + *fLog << err << dbginf << "MPadding::PreProcess : BinningSigmabar not found... aborting." << endl; + return kFALSE; + } + + //-------------------------------------------------------------------- + // plot pedestal sigmas for testing purposes + fHSigmaPedestal = new TH2D("SigPed","padded vs orig. sigma", + 100, 0.0, 5.0, 100, 0.0, 5.0); + fHSigmaPedestal->SetXTitle("orig. Pedestal sigma"); + fHSigmaPedestal->SetYTitle("padded Pedestal sigma"); + + // plot no.of photons (before vs. after padding) for testing purposes + fHPhotons = new TH2D("Photons","Photons after vs.before padding", + 100, -10.0, 90.0, 100, -10, 90); + fHPhotons->SetXTitle("no.of photons before padding"); + fHPhotons->SetYTitle("no.of photons after padding"); + + // plot of added NSB + fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0); + fHNSB->SetXTitle("no.of added NSB photons"); + fHNSB->SetYTitle("no.of pixels"); + + fHSigmaOld = new TH2D(); + fHSigmaOld->SetNameTitle("fHSigmaOld","Sigma before padding"); + MH::SetBinning(fHSigmaOld, binstheta, binssigma); + fHSigmaOld->SetXTitle("Theta"); + fHSigmaOld->SetYTitle("Sigma"); + + fHSigmaNew = new TH2D(); + fHSigmaNew->SetNameTitle("fHSigmaNew","Sigma after padding"); + MH::SetBinning(fHSigmaNew, binstheta, binssigma); + fHSigmaNew->SetXTitle("Theta"); + fHSigmaNew->SetYTitle("Sigma"); + + + //************************************************************************ + // Create fSigmabarMax histogram + // (only if no fixed Sigmabar target value and no histogram have been + // provided) + // + if ( (fFixedSigmabar==0.0) && (fHSigmabarMax==NULL) + && (fHSigmaTheta==NULL) ) + { + *fLog << "MPadding::PreProcess() : fFixedSigmabar, fHSigmabarMax = " + << fFixedSigmabar << ", " << fHSigmabarMax << endl; + *fLog << " create fSigmabarMax histogram" << endl; + + fHSigmabarMax = new TH1D(); + fHSigmabarMax->SetNameTitle("fHSigmabarMax","Sigmabarmax for this analysis"); + TAxis &x = *fHSigmabarMax->GetXaxis(); + 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()); + x.SetTitle("Theta"); + TAxis &y = *fHSigmabarMax->GetYaxis(); + y.SetTitle("SigmabarMax"); + + // ------------------------------------------------- + // 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 MPadding::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 + + } //fFixedSigmabar + //************************************************************************ + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculate Sigmabar for current event +// Then apply padding +// +// 1) have current event +// 2) get sigmabar(theta) +// 3) pad event +// +Bool_t MPadding::Process() +{ + //------------------------------------------- + // Calculate sigmabar of event + // + Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt); + //fSigmabar->Print(""); + + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + mySig = 0.0; + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + + + // Get sigmabar which we have to pad to + // + Double_t otherSig; + Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + + // *fLog << "Theta = " << Theta << endl; + + //------------------------------------------- + // get target sigma for the current Theta from the histogram fHSigmabarMax + // + + if (fHSigmabarMax != NULL) + { + Int_t binNumber = fHSigmabarMax->GetXaxis()->FindBin(Theta); + if (binNumber < 1 || binNumber > fHSigmabarMax->GetNbinsX()) + { + *fLog << err << dbginf + << "Theta of current event is beyond the limits, Theta = " + << kRad2Deg*fMcEvt->GetTelescopeTheta() + << " ...Skipping this event" <GetBinContent(binNumber); + + //*fLog << "Theta, binNumber, otherSig = " + // << kRad2Deg*fMcEvt->GetTelescopeTheta() << ", " + // << binNumber << ", " << otherSig << endl; + } + } + + //------------------------------------------- + // for the current Theta, + // generate a sigma according to the histogram fHSigmaTheta + // + else if (fHSigmaTheta != NULL) + { + Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(Theta); + + if ( binNumber < 1 || binNumber > fHSigmaTheta->GetNbinsX() ) + { + // *fLog << "MPadding::Process(); binNumber out of range, binNumber = " + // << binNumber << ", Skip event " << endl; + return kCONTINUE; + } + else + { + TH1D* fHSigma = + fHSigmaTheta->ProjectionY("", binNumber, binNumber, ""); + if ( fHSigma->GetEntries() == 0.0 ) + { + //*fLog << "MPadding::Process(); projection for Theta bin " << binNumber + // << " has no entries, Skip event" << endl; + return kCONTINUE; + } + else + { + otherSig = fHSigma->GetRandom(); + + //*fLog << "Theta, bin number = " << Theta << ", " << binNumber + // << ", otherSig = " << otherSig << endl; + } + delete fHSigma; + } + } + + //------------------------------------------- + // use a fixed target sigma + // + else if (fFixedSigmabar != 0.0) + { + otherSig = fFixedSigmabar; + } + + //------------------------------------------- + else + { + *fLog << "MPadding::Process(); sigmabar for padding not defined" << endl; + return kFALSE; + } + + //------------------------------------------- + // + + //*fLog << "MPadding::Process(); mySig, otherSig = " << mySig << ", " + // << otherSig << endl; + + // Skip event if target sigma is zero + if (otherSig == 0.0) + { + return kCONTINUE; + } + + // Determine quadratic difference other-mine + Double_t quadraticDiff = otherSig*otherSig - mySig*mySig; + + if (quadraticDiff < 0) { + *fLog << err << dbginf << "Event has higher Sigmabar=" << mySig + <<" than Sigmabarmax=" << otherSig << "; Theta =" + << kRad2Deg*fMcEvt->GetTelescopeTheta() + << " ...Skipping this event" < 0, do the padding; + // + Padding(quadraticDiff); + + // Calculate Sigmabar again and crosscheck + mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt); + + //fSigmabar->Print(""); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Do the padding (mySig ==> otherSig) +// +Bool_t MPadding::Padding(Double_t quadraticDiff) +{ + const UInt_t npix = fEvt->GetNumPixels(); + + // pad only pixels - which are used (before image cleaning) + // - and for which the no.of photons is != 0.0 + // + for (UInt_t i=0; ioperator[](i); + if ( !pix.IsPixelUsed() ) + continue; + + if ( pix.GetNumPhotons() == 0.0) + { + *fLog << "MPadding::Process(); no.of photons is 0 for used pixel" + << endl; + continue; + } + + Int_t j = pix.GetPixId(); + Double_t Area = fCam->GetPixRatio(j); + + // add additional NSB to the number of photons + Double_t oldphotons = pix.GetNumPhotons(); + Double_t NSB = sqrt(quadraticDiff*Area) * fRnd->Gaus(0.0, 1.0); + Double_t newphotons = oldphotons + NSB; + pix.SetNumPhotons( newphotons ); + + fHNSB->Fill( NSB/sqrt(Area) ); + fHPhotons->Fill( newphotons/sqrt(Area), oldphotons/sqrt(Area) ); + + + // error: add sigma of padded noise quadratically + Double_t olderror = pix.GetErrorPhot(); + Double_t newerror = sqrt( olderror*olderror + quadraticDiff*Area ); + pix.SetErrorPhot( newerror ); + + + MPedestalPix &ppix = fPed->operator[](j); + + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + ppix.SetMeanRms(0.0); + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + + Double_t oldsigma = ppix.GetMeanRms(); + Double_t newsigma = sqrt( oldsigma*oldsigma + quadraticDiff*Area ); + ppix.SetMeanRms( newsigma ); + + fHSigmaPedestal->Fill( oldsigma, newsigma ); + fHSigmaOld->Fill( kRad2Deg*fMcEvt->GetTelescopeTheta(), oldsigma ); + fHSigmaNew->Fill( kRad2Deg*fMcEvt->GetTelescopeTheta(), newsigma ); + } //for + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MPadding::PostProcess() +{ + TCanvas &c = *(MH::MakeDefCanvas("Padding", "", 600, 900)); + c.Divide(2,3); + gROOT->SetSelectedPad(NULL); + + + if (fHSigmabarMax != NULL) + { + c.cd(1); + fHSigmabarMax->DrawClone(); + fHSigmabarMax->SetBit(kCanDelete); + } + else if (fHSigmaTheta != NULL) + { + c.cd(1); + fHSigmaTheta->DrawClone(); + fHSigmaTheta->SetBit(kCanDelete); + } + + c.cd(3); + fHSigmaPedestal->DrawClone(); + fHSigmaPedestal->SetBit(kCanDelete); + + c.cd(5); + fHPhotons->DrawClone(); + fHPhotons->SetBit(kCanDelete); + + c.cd(2); + fHNSB->DrawClone(); + fHNSB->SetBit(kCanDelete); + + c.cd(4); + fHSigmaOld->DrawClone(); + fHSigmaOld->SetBit(kCanDelete); + + c.cd(6); + fHSigmaNew->DrawClone(); + fHSigmaNew->SetBit(kCanDelete); + + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadding.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadding.h (revision 1748) +++ trunk/MagicSoft/Mars/manalysis/MPadding.h (revision 1748) @@ -0,0 +1,87 @@ +#ifndef MARS_MPadding +#define MARS_MPadding + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +#include "TRandom3.h" +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TProfile.h" + + +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MMcEvt; +class MPedestalCam; +class MSigmabar; +class MParList; + +class MPadding : public MTask +{ +private: + const MGeomCam *fCam; + MCerPhotEvt *fEvt; + MSigmabar *fSigmabar; + MMcEvt *fMcEvt; + MPedestalCam *fPed; + + TRandom3 *fRnd; + + Int_t fRunType; + Int_t fGroup; + + TH1D *fHSigmabarMax; // histogram (sigmabarmax vs. Theta) + char *fDatabaseFilename; // data file used for generating + // fHSigmabarMax histogram + + TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) + + Double_t fFixedSigmabar; // fixed sigmabar value + + TH2D *fHSigmaPedestal; // for testing : plot of padded vs + // orig. pedestal sigmas + TH2D *fHPhotons; // for testing : no.of photons after + // versus before padding + TH2D *fHSigmaOld; // histogram (sigma vs. Theta) + // before padding + + TH2D *fHSigmaNew ; // histogram (sigma vs. Theta) + // after padding + TH1D *fHNSB; // histogram of added NSB + + +public: + MPadding(const char *name=NULL, const char *title=NULL); + ~MPadding(); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + Bool_t Padding(Double_t quadDiff); + + void SetRunType(Int_t runtype) { fRunType = runtype; } + void SetGroup(Int_t group) { fGroup = group; } + + Bool_t SetDefiningHistogram(TH1D *hist); + void SetDatabaseFile(char *filename) { fDatabaseFilename = filename; } + + Bool_t SetSigmaThetaHist(TH2D *histo); + + void SetTargetLevel(Double_t sigmabar); + + ClassDef(MPadding, 1) // task for the padding +}; + +#endif + + + Index: trunk/MagicSoft/Mars/manalysis/MSigmabar.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSigmabar.cc (revision 1747) +++ trunk/MagicSoft/Mars/manalysis/MSigmabar.cc (revision 1748) @@ -16,5 +16,6 @@ ! ! -! Author(s): Robert Wagner 10/2002 +! Author(s): Robert Wagner 10/2002 +! Wolfgang Wittek 01/2003 ! ! Copyright: MAGIC Software Development, 2002 @@ -27,6 +28,8 @@ // MSigmabar // // // -// This is the storage container to hold information about the mean sigma // -// (aka Sigmabar) of all pedestals // +// This is the storage container to hold information about // +// "average" pedestal sigmas // +// // +// In calculating averages all sigmas are normalized to the area of pixel 0// // // ///////////////////////////////////////////////////////////////////////////// @@ -37,12 +40,16 @@ #include "MLog.h" #include "MLogManip.h" - +#include "MParList.h" #include "MGeomCam.h" #include "MPedestalCam.h" +#include "MPedestalPix.h" #include "MGeomPix.h" -#include "MPedestalPix.h" +#include "MCerPhotEvt.h" +#include "MCerPhotPix.h" ClassImp(MSigmabar); +// -------------------------------------------------------------------------- +// MSigmabar::MSigmabar(const char *name, const char *title) { @@ -50,17 +57,14 @@ 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. } - // -------------------------------------------------------------------------- @@ -68,5 +72,5 @@ // 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 +// outer pixels are treated separately. // // Preliminary! Works for CT1 test, for real MAGIC crosschecks still have @@ -74,77 +78,150 @@ // 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 MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped, + const MCerPhotEvt &evt) +{ + fSigmabar = 0.0; + fSigmabarInner = 0.0; + fSigmabarOuter = 0.0; + + 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 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 0) + fSigmabar=sqrt( (fSigmabarInner + fSigmabarOuter) + / (fInnerPixels + fOuterPixels) ); + + if (fInnerPixels > 0) fSigmabarInner /= fInnerPixels; + if (fOuterPixels > 0) fSigmabarOuter /= fOuterPixels; + + // this is the sqrt of the average sigma**2 + // for the inner and outer pixels respectively + // + fSigmabarInner = sqrt( fSigmabarInner ); + fSigmabarOuter = sqrt( fSigmabarOuter ); + for (UInt_t i=0; i<6; i++) { - fSigmabarInnerSector[i]=innerSquaredSum[i]/innerPixels[i]; - fSigmabarOuterSector[i]=outerSquaredSum[i]/( outerSquaredSum[i]==0 ? 1: outerPixels[i] ); - - fSigmabarSector[i]=sqrt(fSigmabarInnerSector[i] + fSigmabarOuterSector[i]*fRatioA); - fSigmabarInnerSector[i]=sqrt(fSigmabarInnerSector[i]); - fSigmabarOuterSector[i]=sqrt(fSigmabarOuterSector[i]); + fSigmabarSector[i] = 0.0; + fSigmabarInnerSector[i] = 0.0; + fSigmabarOuterSector[i] = 0.0; + + if ( (innerPixels[i]+outerPixels[i]) > 0) + fSigmabarSector[i] = sqrt( (innerSquaredSum[i] + outerSquaredSum[i]) + / (innerPixels[i] + outerPixels[i] ) ); + if ( innerPixels[i] > 0) + fSigmabarInnerSector[i] = innerSquaredSum[i] / innerPixels[i]; + if ( outerPixels[i] > 0) + fSigmabarOuterSector[i] = outerSquaredSum[i] / outerPixels[i]; + + + fSigmabarInnerSector[i] = sqrt( fSigmabarInnerSector[i] ); + fSigmabarOuterSector[i] = sqrt( fSigmabarOuterSector[i] ); } - // Did all our calculations work? fOuterPixels==0 could happen, however. - return (fInnerPixels!=0); -} - - + return (fSigmabar); +} + +// -------------------------------------------------------------------------- +// void MSigmabar::Print(Option_t *) const { - *fLog << all; - *fLog << "Sigmabar Overall " << fSigmabar; + *fLog << endl; + *fLog << "Total number of inner pixels is " << fInnerPixels << endl; + *fLog << "Total number of outer pixels is " << fOuterPixels << endl; + *fLog << endl; + + *fLog << "Sigmabar Overall : " << fSigmabar << " "; *fLog << " Sectors: "; - for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << " "; - *fLog << endl; - *fLog << "Sigmabar Inner " << fSigmabarInner; + for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << ", "; + *fLog << endl; + + + *fLog << "Sigmabar Inner : " << fSigmabarInner << " "; *fLog << " Sectors: "; - for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << " "; - *fLog << endl; - *fLog << "Sigmabar Outer " << fSigmabarOuter; + for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << ", "; + *fLog << endl; + + + *fLog << "Sigmabar Outer : " << fSigmabarOuter << " "; *fLog << " Sectors: "; - for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << " "; - *fLog << endl; - *fLog << "Total number of inner pixels found is " << fInnerPixels << endl; - *fLog << "Total number of outer pixels found is " << fOuterPixels << endl; - *fLog << "Ratio of areas outer/inner pixels found is " << fRatioA << endl; - *fLog << endl; -} + for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << ", "; + *fLog << endl; + +} +// -------------------------------------------------------------------------- + + + + Index: trunk/MagicSoft/Mars/manalysis/MSigmabar.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSigmabar.h (revision 1747) +++ trunk/MagicSoft/Mars/manalysis/MSigmabar.h (revision 1748) @@ -4,4 +4,8 @@ #ifndef MARS_MParContainer #include "MParContainer.h" +#endif + +#ifndef MARS_MParList +#include "MParList.h" #endif @@ -14,18 +18,22 @@ #endif +class MCerPhotEvt; + class MSigmabar : public MParContainer { private: - Float_t fSigmabar; // Sigmabar (mean standard deviation) of pedestal + Float_t fSigmabar; // Sigmabar ("average" RMS) of pedestal + Float_t fSigmabarInner; // --only for inner pixels + Float_t fSigmabarOuter; // --only for outer pixels + Float_t fSigmabarSector[6]; // --for the 6 sectors of the camera Float_t fSigmabarInnerSector[6]; Float_t fSigmabarOuterSector[6]; - Float_t fSigmabarInner; // --only for inner pixels - Float_t fSigmabarOuter; // --only for outer pixels - UInt_t fInnerPixels; // Overall number of inner pixels - UInt_t fOuterPixels; // Overall number of outer pixels - Float_t fRatioA; // Ratio of areas (outer/inner pixels) - Bool_t fCalcPixNum; + UInt_t fInnerPixels; // Overall number of inner pixels + UInt_t fOuterPixels; // Overall number of outer pixels + + Bool_t fCalcPixNum; + public: @@ -46,5 +54,5 @@ // void SetSigmabarOuter(Float_t f) { fSigmabarOuter = f; } - Bool_t Calc(const MGeomCam &geom, const MPedestalCam &ped); + Float_t Calc(const MGeomCam &geom, const MPedestalCam &ped, const MCerPhotEvt &evt); ClassDef(MSigmabar, 1) // Storage Container for Sigmabar @@ -53,2 +61,6 @@ #endif + + + + Index: trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc (revision 1747) +++ trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc (revision 1748) @@ -110,5 +110,12 @@ if (!fMcEvt) { - *fLog << err << dbginf << "MHSigmabarTheta : MMcEvt not found... aborting." << endl; + *fLog << err << dbginf << "MHSigmabaCalc : MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << err << dbginf << "MHSigmabarCalc : MCerPhotEvt not found... aborting." << endl; return kFALSE; } @@ -126,13 +133,13 @@ Bool_t MSigmabarCalc::Process() { - Bool_t rc = fSig->Calc(*fCam, *fPed); - fSigmabarMax = TMath::Max((Double_t)fSig->GetSigmabar(), fSigmabarMax); - fSigmabarMin = TMath::Min((Double_t)fSig->GetSigmabar(), fSigmabarMin); + Float_t rc = fSig->Calc(*fCam, *fPed, *fEvt); + fSigmabarMax = TMath::Max((Double_t)rc, fSigmabarMax); + fSigmabarMin = TMath::Min((Double_t)rc, fSigmabarMin); - if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 90) + if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 120) fThetaMax = TMath::Max(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMax); fThetaMin = TMath::Min(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMin); - return rc; + return kTRUE; } Index: trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.h (revision 1747) +++ trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.h (revision 1748) @@ -43,4 +43,5 @@ MSigmabarParam *fSigParam; MMcEvt *fMcEvt; + MCerPhotEvt *fEvt; void Reset(); @@ -58,2 +59,4 @@ #endif + + Index: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 1748) +++ trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 1748) @@ -0,0 +1,363 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek 1/2003 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHSigmaTheta (extension of Robert's MHSigmabarTheta) // +// // +// calculates - the 2D-histogram sigmabar vs. Theta, and // +// - the 3D-histogram sigma, pixel no., Theta // +// - the 3D-histogram (sigma^2-sigmabar^2), pixel no., Theta // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHSigmaTheta.h" + +#include + +#include "MTime.h" +#include "MMcEvt.hxx" + +#include "MBinning.h" +#include "MParList.h" +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MCerPhotEvt.h" +#include "MCerPhotPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MSigmabar.h" + +ClassImp(MHSigmaTheta); + + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHSigmaTheta::MHSigmaTheta(const char *name, const char *title) + : fSigmaTheta(), fSigmaPixTheta() +{ + fName = name ? name : "MHSigmaTheta"; + fTitle = title ? title : "2D histogram sigmabar vs. Theta"; + + fSigmaTheta.SetDirectory(NULL); + fSigmaTheta.SetName("2D-ThetaSigmabar"); + fSigmaTheta.SetTitle("2D : Sigmabar, \\Theta"); + fSigmaTheta.SetXTitle("\\Theta [\\circ]"); + fSigmaTheta.SetYTitle("Sigmabar"); + + fSigmaPixTheta.SetDirectory(NULL); + fSigmaPixTheta.SetName("3D-ThetaPixSigma"); + fSigmaPixTheta.SetTitle("3D : \\Theta, pixel no., Sigma"); + fSigmaPixTheta.SetXTitle("\\Theta [\\circ]"); + fSigmaPixTheta.SetYTitle("pixel number"); + fSigmaPixTheta.SetZTitle("Sigma"); + + fDiffPixTheta.SetDirectory(NULL); + fDiffPixTheta.SetName("3D-ThetaPixDiff"); + fDiffPixTheta.SetTitle("3D : \\Theta, pixel, Sigma^2-Sigmabar^2"); + fDiffPixTheta.SetXTitle("\\Theta [\\circ]"); + fDiffPixTheta.SetYTitle("pixel number"); + fDiffPixTheta.SetZTitle("Sigma^2-sigmabar^2"); +} + +// -------------------------------------------------------------------------- +// +// Set the binnings and prepare the filling of the histogram +// +Bool_t MHSigmaTheta::SetupFill(const MParList *plist) +{ + + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MHSigmaTheta::SetupFill : MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)plist->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << dbginf << "MHSigmaTheta::SetupFill : MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)plist->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MHSigmaTheta::SetupFill : MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MHSigmaTheta::SetupFill : MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)plist->FindObject("MSigmabar"); + if (!fSigmabar) + { + *fLog << dbginf << "MHSigmaTheta::SetupFill : MSigmabar not found... aborting." << endl; + return kFALSE; + } + + // Get Theta Binning + MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); + if (!binstheta) + { + *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningTheta not found... aborting." << endl; + return kFALSE; + } + + // Get Sigmabar binning + MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar"); + if (!binssigma) + { + *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningSigmabar not found... aborting." << endl; + return kFALSE; + } + + // Get binning for (sigma^2-sigmabar^2) + MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2"); + if (!binsdiff) + { + *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningDiffsigma2 not found... aborting." << endl; + return kFALSE; + } + + + // Set binnings in histograms + SetBinning(&fSigmaTheta, binstheta, binssigma); + + // Get binning for pixel number + UInt_t npix = fPed->GetSize(); + MBinning binspix("BinningPixel"); + MBinning* binspixel = &binspix; + binspixel->SetEdges(npix, -0.5, ((float)npix)-0.5 ); + + SetBinning(&fSigmaPixTheta, binstheta, binspixel, binssigma); + SetBinning(&fDiffPixTheta, binstheta, binspixel, binsdiff); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms +// +Bool_t MHSigmaTheta::Fill(const MParContainer *par) +{ + //*fLog << "entry Fill" << endl; + + Double_t Theta = fMcEvt->GetTelescopeTheta()*kRad2Deg; + Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt); + fSigmaTheta.Fill(Theta, mySig); + + //*fLog << "Theta, mySig = " << Theta << ", " << mySig << endl; + + const UInt_t npix = fEvt->GetNumPixels(); + for (Int_t i=0; ioperator[](i); + if (!cerpix.IsPixelUsed()) + continue; + + if (cerpix.GetNumPhotons() == 0) + continue; + + Int_t j = cerpix.GetPixId(); + const MPedestalPix pix = fPed->operator[](j); + + Double_t Sigma = pix.GetMeanRms(); + Double_t Area = fCam->GetPixRatio(j); + + fSigmaPixTheta.Fill(Theta, (Double_t)j, Sigma); + + Double_t Diff = Sigma*Sigma/Area - mySig*mySig; + fDiffPixTheta.Fill (Theta, (Double_t)j, Diff); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Plot the results +// +Bool_t MHSigmaTheta::Finalize() +{ + DrawClone(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Draw a copy of the histogram +// +TObject *MHSigmaTheta::DrawClone(Option_t *opt) +{ + TCanvas &c = *MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta", + 900, 900); + c.Divide(3, 3); + + gROOT->SetSelectedPad(NULL); + + //-------------------------------------------------------------------- + // draw the 2D histogram Sigmabar versus Theta + TH1D *h; + + c.cd(1); + h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E"); + h->SetTitle("Distribution of \\Theta"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetYTitle("No.of events"); + + h->Draw(opt); + h->SetBit(kCanDelete);; + gPad->SetLogy(); + + c.cd(4); + h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E"); + h->SetTitle("Distribution of Sigmabar"); + h->SetXTitle("Sigmabar"); + h->SetYTitle("No.of events"); + + h->Draw(opt); + h->SetBit(kCanDelete);; + + c.cd(7); + ((TH2*)&fSigmaTheta)->DrawCopy(opt); + + //-------------------------------------------------------------------- + // draw the 3D histogram : Theta, pixel, Sigma^2-Sigmabar^2 + + TH2D *l; + + c.cd(2); + l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zx"); + l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); + l->SetXTitle("\\Theta [\\circ]"); + l->SetYTitle("Sigma^2-Sigmabar^2"); + + l->Draw("box"); + l->SetBit(kCanDelete);; + + c.cd(5); + l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zy"); + l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); + l->SetXTitle("pixel"); + l->SetYTitle("Sigma^2-Sigmabar^2"); + + l->Draw("box"); + l->SetBit(kCanDelete);; + + c.cd(8); + ((TH2*)&fDiffPixTheta)->DrawCopy(opt); + + + //-------------------------------------------------------------------- + // draw the 3D histogram : Theta, pixel, Sigma + + TH2D *k; + + c.cd(3); + k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zx"); + k->SetTitle("Sigma vs. \\Theta (all pixels)"); + k->SetXTitle("\\Theta [\\circ]"); + k->SetYTitle("Sigma"); + + k->Draw("box"); + k->SetBit(kCanDelete);; + + c.cd(6); + k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zy"); + k->SetTitle("Sigma vs. pixel number (all \\Theta)"); + k->SetXTitle("pixel"); + k->SetYTitle("Sigma"); + + k->Draw("box"); + k->SetBit(kCanDelete);; + + c.cd(9); + ((TH2*)&fSigmaPixTheta)->DrawCopy(opt); + + //-------------------------------------------------------------------- + c.Modified(); + c.Update(); + + return &c; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHSigmaTheta::Draw(Option_t *opt) +{ + if (!gPad) + MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta", + 600, 600); + + TH1D *h; + + gPad->Divide(2,2); + + gPad->cd(1); + h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E"); + h->SetTitle("Distribution of \\Theta"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetYTitle("No.of events"); + + h->Draw(opt); + h->SetBit(kCanDelete);; + gPad->SetLogy(); + + gPad->cd(2); + h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E"); + h->SetTitle("Distribution of Sigmabar"); + h->SetXTitle("Sigmabar"); + h->SetYTitle("No.of events"); + + h->Draw(opt); + h->SetBit(kCanDelete);; + + gPad->cd(3); + fSigmaTheta.DrawCopy(opt); + + + + gPad->Modified(); + gPad->Update(); +} +// -------------------------------------------------------------------------- + + + + Index: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 1748) +++ trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 1748) @@ -0,0 +1,73 @@ +#ifndef MARS_MHSigmaTheta +#define MARS_MHSigmaTheta + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +#ifndef ROOT_TH3 +#include "TH3.h" +#endif + +#ifndef ROOT_TProfile2D +#include "TProfile2D.h" +#endif + +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MMcEvt; +class MPedestalCam; +class MSigmabar; +class MParList; + + +class MHSigmaTheta : public MH +{ +private: + const MGeomCam *fCam; + MPedestalCam *fPed; + MCerPhotEvt *fEvt; + MSigmabar *fSigmabar; + MMcEvt *fMcEvt; + + TH2D fSigmaTheta; // 2D-distribution sigmabar versus Theta; + // sigmabar is the average pedestasl sigma + // in an event + TH3D fSigmaPixTheta;// 3D-distr.:Theta, pixel, pedestal sigma + TH3D fDiffPixTheta; // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2 + + +public: + MHSigmaTheta(const char *name=NULL, const char *title=NULL); + + virtual Bool_t SetupFill(const MParList *plist); + virtual Bool_t Fill(const MParContainer *par); + virtual Bool_t Finalize(); + + const TH2D *GetSigmaTheta() { return &fSigmaTheta; } + const TH2D *GetSigmaTheta() const { return &fSigmaTheta; } + TH2D *GetSigmaThetaByName(const TString name) { return &fSigmaTheta; } + + const TH3D *GetSigmaPixTheta() { return &fSigmaPixTheta; } + const TH3D *GetSigmaPixTheta() const { return &fSigmaPixTheta; } + TH3D *GetSigmaPixThetaByName(const TString name) { return &fSigmaPixTheta; } + + const TH3D *GetDiffPixTheta() { return &fDiffPixTheta; } + const TH3D *GetDiffPixTheta() const { return &fDiffPixTheta; } + TH3D *GetDiffPixThetaByName(const TString name) { return &fDiffPixTheta; } + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option=""); + + ClassDef(MHSigmaTheta, 0) //2D-histogram sigmabar vs. Theta +}; + +#endif + + +