Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 9564) +++ trunk/MagicSoft/Mars/Changelog (revision 9565) @@ -40,4 +40,32 @@ * msimreflector/MSimReflector.[h,cc]: - added name of reflector as data member + + * melectronics/MAvalanchePhotoDiode.cc: + - scale the crosstalk probability as the height of the emitted + signal with the recovery time. + + * mhflux/MHEnergyEst.cc: + - added a workaround to get rid of some root-warnings + + * mjtrain/MJTrainEnergy.cc: + - improved axis labels + + * mpointing/MPointingDevCalc.cc: + - added some more comments + + * msimreflector/MMirror.[h,cc]: + - added fShape to allow for parabolic mirrors + + * msimreflector/MReflector.cc: + - implemented the option of parabolic mirrors into the + reflector defintion file + + * msimreflector/MSimReflector.[h,cc]: + - implemented (and fixed) the calculation of the reflection at + parabolic surfaces + - implemented the possibility to switch off the early check for + "photon can hit the camera" (fDetectorMargin<0) + - added setter for fDetectorMargin + Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 9564) +++ trunk/MagicSoft/Mars/NEWS (revision 9565) @@ -7,4 +7,9 @@ * On some systems (version of make?) the __LINUX__ definition was not set when compiling. This lead to MTime(-1) not working properly. Fixed. + + ;ceres: + + * Allow for individual mirrors with parabolic shape (for details see + MReflector::ReadFile) Index: trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc =================================================================== --- trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc (revision 9564) +++ trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc (revision 9565) @@ -123,10 +123,13 @@ // return 0; + // Number of the x/y cell in the one dimensional array // const Int_t cell = fHist.GetBin(x, y); const Int_t cell = x + (fHist.GetNbinsX()+2)*y; - // This is the fastes way to access the bin-contents in fArray + // Getting a reference to the float is the fastes way to + // access the bin-contents in fArray Float_t &cont = fHist.GetArray()[cell]; + // Calculate the time since the last breakdown // const Double_t dt = t-fHist.GetBinContent(x, y)-fDeadTime; // const Float_t dt = t-cont-fDeadTime; @@ -136,7 +139,13 @@ return 0; - // Signal height (in units of one photon) produced after dead time - const Float_t weight = fRecoveryTime<=0 ? 1 : 1.-exp(-dt/fRecoveryTime); - + // The signal height (in units of one photon) produced after dead time + // depends on the recovery of the cell - described by an exponential. + const Float_t weight = fRecoveryTime<=0 ? 1. : 1-TMath::Exp(-dt/fRecoveryTime); + + // The probability that the cell emits a photon causing crosstalk + // scales as the signal height. + const Float_t prob = weight*fCrosstalkProb; + + // Set the contents to the time of the last breakdown (now) cont = t; // fHist.SetBinContent(x, y, t) @@ -161,14 +170,15 @@ */ + //for (int i=0; i<1; i++) while (1) { const Double_t rndm = gRandom->Rndm(); - if (rndm>=fCrosstalkProb) + if (rndm>=prob/*fCrosstalkProb*/) break; - // We can re-use the random number becuase it is uniformely + // We can re-use the random number because it is uniformely // distributed. This saves cpu power - const Int_t dir = TMath::FloorNint(4*rndm/fCrosstalkProb); + const Int_t dir = TMath::FloorNint(4*rndm/prob/*fCrosstalkProb*/); // Get a random neighbor which is hit. Index: trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (revision 9564) +++ trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (revision 9565) @@ -298,9 +298,15 @@ pad->GetPad(1)->cd(1); - if (gPad->FindObject("EnergyEst_ex")) + if ((hx = dynamic_cast(gPad->FindObject("EnergyEst_ex")))) + { + hx->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated hx = fHEnergy.Project3D("ex"); - - if (gPad->FindObject("EnergyEst_ey")) + } + + if ((hy = dynamic_cast(gPad->FindObject("EnergyEst_ey")))) + { + hy->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated hy = fHEnergy.Project3D("ey"); + } if (hx && hy) @@ -317,6 +323,9 @@ { pad->GetPad(1)->GetPad(2)->cd(1); - if (gPad->FindObject("EnergyEst_ez")) + if ((hx = dynamic_cast(gPad->FindObject("EnergyEst_ez")))) + { + hx->GetSumw2()->Set(0); // Workaround. Otherwise we get a warning because it is recreated fHEnergy.Project3D("ez"); + } pad->GetPad(1)->GetPad(2)->cd(2); Index: trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc =================================================================== --- trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc (revision 9564) +++ trunk/MagicSoft/Mars/mjtrain/MJTrainEnergy.cc (revision 9565) @@ -256,12 +256,12 @@ hres2.SetDrawOption("colz profx"); - MHn hres3("Resolution", "Energy Resolution"); + MHn hres3("Resolution", "Energy Resolution squared"); hres3.AddHist("MEnergyEst.fVal", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile); hres3.InitName("ResEest;EnergyEst;"); - hres3.InitTitle(";E_{est} [GeV];Resolution (E_{mc}/E_{est}-1)^{2};"); + hres3.InitTitle(";E_{est} [GeV];Resolution^{2} (E_{mc}/E_{est}-1)^{2};"); hres3.AddHist("MHillas.fSize", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile); hres3.InitName("ResSize;Size;"); - hres3.InitTitle(";S [phe];Resolution (E_{mc}/E_{est}-1)^{2};"); + hres3.InitTitle(";S [phe];Resolution^{2} (E_{mc}/E_{est}-1)^{2};"); /* hres3.AddHist("MMcEvt.fEnergy", "(MEnergyEst.fVal/MMcEvt.fEnergy-1)^2", MH3::kProfile); @@ -271,8 +271,9 @@ hres3.AddHist("MPointingPos.fZd", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile); hres3.InitName("ResTheta;Theta;"); - hres3.InitTitle(";\\Theta [\\circ];Resolution (E_{mc}/E_{est}-1)^{2};"); + hres3.InitTitle(";\\Theta [\\circ];Resolution^{2} (E_{mc}/E_{est}-1)^{2};"); + hres3.AddHist("MMcEvt.fImpact/100", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile); hres3.InitName("ResImpact;Impact;"); - hres3.InitTitle(";I [m];Resolution (E_{mc}/E_{est}-1)^{2};"); + hres3.InitTitle(";I [m];Resolution^{2} (E_{mc}/E_{est}-1)^{2};"); MFillH fillh0(&hist); Index: trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc (revision 9564) +++ trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc (revision 9565) @@ -116,5 +116,10 @@ // ---------------- // -// What we know so far about (maybe) important changes in cosy: +// What we know so far about (maybe) important changes in cosy or to the +// hardware: +// +// 25.01.2010: Changed scaling factors for TPoint camera (M1) +// +// 11.01.2010: New TPoint camera (M1) // // 18.03.2006: The camera holding has been repaired and the camera got @@ -123,5 +128,5 @@ // 16.04.2006: Around this date a roque lamp focussing hass been done // -// 25.04.2006: A missalignment intrduced with the roque adjust has been +// 25.04.2006: A missalignment introduced with the roque adjust has been // corrected // @@ -164,5 +169,6 @@ // 14. May 2009 // M1/M2 after upgrade (from TPoints taken in the days before) // 17. Aug. 2009 New pointing models for both telescopes -// +// 01. Feb. 2010 New pointing models for both telescopes +// 26. Feb. 2010 New pointing models for both telescopes // // From 2.2.2006 beginnig of the night (21:05h, run >=81855) to 24.2.2006 Index: trunk/MagicSoft/Mars/msimreflector/MMirror.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirror.cc (revision 9564) +++ trunk/MagicSoft/Mars/msimreflector/MMirror.cc (revision 9565) @@ -99,4 +99,21 @@ using namespace std; +void MMirror::SetShape(Char_t c) +{ + switch (toupper(c)) + { + case 'S': + fShape = 0; + break; + + case 'P': + fShape = 1; + break; + + default: + fShape = 0; + } +} + // -------------------------------------------------------------------------- // Index: trunk/MagicSoft/Mars/msimreflector/MMirror.h =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirror.h (revision 9564) +++ trunk/MagicSoft/Mars/msimreflector/MMirror.h (revision 9565) @@ -25,8 +25,10 @@ Double_t fSigmaPSF; + Int_t fShape; // Spherical=0, Parabolic=1 + // MMirror *fNeighbors[964]; public: - MMirror() : fSigmaPSF(-1) + MMirror() : fSigmaPSF(-1), fShape(0) { } @@ -48,4 +50,5 @@ fTilt.Rotate(-n.Theta(), TVector3(-n.Y(), n.X(), 0)); } + void SetShape(Char_t c); void SetZ(Double_t z) { fPos.SetZ(z); } @@ -95,5 +98,5 @@ } */ - ClassDef(MMirror, 1) // Base class to describe a mirror + ClassDef(MMirror, 2) // Base class to describe a mirror }; Index: trunk/MagicSoft/Mars/msimreflector/MReflector.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MReflector.cc (revision 9564) +++ trunk/MagicSoft/Mars/msimreflector/MReflector.cc (revision 9565) @@ -162,11 +162,22 @@ atof(arr[5]->GetName())); - const Double_t F = atof(arr[6]->GetName()); - - const TString val = arr[7]->GetName(); + UInt_t n = 6; + + TString six = arr[n]->GetName(); + + Char_t shape = 'S'; + if (!six.IsFloat()) + { + shape = six[0]; + n++; + } + + const Double_t F = atof(arr[n++]->GetName()); + + const TString val = arr[n++]->GetName(); const Double_t psf = val.IsFloat() ? val.Atof() : -1; - const UInt_t n = val.IsFloat() ? 9 : 8; + n += val.IsFloat() ? 1 : 0; TString type = arr[n-1]->GetName(); @@ -201,7 +212,8 @@ } - m->SetFocalLength(F); m->SetPosition(pos); m->SetNorm(norm); + m->SetShape(shape); + m->SetFocalLength(F); m->SetSigmaPSF(psf>=0 ? psf : defpsf); @@ -223,17 +235,18 @@ // x y z nx ny nz F [psf] Type ... // -// x: x-coordinate of the mirror's center -// y: y-coordinate of the mirror's center -// z: z-coordinate of the mirror's center -// nx: x-component of the normal vecor in the mirror center -// ny: y-component of the normal vecor in the mirror center -// nz: z-component of the normal vecor in the mirror center -// F: Focal distance of a spherical mirror -// [psf]: This number is the psf given in the units of x,y,z and -// defined at the focal distance F. It can be used to overwrite -// the second argument given in ReadFile for individual mirrors. -// Type: A instance of a mirrot of the class Type MMirrorType is created -// (Type can be, for example, Hex for for MMirrorHex). -// ...: Additional arguments as defined in MMirrorType::ReadM +// x: x-coordinate of the mirror's center +// y: y-coordinate of the mirror's center +// z: z-coordinate of the mirror's center +// nx: x-component of the normal vecor in the mirror center +// ny: y-component of the normal vecor in the mirror center +// nz: z-component of the normal vecor in the mirror center +// [shape]: S for spherical , P for parabolic +// F: Focal distance of a spherical mirror +// [psf]: This number is the psf given in the units of x,y,z and +// defined at the focal distance F. It can be used to overwrite +// the second argument given in ReadFile for individual mirrors. +// Type: A instance of a mirrot of the class Type MMirrorType is created +// (Type can be, for example, Hex for for MMirrorHex). +// ...: Additional arguments as defined in MMirrorType::ReadM // // Index: trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc (revision 9564) +++ trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc (revision 9565) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 1/2009 ! -! Copyright: CheObs Software Development, 2000-2009 +! Copyright: CheObs Software Development, 2000-2010 ! ! @@ -33,5 +33,5 @@ // obsolete except they can later be "moved" inside the detector, e.g. // if you use MSimPSF to emulate a PSF by moving photons randomly -// on the focal plane. +// on the focal plane. To switch off this check set detector margin to -1. // ////////////////////////////////////////////////////////////////////////////// @@ -206,10 +206,10 @@ // Focal length: F=R/2 | Focal length: F = 1/4p // | -// r^2 + (z-2*F)^2 = (2*F)^2 | z = F/4*r^2 +// r^2 + (z-2*F)^2 = (2*F)^2 | z = r^2/4F // | // z = -sqrt(4*F*F - r*r) + 2*F | // z-2*F = -sqrt(4*F*F - r*r) | // (z-2*F)^2 = 4*F*F - r*r | -// z^2-4*F*z+4*F^2 = 4*F*F - r*r (4F^2-r^2>0) | z - F/4*r^2 = 0 +// z^2-4*F*z+4*F^2 = 4*F*F - r*r (4F^2-r^2>0) | z - r^2/4F = 0 // z^2-4*F*z+r^2 = 0 // @@ -223,10 +223,10 @@ // ------------------ // -// 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 +// 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 // // z = (-b +- sqrt(b*b - 4ac))/(2*a) // // a = 1+|v|^2 | a = |v|^2 -// b = - 2*(2*F+q*v) | b = - 2*(2/F+q*v) +// b = - 2*b' with b' = 2*F+q*v | b = - 2*b' with b' = 2*F+q*v // c = |q|^2 | c = |q|^2 // | @@ -234,11 +234,11 @@ // substitute b := 2*b' // -// z = (-2*b' +- 2*sqrt(b'*b' - ac))/(2*a) -// z = (- b' +- sqrt(b'*b' - ac))/a -// z = (-b'/a +- sqrt(b'*b' - ac))/a +// z = (2*b' +- 2*sqrt(b'*b' - ac))/(2*a) +// z = ( b' +- sqrt(b'*b' - ac))/a +// z = (b'/a +- sqrt(b'*b' - ac))/a // // substitute f := b'/a // -// z = (-f +- sqrt(f^2 - c/a) +// z = f +- sqrt(f^2 - c/a) // // ======================================================================================= @@ -265,5 +265,5 @@ // Find the incident point of the vector to the mirror - // u corresponds to downwaqrd going particles, thus we use -u here + // u corresponds to downward going particles, thus we use -u here const Double_t b = G - q*v; const Double_t a = v.Mod2(); @@ -271,12 +271,13 @@ // Solution for q spherical (a+1) (parabolic mirror (a) instead of (a+1)) - const Double_t f = b/(a+1); - const Double_t g = c/(a+1); + const Double_t A = fShape ? a : a+1; + + const Double_t f = b/A; + const Double_t g = c/A; // Solution of second order polynomial (transformed: a>0) // (The second solution can be omitted, it is the intersection // with the upper part of the sphere) - // const Double_t dz = a==0 ? c/(2*b) : f - TMath::Sqrt(f*f - g); - const Double_t z = f - TMath::Sqrt(f*f - g); + const Double_t z = a==0 ? c/(2*b) : f - TMath::Sqrt(f*f - g); // Move the photon along its trajectory to the x/y plane of the @@ -293,11 +294,12 @@ // Get normal vector for reflection by calculating the derivatives - // of a spherical mirror along x and y - const Double_t d = TMath::Sqrt(G*G - p.R2()); - - // This is a normal vector at the incident point + // of a the mirror's surface along x and y + const Double_t d = fShape ? G : TMath::Sqrt(G*G - p.R2()); + + // The solution for the normal vector is + // TVector3 n(-p.X()/d, -p.Y()/d, 1)); + // Since the normal vector doesn't need to be of normal + // length we can avoid an obsolete division TVector3 n(p.X(), p.Y(), -d); - // This is the obvious solution for the normal vector - // TVector3 n(-p.X()/d, -p.Y()/d, 1)); if (fSigmaPSF>0) @@ -540,5 +542,5 @@ // It is detector specific not reflector specific // Discard all photons which definitly can not hit the detector surface - if (!fGeomCam->HitDetector(p, fDetectorMargin)) + if (fDetectorMargin>=0 && !fGeomCam->HitDetector(p, fDetectorMargin)) continue; Index: trunk/MagicSoft/Mars/msimreflector/MSimReflector.h =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MSimReflector.h (revision 9564) +++ trunk/MagicSoft/Mars/msimreflector/MSimReflector.h (revision 9565) @@ -49,4 +49,6 @@ void SetNameReflector(const char *name="MReflector") { fNameReflector = name; } + void SetDetectorMargin(Double_t m=0) { fDetectorMargin = m; } + ClassDef(MSimReflector, 0) // Task to calculate reflection on a mirror };