Changeset 9565 for trunk


Ignore:
Timestamp:
03/30/10 14:47:27 (15 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r9564 r9565  
    4040   * msimreflector/MSimReflector.[h,cc]:
    4141     - added name of reflector as data member
     42
     43   * melectronics/MAvalanchePhotoDiode.cc:
     44     - scale the crosstalk probability as the height of the emitted
     45       signal with the recovery time.
     46
     47   * mhflux/MHEnergyEst.cc:
     48     - added a workaround to get rid of some root-warnings
     49
     50   * mjtrain/MJTrainEnergy.cc:
     51     - improved axis labels
     52
     53   * mpointing/MPointingDevCalc.cc:
     54     - added some more comments
     55
     56   * msimreflector/MMirror.[h,cc]:
     57     - added fShape to allow for parabolic mirrors
     58
     59   * msimreflector/MReflector.cc:
     60     - implemented the option of parabolic mirrors into the
     61       reflector defintion file
     62
     63   * msimreflector/MSimReflector.[h,cc]:
     64     - implemented (and fixed) the calculation of the reflection at
     65       parabolic surfaces
     66     - implemented the possibility to switch off the early check for
     67       "photon can hit the camera" (fDetectorMargin<0)
     68     - added setter for fDetectorMargin
     69
    4270
    4371

  • trunk/MagicSoft/Mars/NEWS

    r9564 r9565  
    77   * On some systems (version of make?) the __LINUX__ definition was not set
    88     when compiling. This lead to MTime(-1) not working properly. Fixed.
     9
     10 ;ceres:
     11
     12   * Allow for individual mirrors with parabolic shape (for details see
     13     MReflector::ReadFile)
    914
    1015

  • trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc

    r9462 r9565  
    123123    //            return 0;
    124124
     125    // Number of the x/y cell in the one dimensional array
    125126    // const Int_t cell = fHist.GetBin(x, y);
    126127    const Int_t cell = x + (fHist.GetNbinsX()+2)*y;
    127128
    128     // This is the fastes way to access the bin-contents in fArray
     129    // Getting a reference to the float is the fastes way to
     130    // access the bin-contents in fArray
    129131    Float_t &cont = fHist.GetArray()[cell];
    130132
     133    // Calculate the time since the last breakdown
    131134    // const Double_t dt = t-fHist.GetBinContent(x, y)-fDeadTime; //
    132135    const Float_t dt = t-cont-fDeadTime;
    … …  
    136139        return 0;
    137140
    138     // Signal height (in units of one photon) produced after dead time
    139     const Float_t weight = fRecoveryTime<=0 ? 1 : 1.-exp(-dt/fRecoveryTime);
    140 
     141    // The signal height (in units of one photon) produced after dead time
     142    // depends on the recovery of the cell - described by an exponential.
     143    const Float_t weight = fRecoveryTime<=0 ? 1. : 1-TMath::Exp(-dt/fRecoveryTime);
     144
     145    // The probability that the cell emits a photon causing crosstalk
     146    // scales as the signal height.
     147    const Float_t prob = weight*fCrosstalkProb;
     148
     149    // Set the contents to the time of the last breakdown (now)
    141150    cont = t; // fHist.SetBinContent(x, y, t)
    142151
    … …  
    161170     */
    162171
     172
    163173    //for (int i=0; i<1; i++)
    164174    while (1)
    165175    {
    166176        const Double_t rndm = gRandom->Rndm();
    167         if (rndm>=fCrosstalkProb)
     177        if (rndm>=prob/*fCrosstalkProb*/)
    168178            break;
    169179
    170         // We can re-use the random number becuase it is uniformely
     180        // We can re-use the random number because it is uniformely
    171181        // distributed. This saves cpu power
    172         const Int_t dir = TMath::FloorNint(4*rndm/fCrosstalkProb);
     182        const Int_t dir = TMath::FloorNint(4*rndm/prob/*fCrosstalkProb*/);
    173183
    174184        // Get a random neighbor which is hit.

  • trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

    r9301 r9565  
    298298        pad->GetPad(1)->cd(1);
    299299
    300         if (gPad->FindObject("EnergyEst_ex"))
     300        if ((hx = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ex"))))
     301        {
     302            hx->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
    301303            hx = fHEnergy.Project3D("ex");
    302 
    303         if (gPad->FindObject("EnergyEst_ey"))
     304        }
     305
     306        if ((hy = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ey"))))
     307        {
     308            hy->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
    304309            hy = fHEnergy.Project3D("ey");
     310        }
    305311
    306312        if (hx && hy)
    … …  
    317323        {
    318324            pad->GetPad(1)->GetPad(2)->cd(1);
    319             if (gPad->FindObject("EnergyEst_ez"))
     325            if ((hx = dynamic_cast<TH1*>(gPad->FindObject("EnergyEst_ez"))))
     326            {
     327                hx->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated
    320328                fHEnergy.Project3D("ez");
     329            }
    321330
    322331            pad->GetPad(1)->GetPad(2)->cd(2);

  • trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc

    r8888 r9565  
    256256    hres2.SetDrawOption("colz profx");
    257257
    258     MHn hres3("Resolution", "Energy Resolution");
     258    MHn hres3("Resolution", "Energy Resolution squared");
    259259    hres3.AddHist("MEnergyEst.fVal", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
    260260    hres3.InitName("ResEest;EnergyEst;");
    261     hres3.InitTitle(";E_{est} [GeV];Resolution (E_{mc}/E_{est}-1)^{2};");
     261    hres3.InitTitle(";E_{est} [GeV];Resolution^{2} (E_{mc}/E_{est}-1)^{2};");
    262262
    263263    hres3.AddHist("MHillas.fSize", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
    264264    hres3.InitName("ResSize;Size;");
    265     hres3.InitTitle(";S [phe];Resolution (E_{mc}/E_{est}-1)^{2};");
     265    hres3.InitTitle(";S [phe];Resolution^{2} (E_{mc}/E_{est}-1)^{2};");
    266266/*
    267267    hres3.AddHist("MMcEvt.fEnergy", "(MEnergyEst.fVal/MMcEvt.fEnergy-1)^2", MH3::kProfile);
    … …  
    271271    hres3.AddHist("MPointingPos.fZd", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
    272272    hres3.InitName("ResTheta;Theta;");
    273     hres3.InitTitle(";\\Theta [\\circ];Resolution (E_{mc}/E_{est}-1)^{2};");
     273    hres3.InitTitle(";\\Theta [\\circ];Resolution^{2} (E_{mc}/E_{est}-1)^{2};");
     274
    274275    hres3.AddHist("MMcEvt.fImpact/100", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
    275276    hres3.InitName("ResImpact;Impact;");
    276     hres3.InitTitle(";I [m];Resolution (E_{mc}/E_{est}-1)^{2};");
     277    hres3.InitTitle(";I [m];Resolution^{2} (E_{mc}/E_{est}-1)^{2};");
    277278
    278279    MFillH fillh0(&hist);

  • trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc

    r9518 r9565  
    116116// ----------------
    117117//
    118 //  What we know so far about (maybe) important changes in cosy:
     118//  What we know so far about (maybe) important changes in cosy or to the
     119//  hardware:
     120//
     121//   25.01.2010: Changed scaling factors for TPoint camera (M1)
     122//
     123//   11.01.2010: New TPoint camera (M1)
    119124//
    120125//   18.03.2006: The camera holding has been repaired and the camera got
    … …  
    123128//   16.04.2006: Around this date a roque lamp focussing hass been done
    124129//
    125 //   25.04.2006: A missalignment intrduced with the roque adjust has been
     130//   25.04.2006: A missalignment introduced with the roque adjust has been
    126131//               corrected
    127132//
    … …  
    164169//   14. May  2009                // M1/M2 after upgrade (from TPoints taken in the days before)
    165170//   17. Aug. 2009              New pointing models for both telescopes
    166 //
     171//   01. Feb. 2010              New pointing models for both telescopes
     172//   26. Feb. 2010              New pointing models for both telescopes
    167173//
    168174// From 2.2.2006 beginnig of the night (21:05h, run >=81855) to 24.2.2006

  • trunk/MagicSoft/Mars/msimreflector/MMirror.cc

    r9371 r9565  
    9999using namespace std;
    100100
     101void MMirror::SetShape(Char_t c)
     102{
     103    switch (toupper(c))
     104    {
     105    case 'S':
     106        fShape = 0;
     107        break;
     108
     109    case 'P':
     110        fShape = 1;
     111        break;
     112
     113    default:
     114        fShape = 0;
     115    }
     116}
     117
    101118// --------------------------------------------------------------------------
    102119//

  • trunk/MagicSoft/Mars/msimreflector/MMirror.h

    r9564 r9565  
    2525    Double_t  fSigmaPSF;
    2626
     27    Int_t fShape; // Spherical=0, Parabolic=1
     28
    2729    //    MMirror *fNeighbors[964];
    2830
    2931public:
    30     MMirror() : fSigmaPSF(-1)
     32    MMirror() : fSigmaPSF(-1), fShape(0)
    3133    {
    3234    }
    … …  
    4850        fTilt.Rotate(-n.Theta(), TVector3(-n.Y(), n.X(), 0));
    4951    }
     52    void SetShape(Char_t c);
    5053
    5154    void SetZ(Double_t z) { fPos.SetZ(z); }
    … …  
    9598    }
    9699    */
    97     ClassDef(MMirror, 1) // Base class to describe a mirror
     100    ClassDef(MMirror, 2) // Base class to describe a mirror
    98101};
    99102

  • trunk/MagicSoft/Mars/msimreflector/MReflector.cc

    r9518 r9565  
    162162                        atof(arr[5]->GetName()));
    163163
    164     const Double_t F = atof(arr[6]->GetName());
    165 
    166     const TString val = arr[7]->GetName();
     164    UInt_t n = 6;
     165
     166    TString six = arr[n]->GetName();
     167
     168    Char_t shape = 'S';
     169    if (!six.IsFloat())
     170    {
     171        shape = six[0];
     172        n++;
     173    }
     174
     175    const Double_t F = atof(arr[n++]->GetName());
     176
     177    const TString val = arr[n++]->GetName();
    167178
    168179    const Double_t psf = val.IsFloat() ? val.Atof() : -1;
    169180
    170     const UInt_t n = val.IsFloat() ? 9 : 8;
     181    n += val.IsFloat() ? 1 : 0;
    171182
    172183    TString type = arr[n-1]->GetName();
    … …  
    201212    }
    202213
    203     m->SetFocalLength(F);
    204214    m->SetPosition(pos);
    205215    m->SetNorm(norm);
     216    m->SetShape(shape);
     217    m->SetFocalLength(F);
    206218    m->SetSigmaPSF(psf>=0 ? psf : defpsf);
    207219
    … …  
    223235//     x y z nx ny nz F [psf] Type ...
    224236//
    225 //  x:      x-coordinate of the mirror's center
    226 //  y:      y-coordinate of the mirror's center
    227 //  z:      z-coordinate of the mirror's center
    228 //  nx:     x-component of the normal vecor in the mirror center
    229 //  ny:     y-component of the normal vecor in the mirror center
    230 //  nz:     z-component of the normal vecor in the mirror center
    231 //  F:      Focal distance of a spherical mirror
    232 //  [psf]:  This number is the psf given in the units of x,y,z and
    233 //          defined at the focal distance F. It can be used to overwrite
    234 //          the second argument given in ReadFile for individual mirrors.
    235 //  Type:   A instance of a mirrot of the class Type MMirrorType is created
    236 //          (Type can be, for example, Hex for for MMirrorHex).
    237 //  ...:    Additional arguments as defined in MMirrorType::ReadM
     237//  x:         x-coordinate of the mirror's center
     238//  y:         y-coordinate of the mirror's center
     239//  z:         z-coordinate of the mirror's center
     240//  nx:        x-component of the normal vecor in the mirror center
     241//  ny:        y-component of the normal vecor in the mirror center
     242//  nz:        z-component of the normal vecor in the mirror center
     243//  [shape]:   S for spherical <default>, P for parabolic
     244//  F:         Focal distance of a spherical mirror
     245//  [psf]:     This number is the psf given in the units of x,y,z and
     246//             defined at the focal distance F. It can be used to overwrite
     247//             the second argument given in ReadFile for individual mirrors.
     248//  Type:      A instance of a mirrot of the class Type MMirrorType is created
     249//             (Type can be, for example, Hex for for MMirrorHex).
     250//  ...:       Additional arguments as defined in MMirrorType::ReadM
    238251//
    239252//

  • trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc

    r9564 r9565  
    1818!   Author(s): Thomas Bretz,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
    1919!
    20 !   Copyright: CheObs Software Development, 2000-2009
     20!   Copyright: CheObs Software Development, 2000-2010
    2121!
    2222!
    … …  
    3333//  obsolete except they can later be "moved" inside the detector, e.g.
    3434//  if you use MSimPSF to emulate a PSF by moving photons randomly
    35 //  on the focal plane.
     35//  on the focal plane. To switch off this check set detector margin to -1.
    3636//
    3737//////////////////////////////////////////////////////////////////////////////
    … …  
    206206//    Focal length: F=R/2                        |  Focal length: F = 1/4p
    207207//                                               |
    208 //   r^2 + (z-2*F)^2 = (2*F)^2                   |            z = F/4*r^2
     208//   r^2 + (z-2*F)^2 = (2*F)^2                   |           z = r^2/4F
    209209//                                               |
    210210//          z        = -sqrt(4*F*F - r*r) + 2*F  |
    211211//          z-2*F    = -sqrt(4*F*F - r*r)        |
    212212//         (z-2*F)^2 = 4*F*F - r*r               |
    213 //   z^2-4*F*z+4*F^2 = 4*F*F - r*r  (4F^2-r^2>0) |  z - F/4*r^2 = 0
     213//   z^2-4*F*z+4*F^2 = 4*F*F - r*r  (4F^2-r^2>0) |  z - r^2/4F = 0
    214214//   z^2-4*F*z+r^2   = 0
    215215//
    … …  
    223223// ------------------
    224224//
    225 //   z^2*(1+|v|^2) - 2*z*(2*F+q*v) + |q|^2 =  0  |  z^2*|v|^2 - 2*(2/F+q*v)*z + |q|^2 = 0
     225//   z^2*(1+|v|^2) - 2*z*(2*F+q*v) + |q|^2 =  0  |  z^2*|v|^2 - 2*z*(2*F+q*v) + |q|^2 = 0
    226226//
    227227//                                   z = (-b +- sqrt(b*b - 4ac))/(2*a)
    228228//
    229229//              a =   1+|v|^2                    |             a = |v|^2
    230 //              b = - 2*(2*F+q*v)                |             b = - 2*(2/F+q*v)
     230//              b = - 2*b'  with  b' = 2*F+q*v   |             b = - 2*b'  with  b' = 2*F+q*v
    231231//              c =   |q|^2                      |             c = |q|^2
    232232//                                               |
    … …  
    234234//                                        substitute b := 2*b'
    235235//
    236 //                            z = (-2*b' +- 2*sqrt(b'*b' - ac))/(2*a)
    237 //                            z = (-  b' +-   sqrt(b'*b' - ac))/a
    238 //                            z = (-b'/a +-   sqrt(b'*b' - ac))/a
     236//                            z = (2*b' +- 2*sqrt(b'*b' - ac))/(2*a)
     237//                            z = (  b' +-   sqrt(b'*b' - ac))/a
     238//                            z = (b'/a +-   sqrt(b'*b' - ac))/a
    239239//
    240240//                                        substitute f := b'/a
    241241//
    242 //                                      z = (-f +- sqrt(f^2 - c/a)
     242//                                      z = f +- sqrt(f^2 - c/a)
    243243//
    244244// =======================================================================================
    … …  
    265265
    266266    // Find the incident point of the vector to the mirror
    267     // u corresponds to downwaqrd going particles, thus we use -u here
     267    // u corresponds to downward going particles, thus we use -u here
    268268    const Double_t b = G - q*v;
    269269    const Double_t a = v.Mod2();
    … …  
    271271
    272272    // Solution for q spherical (a+1) (parabolic mirror (a) instead of (a+1))
    273     const Double_t f = b/(a+1);
    274     const Double_t g = c/(a+1);
     273    const Double_t A = fShape ? a : a+1;
     274
     275    const Double_t f = b/A;
     276    const Double_t g = c/A;
    275277
    276278    // Solution of second order polynomial (transformed: a>0)
    277279    // (The second solution can be omitted, it is the intersection
    278280    //  with the upper part of the sphere)
    279     //    const Double_t dz = a==0 ? c/(2*b) : f - TMath::Sqrt(f*f - g);
    280     const Double_t z = f - TMath::Sqrt(f*f - g);
     281    const Double_t z = a==0 ? c/(2*b) : f - TMath::Sqrt(f*f - g);
    281282
    282283    // Move the photon along its trajectory to the x/y plane of the
    … …  
    293294
    294295    // Get normal vector for reflection by calculating the derivatives
    295     // of a spherical mirror along x and y
    296     const Double_t d = TMath::Sqrt(G*G - p.R2());
    297 
    298     // This is a normal vector at the incident point
     296    // of a the mirror's surface along x and y
     297    const Double_t d = fShape ? G : TMath::Sqrt(G*G - p.R2());
     298
     299    // The solution for the normal vector is
     300    //    TVector3 n(-p.X()/d, -p.Y()/d, 1));
     301    // Since the normal vector doesn't need to be of normal
     302    // length we can avoid an obsolete division
    299303    TVector3 n(p.X(), p.Y(), -d);
    300     // This is the obvious solution for the normal vector
    301     //  TVector3 n(-p.X()/d, -p.Y()/d, 1));
    302304
    303305    if (fSigmaPSF>0)
    … …  
    540542        // It is detector specific not reflector specific
    541543        // Discard all photons which definitly can not hit the detector surface
    542         if (!fGeomCam->HitDetector(p, fDetectorMargin))
     544        if (fDetectorMargin>=0 && !fGeomCam->HitDetector(p, fDetectorMargin))
    543545            continue;
    544546

  • trunk/MagicSoft/Mars/msimreflector/MSimReflector.h

    r9564 r9565  
    4949    void SetNameReflector(const char *name="MReflector") { fNameReflector = name; }
    5050
     51    void SetDetectorMargin(Double_t m=0) { fDetectorMargin = m; }
     52
    5153    ClassDef(MSimReflector, 0) // Task to calculate reflection on a mirror
    5254};
Note: See TracChangeset for help on using the changeset viewer.