Changeset 9425 for trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc

Timestamp:

04/16/09 12:04:29 (15 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc (modified) (7 diffs)

trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc

@@ -43 +43 @@
 //  Output Containers:
 //   MPhotonEvent
+//   AccidentalPhotonRate [MPedestalCam]
 //
 //////////////////////////////////////////////////////////////////////////////
@@ -62 +63 @@
 #include "MPhotonData.h"
 
+#include "MPedestalCam.h"
+#include "MPedestalPix.h"
+
 #include "MCorsikaRunHeader.h"
+
+#include "MSpline3.h"
+#include "MParSpline.h"
+#include "MReflector.h"
 
 ClassImp(MSimRandomPhotons);
@@ -74 +82 @@
 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 +121 @@
     }
 
+    fRates = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "AccidentalPhotonRates");
+    if (!fRates)
+        return kFALSE;
+
     /*
     fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader");
@@ -132 +144 @@
     }
 
+    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 +233 @@
     for (UInt_t idx=0; idx<npix; idx++)
     {
-        // Scale the rate with the pixel size. The rate must
-        // always be given for the pixel with index 0.
-        const Double_t avglen = 1./(fFreqFixed+fFreqNSB*(*fGeom)[idx].GetA()/100.);
+        // Scale the rate with the pixel size.
+        const Double_t rate = fFreqFixed+fFreqNSB*(*fGeom)[idx].GetA()*fScale;
+
+        (*fRates)[idx].SetPedestal(rate);
+
+        // Calculate the average distance between two consequtive photons
+        const Double_t avglen = 1./rate;
 
         // Start producing photons at time "start"
@@ -215 +292 @@
 //    FrequencyNSB:   0.040
 //
-// The frequency is given in units fitting the units of the time.
-// Usually the time is given in nanoseconds thus, e.g., 40 means 40MHz
+// The fixed frequency is given in units fitting the units of the time.
+// Usually the time is given in nanoseconds thus, e.g., 0.040 means 40MHz.
+//
+// The FrequencyNSB is scaled by the area of the pixel in cm^2. Therefore
+// 0.040 would mean 40MHz/cm^2
 //
 Int_t MSimRandomPhotons::ReadEnv(const TEnv &env, TString prefix, Bool_t print)