Index: trunk/MagicSoft/Mars/mjobs/MJSimulation.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSimulation.cc (revision 9424) +++ trunk/MagicSoft/Mars/mjobs/MJSimulation.cc (revision 9425) @@ -108,5 +108,6 @@ #include "MReflector.h" #include "MParEnv.h" -#include "MPulseShape.h" +#include "MSpline3.h" +#include "MParSpline.h" #include "MGeomCam.h" @@ -152,5 +153,5 @@ //cont.Add(const_cast(&fSequence)); - cont.Add(plist.FindObject("MPulseShape")); + cont.Add(plist.FindObject("PulseShape")); cont.Add(plist.FindObject("Reflector")); cont.Add(plist.FindObject("MGeomCam")); @@ -185,8 +186,15 @@ hist.SetLog(kTRUE, kTRUE, kFALSE); - hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100"); - hist.SetDrawOption("colz"); - hist.InitName("Impact;Impact;Impact"); - hist.InitTitle("Impact;West <--> East [m];South <--> North [m]"); + /* + hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100"); + hist.SetDrawOption("colz"); + hist.InitName("Impact;Impact;Impact"); + hist.InitTitle("Impact;West <--> East [m];South <--> North [m]"); + hist.SetAutoRange(); + */ + + hist.AddHist("MCorsikaEvtHeader.GetImpact/100"); + hist.InitName("Impact"); + hist.InitTitle("Impact;Impact [m]"); hist.SetAutoRange(); @@ -197,5 +205,5 @@ hist.AddHist("(MCorsikaEvtHeader.fAz+MCorsikaRunHeader.fMagneticFieldAz)*TMath::RadToDeg()", "MCorsikaEvtHeader.fZd*TMath::RadToDeg()"); hist.InitName("SkyOrigin;Az;Zd"); - hist.InitTitle("Sky Origin;Az [deg];Zd [deg]"); + hist.InitTitle("Sky Origin;Az [\\deg];Zd [\\deg]"); hist.SetDrawOption("colz"); hist.SetAutoRange(); @@ -279,5 +287,5 @@ plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm"); - MPulseShape shape; + MParSpline shape("PulseShape"); plist.AddToList(&shape); @@ -325,13 +333,23 @@ read.SetForceMode(fForceMode); - for (int i=0; iSetSelectedPad(0); static_cast(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete); - } @@ -749,5 +764,4 @@ c->SetBit(kCanDelete); c->Draw(); - } @@ -767,4 +781,42 @@ } } + + TCanvas &d = fDisplay->AddTab("Info2"); + d.Divide(2,2); + + d.cd(1); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + gROOT->SetSelectedPad(0); + splinepde.DrawClone()->SetBit(kCanDelete); + + d.cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + gROOT->SetSelectedPad(0); + splinemirror.DrawClone()->SetBit(kCanDelete); + + d.cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + gROOT->SetSelectedPad(0); + splinecones.DrawClone()->SetBit(kCanDelete); + + d.cd(4); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + gROOT->SetSelectedPad(0); + MParSpline *all = (MParSpline*)splinepde.DrawClone(); + //all->SetTitle("Combined acceptance"); + all->SetBit(kCanDelete); + all->Multiply(*splinemirror.GetSpline()); } Index: trunk/MagicSoft/Mars/msim/MSimAbsorption.cc =================================================================== --- trunk/MagicSoft/Mars/msim/MSimAbsorption.cc (revision 9424) +++ trunk/MagicSoft/Mars/msim/MSimAbsorption.cc (revision 9425) @@ -30,4 +30,5 @@ // // Input Containers: +// fParName [MParSpline] // MPhotonEvent // MCorsikaEvtHeader @@ -43,5 +44,4 @@ #include -#include #include "MLog.h" @@ -50,8 +50,5 @@ #include "MParList.h" -#include "MBinning.h" -#include "MArrayD.h" - -#include "MSpline3.h" +#include "MParSpline.h" #include "MCorsikaEvtHeader.h" @@ -69,161 +66,8 @@ // MSimAbsorption::MSimAbsorption(const char* name, const char *title) - : fEvt(0), fHeader(0), fSpline(0), fUseTheta(kFALSE) + : fEvt(0), fHeader(0), fSpline(0), fParName("MParSpline"), fUseTheta(kFALSE) { fName = name ? name : "MSimAbsorption"; fTitle = title ? title : "Task to calculate wavelength dependent absorption"; -} - -// -------------------------------------------------------------------------- -// -// Calls Clear() -// -MSimAbsorption::~MSimAbsorption() -{ - Clear(); -} - -// -------------------------------------------------------------------------- -// -// Delete fSpline and set to 0. -// -void MSimAbsorption::Clear(Option_t *) -{ - if (fSpline) - delete fSpline; - fSpline=0; -} - -// -------------------------------------------------------------------------- -// -// Read a TGraph from a file and return a newly allocated MSpline3. -// -MSpline3 *MSimAbsorption::ReadSpline(const char *fname) -{ - *fLog << inf << "Reading spline from " << fname << endl; - - TGraph g(fname); - if (g.GetN()==0) - { - *fLog << err << "ERROR - No data points from " << fname << "." << endl; - return 0; - } - - // option: b1/e1 b2/e2 (first second derivative?) - // option: valbeg/valend (first second derivative?) - - return new MSpline3(g); -} - -// -------------------------------------------------------------------------- -// -// Initializes a spline from min to max with n points with 1. -// -void MSimAbsorption::InitUnity(UInt_t n, Float_t min, Float_t max) -{ - // Delete the existing spline - Clear(); - - // We need at least two points (the edges) - if (n<2) - return; - - // Initialize an array with the x-values - const MBinning bins(n-1, min, max); - - // Initialize an array with all one - MArrayD y(n); - y.Reset(1); - - // Set the spline - fSpline = new MSpline3(bins.GetEdges(), y.GetArray(), n); -} - -// -------------------------------------------------------------------------- -// -// Takes the existing fSpline and multiplies it with the given spline. -// As x-points the values from fSpline are used. -// -void MSimAbsorption::Multiply(const TSpline3 &spline) -{ - if (!fSpline) - { - Clear(); - // WARNING WARNING WARNING: This is a very dangerous cast! - fSpline = (MSpline3*)spline.Clone(); - return; - } - - // Initialize a TGraph with the number of knots from a TSpline - TGraph g(fSpline->GetNp()); - - // Loop over all knots - for (int i=0; iGetNp(); i++) - { - // Get th i-th knot - Double_t x, y; - fSpline->GetKnot(i, x, y); - - // Multiply y by the value from the given spline - y *= spline.Eval(x); - - // push the point "back" - g.SetPoint(i, x, y); - } - - // Clear the spline and initialize a new one from the updated TGraph - Clear(); - fSpline = new MSpline3(g); -} - -// -------------------------------------------------------------------------- -// -// Initializes a TSpline3 with n knots x and y. Call Multiply for it. -// -void MSimAbsorption::Multiply(UInt_t n, const Double_t *x, const Double_t *y) -{ - const TSpline3 spline("Spline", const_cast(x), const_cast(y), n); - Multiply(spline); -} - -// -------------------------------------------------------------------------- -// -// Read a Spline from a file using ReadSpline. -// Call Multiply for it. -// -void MSimAbsorption::Multiply(const char *fname) -{ - const MSpline3 *spline = ReadSpline(fname); - if (!spline) - return; - - fFileName = ""; - - Multiply(*spline); - - delete spline; -} - -// -------------------------------------------------------------------------- -// -// Read a spline from a file and set fSpline accfordingly. -// -Bool_t MSimAbsorption::ReadFile(const char *fname) -{ - if (fname) - fFileName = fname; - - MSpline3 *spline = ReadSpline(fFileName); - if (!spline) - return kFALSE; - - - // option: b1/e1 b2/e2 (first second derivative?) - // option: valbeg/valend (first second derivative?) - - Clear(); - fSpline = spline; - - return kTRUE; } @@ -235,14 +79,21 @@ Int_t MSimAbsorption::PreProcess(MParList *pList) { - if (fFileName.IsNull()) - { - *fLog << inf << "No file given... skipping." << endl; - return kSKIP; - } - fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent"); if (!fEvt) { *fLog << err << "MPhotonEvent not found... aborting." << endl; + return kFALSE; + } + + fSpline = (MParSpline*)pList->FindObject(fParName, "MParSpline"); + if (!fSpline) + { + *fLog << err << fParName << " [MParSpline] not found... aborting." << endl; + return kFALSE; + } + + if (!fSpline->IsValid()) + { + *fLog << err << "Spline in " << fParName << " is inavlid... aborting." << endl; return kFALSE; } @@ -257,5 +108,5 @@ *fLog << inf << "Using " << (fUseTheta?"Theta":"Wavelength") << " for absorption." << endl; - return ReadFile(); + return kTRUE; } @@ -331,8 +182,8 @@ { Bool_t rc = kFALSE; - if (IsEnvDefined(env, prefix, "FileName", print)) + if (IsEnvDefined(env, prefix, "ParName", print)) { rc = kTRUE; - SetFileName(GetEnvValue(env, prefix, "FileName", fFileName)); + SetParName(GetEnvValue(env, prefix, "ParName", fParName)); } Index: trunk/MagicSoft/Mars/msim/MSimAbsorption.h =================================================================== --- trunk/MagicSoft/Mars/msim/MSimAbsorption.h (revision 9424) +++ trunk/MagicSoft/Mars/msim/MSimAbsorption.h (revision 9425) @@ -6,8 +6,6 @@ #endif -class TSpline3; - class MParList; -class MSpline3; +class MParSpline; class MPhotonEvent; class MCorsikaEvtHeader; @@ -19,7 +17,7 @@ MCorsikaEvtHeader *fHeader; //! Header storing event information - MSpline3 *fSpline; //! Spline to interpolate wavelength or incident angle + MParSpline *fSpline; //! Spline to interpolate wavelength or incident angle - TString fFileName; // FileName of the file containing the curve + TString fParName; // Container name of the spline containing the curve Bool_t fUseTheta; // Switches between using wavelength or incident angle @@ -32,25 +30,11 @@ Int_t Process(); - // MSimAbsorption - MSpline3 *ReadSpline(const char *fname); - public: MSimAbsorption(const char *name=NULL, const char *title=NULL); - ~MSimAbsorption(); // MSimAbsorption - Bool_t ReadFile(const char *fname=0); - void SetFileName(const char *fname) { fFileName=fname; } + void SetParName(const char *name) { fParName=name; } void SetUseTheta(Bool_t b=kTRUE) { fUseTheta = b; } - - void InitUnity(UInt_t n, Float_t min, Float_t max); - - void Multiply(const char *fname); - void Multiply(const TSpline3 &spline); - void Multiply(UInt_t n, const Double_t *x, const Double_t *y); - - // TObject - void Clear(Option_t *o=""); ClassDef(MSimAbsorption, 0) // Task to calculate wavelength or incident angle dependent absorption Index: trunk/MagicSoft/Mars/msimcamera/MSimCalibrationSignal.cc =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimCalibrationSignal.cc (revision 9424) +++ trunk/MagicSoft/Mars/msimcamera/MSimCalibrationSignal.cc (revision 9425) @@ -61,5 +61,5 @@ #include "MGeomCam.h" -#include "MPulseShape.h" +#include "MParSpline.h" #include "MTriggerPattern.h" @@ -136,8 +136,8 @@ } - fPulse = (MPulseShape*)pList->FindObject("MPulseShape"); + fPulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline"); if (!fPulse) { - *fLog << err << "MPulsShape not found... aborting." << endl; + *fLog << err << "PulsShape [MParSpline] not found... aborting." << endl; return kFALSE; } @@ -173,4 +173,5 @@ } + // FIXME: Is there a way to write them as LAST entry in the file? fRunHeader->SetReadyToSave(); @@ -195,4 +196,5 @@ for (UInt_t idx=0; idxGetNumPixels(); idx++) { + // FIXME: Scale number of photons with the pixel size! const Int_t num = TMath::Nint(gRandom->Gaus(fNumPhotons, fNumPhotons/10)); @@ -229,5 +231,5 @@ // Length (ns), Pulse position (Units ns) const Float_t pp = fPulsePos->GetVal(); - const Float_t pw = fPulse->GetPulseWidth(); + const Float_t pw = fPulse->GetWidth(); const Float_t first = cnt>0 ? fEvt->GetFirst()->GetTime() : 0; Index: trunk/MagicSoft/Mars/msimcamera/MSimCalibrationSignal.h =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimCalibrationSignal.h (revision 9424) +++ trunk/MagicSoft/Mars/msimcamera/MSimCalibrationSignal.h (revision 9425) @@ -8,5 +8,5 @@ class MGeomCam; class MParList; -class MPulseShape; +class MParSpline; class MPhotonEvent; class MPhotonStatistics; @@ -20,5 +20,5 @@ MParList *fParList; //! Store pointer to MParList for initializing ReInit MGeomCam *fGeom; //! Camera geometry to know the number of expected pixels - MPulseShape *fPulse; //! Pulse Shape to get pulse width from + MParSpline *fPulse; //! Pulse Shape to get pulse width from MParameterD *fPulsePos; //! Expected position at which the pulse should be MParameterD *fTrigger; //! Position in analog channels at which the triggersignal is raised Index: trunk/MagicSoft/Mars/msimcamera/MSimCamera.cc =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimCamera.cc (revision 9424) +++ trunk/MagicSoft/Mars/msimcamera/MSimCamera.cc (revision 9425) @@ -48,8 +48,7 @@ #include "MSpline3.h" -#include "MPulseShape.h" +#include "MParSpline.h" #include "MParList.h" -//#include "MParameters.h" #include "MPhotonEvent.h" @@ -134,8 +133,8 @@ return kFALSE; - MPulseShape *pulse = (MPulseShape*)pList->FindObject("MPulseShape"); + MParSpline *pulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline"); if (!pulse) { - *fLog << err << "MPulseShape not found... aborting." << endl; + *fLog << err << "PulseShape [MParSpline] not found... aborting." << endl; return kFALSE; } Index: trunk/MagicSoft/Mars/msimcamera/MSimCamera.h =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimCamera.h (revision 9424) +++ trunk/MagicSoft/Mars/msimcamera/MSimCamera.h (revision 9425) @@ -28,5 +28,5 @@ MMcEvt *fMcEvt; //! For information stored in MMcEvt - MSpline3 *fSpline; // Pusle Shape + const MSpline3 *fSpline; // Pusle Shape Bool_t fBaselineGain; // Should the gain be applied to baseline and electronic noise? Index: trunk/MagicSoft/Mars/msimcamera/MSimGeomCam.cc =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimGeomCam.cc (revision 9424) +++ trunk/MagicSoft/Mars/msimcamera/MSimGeomCam.cc (revision 9425) @@ -60,7 +60,6 @@ #include "MParameters.h" +#include "MParSpline.h" #include "MRawRunHeader.h" - -#include "MPulseShape.h" ClassImp(MSimGeomCam); @@ -106,5 +105,5 @@ } - fPulse = (MPulseShape*)pList->FindObject("MPulseShape"); + fPulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline"); /* if (!fPulse) @@ -191,6 +190,6 @@ // Length (ns), Pulse position (Units ns) - const Float_t pp = fPulsePos ? fPulsePos->GetVal() : 0; - const Float_t pw = fPulse ? fPulse->GetPulseWidth() : 0; + const Float_t pp = fPulsePos ? fPulsePos->GetVal() : 0; + const Float_t pw = fPulse ? fPulse->GetWidth() : 0; fStat->SetTimeMedDev(fEvt->GetTimeMedianDev()); Index: trunk/MagicSoft/Mars/msimcamera/MSimGeomCam.h =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimGeomCam.h (revision 9424) +++ trunk/MagicSoft/Mars/msimcamera/MSimGeomCam.h (revision 9425) @@ -13,5 +13,5 @@ class MParameterD; class MRawRunHeader; -class MPulseShape; +class MParSpline; class MSimGeomCam : public MTask @@ -23,5 +23,5 @@ MParameterD *fPulsePos; //! Intended pulse position in digitization window [ns] MRawRunHeader *fHeader; //! Length of digitization window - MPulseShape *fPulse; //! + MParSpline *fPulse; //! TString fNameGeomCam; Index: trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc (revision 9424) +++ trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc (revision 9425) @@ -43,4 +43,5 @@ // Output Containers: // MPhotonEvent +// AccidentalPhotonRate [MPedestalCam] // ////////////////////////////////////////////////////////////////////////////// @@ -62,5 +63,12 @@ #include "MPhotonData.h" +#include "MPedestalCam.h" +#include "MPedestalPix.h" + #include "MCorsikaRunHeader.h" + +#include "MSpline3.h" +#include "MParSpline.h" +#include "MReflector.h" ClassImp(MSimRandomPhotons); @@ -74,5 +82,5 @@ MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title) : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0), - fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam") + fRates(0), fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam") { fName = name ? name : "MSimRandomPhotons"; @@ -113,4 +121,8 @@ } + fRates = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "AccidentalPhotonRates"); + if (!fRates) + return kFALSE; + /* fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader"); @@ -132,4 +144,65 @@ } + MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector"); + if (!r) + { + *fLog << err << "Reflector [MReflector] not found... aborting." << endl; + return kFALSE; + } + + const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline"); + const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesAngularAcceptance", "MParSpline"); + const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity", "MParSpline"); + + const Double_t d2 = fGeom->GetCameraDist()*fGeom->GetCameraDist(); + const Double_t pde = s1 && s1->GetSpline() ? s1->GetSpline()->Integral() : 1; + const Double_t sr = s2 && s2->GetSpline() ? s2->GetSpline()->IntegralSolidAngle() : 1; + const Double_t mir = s3 && s3->GetSpline() ? s3->GetSpline()->Integral() : 1; + const Double_t Ar = r->GetA()/1e4; + + // Conversion factor to convert pixel area to steradians (because it + // is a rather small area we can assume it is flat) + const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad(); + + // Multiply all relevant efficiencies + MParSpline *s4 = (MParSpline*)s1->Clone(); + s4->Multiply(*s3->GetSpline()); + + const Double_t nm = s4 && s4->GetSpline() ? s4->GetSpline()->Integral() : 1; + + delete s4; + + // /100 to convert the pixel area from mm^2 to cm^2 + fScale = nm * TMath::Min(Ar, sr*d2) * conv*conv; + + *fLog << inf; + *fLog << "Effective cone acceptance: " << Form("%.2f", sr*d2) << "m^2" << endl; + *fLog << "Reflector area: " << Form("%.2f", Ar) << "m^2" << endl; + *fLog << "Resulting eff. collection area: " << Form("%.1f", TMath::Min(Ar, sr*d2)) << "m^2" << endl; + *fLog << "Eff. wavelength band (PDE): " << Form("%.1f", pde) << "nm" << endl; + *fLog << "Eff. wavelength band (Mirror): " << Form("%.1f", mir) << "nm" << endl; + *fLog << "Eff. wavelength band (PDE+MIR): " << Form("%.1f", nm) << "nm" << endl; + *fLog << "Pixel area of " << fNameGeomCam << "[0]: " << Form("%.2e", (*fGeom)[0].GetA()*conv*conv) << "sr" << endl; + //*fLog << "Effective angular acceptance: " << sr << "sr" << endl; + //*fLog << "Resulting NSB frequency: " << fFreqNSB*nm*Ar*1000 << "MHz/sr" << endl; + *fLog << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", fFreqNSB*(*fGeom)[0].GetA()*fScale*1000) << "MHz" << endl; + + // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader"); + // Set NumPheFromDNSB + + // # Number of photons from the diffuse NSB (nphe / ns 0.1*0.1 deg^2 239 m^2) and + // nsb_mean 0.20 + // Magic pixel: 0.00885361 deg + // dnsbpix = 0.2*50/15 + // ampl = MMcFadcHeader->GetAmplitud() + // sqrt(pedrms*pedrms + dnsbpix*ampl*ampl/ratio) + + return kTRUE; +} + +Bool_t MSimRandomPhotons::ReInit(MParList *pList) +{ + // Overwrite the default set by MGeomApply + fRates->Init(*fGeom); return kTRUE; } @@ -160,7 +233,11 @@ for (UInt_t idx=0; idx