Changeset 9425 for trunk/MagicSoft


Ignore:
Timestamp:
04/16/09 12:04:29 (16 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mjobs/MJSimulation.cc

    r9420 r9425  
    108108#include "MReflector.h"
    109109#include "MParEnv.h"
    110 #include "MPulseShape.h"
     110#include "MSpline3.h"
     111#include "MParSpline.h"
    111112#include "MGeomCam.h"
    112113
    … …  
    152153    //cont.Add(const_cast<MSequence*>(&fSequence));
    153154
    154     cont.Add(plist.FindObject("MPulseShape"));
     155    cont.Add(plist.FindObject("PulseShape"));
    155156    cont.Add(plist.FindObject("Reflector"));
    156157    cont.Add(plist.FindObject("MGeomCam"));
    … …  
    185186    hist.SetLog(kTRUE, kTRUE, kFALSE);
    186187
    187     hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100");
    188     hist.SetDrawOption("colz");
    189     hist.InitName("Impact;Impact;Impact");
    190     hist.InitTitle("Impact;West <--> East [m];South <--> North [m]");
     188    /*
     189     hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100");
     190     hist.SetDrawOption("colz");
     191     hist.InitName("Impact;Impact;Impact");
     192     hist.InitTitle("Impact;West <--> East [m];South <--> North [m]");
     193     hist.SetAutoRange();
     194     */
     195
     196    hist.AddHist("MCorsikaEvtHeader.GetImpact/100");
     197    hist.InitName("Impact");
     198    hist.InitTitle("Impact;Impact [m]");
    191199    hist.SetAutoRange();
    192200
    … …  
    197205    hist.AddHist("(MCorsikaEvtHeader.fAz+MCorsikaRunHeader.fMagneticFieldAz)*TMath::RadToDeg()", "MCorsikaEvtHeader.fZd*TMath::RadToDeg()");
    198206    hist.InitName("SkyOrigin;Az;Zd");
    199     hist.InitTitle("Sky Origin;Az [deg];Zd [deg]");
     207    hist.InitTitle("Sky Origin;Az [\\deg];Zd [\\deg]");
    200208    hist.SetDrawOption("colz");
    201209    hist.SetAutoRange();
    … …  
    279287    plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm");
    280288
    281     MPulseShape shape;
     289    MParSpline shape("PulseShape");
    282290    plist.AddToList(&shape);
    283291
    … …  
    325333    read.SetForceMode(fForceMode);
    326334
    327     for (int i=0; i<args.GetNumArguments();i ++)
     335    for (int i=0; i<args.GetNumArguments(); i++)
    328336        read.AddFile(args.GetArgumentStr(i));
    329337
    330338    MSimMMCS simmmcs;
    331339
     340    MParSpline splinepde("PhotonDetectionEfficiency");
     341    MParSpline splinemirror("MirrorReflectivity");
     342    MParSpline splinecones("ConesAngularAcceptance");
     343    plist.AddToList(&splinepde);
     344    plist.AddToList(&splinemirror);
     345    plist.AddToList(&splinecones);
     346
    332347    MSimAtmosphere simatm;
    333     MSimAbsorption absapd("PhotonDetectionEfficiency");
    334     MSimAbsorption absmir("MirrorReflectivity");
    335     MSimAbsorption cones("ConesAngularAcceptance");
     348    MSimAbsorption absapd("SimPhotonDetectionEfficiency");
     349    MSimAbsorption absmir("SimMirrorReflectivity");
     350    MSimAbsorption cones("SimConesAngularAcceptance");
     351    absapd.SetParName("PhotonDetectionEfficiency");
     352    absmir.SetParName("MirrorReflectivity");
     353    cones.SetParName("ConesAngularAcceptance");
    336354    cones.SetUseTheta();
    337355
    … …  
    344362    MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
    345363    MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");
     364    MContinue cont3("MPhotonEvent.GetNumPhotons<2", "ContEmpty3");
    346365
    347366    // -------------------------------------------------------------------
    348367
    349368    MBinning binse( 100,     1,   100000, "BinningEnergy",    "log");
    350     MBinning binsth( 35,   0.9,   900000, "BinningThreshold", "log");
     369    MBinning binsth( 70,   0.9,   900000, "BinningThreshold", "log");
    351370    MBinning binsee( 35,   0.9,   900000, "BinningEnergyEst", "log");
    352371    MBinning binss( 100,     1, 10000000, "BinningSize",      "log");
    353     MBinning binsi( 100,  -500,      500, "BinningImpact");
     372//    MBinning binsi( 100,  -500,      500, "BinningImpact");
     373    MBinning binsi(  55,     0,     1100, "BinningImpact");
    354374    MBinning binsh( 150,     0,       50, "BinningHeight");
    355375    MBinning binsaz(720,  -360,      360, "BinningAz");
    … …  
    377397    plist.AddToList(&binsd0);
    378398
    379     MHn mhn1, mhn2, mhn3;
     399    MHn mhn1, mhn2, mhn3, mhn4;
    380400    SetupHist(mhn1);
    381401    SetupHist(mhn2);
    382402    SetupHist(mhn3);
    383 
    384     MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-MPulseShape.GetPulseWidth");
     403    SetupHist(mhn4);
     404
     405    MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-PulseShape.GetWidth");
    385406    mhtp.SetName("TrigPos");
    386407    mhtp.SetTitle("Trigger position w.r.t. the first photon hitting a detector");
    … …  
    397418    MFillH fillh2(&mhn2, "", "FillH2");
    398419    MFillH fillh3(&mhn3, "", "FillH3");
     420    MFillH fillh4(&mhn4, "", "FillH4");
    399421    MFillH filltp(&mhtp, "", "FillTriggerPos");
    400422    MFillH fillew(&mhew, "", "FillEvtWidth");
    … …  
    403425    fillh2.SetNameTab("Detectable", "Distribution of events hit the detector");
    404426    fillh3.SetNameTab("Triggered",  "Distribution of triggered events");
     427    fillh4.SetNameTab("Cleaned",    "Distribution after cleaning and cuts");
    405428    filltp.SetNameTab("TrigPos",    "TriggerPosition w.r.t the first photon");
    406429    fillew.SetNameTab("EvtWidth",   "Time between first and last photon hitting a detector");
    … …  
    491514    simcal.SetNameGeomCam("GeomCones");
    492515
    493     //  Dark Current: 4MHz = 0.004 GHz
    494     //  Light of night sky at La Palma: ~ 0.2 ph / cm^2 / sr / ns
    495     //  5deg camera: 6e-3 sr
    496     //  NSB = 0.2*6e-3
    497 
    498     // MAGIC: 84MHz       / Pixel
    499     // DWARF: 14MHz-25Mhz / APD
    500 
    501516    // FIXME: Simulate photons before cones and QE!
    502517
    503518    MSimRandomPhotons  simnsb;
    504     simnsb.SetFreq(0.0015, 0.04);  // 1.5MHz dark count rate, 40MHZ/cm^2 NSB
     519    simnsb.SetFreq(5.8, 0.004);  // ph/m^2/nm/sr/ns NSB, 4MHz dark count rate
    505520    simnsb.SetNameGeomCam("GeomCones");
    506521
    … …  
    509524
    510525    MSimAPD simapd;
    511     //simapd.SetFreq(0.04); // Must be identical to MSimRandomPhotons!!
    512526    simapd.SetNameGeomCam("GeomCones");
    513527
    … …  
    555569    MFillH fillx0d(&evt0d,             "MSignalCam",      "FillArrTm");
    556570    MFillH fillx1("MHHillas",          "MHillas",         "FillHillas");
    557     //MFillH fillx2("MHHillasExt",       "",                "FillHillasExt");
    558571    MFillH fillx3("MHHillasSrc",       "MHillasSrc",      "FillHillasSrc");
     572    MFillH fillx4("MHNewImagePar",     "MNewImagePar",    "FillNewImagePar");
    559573    MFillH fillth("MHThreshold",       "",                "FillThreshold");
    560574    MFillH fillca("MHCollectionArea",  "",                "FillTrigArea");
    561     //MFillH fillx4("MHImagePar",        "MImagePar",       "FillImagePar");
    562     //MFillH fillx5("MHNewImagePar",     "MNewImagePar",    "FillNewImagePar");
    563575
    564576    fillth.SetNameTab("Threshold");
    … …  
    603615        }
    604616        tasks.AddToList(&cones);
    605         tasks.AddToList(&fillF2);
    606617        //if (header.IsPointRun())
    607618        //    tasks.AddToList(&fillP);
    608619        tasks.AddToList(&cont1);
    609         if (!header.IsPointRun())
    610             tasks.AddToList(&fillh2);
    611620    }
    612621    // -------------------------------
    … …  
    617626        tasks.AddToList(&simgeom);
    618627        tasks.AddToList(&cont2);
     628        if (!header.IsPointRun())
     629        {
     630            tasks.AddToList(&fillF2);
     631            tasks.AddToList(&fillh2);
     632        }
     633        tasks.AddToList(&cont3);
    619634    }
    620635    if (fCamera)
    … …  
    669684        //tasks.AddToList(&fillx2);
    670685        tasks.AddToList(&fillx3);
    671         //tasks.AddToList(&fillx4);
     686        tasks.AddToList(&fillx4);
    672687        //tasks.AddToList(&fillx5);
    673688    }
    674689    if (header.IsDataRun())
    675690    {
     691        tasks.AddToList(&fillh4);
    676692        tasks.AddToList(&fillth);
    677693        tasks.AddToList(&fillca);
    … …  
    695711
    696712    if (binstr.IsDefault())
    697         binstr.SetEdgesLin(150, -shape.GetPulseWidth(),
    698                            header.GetFreqSampling()/1000.*header.GetNumSamples()+shape.GetPulseWidth());
     713        binstr.SetEdgesLin(150, -shape.GetWidth(),
     714                           header.GetFreqSampling()/1000.*header.GetNumSamples()+shape.GetWidth());
    699715
    700716    if (binsd.IsDefault() && cam)
    … …  
    731747            gROOT->SetSelectedPad(0);
    732748            static_cast<MReflector*>(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete);
    733 
    734749        }
    735750
    … …  
    749764                c->SetBit(kCanDelete);
    750765                c->Draw();
    751 
    752766            }
    753767
    … …  
    767781            }
    768782        }
     783
     784        TCanvas &d = fDisplay->AddTab("Info2");
     785        d.Divide(2,2);
     786
     787        d.cd(1);
     788        gPad->SetBorderMode(0);
     789        gPad->SetFrameBorderMode(0);
     790        gPad->SetGridx();
     791        gPad->SetGridy();
     792        gROOT->SetSelectedPad(0);
     793        splinepde.DrawClone()->SetBit(kCanDelete);
     794
     795        d.cd(2);
     796        gPad->SetBorderMode(0);
     797        gPad->SetFrameBorderMode(0);
     798        gPad->SetGridx();
     799        gPad->SetGridy();
     800        gROOT->SetSelectedPad(0);
     801        splinemirror.DrawClone()->SetBit(kCanDelete);
     802
     803        d.cd(3);
     804        gPad->SetBorderMode(0);
     805        gPad->SetFrameBorderMode(0);
     806        gPad->SetGridx();
     807        gPad->SetGridy();
     808        gROOT->SetSelectedPad(0);
     809        splinecones.DrawClone()->SetBit(kCanDelete);
     810
     811        d.cd(4);
     812        gPad->SetBorderMode(0);
     813        gPad->SetFrameBorderMode(0);
     814        gPad->SetGridx();
     815        gPad->SetGridy();
     816        gROOT->SetSelectedPad(0);
     817        MParSpline *all = (MParSpline*)splinepde.DrawClone();
     818        //all->SetTitle("Combined acceptance");
     819        all->SetBit(kCanDelete);
     820        all->Multiply(*splinemirror.GetSpline());
    769821    }
    770822

  • trunk/MagicSoft/Mars/msim/MSimAbsorption.cc

    r9348 r9425  
    3030//
    3131//  Input Containers:
     32//   fParName [MParSpline]
    3233//   MPhotonEvent
    3334//   MCorsikaEvtHeader
    … …  
    4344
    4445#include <TRandom.h>
    45 #include <TSpline.h>
    4646
    4747#include "MLog.h"
    … …  
    5050#include "MParList.h"
    5151
    52 #include "MBinning.h"
    53 #include "MArrayD.h"
    54 
    55 #include "MSpline3.h"
     52#include "MParSpline.h"
    5653
    5754#include "MCorsikaEvtHeader.h"
    … …  
    6966//
    7067MSimAbsorption::MSimAbsorption(const char* name, const char *title)
    71     : fEvt(0), fHeader(0), fSpline(0), fUseTheta(kFALSE)
     68    : fEvt(0), fHeader(0), fSpline(0), fParName("MParSpline"), fUseTheta(kFALSE)
    7269{
    7370    fName  = name  ? name  : "MSimAbsorption";
    7471    fTitle = title ? title : "Task to calculate wavelength dependent absorption";
    75 }
    76 
    77 // --------------------------------------------------------------------------
    78 //
    79 //  Calls Clear()
    80 //
    81 MSimAbsorption::~MSimAbsorption()
    82 {
    83     Clear();
    84 }
    85 
    86 // --------------------------------------------------------------------------
    87 //
    88 // Delete fSpline and set to 0.
    89 //
    90 void MSimAbsorption::Clear(Option_t *)
    91 {
    92     if (fSpline)
    93         delete fSpline;
    94     fSpline=0;
    95 }
    96 
    97 // --------------------------------------------------------------------------
    98 //
    99 // Read a TGraph from a file and return a newly allocated MSpline3.
    100 //
    101 MSpline3 *MSimAbsorption::ReadSpline(const char *fname)
    102 {
    103     *fLog << inf << "Reading spline from " << fname << endl;
    104 
    105     TGraph g(fname);
    106     if (g.GetN()==0)
    107     {
    108         *fLog << err << "ERROR - No data points from " << fname << "." << endl;
    109         return 0;
    110     }
    111 
    112     // option: b1/e1 b2/e2   (first second derivative?)
    113     // option: valbeg/valend (first second derivative?)
    114 
    115     return new MSpline3(g);
    116 }
    117 
    118 // --------------------------------------------------------------------------
    119 //
    120 // Initializes a spline from min to max with n points with 1.
    121 //
    122 void MSimAbsorption::InitUnity(UInt_t n, Float_t min, Float_t max)
    123 {
    124     // Delete the existing spline
    125     Clear();
    126 
    127     // We need at least two points (the edges)
    128     if (n<2)
    129         return;
    130 
    131     // Initialize an array with the x-values
    132     const MBinning bins(n-1, min, max);
    133 
    134     // Initialize an array with all one
    135     MArrayD y(n);
    136     y.Reset(1);
    137 
    138     // Set the spline
    139     fSpline = new MSpline3(bins.GetEdges(), y.GetArray(), n);
    140 }
    141 
    142 // --------------------------------------------------------------------------
    143 //
    144 // Takes the existing fSpline and multiplies it with the given spline.
    145 // As x-points the values from fSpline are used.
    146 //
    147 void MSimAbsorption::Multiply(const TSpline3 &spline)
    148 {
    149     if (!fSpline)
    150     {
    151         Clear();
    152         // WARNING WARNING WARNING: This is a very dangerous cast!
    153         fSpline = (MSpline3*)spline.Clone();
    154         return;
    155     }
    156 
    157     // Initialize a TGraph with the number of knots from a TSpline
    158     TGraph g(fSpline->GetNp());
    159 
    160     // Loop over all knots
    161     for (int i=0; i<fSpline->GetNp(); i++)
    162     {
    163         // Get th i-th knot
    164         Double_t x, y;
    165         fSpline->GetKnot(i, x, y);
    166 
    167         // Multiply y by the value from the given spline
    168         y *= spline.Eval(x);
    169 
    170         // push the point "back"
    171         g.SetPoint(i, x, y);
    172     }
    173 
    174     // Clear the spline and initialize a new one from the updated TGraph
    175     Clear();
    176     fSpline = new MSpline3(g);
    177 }
    178 
    179 // --------------------------------------------------------------------------
    180 //
    181 // Initializes a TSpline3 with n knots x and y. Call Multiply for it.
    182 //
    183 void MSimAbsorption::Multiply(UInt_t n, const Double_t *x, const Double_t *y)
    184 {
    185     const TSpline3 spline("Spline", const_cast<Double_t*>(x), const_cast<Double_t*>(y), n);
    186     Multiply(spline);
    187 }
    188 
    189 // --------------------------------------------------------------------------
    190 //
    191 // Read a Spline from a file using ReadSpline.
    192 // Call Multiply for it.
    193 //
    194 void MSimAbsorption::Multiply(const char *fname)
    195 {
    196     const MSpline3 *spline = ReadSpline(fname);
    197     if (!spline)
    198         return;
    199 
    200     fFileName = "";
    201 
    202     Multiply(*spline);
    203 
    204     delete spline;
    205 }
    206 
    207 // --------------------------------------------------------------------------
    208 //
    209 // Read a spline from a file and set fSpline accfordingly.
    210 //
    211 Bool_t MSimAbsorption::ReadFile(const char *fname)
    212 {
    213     if (fname)
    214         fFileName = fname;
    215 
    216     MSpline3 *spline = ReadSpline(fFileName);
    217     if (!spline)
    218         return kFALSE;
    219 
    220 
    221     // option: b1/e1 b2/e2   (first second derivative?)
    222     // option: valbeg/valend (first second derivative?)
    223 
    224     Clear();
    225     fSpline = spline;
    226 
    227     return kTRUE;
    22872}
    22973
    … …  
    23579Int_t MSimAbsorption::PreProcess(MParList *pList)
    23680{
    237     if (fFileName.IsNull())
    238     {
    239         *fLog << inf << "No file given... skipping." << endl;
    240         return kSKIP;
    241     }
    242 
    24381    fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
    24482    if (!fEvt)
    24583    {
    24684        *fLog << err << "MPhotonEvent not found... aborting." << endl;
     85        return kFALSE;
     86    }
     87
     88    fSpline = (MParSpline*)pList->FindObject(fParName, "MParSpline");
     89    if (!fSpline)
     90    {
     91        *fLog << err << fParName << " [MParSpline] not found... aborting." << endl;
     92        return kFALSE;
     93    }
     94
     95    if (!fSpline->IsValid())
     96    {
     97        *fLog << err << "Spline in " << fParName << " is inavlid... aborting." << endl;
    24798        return kFALSE;
    24899    }
    … …  
    257108    *fLog << inf << "Using " << (fUseTheta?"Theta":"Wavelength") << " for absorption." << endl;
    258109
    259     return ReadFile();
     110    return kTRUE;
    260111}
    261112
    … …  
    331182{
    332183    Bool_t rc = kFALSE;
    333     if (IsEnvDefined(env, prefix, "FileName", print))
     184    if (IsEnvDefined(env, prefix, "ParName", print))
    334185    {
    335186        rc = kTRUE;
    336         SetFileName(GetEnvValue(env, prefix, "FileName", fFileName));
     187        SetParName(GetEnvValue(env, prefix, "ParName", fParName));
    337188    }
    338189

  • trunk/MagicSoft/Mars/msim/MSimAbsorption.h

    r9232 r9425  
    66#endif
    77
    8 class TSpline3;
    9 
    108class MParList;
    11 class MSpline3;
     9class MParSpline;
    1210class MPhotonEvent;
    1311class MCorsikaEvtHeader;
    … …  
    1917    MCorsikaEvtHeader *fHeader;  //! Header storing event information
    2018
    21     MSpline3          *fSpline;  //! Spline to interpolate wavelength or incident angle
     19    MParSpline        *fSpline;  //! Spline to interpolate wavelength or incident angle
    2220
    23     TString fFileName;           // FileName of the file containing the curve
     21    TString fParName;            // Container name of the spline containing the curve
    2422    Bool_t  fUseTheta;           // Switches between using wavelength or incident angle
    2523
    … …  
    3230    Int_t  Process();
    3331
    34     // MSimAbsorption
    35     MSpline3 *ReadSpline(const char *fname);
    36 
    3732public:
    3833    MSimAbsorption(const char *name=NULL, const char *title=NULL);
    39     ~MSimAbsorption();
    4034
    4135    // MSimAbsorption
    42     Bool_t ReadFile(const char *fname=0);
    43     void SetFileName(const char *fname) { fFileName=fname; }
     36    void SetParName(const char *name) { fParName=name; }
    4437
    4538    void SetUseTheta(Bool_t b=kTRUE) { fUseTheta = b; }
    46 
    47     void InitUnity(UInt_t n, Float_t min, Float_t max);
    48 
    49     void Multiply(const char *fname);
    50     void Multiply(const TSpline3 &spline);
    51     void Multiply(UInt_t n, const Double_t *x, const Double_t *y);
    52 
    53     // TObject
    54     void Clear(Option_t *o="");
    5539
    5640    ClassDef(MSimAbsorption, 0) // Task to calculate wavelength or incident angle dependent absorption

  • trunk/MagicSoft/Mars/msimcamera/MSimCalibrationSignal.cc

    r9378 r9425  
    6161
    6262#include "MGeomCam.h"
    63 #include "MPulseShape.h"
     63#include "MParSpline.h"
    6464#include "MTriggerPattern.h"
    6565
    … …  
    136136    }
    137137
    138     fPulse = (MPulseShape*)pList->FindObject("MPulseShape");
     138    fPulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline");
    139139    if (!fPulse)
    140140    {
    141         *fLog << err << "MPulsShape not found... aborting." << endl;
     141        *fLog << err << "PulsShape [MParSpline] not found... aborting." << endl;
    142142        return kFALSE;
    143143    }
    … …  
    173173    }
    174174
     175    // FIXME: Is there a way to write them as LAST entry in the file?
    175176    fRunHeader->SetReadyToSave();
    176177
    … …  
    195196        for (UInt_t idx=0; idx<fGeom->GetNumPixels(); idx++)
    196197        {
     198            // FIXME: Scale number of photons with the pixel size!
    197199            const Int_t num = TMath::Nint(gRandom->Gaus(fNumPhotons, fNumPhotons/10));
    198200
    … …  
    229231    // Length (ns), Pulse position (Units ns)
    230232    const Float_t pp = fPulsePos->GetVal();
    231     const Float_t pw = fPulse->GetPulseWidth();
     233    const Float_t pw = fPulse->GetWidth();
    232234
    233235    const Float_t first = cnt>0 ? fEvt->GetFirst()->GetTime() : 0;

  • trunk/MagicSoft/Mars/msimcamera/MSimCalibrationSignal.h

    r9306 r9425  
    88class MGeomCam;
    99class MParList;
    10 class MPulseShape;
     10class MParSpline;
    1111class MPhotonEvent;
    1212class MPhotonStatistics;
    … …  
    2020    MParList          *fParList;    //! Store pointer to MParList for initializing ReInit
    2121    MGeomCam          *fGeom;       //! Camera geometry to know the number of expected pixels
    22     MPulseShape       *fPulse;      //! Pulse Shape to get pulse width from
     22    MParSpline        *fPulse;      //! Pulse Shape to get pulse width from
    2323    MParameterD       *fPulsePos;   //! Expected position at which the pulse should be
    2424    MParameterD       *fTrigger;    //! Position in analog channels at which the triggersignal  is raised

  • trunk/MagicSoft/Mars/msimcamera/MSimCamera.cc

    r9413 r9425  
    4848
    4949#include "MSpline3.h"
    50 #include "MPulseShape.h"
     50#include "MParSpline.h"
    5151
    5252#include "MParList.h"
    53 //#include "MParameters.h"
    5453
    5554#include "MPhotonEvent.h"
    … …  
    134133        return kFALSE;
    135134
    136     MPulseShape *pulse = (MPulseShape*)pList->FindObject("MPulseShape");
     135    MParSpline *pulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline");
    137136    if (!pulse)
    138137    {
    139         *fLog << err << "MPulseShape not found... aborting." << endl;
     138        *fLog << err << "PulseShape [MParSpline] not found... aborting." << endl;
    140139        return kFALSE;
    141140    }

  • trunk/MagicSoft/Mars/msimcamera/MSimCamera.h

    r9328 r9425  
    2828    MMcEvt            *fMcEvt;           //! For information stored in MMcEvt
    2929
    30     MSpline3 *fSpline;     // Pusle Shape
     30    const MSpline3    *fSpline;          // Pusle Shape
    3131
    3232    Bool_t fBaselineGain;  // Should the gain be applied to baseline and electronic noise?

  • trunk/MagicSoft/Mars/msimcamera/MSimGeomCam.cc

    r9374 r9425  
    6060
    6161#include "MParameters.h"
     62#include "MParSpline.h"
    6263#include "MRawRunHeader.h"
    63 
    64 #include "MPulseShape.h"
    6564
    6665ClassImp(MSimGeomCam);
    … …  
    106105    }
    107106
    108     fPulse = (MPulseShape*)pList->FindObject("MPulseShape");
     107    fPulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline");
    109108/*
    110109    if (!fPulse)
    … …  
    191190
    192191    // Length (ns), Pulse position (Units ns)
    193     const Float_t pp   = fPulsePos ? fPulsePos->GetVal()     : 0;
    194     const Float_t pw   = fPulse    ? fPulse->GetPulseWidth() : 0;
     192    const Float_t pp   = fPulsePos ? fPulsePos->GetVal() : 0;
     193    const Float_t pw   = fPulse    ? fPulse->GetWidth() : 0;
    195194
    196195    fStat->SetTimeMedDev(fEvt->GetTimeMedianDev());

  • trunk/MagicSoft/Mars/msimcamera/MSimGeomCam.h

    r9274 r9425  
    1313class MParameterD;
    1414class MRawRunHeader;
    15 class MPulseShape;
     15class MParSpline;
    1616
    1717class MSimGeomCam : public MTask
    … …  
    2323    MParameterD       *fPulsePos; //! Intended pulse position in digitization window [ns]
    2424    MRawRunHeader     *fHeader;   //! Length of digitization window
    25     MPulseShape       *fPulse;    //!
     25    MParSpline        *fPulse;    //!
    2626
    2727    TString fNameGeomCam;

  • trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc

    r9375 r9425  
    4343//  Output Containers:
    4444//   MPhotonEvent
     45//   AccidentalPhotonRate [MPedestalCam]
    4546//
    4647//////////////////////////////////////////////////////////////////////////////
    … …  
    6263#include "MPhotonData.h"
    6364
     65#include "MPedestalCam.h"
     66#include "MPedestalPix.h"
     67
    6468#include "MCorsikaRunHeader.h"
     69
     70#include "MSpline3.h"
     71#include "MParSpline.h"
     72#include "MReflector.h"
    6573
    6674ClassImp(MSimRandomPhotons);
    … …  
    7482MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title)
    7583    : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0),
    76     fSimulateWavelength(kFALSE),  fNameGeomCam("MGeomCam")
     84    fRates(0), fSimulateWavelength(kFALSE),  fNameGeomCam("MGeomCam")
    7785{
    7886    fName  = name  ? name  : "MSimRandomPhotons";
    … …  
    113121    }
    114122
     123    fRates = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "AccidentalPhotonRates");
     124    if (!fRates)
     125        return kFALSE;
     126
    115127    /*
    116128    fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader");
    … …  
    132144    }
    133145
     146    MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector");
     147    if (!r)
     148    {
     149        *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
     150        return kFALSE;
     151    }
     152
     153    const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline");
     154    const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesAngularAcceptance",    "MParSpline");
     155    const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity",        "MParSpline");
     156
     157    const Double_t d2  = fGeom->GetCameraDist()*fGeom->GetCameraDist();
     158    const Double_t pde = s1 && s1->GetSpline() ? s1->GetSpline()->Integral() : 1;
     159    const Double_t sr  = s2 && s2->GetSpline() ? s2->GetSpline()->IntegralSolidAngle() : 1;
     160    const Double_t mir = s3 && s3->GetSpline() ? s3->GetSpline()->Integral() : 1;
     161    const Double_t Ar  = r->GetA()/1e4;
     162
     163    // Conversion factor to convert pixel area to steradians (because it
     164    // is a rather small area we can assume it is flat)
     165    const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad();
     166
     167    // Multiply all relevant efficiencies
     168    MParSpline *s4 = (MParSpline*)s1->Clone();
     169    s4->Multiply(*s3->GetSpline());
     170
     171    const Double_t nm = s4 && s4->GetSpline() ? s4->GetSpline()->Integral() : 1;
     172
     173    delete s4;
     174
     175    // /100 to convert the pixel area from mm^2 to cm^2
     176    fScale =  nm * TMath::Min(Ar, sr*d2) * conv*conv;
     177
     178    *fLog << inf;
     179    *fLog << "Effective cone acceptance:      " << Form("%.2f", sr*d2) << "m^2" << endl;
     180    *fLog << "Reflector area:                 " << Form("%.2f", Ar) << "m^2" << endl;
     181    *fLog << "Resulting eff. collection area: " << Form("%.1f", TMath::Min(Ar, sr*d2)) << "m^2" << endl;
     182    *fLog << "Eff. wavelength band (PDE):     " << Form("%.1f", pde) << "nm" << endl;
     183    *fLog << "Eff. wavelength band (Mirror):  " << Form("%.1f", mir) << "nm" << endl;
     184    *fLog << "Eff. wavelength band (PDE+MIR): " << Form("%.1f", nm) << "nm" << endl;
     185    *fLog << "Pixel area of " << fNameGeomCam << "[0]: " << Form("%.2e", (*fGeom)[0].GetA()*conv*conv) << "sr" << endl;
     186    //*fLog << "Effective angular acceptance:     " << sr << "sr" << endl;
     187    //*fLog << "Resulting NSB frequency:           " << fFreqNSB*nm*Ar*1000 << "MHz/sr" << endl;
     188    *fLog << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", fFreqNSB*(*fGeom)[0].GetA()*fScale*1000) << "MHz" << endl;
     189
     190    // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
     191    // Set NumPheFromDNSB
     192
     193    // # Number of photons from the diffuse NSB (nphe / ns 0.1*0.1 deg^2 239 m^2) and
     194    // nsb_mean 0.20
     195    // Magic pixel: 0.00885361 deg
     196    // dnsbpix = 0.2*50/15
     197    // ampl = MMcFadcHeader->GetAmplitud()
     198    // sqrt(pedrms*pedrms + dnsbpix*ampl*ampl/ratio)
     199
     200    return kTRUE;
     201}
     202
     203Bool_t MSimRandomPhotons::ReInit(MParList *pList)
     204{
     205    // Overwrite the default set by MGeomApply
     206    fRates->Init(*fGeom);
    134207    return kTRUE;
    135208}
    … …  
    160233    for (UInt_t idx=0; idx<npix; idx++)
    161234    {
    162         // Scale the rate with the pixel size. The rate must
    163         // always be given for the pixel with index 0.
    164         const Double_t avglen = 1./(fFreqFixed+fFreqNSB*(*fGeom)[idx].GetA()/100.);
     235        // Scale the rate with the pixel size.
     236        const Double_t rate = fFreqFixed+fFreqNSB*(*fGeom)[idx].GetA()*fScale;
     237
     238        (*fRates)[idx].SetPedestal(rate);
     239
     240        // Calculate the average distance between two consequtive photons
     241        const Double_t avglen = 1./rate;
    165242
    166243        // Start producing photons at time "start"
    … …  
    215292//    FrequencyNSB:   0.040
    216293//
    217 // The frequency is given in units fitting the units of the time.
    218 // Usually the time is given in nanoseconds thus, e.g., 40 means 40MHz
     294// The fixed frequency is given in units fitting the units of the time.
     295// Usually the time is given in nanoseconds thus, e.g., 0.040 means 40MHz.
     296//
     297// The FrequencyNSB is scaled by the area of the pixel in cm^2. Therefore
     298// 0.040 would mean 40MHz/cm^2
    219299//
    220300Int_t MSimRandomPhotons::ReadEnv(const TEnv &env, TString prefix, Bool_t print)

  • trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h

    r9241 r9425  
    1212//class MCorsikaEvtHeader;
    1313class MCorsikaRunHeader;
     14class MPedestalCam;
    1415
    1516class MSimRandomPhotons : public MTask
    … …  
    2122//    MCorsikaEvtHeader *fEvtHeader;  //! Header storing event information
    2223    MCorsikaRunHeader *fRunHeader;  //! Header storing run information
    23 
     24    MPedestalCam      *fRates;   // Random count rate per pixel
    2425
    2526    // FIXME: Make this a single number per Pixel/APD
    2627    Double_t fFreqFixed; // [1/ns]      A fixed frequency per pixel
    2728    Double_t fFreqNSB;   // [1/ns/cm^2] A frequency depending on area
     29
     30    Double_t fScale;
    2831
    2932    Bool_t fSimulateWavelength;
    … …  
    3235
    3336    // MTask
    34     Int_t PreProcess(MParList *pList);
    35     Int_t Process();
     37    Int_t  PreProcess(MParList *pList);
     38    Bool_t ReInit(MParList *pList);
     39    Int_t  Process();
    3640
    3741    // MParContainer
Note: See TracChangeset for help on using the changeset viewer.