Index: old/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- old/WuerzburgSoft/Thomas/mphys/Changelog (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/Changelog (revision 15424) @@ -0,0 +1,323 @@ + -*-*- END -*-*- + 2002/08/19: Thomas Bretz + + * MCascade.[h,cc]: + - removed argument of DoInvCompton + - delete GRandom only if set + - set gRandom to 0 after deletion + + * MElectron.[h,cc]: + - added Sigma_ge + - removed argument from DoInvCompton + - Implemented Omega_sigmae and RandomThetaE + - Changed old wrong dertermination of interaction angle to + correct simulation. + - Reengenieered some formulas to be more correct in terms of numerics + + + + 2002/08/02: Thomas Bretz + + * MCascade.[h,cc]: + - doesn't overwrite existing files anymore + - removed redefinition of default argument + + * MFit.cc: + - removed redefinition of default argument + + * analp.C: + - added smearing + + + + 2002/07/29: Thomas Bretz + + * MCascade.[h,cc]: + - replaces draw by fIsBatch, using gROOT->IsBatch in addition. + + * MFit.h: + - added default=0 to SetParameter for min and max + + + + 2002/07/29: Thomas Bretz + + * MElectron.[h,cc]: + - if the photon energy is to small return a straight line in + Compton + - changed lolimit of p_e integration -0.5 + - changed number of calculated points to 50 for fP, fQ + - changed usage of TRandom to gRandom + - changed some TVector and TMatrix objects to static + - added new copy constructors + + * MPairProduction.cc: + - use the new copy constructor of MElectron + + * MParticle.[h,cc]: + - added new copy constructors + - Added a TAxis object in Planck + - replaced usage of TRandom by gRandom + - replaced TVector, TMatrix by TVector3, TRotation + - fixed a bug in SetNewPosition (fX was not incremented if the + particle reached the observers plain) + - incremented version-number + + * MPhoton.[h,cc]: + - added Fill member function + - added new copy constructors + + * MFit.[h,cc], MCascade.[h,cc]: + - added + + * phys.C: + - moved code to MCascade + + * Makefile, PhysLinkDef.h: + - added MFit, MCascade + + * analp.C: + - simplified + - added fit for cutoff + + + + 2002/07/26: Thomas Bretz + + * MParticle.[h,cc]: + - added fX data member + - corrected a small bug in the SetNewPosition calculation + (dr*cos(fTheta)>R*cos(fTheta) --> x(2)>R*cos(fTheta)) + - added usage of fX to the =-operator and SetNewPosition + - added fSrcR (unsaved) + + * MElectron.cc: + - added usage of fX to the =-operator and SetNewPosition + - changed to use the more accurate SetNewPosition in case B==0 + + * MPhoton.cc: + - return 1e50 as InteractionLength in case Eg<100GeV + + * phys.C: + - added support for a B-Field-Bubble + + + + 2002/07/21: Thomas Bretz + + * MElectron.cc: + - changed upper limit in InteractionLength from 1e-11 to 2e-11 + - changed upper limit of p_e in EnergyLoss from -10.8 to -10.6 + + * MParticle.[h,cc]: + - replaced Planck by tanjas background + + * background.txt: + - added + + * MPhoton.[h,cc]: + - changed lower limit of Int2-Integration for Interaction Length + in case Eg<5e4 to fit better its needs + - changed max/min of drawing interaction length + - changed lower limit of interaction length plot + + * analp.C: + - added plotting 1h1426 data + - changed alpha/plot from -2/2 to -1/1 (which means I do the well + known E^2 plot of an E^-2 spectrum) + - added a cutoff to simulate cutoff spectra + - some small corrections for the layout of the plots + + * phys.C: + - added a maximum number of inv. Compton processes + - don't throw away photons which were reemitted lower than the lower + limit of the plot (so they are still written to the file) + - changed the lower limit of the PhotonSpect simulation function + to a physically more realistic value. No photon which don't have + an corresponding angle are produced anymore. + - stop the electron only if its energy is less than a thousand + of its primary energy + + + + 2002/06/21: Thomas Bretz + + * analp.C: + - clean up the code + - added the point-spread histograms from different directions + + * mphys.C: + - removed the angular histograms + - added the plot spectral index + - added the timer to get a usable ubdated window + + + + 2002/06/20: Thomas Bretz + + * MElectron.cc: + - Some changes to the use of beta (simplifications) + - changed the default of rho (B==0) to 0 + + * MParticle.cc: + - made TRandom static (InitRandom) + + * mphys.C: + - added support for Emis and Esum + + + + 2002/06/17: Thomas Bretz + + * MElectron.[h,cc]: + - Added SetNewPositionB + - changed the constructor to support positrons + + * MPairProduction.cc: + - use new constructor of MElectron + + * MParticle.[h,cc]: + - added InitRandom + - added GetType + + * anale.C: + - added SetDirectory(NULL) + + * analp.C: + - changed to an automatic scaling + + * phys.C: + - changed limit to 60min + - added lower limit for electrons 5*lo + - Changed to use a B-field + + + + 2002/06/17: Thomas Bretz + + * MElectron.cc: + - Added SetNewDirection for the electron itself + + * MParticle.cc: + - changed sqrt() to asin(sqrt()) - no effect (but correct) + + + + 2002/06/14: Thomas Bretz + + * MElectron.[h,cc]: + - Moved Planck function to MParticle + - Added the DiSum function + - changed Li2 to use the DiSum function (speed reasons) + - changed the eps of all integrals to 1e-2 + - changed the p_e function to use the Compton function + - made all function which are used for integration only static + and changed the function definition from=0, to=0 + + * MPhoton.[h,cc]: + - removed the planck function + - changed the eps of all integrals to 1e-2 + + * MParticle.[h,cc]: + - added the planck function + + * Makefile: + - removed unused source files + + * analp.C: + - did some small fixes + - changed the sclae of the radius- and phi-plots + + * phys.C: + - changed the integral eps to 1e-2 + - small changed to the output + - removed integral==0, replaced by Ep==pow(10,0) + + + + + 2002/06/13: Thomas Bretz + + * MParticle.h: + - Implemented Set/IsPrimary + - removed pure abstract method (root can't read this) + - added default constructor + + * MElectron.cc: + - changed many variables to const + - changed the output + + * MPairProduction.cc: + - changed many variables to const + + * MPhoton.[h,cc]: + - changed the upper and lower limit of the integration in Int2 + - added upper and lower limit of draw to DrawInteractionLength + + * energyloss.C: + - changed output range for ZofR + + * phys.C: + - changed energy randomizer + - added root file output + - removed log scale from point spread plots + + * analp.C: + - added macro to analize the photons from the root file + + * anale.C: + - added macro to analize the electrons from the root file + + + + 2002/06/12: Thomas Bretz + + * MElectron.[h,cc]: + - added a primitive theta dependancy to DoInvCompton + - added DrawInteractionLength + - changed algorithm in ZofR and RofZ + + * MParticle.[h,cc]: + - added RofZ and ZofR as static functions + - changed from MParContainer to TObject + - removed unnecessary or commented inline functions + + * MPhoton.[h,cc]: + - added DrawInteractionLength + - fixed the bug causing the 'strange factor': sigma_gg needs + the sqrt of s. + + * energyloss.C: + - added + + * MPairProduction.[h,cc]: + - simplified formula + - changed arguments to Energy of Photon and interaction angle + + * phys.C: + - added correct photon energy for pair production + - added correct angle for pair production + - changed output to E^2*Counts + - changed gamma production from E^-2 to Uniform by weighting + - removed all unecessary functions and code (s.energyloss.C) + + + + 2002/06/10: Thomas Bretz + + * MElectron.[h,cc], MPhoton.[h,cc], MParticle.[h,cc]: + - implemented SetNewDirection, SetNewPosition + - implemented DoInvCompton + - =operator + + * MPairProduction.cc: + - changed to use SetNewDirection Method + + * phys.C: + - changed to new style functions + + * MGeinIRPhoton.cc, MGenPrimaryParticle.[h,cc]: + - make compile (not used) + - changed Process to get angle + + Index: old/WuerzburgSoft/Thomas/mphys/MCascade.cc =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MCascade.cc (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MCascade.cc (revision 15424) @@ -0,0 +1,491 @@ +#include "MCascade.h" + +#include // fabs, for alpha +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "MPhoton.h" +#include "MElectron.h" + +#include "MH.h" +#include "MBinning.h" + +ClassImp(MCascade); + +Double_t PrimSpect(Double_t *x, Double_t *k) +{ + return pow(pow(10, x[0]), k[0]); +} + +Double_t PhotonSpect(Double_t *x, Double_t *k=NULL) +{ + Double_t Ep = pow(10, x[0]); + + Double_t res = MPhoton::Int2(&Ep, k); + return res*1e55; //65/k[0]; + + //return MPhoton::Planck(&Ep, &k[1]); +} + +Double_t Sbar_sigmas(Double_t *x, Double_t *k) +{ + Double_t sbar = pow(10, x[0]); + + Double_t s = 1./(sbar*4); + + Double_t sigma = MPhoton::Sigma_gg(&s); + + return sigma*sbar*1e28; +} + +Double_t RandomThetaG(Double_t Eg, Double_t Ep) +{ + Double_t E0 = 511e-6; // [GeV] + + Double_t f = Eg/E0*Ep/E0; + + if (f<1) + return 0; + + static TF1 func("RndThetaG", Sbar_sigmas, 0, 0, 0); + + func.SetRange(0, log10(f)); + func.SetNpx(50); + + Double_t sbar = pow(10, func.GetRandom()); + Double_t theta = acos(1.-sbar*2/f); + + return theta; +} + +Double_t MCascade::GetEnergy() +{ + static int bin=0; + Double_t w = log10(fEHi/fELo)/fNumBins; + + Double_t E = fELo*pow(10, gRandom->Uniform(w) + w*(fNumBins-bin-1)); + + if (++bin==fNumBins) + bin=0; + + return E; +} + +TFile *MCascade::OpenFile(TString fFilename) +{ + TFile *fFile = new TFile(fFilename, "CREATE", "Intergalactic cascade", 9); + + if (fFile->IsZombie()) + { + delete fFile; + return NULL; + } + + Write("Setup"); + + cout << "Trees... " << flush; + + TTree *T1 = new TTree ("Photons", "Photons from Cascade"); + TTree *T2 = new TTree ("Electrons", "Electrons in the Cascade"); + + cout << "Branches... " << flush; + + MPhoton dummyp; + void *ptr = &dummyp; + fBranchGammas = T1->Branch("MPhoton.", "MPhoton", &ptr); + + MElectron dummye; + ptr = &dummye; + fBranchElectrons = T2->Branch("MElectron.", "MElectron", &ptr); + + return fFile; +} + +void MCascade::CloseFile(TFile *fFile) +{ + fFile->Write(); + cout << "Wrote: " << fFile->GetName() << endl; + delete fFile; +} + +void MCascade::ProcessElectron(MElectron &e, TList &fListGammas) +{ + Double_t Ee = e.GetEnergy(); + + cout << ":" << flush; + + int test = fNumMaxInvCompton<0 ? -1 : 0; + int n=0; + while ((test<0 ? true : (test++ fRatioInvCompton*Ee)) + { + n++; + + if (!e.SetNewPositionB(e.GetZ()>fBubbleZ ? fB : 0)) + { + cout << "!" << flush; + return; + } + + MPhoton *p = e.DoInvCompton(); + + fBranchElectrons->GetTree()->Fill(); + + //cout << "." << flush; + fListGammas.Add(p); + + if (fabs(e.GetTheta()*3437)>60) // < 60min + { + cout << "T" << flush; + return; + } + + if (e.GetEnergy()SetAddress(&p); + + cout << ":" << fListGammas.GetSize() << ":" << flush; + + TIter NextP(&fListGammas); + while ((p=(MPhoton*)NextP())) + { + if (!ProcessGamma(*p, weight, fListElectrons)) + { + delete fListGammas.Remove(p); + cout << "." << flush; + continue; + } + + fBranchGammas->GetTree()->Fill(); + Esum += p->GetEnergy(); + //cout << "!" << flush; + } + + return Esum; +} + +Double_t MCascade::ProcessElectrons(TList &fListElectrons, TList &fListGammas) +{ + Double_t E = 0; + + cout << ":" << fListElectrons.GetSize() << flush; + + TIter Next(&fListElectrons); + MElectron *e = NULL; + fBranchElectrons->SetAddress(&e); + while ((e=(MElectron*)Next())) + { + e->SetIsPrimary(kTRUE); + fBranchElectrons->GetTree()->Fill(); + e->SetIsPrimary(kFALSE); + + ProcessElectron(*e, fListGammas); + + E += e->GetEnergy(); + } + fListElectrons.Delete(); + + return E; +} + +void MCascade::ProcessPrimaryGamma(Double_t E, Double_t weight) +{ + TList fListGammas; + TList fListElectrons; + + fListGammas.SetOwner(); + fListElectrons.SetOwner(); + + MPhoton *gamma=new MPhoton(E, fSrcZ); + gamma->SetSrcR(fSrcR); + gamma->InitRandom(); + fListGammas.Add(gamma); + + gamma->SetIsPrimary(kTRUE); + fBranchGammas->SetAddress(&gamma); + fBranchGammas->GetTree()->Fill(); + gamma->SetIsPrimary(kFALSE); + + Double_t Esum=0; // sum of all energies + Double_t Emis=0; // sum of the energies thrown away + while (1) + { + if (fListGammas.GetSize()==0) + break; + + cout << " |P" << flush; + + Esum += ProcessGammas(fListGammas, fListElectrons, weight); + + if (!fIsBatch) + fListGammas.ForEach(MPhoton, Fill)(fHist, fDisplayIndex, weight); + fListGammas.Delete(); + + if (fListElectrons.GetSize()==0) + break; + + cout << " |E" << flush; + + Emis += ProcessElectrons(fListElectrons, fListGammas); + } + Esum += Emis; + + cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl; +} + +MCascade::MCascade() +{ + if (gRandom) + delete gRandom; + + TRandom r(0); + gRandom = new TRandom3(r.GetSeed()); + + fHist.SetName("Spectrum"); + fHist.SetXTitle("E [GeV]"); + fHist.SetYTitle(Form("E^{%.1f} Counts", fDisplayIndex)); + fHist.GetXaxis()->SetLabelOffset(-0.015); + fHist.GetXaxis()->SetTitleOffset(1.1); + fHist.SetFillStyle(0); + fHist.SetMarkerStyle(kPlus); +} + +MCascade::~MCascade() +{ + delete gRandom; + gRandom = 0; +} + +void MCascade::SetSourceZ(Double_t z) +{ + fSrcZ = z; + fSrcR = MParticle::RofZ(&z); +} + +void MCascade::SetSourceRZ(Double_t r) // [kpc] +{ + fSrcZ = MParticle::ZofR(&r); + fSrcR = r; +} + +void MCascade::SetEnergyBins(Int_t n, Double_t lo, Double_t hi) +{ + fNumBins = n; // number of bins produced in energy spectrum + + fELo = lo; // lower limit of displayed spectrum + fEHi = hi; // upper limit of spectrum (cutoff) +} + +void MCascade::SetBradius(Double_t r) // [Mpc] +{ + Double_t bubbler = fSrcR-1e3*r; + fBubbleZ = MParticle::ZofR(&bubbler); +} + +void MCascade::Run(TString filename, Bool_t draw) +{ + fIsBatch = gROOT->IsBatch() ? kFALSE : draw; + + // ------------------------------ + + cout << "Output File '" << filename << "'... " << flush; + + TFile *file=OpenFile(filename); + if (!file) + return; + + // ------------------------------ + + cout << endl; + + cout << "R = " << fSrcR << "kpc" << endl; + cout << "Z = " << fSrcZ << endl; + + cout << "Setting up: Histograms... " << flush; + + fHist.Reset(); + TH1D histsrc; + + MBinning bins; + bins.SetEdgesLog(fNumBins, fELo, fEHi); + + MH::SetBinning(&fHist, &bins); + MH::SetBinning(&histsrc, &bins); + + TCanvas *c=NULL; + + if (!fIsBatch) + { + fHist.SetMinimum(pow(fELo, fSpectralIndex+fDisplayIndex)/100); + histsrc.SetMinimum(pow(fELo, fSpectralIndex+fDisplayIndex)/100); + + gStyle->SetOptStat(10); + + // + // Don't change the order!!! + // + histsrc.SetFillStyle(0); + histsrc.SetMarkerStyle(kMultiply); + histsrc.SetMarkerColor(kRed); + histsrc.SetLineColor(kRed); + + c=MH::MakeDefCanvas("Cascade", "Cascade"); + + c->SetGrid(); + c->SetLogx(); + c->SetLogy(); + + fHist.Draw("P"); + histsrc.Draw("Psame"); + histsrc.Draw("same"); + fHist.Draw("same"); + } + + // ------------------------------ + + cout << "Timers... " << flush; + + TTimer timer("gSystem->ProcessEvents();", 333, kFALSE); + if (!fIsBatch) + timer.TurnOn(); + + TStopwatch clock; + clock.Start(); + + cout << "Done. " << endl; + + Int_t n=0; + Double_t starttime = TStopwatch::GetRealTime(); // s + while (TStopwatch::GetRealTime() " << n << ". " << Form("%d: %3.1e", (int)(starttime+fRuntime-TStopwatch::GetRealTime()), E) << flush; + + ProcessPrimaryGamma(E, weight); + + if (fIsBatch) + continue; + + fHist.SetTitle(Form("E^{%.1f} z=%f T=%d'%d\" N=%d", fSpectralIndex, fSrcZ, + (int)fRuntime/60, (int)fRuntime%60, n)); + + timer.Stop(); + c->Update(); + timer.Start(250); + } + + cout << endl; + + clock.Stop(); + clock.Print(); + + timer.Stop(); + + CloseFile(file); + + cout << "Created " << n << " gammas (" << n/fNumBins << "rows) from source with E^" << fSpectralIndex << endl; + cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma ("; + cout << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n*fNumBins/60) << " min/row)" << endl; + + // ------------------------------ + + if (!fIsBatch) + { + c->Clear(); + + fHist.SetTitle(Form("E^{%.1f} z=%f T=%d'%d\" N=%d", fSpectralIndex, fSrcZ, (int)fRuntime/60, (int)fRuntime%60, n)); + + TH1 &h1 = *fHist.DrawCopy("P"); + TH1 &h2 = *histsrc.DrawCopy("Psame"); + h2.Draw("Csame"); + h1.Draw("Csame"); + } +} + Index: old/WuerzburgSoft/Thomas/mphys/MCascade.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MCascade.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MCascade.h (revision 15424) @@ -0,0 +1,79 @@ +#ifndef MARS_MCascade +#define MARS_MCascade + +#ifndef ROOT_TH1 +#include +#endif + +#ifndef MARS_MPairProduction +#include "MPairProduction.h" +#endif + +class TList; +class TFile; +class TBranch; + +class MPhoton; +class MElectron; + +class MCascade : public TObject +{ +private: + Double_t fSrcZ; + Double_t fSrcR; + + Double_t fB; + Double_t fBubbleZ; + + Double_t fRuntime; + + Double_t fEHi; + Double_t fELo; + + Double_t fSpectralIndex; + Double_t fDisplayIndex; + + Int_t fNumBins; + Int_t fNumMaxInvCompton; + Double_t fRatioInvCompton; + + // ------------------------------- + + MPairProduction fPair; //! + TBranch *fBranchGammas; //! + TBranch *fBranchElectrons; //! + + TH1D fHist; //! + + Bool_t fIsBatch; //! + + Double_t GetEnergy(); + TFile *OpenFile(TString fFilename); + void CloseFile(TFile *fFile); + void ProcessElectron(MElectron &e, TList &fListGammas); + Bool_t ProcessGamma(MPhoton &p, Double_t weight, TList &fListElectrons); + Double_t ProcessGammas(TList &fListGammas, TList &fListElectrons, Double_t weight); + Double_t ProcessElectrons(TList &fListElectrons, TList &fListGammas); + void ProcessPrimaryGamma(Double_t E, Double_t weight); + +public: + MCascade(); + ~MCascade(); + + void SetSourceZ(Double_t z); + void SetSourceRZ(Double_t r); // [kpc] + void SetB(Double_t B) { fB = B*1e-4; } // [G=1e-4T] mean magnetic field + void SetRuntime(Double_t min) { fRuntime = min*60; } // [s] maximum time to run the simulation + void SetEnergyBins(Int_t n, Double_t lo, Double_t hi); + void SetMaxInvCompton(Int_t max) { fNumMaxInvCompton = max; } // maximum number of inv. Compton (-1 means infinite) + void SetRatioInvCompton(Double_t r) { fRatioInvCompton=r; } + void SetSpectralIndex(Double_t idx) { fSpectralIndex=idx; } + void SetDisplayIndex(Double_t idx) { fDisplayIndex=idx; } + void SetBradius(Double_t r); // [Mpc] + + void Run(TString filename, Bool_t draw); + + ClassDef(MCascade, 1) +}; + +#endif Index: old/WuerzburgSoft/Thomas/mphys/MElectron.cc =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 15424) @@ -0,0 +1,607 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MElectron.h" + +#include // for aqlphas +#include + +#include +#include +#include +#include + +#include + +#include +#include + +#include "MPhoton.h" + +ClassImp(MElectron); + +Double_t MElectron::Sigma_ge(Double_t *x, Double_t *k) +{ + const Double_t E0 = 511e-6; // [GeV] + const Double_t c = 299792458; // [m/s] + const Double_t e = 1.602176462e-19; // [C] + const Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + const Double_t a = 1./137; // [1] + + const Double_t o = x[0]; + + const Double_t s1 = TMath::Pi()*2; + const Double_t s2 = a*h*c/E0; // [m] + + if (o<1e-4) + return s2*s2/s1 * 4/3; + + const Double_t o1 = o+1; + const Double_t o21 = o*2+1; + + const Double_t s3 = o1/(o*o*o); + const Double_t s4 = o*2*o1/o21; + const Double_t s5 = log(o21); + const Double_t s6 = s5/o/2; + const Double_t s7 = (o*3+1)/(o21*o21); + + return s2*s2/s1*(s3*(s4-s5)+s6-s7); // [m^2] +} + +Double_t MElectron::Li(Double_t *x, Double_t *k) +{ + const Double_t t = x[0]; + return log(1.-t)/t; +} + +Double_t MElectron::DiSum(Double_t *x, Double_t *k) +{ + Double_t t = x[0]; + + const Double_t eps = fabs(t*1e-2); + + Double_t disum = t; + Double_t add = 0; + + Int_t n = 2; + Double_t pow = t*t; // t^2 + + do + { + add = pow/n/n; + + pow *= t; // pow = t^n + n++; + + disum += add; + + } while (fabs(add)>eps); + + return disum; +} + +Double_t MElectron::Li2(Double_t *x, Double_t *k) +{ + // + // Dilog, Li2 + // ---------- + // + // Integral(0, 1) = konst; + // Double_t konst = 1./6*TMath::Pi()*TMath::Pi(); + // + // x[0]: z + // + const Double_t z = x[0]; + + if (fabs(z)<1) + return DiSum(x); + + // TF1 IntLi("Li", Li, 0, z, 0); + static TF1 IntLi("Li", Li, 0, 0, 0); + const Double_t integ = IntLi.Integral(0, z, (Double_t*)NULL, 1e-2); + return -integ; +} + +Double_t MElectron::Flim(Double_t *x, Double_t *k) // F(omegap)-F(omegam) mit b-->1 (Maple) +{ + const Double_t w = x[0]; + + const Double_t w4 = w*4; + const Double_t wsqr = w*w; + + const Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1); + const Double_t u2 = -w4*(wsqr*2 + w*9 + 2); + const Double_t d = w4*(w4 + 1); + + Double_t s = -w*2*(1+1); // -2*omega*(1+beta) + const Double_t li2 = Li2(&s); + + const Double_t res = (u1+u2)/d + li2; + + return res; //<1e-10? 0 : res; +} + +Double_t MElectron::Compton(Double_t *x, Double_t *k) +{ + const Double_t E0 = 511e-6; //[GeV] + + Double_t epsilon = x[0]; + Double_t z = k[1]; + + const Double_t E = k[0]; + + Double_t flim; + if (epsilon<1e-14) + { + const Double_t d = E/(E0*E0); + + Double_t omega1 = 1e-13*d; + Double_t omega2 = 1e-12*d; + + const Double_t f1 = Flim(&omega1); + const Double_t f2 = Flim(&omega2); + + const Double_t m = log10(f2/f1); + const Double_t t = pow(f2, 13)/pow(f1, 12); + + flim = pow(epsilon, m) * t; + } + else + { + Double_t omega = epsilon*E/(E0*E0); + flim = Flim(&omega); + } + + const Double_t n = MParticle::Planck(&epsilon, &z)/epsilon/epsilon; // [1] + return flim*n; +} + +Double_t MElectron::InteractionLength(Double_t *E, Double_t *k) +{ + // E = electron energy, ~ TeV(?) 1e12 + // e = photon energy, ~ meV(?) 1e-3 + // mc^2 = electron rest mass energy ~.5keV(?) .5e3 + // + // x^-1 = int( n(epsilon)/2beta * ((mc^2)^2/eE)^2 * int ( omega*sigma(omega), omega=o-..o+), epsilon=0..inf) + // + // o+/- = omage_0 (1 +- beta) + // + // omega_0 = eE/(mc^2)^2 ~1e12*1e-3/.25e6=4e3 + // + // --> x^-1 = (alpha*hc)^2/4pibetaE^2 * int(n(epsilon)/epsilon^2 *( F(o+)-F(o-)), epsilon=0..inf) + // + // F(o) = -o/4 + (9/4 + 1/o + o/2) * ln(1+2o) + 1/8(1+2o) - 3/8 + Li2(-2o) + // + // Li2(x) = int(ln(1-t)/t, t=0..x) + // + // F(o+)~F(2o) = -o/2 + (9/4 + 1/2o + o) * ln(1+4o) + 1/8(1+4o) - 3/8 + Li2(-4o) + // F(o-)~F(0) = 14/8 = 1.75 + + const Double_t E0 = 511e-6; // [GeV] + const Double_t E02 = E0*E0; // [GeV^2] + const Double_t c = 299792458; // [m/s] + const Double_t e = 1.602176462e-19; // [C] + const Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + const Double_t hc = h*c; // [GeVm] + const Double_t alpha = 1./137.; // [1] + + const Double_t z = k ? k[0] : 0; + + /* -------------- old ---------------- + Double_t from = 1e-15; + Double_t to = 1e-11; + eps = [default]; + ----------------------------------- + */ + static TF1 func("Compton", Compton, 0, 0, 2); // [0, inf] + + const Double_t from = 1e-17; + const Double_t to = 2e-11; + + Double_t val[3] = { E[0], z }; // E[GeV] + + Double_t integ = func.Integral(from, to, val, 1e-2); // [Gev] [0, inf] + + const Double_t aE = alpha/E[0]; // [1/GeV] + + const Double_t beta = 1; + + const Double_t konst = 2.*E02/hc/beta; // [1 / GeV m] + const Double_t ret = konst * (aE*aE) * integ; // [1 / m] + + const Double_t ly = 3600.*24.*365.*c; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + + return (1./ret)/ly*pc/1000; // [kpc] +} + +Double_t MElectron::GetInteractionLength(Double_t energy, Double_t z) +{ + return InteractionLength(&energy, &z); +} + +Double_t MElectron::GetInteractionLength() const +{ + return InteractionLength((Double_t*)&fEnergy, (Double_t*)&fZ); +} + +// -------------------------------------------------------------------------- + +/*inline*/ Double_t MElectron::p_e(Double_t *x, Double_t *k) +{ + Double_t e = pow(10, x[0]); + return Compton(&e, k); + /* + Double_t z = k[1]; + + const Double_t E = k[0]; + + const Double_t E0 = 511e-6; //[GeV] + const Double_t E02 = E0*E0; + + Double_t omega = e*E/E02; + + const Double_t n = Planck(&e, &z); + + const Double_t F = Flim(&omega)/omega/omega; + + return n*F*1e26; + */ +} + +Double_t MElectron::G_q(Double_t *x, Double_t *k) +{ + const Double_t q = x[0]; + const Double_t Gamma = k[0]; + + const Double_t Gq = Gamma*q; + + const Double_t s1 = 2.*q*log(q); + const Double_t s2 = (1.+2.*q); + const Double_t s3 = (Gq*Gq)/(1.+Gq)/2.; + + return s1+(s2+s3)*(1.-q); +} + + +Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep) +{ + const Double_t E = x[0]; + const Double_t z = k ? k[0] : 0; + + const Double_t E0 = 511e-6; //[GeV] + + const Double_t lolim = -log10(E)/7*4-13.5; + + static TF1 fP("p_e", p_e, lolim, -10.6, 2); + static TF1 fQ("G", G_q, 0, 1., 1); + + fP.SetNpx(50); + fQ.SetNpx(50); + + fP.SetRange(lolim, -10.6); + fP.SetParameter(0, E); + fP.SetParameter(1, z); + + const Double_t e = pow(10, fP.GetRandom()); + + if (ep) + *ep = e; + + const Double_t omega = (e/E0)*(E/E0); + const Double_t Gamma = 4.*omega; + + fQ.SetParameter(0, Gamma); + + const Double_t q = fQ.GetRandom(); + const Double_t Gq = Gamma*q; + + const Double_t e1 = Gq*E/(1.+Gq); + + return e1; +} + +Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z, Double_t *ep) +{ + return EnergyLoss(&E, &z); +} + +Double_t MElectron::GetEnergyLoss(Double_t *ep) const +{ + return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ, ep); +} + + +Double_t Omega_sigmae(Double_t *x, Double_t *k) +{ + Double_t sbar = pow(10,x[0]); + + Double_t omega = (sbar-1)/2; + + Double_t sigma = MElectron::Sigma_ge(&omega); + + return (sbar-1)*sigma*1e28; +} + +Double_t RandomThetaE(Double_t Ee, Double_t Ep) +{ + Double_t E0 = 511e-6; // [GeV] + + Double_t f = 2*Ee/E0*Ep/E0; + + static TF1 func("RndThetaE", Omega_sigmae, 0, 0, 0); + + Double_t beta = sqrt(1-E0/Ee*E0/Ee); + + //func.SetRange(0, log10(1+f*(1+beta))); + func.SetRange(log10(1+f*(1-beta)), log10(1+f*(1+beta))); + func.SetNpx(50); + + Double_t sbar = pow(10, func.GetRandom()); + + Double_t bcost = 1 - (sbar-1)/f; + return bcost; + + /* + Double_t theta = acos(bcost/beta); + return theta; + */ +} + +MPhoton *MElectron::DoInvCompton() +{ + const Double_t E0 = 511e-6; //[GeV] + + Double_t epsilon; + const Double_t e = GetEnergyLoss(&epsilon); + + // epsilon: photon energy befor interaction, lab + // e: photon energy after interaction, lab + + const Double_t gamma = fEnergy/E0; + const Double_t gammabeta = sqrt(gamma*gamma-1); + const Double_t beta = gammabeta/gamma; + + const Double_t bcost = RandomThetaE(fEnergy, epsilon); + const Double_t t = acos(bcost/beta); + + const Double_t f = epsilon/fEnergy; + const Double_t r = gamma*(1-bcost); + + Double_t arg = (1 - 1/(gamma*r) - f)/(beta-f); + + if (arg<-1 || arg>1) + cout << "<" << (int)(t*180/TMath::Pi()) << "°>" << flush; + + // + // Due to numerical uncertanties arg can be something like: + // 1+1e-15 which we cannot allow + // + if (arg<-1) + arg = -1; + if (arg>1) + arg = 1; + + const Double_t theta1s = acos(arg); + const Double_t thetas = atan(-sin(t)/(beta*r)/*(1-bcost)/gammabeta*/); + + const Double_t thetastar = thetas-theta1s; + + // const Double_t theta1 = atan(sin(thetastar)/(gamma*cos(thetastar)+gammabeta)); + const Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta))); + + fEnergy -= e; + + const Double_t phi = gRandom->Uniform(TMath::Pi()*2); + + MPhoton &p = *new MPhoton(*this, e); + p.SetNewDirection(theta1, phi); + + /* + const Double_t beta2 = sqrt(1.-E0/fEnergy*E0/fEnergy); + const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/fEnergy/beta2); + */ + + const Double_t div = sqrt(fEnergy*fEnergy-E0*E0); + const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/div); + + SetNewDirection(theta2, phi); + + return &p; +} + +void MElectron::DrawInteractionLength(Double_t z) +{ + if (!gPad) + new TCanvas("IL_Electron", "Mean Interaction Length Electron"); + else + gPad->GetVirtCanvas()->cd(3); + + TF1 f2("length", MElectron::InteractionLength, 1e3, 1e10, 0); + f2.SetParameter(0, z); + + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + f2.SetLineWidth(1); + + TH1 &h=*f2.DrawCopy()->GetHistogram(); + + h.SetTitle("Mean Interaction Length (Electron)"); + h.SetXTitle("E [GeV]"); + h.SetYTitle("x [kpc]"); + + gPad->Modified(); + gPad->Update(); +} + +void MElectron::DrawInteractionLength() const +{ + DrawInteractionLength(fZ); +} + +Bool_t MElectron::SetNewPositionB(Double_t B) +{ + if (B==0) + return SetNewPosition(); + + Double_t x = gRandom->Exp(GetInteractionLength()); + + // ----------------------------- + + const Double_t E0 = 511e-6; //[GeV] + + Double_t B_theta = gRandom->Uniform(TMath::Pi()); + Double_t B_phi = gRandom->Uniform(TMath::Pi()*2); + + static TMatrixD M(3,3); + + M(0, 0) = sin(B_theta)*cos(B_phi); + M(1, 0) = cos(B_theta)*cos(B_phi); + M(2, 0) = -sin(B_phi); + + M(0, 1) = sin(B_theta)*sin(B_phi); + M(1, 1) = cos(B_theta)*sin(B_phi); + M(2, 1) = cos(B_phi); + + M(0, 2) = cos(B_theta); + M(1, 2) = -sin(B_theta); + M(2, 2) = 0; + + const Double_t beta = sqrt(1.-E0/fEnergy*E0/fEnergy); + + // + // rotate vector of velocity into a system in + // which the x-axis is identical with the B-Field + // + static TVectorD v(3); + v(0) = beta*sin(fTheta)*cos(fPsi); + v(1) = beta*sin(fTheta)*sin(fPsi); + v(2) = beta*cos(fTheta); + v *= M; + + // + // Now rotate the system so, that the velocity vector + // lays in the y-z-plain + // + Double_t chi = atan(-v(2)/v(1)); + + // ----------------------------- + + static TMatrixD N(3,3); + N(0, 0) = 1; + N(1, 0) = 0; + N(2, 0) = 0; + + N(0, 1) = 0; + N(1, 1) = -cos(chi); + N(2, 1) = -sin(chi); + + N(0, 2) = 0; + N(1, 2) = sin(chi); + N(2, 2) = -cos(chi); + v *= N; + + Double_t beta_p = v(0); // beta parallel + Double_t beta_o = v(1); // beta orthogonal + // v(2) = 0 + // ----------------------------- + static TVectorD p(3); + + Double_t rho = 0; + if (B>0) + { + const Double_t c = 299792458; // [m/s] + const Double_t ly = 3600.*24.*365.*c; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + + Double_t r = (fEnergy*1e9)*beta_o/(c*B)*(pc/ly/1e3); // [kpc] + + rho = beta_o/beta*x/r; // [2pi] + + // ----------------------------- + + if (fabs(rho*3437)>5*60) // > 5*60min + { + cout << "r" << flush; + return kFALSE; + } + + p(0) = beta_p/beta*x; + p(1) = GetPType()==kEElectron?r*sin(rho):-r*sin(rho); + p(2) = r*(1.-cos(rho)); + } + else + { + p(0) = beta_p/beta*x; + p(1) = beta_o/beta*x; + p(2) = 0; + cout << "------------- HEY! --------------" << endl; + } + + static TVectorD s(3); + s(0) = fR*cos(fPhi); + s(1) = fR*sin(fPhi); + s(2) = RofZ(&fZ); + + TMatrixD N2(TMatrixD::kTransposed, N); + TMatrixD M2(TMatrixD::kTransposed, M); + p *= N2; + p *= M2; + + if (p(2) + +#include +#include + +#include + +ClassImp(MFit); + +TF1 *MFit::fgFunc=NULL; +TH1 *MFit::fgHist=NULL; + +Int_t MFit::fgBinMin; +Int_t MFit::fgBinMax; + +void MFit::Func(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag) +{ + Double_t chisq = 0; + + for (int i=fgBinMin+1; iGetBinCenter(i); + + Double_t y = fgFunc->EvalPar(&x, par); + Double_t d = (y - fgHist->GetBinContent(i))/fgHist->GetBinError(i); + + chisq += d*d; + } + + f = chisq / fgFunc->GetNDF(); + fgFunc->SetChisquare(f); +} + +void MFit::FuncLog(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag) +{ + Double_t chisq = 0; + + for (int i=fgBinMin+1; iGetBinLowEdge(i+1); + Double_t x0 = fgHist->GetBinLowEdge(i); + + Double_t x = sqrt(x1*x0); + + Double_t y = fgFunc->EvalPar(&x, par); + Double_t d = (y - fgHist->GetBinContent(i))/fgHist->GetBinError(i); + + chisq += d*d; + } + + f = chisq / fgFunc->GetNDF(); + fgFunc->SetChisquare(f); +} + +void MFit::DeleteF() +{ + if (fFunc) + if (fFunc->TestBit(kIsOwner)) + delete fFunc; + + if (fMinuit) + delete fMinuit; +} + +void MFit::DeleteH() +{ + if (fHist) + if (fHist->TestBit(kIsOwner)) + delete fHist; +} + +void MFit::SetFunc(TF1 *f) +{ + DeleteF(); + + fFunc = f; + + fMinuit = new TMinuit(fFunc->GetNpar()); + fMinuit->SetPrintLevel(-1); +} + +void MFit::SetFunc(TString formula, Axis_t min, Axis_t max) +{ + TF1 *f = new TF1("", formula, min, max); + f->SetBit(kIsOwner); + + SetFunc(f); +} + +void MFit::SetRange(Axis_t min, Axis_t max) +{ + if (fFunc) + fFunc->SetRange(min, max); +} + +void MFit::SetHist(TH1 *h, Bool_t candelete) +{ + DeleteH(); + + fHist = h; + + if (candelete) + fHist->SetBit(kIsOwner); +} + +void MFit::Fit(Bool_t log) +{ + fgBinMin = fHist->FindBin(fFunc->GetXmin()); // excluded + fgBinMax = fHist->FindBin(fFunc->GetXmax()); // excluded + + fFunc->SetNDF(fgBinMax-fgBinMin-1); + + fgHist = fHist; + fgFunc = fFunc; + + fMinuit->SetFCN(log ? FuncLog : Func); + fMinuit->Migrad(); + + for (int i=0; iGetNpar(); i++) + { + Double_t par, err; + fMinuit->GetParameter(i, par, err); + fFunc->SetParameter(i, par); + fFunc->SetParError(i, err); + } +} + +void MFit::Print(Option_t *opt) const +{ + for (int i=0; iGetNpar(); i++) + cout << "Par["<GetChisquare() << endl; + cout << "NDF: " << fFunc->GetNDF() << endl; +} + +void MFit::SetParameter(Int_t n, TString name, Axis_t start, Axis_t min, Axis_t max) +{ + if (n>=fFunc->GetNpar()) + return; + + fMinuit->DefineParameter(n, name, start, min/10, min, max); +} + +Double_t MFit::GetParameter(Int_t n) const +{ + return fFunc->GetParameter(n); +} + +Double_t MFit::GetParError(Int_t n) const +{ + return fFunc->GetParError(n); +} + +Double_t MFit::operator[](Int_t n) const +{ + return fFunc->GetParameter(n); +} + +TF1 *MFit::DrawCopy(Option_t *o) const +{ + return fFunc->DrawCopy(o); +} + Index: old/WuerzburgSoft/Thomas/mphys/MFit.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MFit.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MFit.h (revision 15424) @@ -0,0 +1,103 @@ +#ifndef MARS_MFit +#define MARS_Mfit + +/* +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "MH.h" +#include "MBinning.h" + +#include "MPhoton.h" +*/ +#ifndef ROOT_TString +#include +#endif +#ifndef ROOT_TObject +#include +#endif + +class TH1; +class TF1; + +class TMinuit; + +class MFit : public TObject +{ +private: + static void Func(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag); + static void FuncLog(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag); + + static TF1 *fgFunc; + static TH1 *fgHist; + + static Int_t fgBinMin; + static Int_t fgBinMax; + +private: + enum { + kIsOwner = BIT(14) + }; + + TH1 *fHist; + TF1 *fFunc; + + TMinuit *fMinuit; + + void DeleteF(); + void DeleteH(); + +public: + MFit() : fHist(NULL), fFunc(NULL), fMinuit(NULL) + { + } + + MFit(TString formula, Double_t min=0, Double_t max=1) : fHist(NULL), fFunc(NULL), fMinuit(NULL) + { + SetFunc(formula, min, max); + } + + ~MFit() + { + DeleteF(); + DeleteH(); + } + + void SetFunc(TF1 *f); + void SetFunc(TString formula, Double_t min, Double_t max); + void SetHist(TH1 *h, Bool_t candelete=kFALSE); + void SetRange(Double_t lo, Double_t hi); + + void Fit(Bool_t log=kFALSE); + void FitLog() { Fit(kTRUE); } + + void Fit(TH1 *h) { SetHist(h); Fit(); } + void FitLog(TH1 *h) { SetHist(h); FitLog(); } + + void Print(Option_t *opt="") const; + void SetParameter(Int_t n, TString name, Double_t start, Double_t min=0, Double_t max=0); + + Double_t GetParameter(Int_t n) const; + Double_t GetParError(Int_t n) const; + + Double_t operator[](Int_t n) const; + + TF1 *DrawCopy(Option_t *o="") const; + + ClassDef(MFit, 0) +}; + +#endif Index: old/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc (revision 15424) @@ -0,0 +1,116 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MGenIRPhoton.h" + +#include +#include + +#include "MParList.h" +#include "MPhoton.h" + +ClassImp(MGenIRPhoton); + +/* +Double_t Planck(Double_t *x, Double_t *k) +{ + const Double_t E = x[0]; //[eV] + + const Double_t T = 2.87; //[K] + + const Double_t e = 1.602176462e-19; //[C] + const Double_t kB = 1.3806503e-23/e; //[eV/K] + + const Double_t E3 = E*E*E; + const Double_t EkT = E/kB/T; + + //const Double_t h = 1e-9/e*6.62606876e-34; //[GeVs] + //const Double_t h2 = h*h; + //const Double_t c = 299792458; //[m/s] + //const Double_t c3 = c*c*c; + //const Double_t konst = 4/c3/h2; + + return E3/(exp(EkT)-1); +} +*/ + +Double_t MGenIRPhoton::Planck(Double_t *x, Double_t *k=NULL) +{ + // + // Planck, per unit volume, per unit energy + // + // constants moved out of function + // + Double_t E = x[0]; // [GeV] + Double_t z = k ? k[0] : 0; + + Double_t T = 2.96*(z+1); // [K] + Double_t e = 1.602176462e-19; // [C] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + + Double_t EkT = E/kB/T; + + /* + //Double_t c = 299792458; // [m/s] + //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + //Double_t hc = h*c; // [GeVm] + Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ] + */ + + return E*E / (exp(EkT)-1.); // [GeV^2] +} + + +// -------------------------------------------------------------------------- +MGenIRPhoton::MGenIRPhoton() +{ + fSrc = new TF1("Planck", Planck, 0, .5e-2, 0); +} + +MGenIRPhoton::~MGenIRPhoton() +{ + delete fSrc; +} + +// -------------------------------------------------------------------------- +MParticle *MGenIRPhoton::GetRandom() +{ + const Double_t pi2 = TMath::Pi() * 2; + + MPhoton *p = new MPhoton; + /* + MParticle *p = new MParticle(MParticle::kEGamma); + p->SetAngle(fRand.Uniform(pi2)); + p->SetEnergy(fSrc->GetRandom()*1e-9); + */ + + return p; +} + Index: old/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.h (revision 15424) @@ -0,0 +1,34 @@ +#ifndef MARS_MGenIRPhoton +#define MARS_MGenIRPhoton + +#ifndef ROOT_TRandom +#include +#endif + +class TF1; +class MParticle; + +class MGenIRPhoton +{ +private: + TF1 *fSrc; + + TRandom fRand; + + Double_t fZ; + + static Double_t Planck(Double_t *x, Double_t *k=NULL); + +public: + MGenIRPhoton(); + virtual ~MGenIRPhoton(); + + void SetZ(Double_t z) { fZ = z; } + + MParticle *GetRandom(); + + ClassDef(MGenIRPhoton, 0) // +}; + +#endif + Index: old/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc (revision 15424) @@ -0,0 +1,81 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MGenPrimaryParticle.h" + +#include +#include + +#include "MParList.h" +#include "MParticle.h" + +ClassImp(MGenPrimaryParticle); + +Double_t PowerLaw(Double_t *x, Double_t *k) +{ + return 1/x[0]; +} + +// -------------------------------------------------------------------------- +MGenPrimaryParticle::MGenPrimaryParticle(const char *name, const char *title) + : fPi2(TMath::Pi()*2) +{ + fName = name ? name : "MGenPrimaryParticle"; + fTitle = title ? title : ""; + + fSrc = new TF1("Power-Law", PowerLaw, 10, 1e8, 0); // 10GeV-100PeV +} + +MGenPrimaryParticle::~MGenPrimaryParticle() +{ + delete fSrc; +} + +// -------------------------------------------------------------------------- +Bool_t MGenPrimaryParticle::PreProcess(MParList *pList) +{ + fList = (MParList*)pList->FindCreateObj("MParticleList"); + if (!fList) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +Bool_t MGenPrimaryParticle::Process() +{ +// MParticle *p = new MParticle(MParticle::kEGamma); +// p->SetEnergy(fSrc->GetRandom()); + return kTRUE; +} + +Bool_t MGenPrimaryParticle::PostProcess() +{ + return kTRUE; +} Index: old/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.h (revision 15424) @@ -0,0 +1,39 @@ +#ifndef MARS_MGenPrimaryParticle +#define MARS_MGenPrimaryParticle + +#ifndef ROOT_TRandom +#include +#endif + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class TF1; +class MParList; +class MSinglePairInput; + +class MGenPrimaryParticle : public MTask +{ +private: + MParList *fList; + + TF1 *fSrc; + + TRandom fRand; + + Double_t fPi2; + +public: + MGenPrimaryParticle(const char *name=NULL, const char *title=NULL); + ~MGenPrimaryParticle(); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + ClassDef(MGenPrimaryParticle, 0) // +}; + +#endif + Index: old/WuerzburgSoft/Thomas/mphys/MPairProduction.cc =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 15424) @@ -0,0 +1,195 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 12/2000 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MPairProduction.h" + +#include +#include +#include + +#include "TList.h" +#include "MElectron.h" + +ClassImp(MPairProduction); + +Double_t AngleDistrib(Double_t *x, Double_t *k) +{ + + const Double_t c = x[0]; // cos(alpha) + const Double_t b = k[0]; // sqrt(1-1/s) + + const Double_t b2 = b*b; + const Double_t b4 = b2*b2; + + const Double_t c2 = c*c; + const Double_t c4 = c2*c2; + + const Double_t u = 1 - b4*c4 +2*b2*(1-b2)*(1-c2); + const Double_t d = 1-b2*c2; + + return u/(d*d); +} + +// -------------------------------------------------------------------------- +MPairProduction::MPairProduction() +{ + fAngle = new TF1("AngleDistrib", AngleDistrib, -1, 1, 1); +} + +MPairProduction::~MPairProduction() +{ + delete fAngle; +} + +#include + +// -------------------------------------------------------------------------- +Bool_t MPairProduction::Process(MParticle *gamma, const Double_t Ep, const Double_t theta, TList *list) +{ + // + // gamma: primary particle from source + // phot: infrared photon from background. (angle = interaction angle) + // Ep: Energy photon [GeV] + // theta: interaction angle [rad] + // + const Double_t E0 = 511e-6; // [GeV] + const Double_t Eg = gamma->GetEnergy(); // [GeV] + + const Double_t ctheta = cos(theta); + + const Double_t s = Ep*Eg*2*(1-ctheta); //[1] + const Double_t S = s/(E0*E0*4); //[1] + if (S<1) + return kFALSE; + + const Double_t stheta = sin(theta); + + const Double_t sqrbetah = s/((Eg+Ep)*(Eg+Ep)) + 1; + const Double_t sqrbetae = 1.-1./S; + + const Double_t GammaH = (Eg+Ep)/sqrt(s); + + const Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah))); + + fAngle->SetParameter(0, sqrt(sqrbetae)); + + const Double_t alpha = psi-acos(fAngle->GetRandom()); + + const Double_t salpha = sin(alpha); + const Double_t calpha = cos(alpha); + + const Double_t tphi = stheta/(Eg/Ep+ctheta); // tan(phi) + + const Double_t bb = sqrt(sqrbetah/sqrbetae); + + const Double_t s1 = calpha/GammaH; + const Double_t s2 = tphi*s1 - salpha - bb; + + const Double_t tan1 = ((salpha+bb)*tphi+s1)/s2; + const Double_t tan2 = ((salpha-bb)*tphi+s1)/s2; + + const Double_t E = (Eg+Ep)/2;; + const Double_t f = sqrt(sqrbetah*sqrbetae)*salpha; + + // cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush; + + MElectron &p0 = *new MElectron(*gamma, E*(1.-f), kTRUE); + MElectron &p1 = *new MElectron(*gamma, E*(1.+f), kFALSE); + + Double_t rnd = gRandom->Uniform(TMath::Pi() * 2); + p0.SetNewDirection(atan(+tan1), rnd); + p1.SetNewDirection(atan(-tan2), rnd+TMath::Pi()); + + list->Add(&p0); + list->Add(&p1); + + /* + const Double_t Epg = Ep/Eg; // 1+Epg ~ 1 + const Double_t Egp = Eg/Ep; // 1+Egp ~ Egp + + const Double_t phi = atan(sin(theta)/(Egp+ctheta)); + + const Double_t sphi = sin(phi); + const Double_t cphi = cos(phi); + + const Double_t alpha = theta-phi; + const Double_t calpha = cos(alpha); + const Double_t salpha = sin(alpha); + + const Double_t beta = (Eg*cphi+Ep*calpha)/(Ep+Eg); + + // + // gamma1 = 1/gamma = sqrt(1-beta*beta) + // + const Double_t gamma1 = sqrt((sphi*phi+Epg*salpha*Epg*salpha+2*Epg*(1-cphi*calpha))); + + const Double_t Beta = sqrt(1-1/s); //[1] + + fAngle->SetParameter(0, Beta); + + const Double_t psi = atan(gamma1*sphi/(cphi-beta)); + const Double_t delta = acos(fAngle->GetRandom()) - psi; + + const Double_t Bcosd = Beta*cos(delta); + const Double_t Bsind = Beta*sin(delta); + + const Double_t E = sqrt(s)*E0/gamma1; + const Double_t dE = E*Bcosd; + + const Double_t E1 = E0/(E+dE); + const Double_t E2 = E0/(E-dE); + + const Double_t beta1 = sqrt(1.-E1*E1); + const Double_t beta2 = sqrt(1.-E2*E2); + + const Double_t Bscp = Bsind*cphi; + const Double_t spg = sphi/gamma1; + const Double_t cpg = cphi/gamma1; + + const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp)); + const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp)); + + MElectron &p0 = *new MElectron(E+dE, gamma->GetZ()); + MElectron &p1 = *new MElectron(E-dE, gamma->GetZ()); + p0 = *gamma; // Set Position and direction + p1 = *gamma; // Set Position and direction + + p0.SetBeta(beta1); + p1.SetBeta(beta2); + + static TRandom rand(0); + Double_t rnd = rand.Uniform(TMath::Pi() * 2); + p0.SetNewDirection(atan(tan1), rnd); + p1.SetNewDirection(atan(tan2), rnd+TMath::Pi()); + + list->Add(&p0); + list->Add(&p1); + */ + return kTRUE; +} + Index: old/WuerzburgSoft/Thomas/mphys/MPairProduction.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MPairProduction.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MPairProduction.h (revision 15424) @@ -0,0 +1,27 @@ +#ifndef MARS_MPairProduction +#define MARS_MPairProduction + +#ifndef ROOT_TObject +#include +#endif + +class TF1; +class MParticle; +class TList; + +class MPairProduction +{ +private: + TF1 *fAngle; + +public: + MPairProduction(); + virtual ~MPairProduction(); + + Bool_t Process(MParticle *g, const Double_t Ep, const Double_t theta, TList *l); + + ClassDef(MPairProduction, 0) // +}; + +#endif + Index: old/WuerzburgSoft/Thomas/mphys/MParticle.cc =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 15424) @@ -0,0 +1,315 @@ +/////////////////////////////////////////////////////////////////////// +// +// MParticle +// +/////////////////////////////////////////////////////////////////////// + +#include "MParticle.h" + +#include +#include +#include + +ClassImp(MParticle); + +/************************************************** + * + * H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + * + **************************************************/ + +Double_t MParticle::ZofR(Double_t *x, Double_t *k) +{ + /* + const Double_t c = 299792458; // [m/s] + const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + + const Double_t ly = 3600.*24.*365.*c; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + const Double_t r = x[0] /pc*ly*1e3; // [m] + + const Double_t R = r*H0/c; // [1] + + return (R+1+sqrt(R*2+1))/2 - 1; + */ + const Double_t c = 299792458; // [m/s] + const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + + const Double_t ly = 3600.*24.*365.*c; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + const Double_t r = x[0] /pc*ly*1e3; // [m] + + const Double_t R = 1./(1.-r*H0/c/2); // [1] + + return R*R - 1; +} + +Double_t MParticle::RofZ(Double_t *x, Double_t *k) +{ + /* + Double_t z1 = x[0] + 1; + + const Double_t c = 299792458; // [m/s] + const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + + const Double_t ly = 3600.*24.*365.*c; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + + const Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m] + + return R * pc/ly/1e3; // [kpc] + */ + Double_t z1 = x[0] + 1; + + const Double_t c = 299792458; // [m/s] + const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + + const Double_t ly = 3600.*24.*365.*c; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + + const Double_t R = c/H0 * 2 * (1. - 1./sqrt(z1)); // [m] + + return R * pc/ly/1e3; // [kpc] +} + +#include +#include +#include "../mhist/MBinning.h" +#include "../mhist/MH.h" + +TH2D *hist2; + +Double_t MParticle::Planck(Double_t *x, Double_t *k) +{ + static Bool_t isloaded = kFALSE; + + if (!isloaded) + { + Double_t c = 299792458; // [m/s] + Double_t e = 1.602176462e-19; // [C] + Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + Double_t hc = h*c; // [GeVm] + Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + + ifstream fin("mphys/background.txt"); + + hist2 = new TH2D; + + MBinning binsz; + MBinning binse; + binsz.SetEdgesLog(100, 1e-6, 1); // --> 101 Edges / 100 bins + binse.SetEdgesLog(100, 7e-15, 3e-8); // --> 101 Edges / 100 bins + + MH::SetBinning(hist2, &binsz, &binse); + + for (int y=0; y<101; y++) + { + Double_t val; + fin >> val; + for (int x=0; x<101; x++) + { + fin >> val; + + val += 9; + + Double_t z = binsz.GetEdges()[x]; + Double_t f = z+1; + + Double_t newval = pow(10, val)/konst; + hist2->SetBinContent(x, y, newval*f*f*f); + + } + } + isloaded = kTRUE; + } + + static TAxis &axez = *hist2->GetXaxis(); + static TAxis &axee = *hist2->GetYaxis(); + + // + // y = (y1-y0)/(x1-x0) * (x-x0) + y0 + // + Double_t z = k ? k[0] : 0; + Double_t E = x[0]; + + Int_t i = axez.FindFixBin(z); + Int_t j = axee.FindFixBin(E); + + Double_t z1 = axez.GetBinLowEdge(i+1); + Double_t z0 = axez.GetBinLowEdge(i); + + Double_t E1 = axee.GetBinLowEdge(j+1); + Double_t E0 = axee.GetBinLowEdge(j); + + Double_t n00 = hist2->GetBinContent(i, j); + Double_t n01 = hist2->GetBinContent(i+1, j); + Double_t n10 = hist2->GetBinContent(i, j+1); + Double_t n11 = hist2->GetBinContent(i+1, j+1); + + Double_t dz = (z-z0)/(z1-z0); + Double_t dE = (E-E0)/(E1-E0); + + Double_t n0 = dz*(n01-n00)+n00; + Double_t n1 = dz*(n11-n10)+n10; + + Double_t n = dE*(n1-n0)+n0; + + return n; + /* + // + // TANJA2 + // + Double_t E1 = x[0]; + Double_t E2 = x[0]/8; + return (MParticle::Planck0(&E1, k)+MParticle::Planck0(&E2, k)/40e3)*0.7/0.4; + */ + /* + // + // TANJA + // + Double_t E1 = x[0]; + Double_t E2 = x[0]/8; + return Planck0(&E1, k)+Planck0(&E2, k)/5e3; + */ +} + +Double_t MParticle::Planck0(Double_t *x, Double_t *k) +{ + // + // Planck, per unit volume, per unit energy + // + // constants (see below) moved out of function + // + const Double_t E = x[0]; // [GeV] + const Double_t z = k ? k[0] : 0; + + const Double_t T = 2.96*(z+1); // [K] + const Double_t e = 1.602176462e-19; // [C] + const Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + + const Double_t EkT = E/kB/T; + + /* + //Double_t c = 299792458; // [m/s] + //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + //Double_t hc = h*c; // [GeVm] + Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ] + */ + + return E*E / (exp(EkT)-1.); // [GeV^2] +} + +MParticle::MParticle(ParticleType_t t, const char *name, const char *title) + : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0), fX(0) +{ + // + // default constructor + // +} + +void MParticle::InitRandom() +{ + fPhi = gRandom->Uniform(TMath::Pi()*2); + fPsi = gRandom->Uniform(TMath::Pi()*2); +} + +void MParticle::SetNewDirection(Double_t theta, Double_t phi) +{ + static TMatrixD B(3, 3); + + const Double_t ct = cos(fTheta); + const Double_t st = sin(fTheta); + const Double_t cp = cos(fPsi); + const Double_t sp = sin(fPsi); + + /* + TRotation B( ct*cp, ct*sp, -st, + -sp, cp, 0, + st*cp, st*sp, ct); + */ + + // first row + B(0, 0) = ct*cp; + B(1, 0) = ct*sp; + B(2, 0) = -st; + + // second row + B(0, 1) = -sp; + B(1, 1) = cp; + B(2, 1) = 0; + + // third row + B(0, 2) = st*cp; + B(1, 2) = st*sp; + B(2, 2) = ct; + + // ------------------------------ + + static TVectorD r(3); + + const Double_t sint = sin(theta); + + /* + TVector3 r(sint*cos(phi), sint*sin(phi), cos(theta)); + */ + + r(0) = sint*cos(phi); + r(1) = sint*sin(phi); + r(2) = cos(theta); + + // ------------------------------ + + r *= B; + + fTheta = asin(sqrt(r(0)*r(0)+r(1)*r(1))); // Numerically bad: acos(r(2)); + fPsi = atan2(r(1), r(0)); + + /* + if (fTheta*2 > TMath::Pi()) + fTheta = fabs(fTheta-TMath::Pi()); + */ +} + +Bool_t MParticle::SetNewPosition(Double_t dr) +{ + Bool_t rc=kTRUE; + + const Double_t st = sin(fTheta); + + TVector3 d(st*cos(fPsi), st*sin(fPsi), cos(fTheta)); + + // ------------------------------ + + const Double_t R = RofZ(&fZ); + const Double_t dx = R - dr*d.z(); + if (dx < 0) + { + dr = R/d.z(); // R>0 --> x(2)>0 + rc = kFALSE; + } + + fX += dr*(1.-d(2)); + d *= dr; + + // ------------------------------ + + TVector3 r(fR*cos(fPhi), fR*sin(fPhi), R); + + r -= d; + + // ------------------------------ + + if (fR!=0) + fPhi = atan2(r.y(), r.x()); + fR = sqrt(r.x()*r.x()+r.y()*r.y()); + fZ = ZofR(&r(2)); + + return rc; +} + +Bool_t MParticle::SetNewPosition() +{ + Double_t r = gRandom->Exp(GetInteractionLength()); + return SetNewPosition(r); +} Index: old/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 15424) @@ -0,0 +1,97 @@ +#ifndef MARS_MParticle +#define MARS_MParticle + +#ifndef ROOT_TObject +#include +#endif + +/* + #ifndef MARS_MParContainer + #include "MParContainer.h" + #endif + */ + +class MParticle : public TObject +{ +public: + typedef enum { kEGamma, kEElectron, kEPositron, kENone } ParticleType_t; + enum { kEIsPrimary = BIT(14) }; + +private: + const ParticleType_t fPType; //! Particle type + +protected: + Double_t fSrcR; //! [kpc] Distance from observer to source + + Double_t fEnergy; // [GeV] Energy + + Double_t fZ; // [1] Red shift + Double_t fR; // [kpc] Radius from line of view + Double_t fPhi; // [rad] Phi@z from line of view + + Double_t fTheta; // [rad] Direction of momentum, angle || line of view + Double_t fPsi; // [rad] Direction of momentum, az. angle + + Double_t fX; // [kpc] real traveled distance minus observers distance (added in version 2) + +public: + MParticle(ParticleType_t t=kENone, const char *name=NULL, const char *title=NULL); + MParticle(MParticle &p, ParticleType_t t=kENone) : fPType(t) { *this = p; } + MParticle(MParticle &p, Double_t e, ParticleType_t t=kENone) : fPType(t) { *this = p; fEnergy = e; } + + void InitRandom(); + + static Double_t ZofR(Double_t *x, Double_t *k=NULL); + static Double_t RofZ(Double_t *x, Double_t *k=NULL); + + void operator=(MParticle &p) + { + fSrcR = p.fSrcR; + fZ = p.fZ; + fR = p.fR; + fPhi = p.fPhi; + fTheta = p.fTheta; + fPsi = p.fPsi; + fX = p.fX; + } + + // ---------------------------------------------------------------- + + static Double_t Planck0(Double_t *x, Double_t *k=NULL); + static Double_t Planck(Double_t *x, Double_t *k=NULL); + + // ---------------------------------------------------------------- + + void SetIsPrimary(Bool_t is=kTRUE) { is ? SetBit(kEIsPrimary) : ResetBit(kEIsPrimary); } + Bool_t IsPrimary() const { return TestBit(kEIsPrimary); } + + void SetSrcR(Double_t r) { fSrcR = r; } + + // ---------------------------------------------------------------- + + virtual Double_t GetInteractionLength() const { AbstractMethod("GetInteractionLength"); return 0; } + + void SetEnergy(Double_t e) { fEnergy = e; } + void SetZ(Double_t z) { fZ = z; } + + void SetNewDirection(Double_t theta, Double_t phi); + Bool_t SetNewPosition(Double_t dr); + Bool_t SetNewPosition(); + + ParticleType_t GetPType() const { return fPType; } + + Double_t GetEnergy() const { return fEnergy; } + + Double_t GetZ() const { return fZ; } + Double_t GetPhi() const { return fPhi; } + Double_t GetR() const { return fR; } + Double_t GetX() const { return fX; } + + Double_t GetTheta() const { return fTheta; } + Double_t GetPsi() const { return fPsi; } + + ClassDef(MParticle, 2) // Container which holds hostograms for the Hillas parameters +}; + +#endif + Index: old/WuerzburgSoft/Thomas/mphys/MPhoton.cc =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 15424) @@ -0,0 +1,206 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MPhoton.h" + +#include + +#include +#include +#include +#include + +ClassImp(MPhoton); + +Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k) +{ + const Double_t m2 = x[0]; // m2: (E0/sqrt(s))^2 + + const Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2 + + const Double_t beta2 = 1.-m2; + const Double_t beta = sqrt(beta2); + + const Double_t p1 = r0*r0*TMath::Pi()/2; + + // ----- Extreme Relativistic ------- + // return p1*2 * m*m*m* (log(2./m)-1); + // ---------------------------------- + + const Double_t p2 = 3.-beta2*beta2; + const Double_t p3 = log((1.+beta)/(1.-beta)); + const Double_t p4 = beta*2*(1.+m2); + + const Double_t sigma = p1*m2*(p2*p3-p4); // [m^2] + + return sigma; +} + +Double_t MPhoton::Int1(Double_t *x, Double_t *k) +{ + const Double_t costheta = x[0]; + + const Double_t Eg = k[0]; + const Double_t Ep = k[1]; + + const Double_t E0 = 511e-6; // [GeV] + + Double_t s = E0/Eg*E0/Ep/(1.-costheta)/2; + if (s>1) // Why is this necessary??? + return 0; + + const Double_t sigma = Sigma_gg(&s); // [m^2] + + return sigma/2 * (1.-costheta); // [m^2] +} + +Double_t MPhoton::Int2(Double_t *x, Double_t *k) +{ + const Double_t E0 = 511e-6; // [GeV] + + Double_t Ep = x[0]; + Double_t z = k[1]; + + const Double_t Eg = k[0]; + + Double_t val[2] = { Eg, Ep }; + + static TF1 f("int1", Int1, 0, 0, 2); + + const Double_t from = -1.0; + const Double_t to = 1.-E0*E0/(2.*Eg*Ep); // Originally Was: 1. + const Double_t int1 = f.Integral(from, to, val, 1e-2); // [m^2] + const Double_t planck = MParticle::Planck(&Ep, &z); // [GeV^2] + + const Double_t res = planck * int1; + + return res; // [GeV^2 m^2] +} + +// -------------------------------------------------------------------------- +// +// Returns 0 in case IL becomes (numerically) infinite. +// +Double_t MPhoton::InteractionLength(Double_t *x, Double_t *k) +{ + Double_t E0 = 511e-6; // [GeV] + Double_t c = 299792458; // [m/s] + Double_t e = 1.602176462e-19; // [C] + Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + Double_t hc = h*c; // [GeVm] + Double_t pc = 1./3.258; // [pc/ly] + Double_t ly = 3600.*24.*365.*c; // [m/ly] + + Double_t Eg = x[0]; + Double_t z = k ? k[0] : 0; + + if (Eg<100) + return 1e50; + + Double_t val[2] = { Eg, z }; + + static TF1 f("int2", Int2, 0, 0, 2); + + Double_t lolim = E0*E0/Eg; + Double_t inf = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1)); + + if (Eg<5e4) + //inf = 3e-11*(z+1)*pow(10, 4.7*0.5-log10(Eg)*0.5); + inf = 3e-11*(z+1)*pow(10, 4.7-log10(Eg)); + //inf = 3e-11*(z+1)*pow(10, 7.0-log10(Eg)*1.5); + //inf = 3e-11*(z+1)*pow(10, 8.2-log10(Eg)*1.75); + //inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2); + + Double_t int2 = f.Integral(lolim, inf, val, 1e-2); //[GeV^3 m^2] + + if (int2==0) + { + //cout << "---> Int2==0 <---" << endl; + return 0; + } + + /* Planck constants: konst */ + const Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + int2 *= konst; // [1 / m] + + Double_t res = 1./ int2; // [m] + res *= pc/ly * 1e-3; // [kpc] + + if (res > 1e50) return 1e50; + if (res < 0) return 1e35; + + return res; //[kpc] +} + +Double_t MPhoton::GetInteractionLength(Double_t energy, Double_t z) +{ + return InteractionLength(&energy, &z); +} + +Double_t MPhoton::GetInteractionLength() const +{ + return InteractionLength((Double_t*)&fEnergy, (Double_t*)&fZ); +} + +void MPhoton::DrawInteractionLength(Double_t z, Double_t from, Double_t to, Option_t *opt) +{ + if (!gPad) + new TCanvas("ILPhoton", "Mean Interaction Length Photon"); + else + gPad->GetVirtCanvas()->cd(4); + + TF1 f1("length", InteractionLength, from, to, 1); + f1.SetParameter(0, z); + + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + f1.SetMinimum(1); + f1.SetMaximum(1e9); + f1.SetLineWidth(1); + + TH1 &h=*f1.DrawCopy(opt)->GetHistogram(); + + h.SetTitle("Mean Interaction Length (Photon)"); + h.SetXTitle("E [GeV]"); + h.SetYTitle("x [kpc]"); + + gPad->Modified(); + gPad->Update(); +} + +void MPhoton::DrawInteractionLength() const +{ + DrawInteractionLength(fZ); +} + +void MPhoton::Fill(TH1 &h, Double_t idx, Double_t w) const +{ + h.Fill(fEnergy, pow(fEnergy, idx)*w); +} Index: old/WuerzburgSoft/Thomas/mphys/MPhoton.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 15424) @@ -0,0 +1,46 @@ +#ifndef MARS_MPhoton +#define MARS_MPhoton + +#ifndef MARS_MParticle +#include "MParticle.h" +#endif + +class TH1; + +class MPhoton : public MParticle +{ +public: + + MPhoton(Double_t e=0, Double_t z=0) + : MParticle(MParticle::kEGamma) + { + fEnergy = e; + fZ = z; + } + + MPhoton(MParticle &p) : MParticle(p, MParticle::kEGamma) { } + MPhoton(MParticle &p, Double_t e) : MParticle(p, e, MParticle::kEGamma) { } + + void operator=(MParticle &p) { MParticle::operator=(p); } + + void Fill(TH1 &h, Double_t idx, Double_t w) const; + + // ---------------------------------------------------------------- + + static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL); + static Double_t Int1(Double_t *x, Double_t *k=NULL); + static Double_t Int2(Double_t *x, Double_t *k); + static Double_t InteractionLength(Double_t *x, Double_t *k=NULL); + static Double_t GetInteractionLength(Double_t energy, Double_t z=0); + + Double_t GetInteractionLength() const; + + // ---------------------------------------------------------------- + + static void DrawInteractionLength(Double_t z, Double_t from=2e2, Double_t to=1e11, Option_t*opt=""); + void DrawInteractionLength() const; + + ClassDef(MPhoton, 1) +}; + +#endif Index: old/WuerzburgSoft/Thomas/mphys/Makefile =================================================================== --- old/WuerzburgSoft/Thomas/mphys/Makefile (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/Makefile (revision 15424) @@ -0,0 +1,52 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +# +# Handling name of the Root Dictionary Files +# +CINT = Phys + +# +# Library name to creatre +# +LIB = mphys.a + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mhist + +#------------------------------------------------------------------------------ + +.SUFFIXES: .c .cc .cxx .h .hxx .o + +SRCFILES = MParticle.cc \ + MPairProduction.cc \ + MPhoton.cc \ + MFit.cc \ + MCascade.cc \ + MElectron.cc + +SRCS = $(SRCFILES) +HEADERS = $(SRCFILES:.cc=.h) +OBJS = $(SRCFILES:.cc=.o) + +############################################################ + +all: $(LIB) + +include ../Makefile.rules + +clean: rmcint rmobjs rmcore rmlib + +mrproper: clean rmbak + +# @endcode + Index: old/WuerzburgSoft/Thomas/mphys/PhysIncl.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/PhysIncl.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/PhysIncl.h (revision 15424) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: old/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h =================================================================== --- old/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h (revision 15424) @@ -0,0 +1,14 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MParticle+; +#pragma link C++ class MElectron+; +#pragma link C++ class MPhoton+; +#pragma link C++ class MCascade+; +#pragma link C++ class MFit+; +#pragma link C++ class MPairProduction+; + +#endif Index: old/WuerzburgSoft/Thomas/mphys/anale.C =================================================================== --- old/WuerzburgSoft/Thomas/mphys/anale.C (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/anale.C (revision 15424) @@ -0,0 +1,139 @@ +//Double_t kRad2Deg = 180./TMath::Pi(); + +void anale() +{ + TChain chain("Electrons"); +// chain.Add("cascade_100kpc_0*.root"); +// chain.Add("cascade_100kpc_14*.root"); +// chain.Add("cascade_100kpc_4*.root"); + chain.Add("cascade_100kpc_21d.root"); +/* + chain.Add("cascade_100kpc_01.root"); + chain.Add("cascade_100kpc_02.root"); + chain.Add("cascade_100kpc_03.root"); + */ + + + MElectron *e = new MElectron; + chain.SetBranchAddress("MElectron.", &e); + + Int_t n = chain.GetEntries(); + + cout << "Found " << n << " entries." << endl; + + MBinning binsx; + MBinning binsy; + MBinning binsr; + binsx.SetEdgesLog(80, 1e1, 1e9); + binsy.SetEdgesLog(80, 1e1, 1e9); + binsr.SetEdgesLog(84, 1e-5, 1.8); + + TH2D h; + h.SetName("Photons"); + h.SetTitle(" Photons from inv.Compton "); + h.SetXTitle("E_{e} [GeV]"); + h.SetYTitle("E\\gamma [GeV]"); + h.SetZTitle("Counts"); + MH::SetBinning(&h, &binsx, &binsy); + + TH2D r; + r.SetName("ratio"); + r.SetTitle("Ratio E\\gamma / E_{e}"); + r.SetXTitle("E_{e} [GeV]"); + r.SetYTitle("Ratio"); + r.SetZTitle("Counts"); + MH::SetBinning(&r, &binsx, &binsr); + + Double_t E = -1; + for (int i=0; iIsPrimary()) + { + E = e->GetEnergy(); + continue; + } + + Double_t dE = E - e->GetEnergy(); + + h.Fill(E, dE); + r.Fill(E, dE/E); + + E = e->GetEnergy(); + } + + delete e; + + gStyle->SetOptStat(10); + + TLine line; + line.SetLineColor(kBlack); + line.SetLineWidth(1); + + TCanvas *c = new TCanvas("c1", "name"); + c->Divide(2,2); + + c->cd(1); + gPad->SetLogx(); + gPad->SetLogy(); + h.GetXaxis()->SetRangeUser(5e3, 1e9); + h.GetXaxis()->SetLabelOffset(-0.015); + h.GetXaxis()->SetTitleOffset(1.1); + h.GetYaxis()->SetTitleOffset(1.1); +// h.SetMinimum(1e1); +// h.SetMaximum(1e9); + h.DrawCopy(); + TH1 *p=h.ProfileX(); + p->SetDirectory(NULL); + p->SetBit(kCanDelete); + p->SetLineColor(kRed); + p->Draw("same"); + line.DrawLine(4.2, 4.2, 8.8, 8.8); + + c->cd(2); + gPad->SetLogx(); + gPad->SetLogy(); + r.GetXaxis()->SetLabelOffset(-0.015); + r.GetXaxis()->SetTitleOffset(1.1); + r.GetXaxis()->SetRangeUser(1e4, 1e9); + r.DrawCopy(); + p=r.ProfileX(); + p->SetDirectory(NULL); + p->SetBit(kCanDelete); + p->SetLineColor(kRed); + p->Draw("same"); + + c->cd(3); + gPad->SetLogx(); + p=h.ProjectionX(); + p->SetDirectory(NULL); + p->SetBit(kCanDelete); + p->SetTitle("e^{-} / \\gamma Distribution"); + p->GetXaxis()->SetLabelOffset(-0.015); + p->GetXaxis()->SetTitleOffset(1.1); + p->GetYaxis()->SetTitleOffset(1.3); + p->SetXTitle("E [GeV]"); + p->SetYTitle("Counts"); + p->Draw(); + p=h.ProjectionY(); + p->SetDirectory(NULL); + p->SetBit(kCanDelete); + p->SetTitle("Projection Y"); + p->SetLineColor(kBlue); + p->Draw("same"); + + c->cd(4); + gPad->SetLogx(); + p=r.ProjectionY(); + p->SetDirectory(NULL); + p->SetBit(kCanDelete); + p->SetTitle("Ratio E\\gamma / E_{e}"); + p->GetXaxis()->SetLabelOffset(-0.015); + p->GetXaxis()->SetTitleOffset(1.1); + p->GetYaxis()->SetTitleOffset(1.3); + p->SetXTitle("Ratio"); + p->SetYTitle("Counts"); + p->SetLineColor(kBlue); + p->Draw(); +} Index: old/WuerzburgSoft/Thomas/mphys/analp.C =================================================================== --- old/WuerzburgSoft/Thomas/mphys/analp.C (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/analp.C (revision 15424) @@ -0,0 +1,600 @@ +//Double_t kRad2Deg = 180./TMath::Pi(); + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "MH.h" +#include "MBinning.h" + +#include "MFit.h" +#include "MPhoton.h" + +TH1D Smear(const TH1 &h, Double_t sigma) +{ + TH1D ret; + + MH::SetBinning((TH1*)&ret, (TH1*)&h); + + Int_t n = h.GetNbinsX(); + + for (int i=1; i<=n; i++) + { + Double_t xi = sqrt(h.GetBinLowEdge(i+1)*h.GetBinLowEdge(i)); + for (int j=1; j<=n; j++) + { + Double_t xj = sqrt(h.GetBinLowEdge(j+1)*h.GetBinLowEdge(j)); + + Double_t dx = log10(xj/xi)/sigma; + + Double_t gaus = exp(-dx*dx/2)/(sigma*sqrt(TMath::Pi()*2)); + + gaus *= h.GetBinContent(i); + + ret.Fill(xj, gaus/sqrt(n)); + } + } + return ret; +} + +void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1, Bool_t zero=kTRUE) +{ + TString str("MPhoton.MParticle."); + str += name; + + MPhoton *p = new MPhoton; + + chain->SetBranchAddress("MPhoton.", &p); + chain->SetBranchStatus("*", 0); + chain->SetBranchStatus(str, 1); + + min = FLT_MAX; + max = -FLT_MAX; + + Int_t numt = chain->GetTreeNumber(); + TLeaf *leaf = chain->GetLeaf(str); + + if (!leaf && numt>=0) + return; + + Stat_t n = chain->GetEntries(); + for (int i=0; iGetEntry(i); + + if (numt != chain->GetTreeNumber()) + { + if (!(leaf = chain->GetLeaf(str))) + return; + + numt = chain->GetTreeNumber(); + } + + Double_t val = leaf->GetValue(); + + if (p->IsPrimary()) + continue; + + if (val < min && (zero || val!=0)) + min = val; + if (val > max) + max = val; + } + + min *= conv; + max *= conv; + + Int_t newn=0; + MH::FindGoodLimits(nbins, newn, min, max, kFALSE); + nbins = newn; +} + +void draw1h1426(Double_t E, Double_t plot=1) +{ + plot += 1; + + const Double_t level = E / (5.73e-0*pow( 0.39e3, plot)); + + TPolyLine *l = new TPolyLine(6); + TPolyMarker *m = new TPolyMarker(6); + + // m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level)); + m->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level)); + m->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level)); + m->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level)); + m->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level)); + m->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level)); + m->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level)); + m->SetMarkerStyle(kOpenStar); + m->SetMarkerColor(kMagenta); + m->SetBit(kCanDelete); + m->Draw(); + + // l->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level)); + l->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level)); + l->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level)); + l->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level)); + l->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level)); + l->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level)); + l->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level)); + l->SetLineColor(kMagenta); + l->SetBit(kCanDelete); + l->Draw(); +} + +void analp() +{ + // ------------------------------------------------------------------- + + TString dir = "~/data/"; + + TChain chain("Photons"); + /* + chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_01.root"); + chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_02.root"); + chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_03.root"); + chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_04.root"); + + chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_01.root"); + chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_02.root"); + chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_03.root"); + chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_05.root"); + */ + + //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B0_512_01.root"); + //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_10Mpc_512_01.root"); + //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_50Mpc_512_01.root"); + chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_100Mpc_512_01.root"); + + //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-6_10Mpc_512_01.root"); + //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-6_50Mpc_512_01.root"); + //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-7_100Mpc_512_01.root"); + + /* + chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_01.root"); + chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_02.root"); + chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_03.root"); + chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_04.root"); + */ + + //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-6_10Mpc_512_01.root"); + //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-6_50Mpc_512_01.root"); + //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-7_100Mpc_512_01.root"); + + Int_t nbinse = 18; // number of bins in your histogram + Double_t lo = 1e2/1e3; // lower limit of the energy histogram + Double_t hi = 1e5; // upper limit of the energy histogram + + Double_t minE = 1e2; // minimum value produced in the simulation + + Double_t cutoff = 1e12; // arbitrary energy cutoff (throw away the events) + + Double_t alpha = -1; // alpha of the energy spectrum (-1 is E^-2) + Double_t plot = 1; // plot of the enrgy spectrum (1 is E^2) + + Int_t nbins = 16; // number of bins you want to have in psi/phi + + // ------------------------------------------------------------------- + + MBinning binspolx; + binspolx.SetEdges(16, -180, 180); + + // ------------------------------------------------------------------- + + Stat_t n = chain.GetEntries(); + + cout << endl << "Found " << n << " entries." << endl; + + if (n==0) + return; + + // ----------------------------------- + + MBinning bins; + + cout << "Determin min/max..." << endl; + + Double_t min, max; + GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60); + cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl; + + bins.SetEdges(nbins, min, max*1.1); + + TH2D a, a2, a3; + MH::SetBinning(&a, &binspolx, &bins); + MH::SetBinning(&a2, &binspolx, &bins); + MH::SetBinning(&a3, &binspolx, &bins); + + // ----------------------------------- + + GetRange(&chain, "fR", nbins, min, max); + cout << "fR, " << nbins << ": " << min << " ... " << max << " [kpc]" << endl; + + bins.SetEdges(nbins, min, max*1.1); + + TH2D r, r2, r3; + MH::SetBinning(&r, &binspolx, &bins); + MH::SetBinning(&r2, &binspolx, &bins); + MH::SetBinning(&r3, &binspolx, &bins); + + // ----------------------------------- + + MBinning binsx; + binsx.SetEdgesLog(nbinse, lo, hi); + + TH1D prim, h, h2, h3; + MH::SetBinning(&h, &binsx); + MH::SetBinning(&h2, &binsx); + MH::SetBinning(&h3, &binsx); + MH::SetBinning(&prim, &binsx); +// prim.Sumw2(); + h.Sumw2(); + h2.Sumw2(); + h3.Sumw2(); + + // ----------------------------------- + + MPhoton *p = new MPhoton; + chain.SetBranchAddress("MPhoton.", &p); + chain.SetBranchStatus("*", 1); + + chain.GetEntry(0); + if (!p->IsPrimary()) + return; + + Double_t z = p->GetZ(); + Double_t R = MParticle::RofZ(&z); + + cout << "Z = " << z << endl; + cout << "R = " << R << "kpc" << endl; + + // ----------------------------------- + + const Double_t skpc = 3600.*24.*365.*3.258*1e3; // [s/kpc] + GetRange(&chain, "fX", nbins, min, max, 1, kFALSE); + min *= skpc; + max *= skpc; + if (min<1e-2) + min=1e-2; + bins.SetEdgesLog(nbins, min, max); + cout << "dX, " << nbins << ": " << min << " ... " << max << " [s]" << endl; + + TH1D t; + MH::SetBinning(&t, &bins); + t.Sumw2(); + + // ---------------------------------------------------------------------- + + chain.SetBranchAddress("MPhoton.", &p); + chain.SetBranchStatus("*", 1); + + cout << "Filling: " << nbinse << " bins @ " << minE << " - " << hi << "... " << flush; + + Double_t weight = 0; + + Bool_t isprim = kFALSE; + Bool_t ignore = kFALSE; + for (int i=0; iGetEnergy(); + + if (p->IsPrimary()) + { + if (Ep>cutoff) + { + ignore = kTRUE; + continue; + } + ignore = kFALSE; + weight = pow(Ep, alpha); + prim.Fill(Ep, pow(Ep, plot+alpha)); + isprim = kTRUE; + continue; + } + + if (ignore) + continue; + + Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180; + Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180; + + if (isprim) + { + h3.Fill(Ep, pow(Ep,plot) * weight); + r3.Fill(phi, 0.0, weight); + a3.Fill(psi, 0.0, weight); + isprim = kFALSE; + continue; + } + + h2.Fill(Ep, pow(Ep,plot) * weight); + r2.Fill(phi, p->GetR(), weight); + a2.Fill(psi, p->GetTheta()*kRad2Deg*60, weight); + + Double_t tm1 = sqrt(R*R+p->GetR()*p->GetR()); + Double_t tm2 = (R+p->GetX())/tm1 - 1; + t.Fill(p->GetX()<=0 ? 0 : R*tm2*skpc, weight); + } + + // ---------------------------------------------------------------------- + + delete p; + + cout << "Done." << endl; + + Double_t N = prim.GetBinContent(nbinse); + Double_t Nerr = prim.GetBinError(nbinse); + for (int i=1; i<=nbinse; i++) + { + Double_t E = prim.GetXaxis()->GetBinCenter(i); + if (E>minE) + continue; + + Double_t E1 = prim.GetXaxis()->GetBinLowEdge(i+1); + Double_t E0 = prim.GetXaxis()->GetBinLowEdge(i); + + E = sqrt(E1*E0); + + prim.SetBinContent(i, N); + h3.SetBinContent(i, N); + prim.SetBinError(i, Nerr); + h3.SetBinError(i, Nerr); + } + + h.Add(&h2, &h3); + r.Add(&r2, &r3); + a.Add(&a2, &a3); + + if (alpha==-plot) + { + cout << "Observed Energy: " << h.Integral() << "GeV" << endl; + cout << "Emitted Energy: " << prim.Integral() << "GeV" << endl; + cout << "Ratio: " << h.Integral()/prim.Integral() << endl; + } + + cout << "Mean fPhi: " << r.GetMean() << " +- " << r.GetRMS() << endl; + cout << "Mean fPsi: " << a.GetMean() << " +- " << a.GetRMS() << endl; + + gStyle->SetOptStat(10); + gStyle->SetPalette(1, 0); + + TLine line; + + // ---------------------------------------------------------------------- + +// delete gROOT->FindObject("Analysis Arrival"); + TCanvas *c = MH::MakeDefCanvas("Analysis Arrival", ""); + c->Divide(2,2); + + c->cd(1); + r.SetTitle(" Distance from Observer "); + r.GetXaxis()->SetLabelOffset(-0.015); + r.GetXaxis()->SetTitleOffset(1.1); + r.GetXaxis()->SetRangeUser(1e4, 1e9); + r.SetXTitle("\\Phi [\\circ]"); + r.SetYTitle("R [kpc]"); + TVirtualPad &pad = *gPad; + pad.Divide(2,2); + pad.cd(1); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(0); + gPad->SetPhi(90); + TH1 *g = r.DrawCopy("surf1pol"); + pad.cd(2); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(70); + gPad->SetPhi(90); + g->Draw("surf1pol"); + pad.cd(3); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(90); + gPad->SetPhi(90); + g->Draw("surf1pol"); + pad.cd(4); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(20); + gPad->SetPhi(90); + g->Draw("surf1pol"); + + c->cd(2); + a.SetTitle(" Hit Angle for Observer "); + a.GetXaxis()->SetLabelOffset(-0.015); + a.GetXaxis()->SetTitleOffset(1.1); + a.GetXaxis()->SetRangeUser(1e4, 1e9); + a.SetXTitle("\\Psi [\\circ]"); + a.SetYTitle("\\Theta [']"); + TVirtualPad &pad2 = *gPad; + pad2.Divide(2,2); + pad2.cd(1); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(0); + gPad->SetPhi(90); + g = a.DrawCopy("surf1pol"); + pad2.cd(2); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(70); + gPad->SetPhi(90); + g->Draw("surf1pol"); + pad2.cd(3); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(90); + gPad->SetPhi(90); + g->Draw("surf1pol"); + pad2.cd(4); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(20); + gPad->SetPhi(90); + g->Draw("surf1pol"); + + c->cd(3); + gPad->SetLogy(); + g=r.ProjectionY(); + g->SetDirectory(NULL); + TH1 *g2=r2.ProjectionY(); + g->SetMinimum(MH::GetMinimumGT(*g2)/2); + g->SetBit(kCanDelete); + g->SetTitle(" Hit Observers Plain "); + g->GetXaxis()->SetLabelOffset(0); + g->GetXaxis()->SetTitleOffset(1.1); + g->GetYaxis()->SetTitleOffset(1.3); + g->SetXTitle("R [kpc]"); + g->SetYTitle("Counts"); + g->Draw(); + g2->SetDirectory(NULL); + g2->SetBit(kCanDelete); + g2->SetLineColor(kBlue); + g2->Draw("same"); + g=r3.ProjectionY(); + g->SetDirectory(NULL); + g->SetBit(kCanDelete); + g->SetLineColor(kGreen); + g->Draw("same"); + + c->cd(4); + gPad->SetLogy(); + g=a.ProjectionY(); + g->SetMinimum(MH::GetMinimumGT(*g)/2); + g->SetDirectory(NULL); + g->SetBit(kCanDelete); + g->SetTitle("Hit Angle"); + g->GetXaxis()->SetLabelOffset(0); + g->GetXaxis()->SetTitleOffset(1.1); + g->GetYaxis()->SetTitleOffset(1.3); + g->SetXTitle("\\Phi [']"); + g->SetYTitle("Counts"); + g->Draw(); + g=a2.ProjectionY(); + g->SetDirectory(NULL); + g->SetBit(kCanDelete); + g->SetLineColor(kBlue); + g->Draw("same"); + g=a3.ProjectionY(); + g->SetDirectory(NULL); + g->SetBit(kCanDelete); + g->SetLineColor(kGreen); + g->Draw("same"); + + // ---------------------------------------------------------------------- + + // delete gROOT->FindObject("Analysis Photons"); + c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870); + c->Divide(1,2); + + c->cd(1); + gPad->SetLogx(); + gPad->SetLogy(); + h3.SetTitle(Form(" E^{%.1f} Spectra @ z=%.1e (R=%.0fkpc) ", alpha, z, R)); + h3.SetXTitle("E\\gamma [GeV]"); + h3.SetYTitle(Form("E^{%.1f} * Counts", plot)); + h3.GetXaxis()->SetLabelOffset(-0.015); + h3.GetXaxis()->SetTitleOffset(1.1); + h3.SetFillStyle(0); + h3.SetName("PrimPhotons"); + h3.SetMarkerStyle(kFullCircle); + h3.SetMarkerColor(kRed); + h3.SetMarkerSize(0.8); + h3.SetLineColor(kRed); + h3.DrawCopy("C"); //E1 + + TH1D hs = Smear(h, 0.2); + hs.SetLineColor(kGreen); + hs.SetFillStyle(0); + hs.DrawCopy("Csame"); + + h.SetFillStyle(0); + h.SetName("Photons"); + h.SetMarkerStyle(kFullCircle); + h.SetMarkerSize(0.8); + h.DrawCopy("Csame E4"); + prim.SetFillStyle(0); + prim.SetMarkerStyle(kPlus); + prim.SetMarkerColor(kRed); + prim.SetLineColor(kRed); + prim.SetMarkerSize(0.8); + prim.DrawCopy("Csame"); + h2.SetFillStyle(0); + h2.SetName("SecPhotons"); + h2.SetMarkerStyle(kFullCircle); + h2.SetMarkerColor(kBlue); + h2.SetMarkerSize(0.8); + h2.SetLineColor(kBlue); + h2.DrawCopy("Csame E4"); //E2 + + draw1h1426(prim.GetBinContent(2),plot); + + c->cd(2); + gPad->SetLogx(); + TH1D div; + MH::SetBinning(&div, &binsx); + div.Sumw2(); + div.Divide(&h, &prim); + div.SetTitle(" Spectrum / Primary Spectrum "); + div.GetXaxis()->SetLabelOffset(-0.015); + div.GetXaxis()->SetTitleOffset(1.1); + div.SetXTitle("E\\gamma [GeV]"); + div.SetYTitle("Ratio"); + div.SetMarkerStyle(kPlus); + TH1 *gHist = div.DrawCopy("E4"); + + line.SetLineWidth(1); + line.SetLineColor(kGreen); + line.DrawLine(log10(lo), exp(-1), log10(hi), exp(-1)); + line.SetLineColor(kBlue); + line.DrawLine(log10(lo), 1, log10(hi), 1); + + MFit fit("exp(-1)-[1]*log10(x/[0])"); + for (int i=0; iGetEntries(); i++) + { + if (gHist->GetBinContent(i+1)GetBinCenter(i-2), gHist->GetBinCenter(i+2)); + cout << "Fitting from " << gHist->GetBinLowEdge(i-1) << " to "; + cout << gHist->GetBinLowEdge(i+1) << endl; + break; + } + } + fit.SetParameter(0, "t", 750, 10, 1e5); + fit.SetParameter(1, "m", 0.5, 0.1, 10); + fit.FitLog(gHist); + fit.Print(); + fit.DrawCopy("same"); + + cout << "Cutoff: " << setprecision(2) << fit[0]/1e3 << " +- "; + cout << fit.GetParError(0)/1e3 << " TeV" << endl; + + // ---------------------------------------------------------------------- + + c = MH::MakeDefCanvas("Time Analysis", "", 580, 435); + gPad->SetLogx(); + // gPad->SetLogy(); + t.SetName("Times"); + t.SetTitle(" Arrival time distribution "); + t.SetXTitle("\\Delta t [s]"); + t.GetXaxis()->SetLabelOffset(-0.015); + t.GetXaxis()->SetTitleOffset(1.1); + t.DrawCopy("E4"); + line.DrawLine(log10(1), 0, log10(1), t.GetMaximum()*1.05); + line.DrawLine(log10(3600), 0, log10(3600), t.GetMaximum()*1.05); + line.DrawLine(log10(3600*24), 0, log10(3600*24), t.GetMaximum()*1.05); + line.DrawLine(log10(3600*24*365), 0, log10(3600*24*365), t.GetMaximum()*1.05); +} Index: old/WuerzburgSoft/Thomas/mphys/background.txt =================================================================== --- old/WuerzburgSoft/Thomas/mphys/background.txt (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/background.txt (revision 15424) @@ -0,0 +1,101 @@ + -5.155 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.632 10.631 10.631 10.631 10.631 10.631 10.631 10.631 10.630 10.630 10.630 10.630 10.629 10.629 10.628 10.628 10.627 10.626 10.625 10.624 10.623 10.622 10.620 10.619 10.617 10.614 10.612 10.609 10.605 10.601 10.597 10.592 10.586 10.579 10.572 10.563 10.553 10.542 10.530 10.516 10.500 10.482 10.461 10.439 10.413 10.385 10.354 10.319 10.281 10.240 10.194 10.145 10.091 10.033 + -5.089 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.697 10.696 10.696 10.696 10.696 10.696 10.695 10.695 10.695 10.694 10.694 10.693 10.693 10.692 10.691 10.691 10.690 10.688 10.687 10.686 10.684 10.682 10.680 10.677 10.674 10.671 10.667 10.662 10.657 10.651 10.645 10.637 10.629 10.619 10.608 10.595 10.581 10.565 10.547 10.527 10.504 10.479 10.451 10.420 10.385 10.347 10.305 10.260 10.210 10.157 10.099 + -5.022 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.762 10.761 10.761 10.761 10.761 10.761 10.760 10.760 10.760 10.759 10.759 10.759 10.758 10.757 10.757 10.756 10.755 10.754 10.752 10.751 10.749 10.747 10.745 10.742 10.739 10.736 10.732 10.727 10.722 10.716 10.710 10.702 10.694 10.684 10.673 10.660 10.646 10.630 10.612 10.592 10.570 10.544 10.516 10.485 10.451 10.413 10.371 10.325 10.276 10.222 10.165 + -4.956 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.827 10.826 10.826 10.826 10.826 10.826 10.826 10.825 10.825 10.825 10.824 10.824 10.823 10.823 10.822 10.821 10.821 10.820 10.818 10.817 10.816 10.814 10.812 10.810 10.807 10.804 10.801 10.797 10.792 10.787 10.781 10.775 10.767 10.759 10.749 10.738 10.725 10.711 10.695 10.677 10.657 10.635 10.610 10.581 10.550 10.516 10.478 10.436 10.391 10.341 10.288 10.230 + -4.890 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.892 10.891 10.891 10.891 10.891 10.891 10.891 10.891 10.891 10.890 10.890 10.890 10.890 10.889 10.889 10.889 10.888 10.887 10.887 10.886 10.885 10.884 10.883 10.882 10.880 10.878 10.876 10.874 10.872 10.869 10.865 10.861 10.857 10.852 10.846 10.840 10.832 10.823 10.814 10.803 10.790 10.776 10.760 10.742 10.722 10.700 10.675 10.646 10.615 10.581 10.543 10.502 10.456 10.407 10.353 10.296 + -4.823 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.956 10.955 10.955 10.955 10.955 10.955 10.955 10.954 10.954 10.954 10.953 10.953 10.952 10.952 10.951 10.950 10.950 10.949 10.947 10.946 10.945 10.943 10.941 10.939 10.936 10.933 10.930 10.926 10.921 10.916 10.911 10.904 10.896 10.888 10.878 10.867 10.855 10.841 10.825 10.807 10.787 10.764 10.739 10.711 10.680 10.646 10.608 10.567 10.521 10.472 10.419 10.361 + -4.757 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.020 11.019 11.019 11.019 11.019 11.019 11.019 11.018 11.018 11.018 11.017 11.017 11.016 11.016 11.015 11.014 11.014 11.013 11.011 11.010 11.009 11.007 11.005 11.003 11.000 10.997 10.994 10.990 10.985 10.980 10.975 10.968 10.961 10.952 10.942 10.932 10.919 10.905 10.889 10.871 10.851 10.829 10.804 10.776 10.745 10.711 10.673 10.631 10.586 10.537 10.484 10.426 + -4.691 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.084 11.083 11.083 11.083 11.083 11.083 11.083 11.083 11.083 11.083 11.082 11.082 11.082 11.082 11.081 11.081 11.081 11.080 11.079 11.079 11.078 11.077 11.076 11.075 11.074 11.072 11.071 11.069 11.066 11.064 11.061 11.057 11.054 11.049 11.044 11.038 11.032 11.024 11.016 11.006 10.995 10.983 10.969 10.953 10.935 10.915 10.893 10.868 10.840 10.809 10.775 10.737 10.696 10.651 10.602 10.549 10.491 + -4.624 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.147 11.146 11.146 11.146 11.146 11.146 11.146 11.146 11.145 11.145 11.145 11.145 11.144 11.144 11.143 11.143 11.142 11.141 11.140 11.139 11.138 11.137 11.135 11.134 11.132 11.130 11.127 11.124 11.121 11.117 11.113 11.108 11.102 11.095 11.088 11.079 11.070 11.059 11.047 11.033 11.017 10.999 10.979 10.957 10.932 10.904 10.873 10.839 10.802 10.760 10.715 10.666 10.613 10.556 + -4.558 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.210 11.209 11.209 11.209 11.209 11.209 11.209 11.209 11.209 11.209 11.209 11.209 11.209 11.208 11.208 11.208 11.208 11.208 11.207 11.207 11.206 11.206 11.205 11.205 11.204 11.203 11.202 11.201 11.200 11.198 11.196 11.195 11.192 11.190 11.187 11.183 11.180 11.175 11.170 11.165 11.158 11.151 11.142 11.133 11.122 11.110 11.096 11.080 11.062 11.043 11.020 10.995 10.968 10.937 10.903 10.866 10.824 10.779 10.731 10.678 10.621 + -4.492 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.272 11.271 11.271 11.271 11.271 11.271 11.271 11.271 11.271 11.271 11.271 11.271 11.271 11.270 11.270 11.270 11.270 11.269 11.269 11.268 11.268 11.267 11.267 11.266 11.265 11.264 11.263 11.262 11.260 11.259 11.257 11.254 11.252 11.249 11.246 11.242 11.238 11.233 11.227 11.220 11.213 11.205 11.195 11.184 11.172 11.158 11.143 11.125 11.105 11.083 11.058 11.031 11.000 10.966 10.929 10.888 10.843 10.794 10.742 10.685 + -4.425 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.333 11.332 11.332 11.332 11.332 11.332 11.332 11.332 11.331 11.331 11.331 11.331 11.330 11.330 11.329 11.329 11.328 11.327 11.326 11.326 11.324 11.323 11.322 11.320 11.318 11.316 11.313 11.310 11.307 11.303 11.299 11.294 11.288 11.282 11.275 11.266 11.257 11.246 11.234 11.220 11.205 11.187 11.168 11.146 11.121 11.093 11.063 11.029 10.992 10.951 10.907 10.858 10.805 10.749 + -4.359 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.394 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.393 11.392 11.392 11.392 11.392 11.392 11.391 11.391 11.391 11.390 11.390 11.389 11.389 11.388 11.387 11.386 11.385 11.384 11.382 11.381 11.379 11.376 11.374 11.371 11.368 11.364 11.360 11.355 11.349 11.343 11.336 11.327 11.318 11.307 11.295 11.281 11.266 11.249 11.229 11.207 11.183 11.155 11.125 11.091 11.054 11.014 10.969 10.921 10.869 10.812 + -4.293 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.453 11.452 11.452 11.452 11.452 11.452 11.452 11.452 11.452 11.452 11.451 11.451 11.451 11.450 11.450 11.450 11.449 11.449 11.448 11.447 11.446 11.446 11.444 11.443 11.442 11.440 11.438 11.436 11.434 11.431 11.427 11.424 11.419 11.415 11.409 11.403 11.395 11.387 11.378 11.367 11.355 11.342 11.327 11.309 11.290 11.268 11.244 11.217 11.186 11.153 11.116 11.076 11.032 10.983 10.931 10.875 + -4.226 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.511 11.510 11.510 11.510 11.510 11.510 11.510 11.509 11.509 11.509 11.508 11.508 11.508 11.507 11.506 11.506 11.505 11.504 11.503 11.502 11.500 11.498 11.497 11.494 11.492 11.489 11.486 11.482 11.478 11.473 11.468 11.461 11.454 11.446 11.437 11.426 11.415 11.401 11.386 11.369 11.350 11.328 11.304 11.277 11.247 11.214 11.177 11.137 11.093 11.045 10.993 10.937 + -4.160 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.568 11.567 11.567 11.567 11.567 11.567 11.567 11.567 11.566 11.566 11.566 11.565 11.565 11.564 11.564 11.563 11.563 11.562 11.561 11.560 11.559 11.557 11.556 11.554 11.552 11.549 11.546 11.543 11.539 11.535 11.531 11.525 11.519 11.512 11.504 11.495 11.484 11.473 11.459 11.444 11.427 11.408 11.387 11.363 11.336 11.307 11.274 11.237 11.198 11.154 11.106 11.055 10.999 + -4.094 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.623 11.622 11.622 11.622 11.622 11.622 11.622 11.621 11.621 11.621 11.620 11.620 11.619 11.619 11.618 11.617 11.616 11.615 11.614 11.613 11.611 11.609 11.607 11.605 11.602 11.599 11.595 11.591 11.586 11.581 11.575 11.568 11.560 11.551 11.541 11.529 11.516 11.501 11.484 11.466 11.444 11.421 11.394 11.365 11.332 11.296 11.257 11.214 11.167 11.115 11.060 + -4.027 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.677 11.676 11.676 11.676 11.676 11.676 11.676 11.676 11.676 11.676 11.676 11.676 11.676 11.676 11.676 11.675 11.675 11.675 11.675 11.674 11.674 11.674 11.673 11.673 11.672 11.671 11.671 11.670 11.669 11.667 11.666 11.664 11.663 11.661 11.658 11.656 11.652 11.649 11.645 11.640 11.635 11.629 11.622 11.614 11.605 11.595 11.584 11.571 11.556 11.540 11.521 11.500 11.477 11.451 11.422 11.390 11.354 11.315 11.272 11.226 11.175 11.120 + -3.961 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.728 11.727 11.727 11.727 11.727 11.727 11.727 11.727 11.727 11.727 11.727 11.727 11.727 11.727 11.726 11.726 11.726 11.726 11.725 11.725 11.725 11.724 11.724 11.723 11.722 11.722 11.721 11.720 11.719 11.717 11.716 11.714 11.712 11.710 11.707 11.704 11.701 11.697 11.692 11.687 11.681 11.674 11.667 11.658 11.648 11.637 11.624 11.610 11.593 11.575 11.555 11.532 11.506 11.477 11.446 11.411 11.372 11.330 11.283 11.233 11.179 + -3.895 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.776 11.775 11.775 11.775 11.775 11.775 11.775 11.775 11.775 11.774 11.774 11.774 11.773 11.773 11.773 11.772 11.772 11.771 11.770 11.769 11.768 11.767 11.766 11.764 11.763 11.761 11.758 11.756 11.753 11.750 11.746 11.741 11.736 11.730 11.724 11.716 11.708 11.698 11.687 11.675 11.661 11.645 11.627 11.607 11.584 11.559 11.531 11.500 11.465 11.427 11.386 11.340 11.290 11.237 + -3.828 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.821 11.820 11.820 11.820 11.820 11.820 11.820 11.820 11.819 11.819 11.819 11.819 11.818 11.818 11.817 11.817 11.816 11.815 11.815 11.814 11.813 11.811 11.810 11.808 11.806 11.804 11.802 11.799 11.795 11.792 11.787 11.782 11.777 11.771 11.763 11.755 11.746 11.735 11.723 11.709 11.694 11.676 11.656 11.634 11.610 11.582 11.552 11.518 11.481 11.439 11.395 11.346 11.293 + -3.762 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.862 11.861 11.861 11.861 11.861 11.861 11.861 11.860 11.860 11.860 11.859 11.859 11.858 11.858 11.857 11.856 11.855 11.854 11.853 11.851 11.850 11.848 11.846 11.844 11.841 11.838 11.834 11.830 11.825 11.820 11.814 11.807 11.799 11.790 11.779 11.767 11.754 11.739 11.722 11.703 11.682 11.658 11.631 11.601 11.568 11.531 11.491 11.447 11.399 11.347 + -3.696 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.899 11.898 11.898 11.898 11.898 11.898 11.898 11.898 11.898 11.898 11.897 11.897 11.897 11.897 11.896 11.896 11.895 11.895 11.894 11.894 11.893 11.892 11.891 11.890 11.889 11.887 11.885 11.883 11.881 11.878 11.875 11.872 11.868 11.863 11.858 11.852 11.846 11.838 11.829 11.819 11.808 11.795 11.781 11.765 11.746 11.726 11.702 11.676 11.647 11.615 11.580 11.541 11.498 11.451 11.399 + -3.630 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.930 11.929 11.929 11.929 11.929 11.929 11.929 11.929 11.929 11.929 11.928 11.928 11.928 11.928 11.927 11.927 11.926 11.926 11.925 11.924 11.924 11.923 11.921 11.920 11.919 11.917 11.915 11.913 11.911 11.908 11.905 11.901 11.897 11.892 11.886 11.880 11.873 11.864 11.855 11.844 11.832 11.818 11.803 11.785 11.765 11.743 11.718 11.690 11.659 11.625 11.587 11.545 11.499 11.449 + -3.563 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.954 11.953 11.953 11.953 11.953 11.953 11.953 11.952 11.952 11.952 11.951 11.951 11.950 11.950 11.949 11.948 11.947 11.946 11.945 11.944 11.942 11.941 11.939 11.936 11.934 11.931 11.927 11.923 11.919 11.914 11.908 11.901 11.893 11.885 11.875 11.863 11.850 11.835 11.819 11.800 11.779 11.755 11.729 11.699 11.666 11.629 11.589 11.545 11.496 + -3.497 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.970 11.969 11.969 11.969 11.969 11.969 11.968 11.968 11.968 11.967 11.967 11.966 11.965 11.964 11.964 11.963 11.961 11.960 11.958 11.957 11.954 11.952 11.949 11.946 11.943 11.939 11.934 11.928 11.922 11.915 11.907 11.898 11.888 11.876 11.862 11.847 11.829 11.809 11.787 11.762 11.734 11.702 11.667 11.629 11.586 11.540 + -3.431 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.977 11.976 11.976 11.976 11.976 11.976 11.976 11.976 11.975 11.975 11.975 11.974 11.974 11.974 11.973 11.972 11.972 11.971 11.970 11.969 11.968 11.967 11.965 11.963 11.961 11.959 11.956 11.953 11.949 11.945 11.940 11.935 11.929 11.921 11.913 11.904 11.893 11.881 11.867 11.851 11.833 11.812 11.789 11.763 11.733 11.700 11.664 11.623 11.579 + -3.364 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.972 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.971 11.970 11.970 11.970 11.970 11.969 11.969 11.968 11.968 11.967 11.967 11.966 11.965 11.964 11.963 11.962 11.960 11.959 11.957 11.954 11.952 11.949 11.945 11.941 11.937 11.932 11.926 11.919 11.910 11.901 11.890 11.878 11.864 11.848 11.829 11.808 11.784 11.757 11.727 11.693 11.655 11.613 + -3.298 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.952 11.951 11.951 11.951 11.951 11.951 11.951 11.951 11.951 11.951 11.950 11.950 11.950 11.950 11.949 11.949 11.948 11.948 11.947 11.947 11.946 11.945 11.944 11.943 11.941 11.939 11.937 11.935 11.933 11.930 11.926 11.922 11.917 11.912 11.905 11.898 11.889 11.879 11.867 11.853 11.837 11.819 11.798 11.773 11.746 11.715 11.680 11.641 + -3.232 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.915 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.914 11.913 11.913 11.913 11.913 11.912 11.912 11.912 11.911 11.911 11.910 11.910 11.909 11.908 11.907 11.905 11.904 11.902 11.900 11.897 11.894 11.890 11.885 11.880 11.874 11.866 11.857 11.846 11.833 11.818 11.800 11.780 11.756 11.729 11.697 11.662 + -3.165 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.856 11.855 11.855 11.854 11.853 11.852 11.850 11.848 11.845 11.842 11.837 11.831 11.824 11.815 11.804 11.790 11.774 11.755 11.732 11.705 11.674 + -3.099 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.771 11.772 11.772 11.772 11.772 11.772 11.772 11.772 11.772 11.772 11.772 11.773 11.773 11.773 11.773 11.774 11.774 11.774 11.775 11.775 11.776 11.777 11.777 11.778 11.779 11.780 11.781 11.782 11.783 11.784 11.785 11.787 11.788 11.789 11.790 11.790 11.790 11.789 11.787 11.784 11.780 11.773 11.765 11.754 11.740 11.722 11.701 11.675 + -3.033 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.655 11.656 11.656 11.656 11.656 11.656 11.656 11.656 11.656 11.656 11.656 11.656 11.656 11.657 11.657 11.657 11.657 11.658 11.658 11.658 11.659 11.659 11.660 11.660 11.661 11.662 11.663 11.664 11.665 11.667 11.668 11.670 11.672 11.674 11.677 11.679 11.682 11.685 11.689 11.693 11.697 11.701 11.705 11.709 11.713 11.717 11.720 11.722 11.724 11.723 11.721 11.716 11.708 11.697 11.683 11.664 + -2.966 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.502 11.503 11.503 11.503 11.503 11.503 11.503 11.503 11.504 11.504 11.504 11.505 11.505 11.506 11.506 11.507 11.507 11.508 11.509 11.510 11.511 11.513 11.514 11.516 11.518 11.520 11.523 11.526 11.529 11.533 11.537 11.542 11.547 11.553 11.559 11.566 11.573 11.581 11.589 11.598 11.607 11.616 11.625 11.634 11.642 11.649 11.654 11.657 11.657 11.654 11.647 11.636 + -2.900 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.304 11.305 11.305 11.305 11.305 11.305 11.305 11.306 11.306 11.306 11.307 11.307 11.308 11.308 11.309 11.310 11.311 11.312 11.313 11.314 11.316 11.318 11.320 11.322 11.324 11.327 11.331 11.335 11.339 11.344 11.349 11.355 11.362 11.370 11.379 11.388 11.399 11.410 11.423 11.436 11.450 11.466 11.481 11.498 11.514 11.530 11.545 11.559 11.571 11.581 11.588 11.591 11.589 + -2.834 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.053 11.054 11.054 11.054 11.054 11.054 11.055 11.055 11.055 11.056 11.056 11.057 11.057 11.058 11.059 11.060 11.061 11.062 11.064 11.065 11.067 11.069 11.072 11.074 11.077 11.081 11.085 11.090 11.095 11.101 11.108 11.115 11.124 11.134 11.144 11.156 11.170 11.185 11.201 11.219 11.238 11.259 11.281 11.305 11.330 11.355 11.381 11.406 11.431 11.455 11.476 11.494 11.509 11.519 + -2.767 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.739 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.740 10.741 10.741 10.741 10.741 10.741 10.742 10.742 10.743 10.743 10.744 10.744 10.745 10.746 10.747 10.748 10.749 10.750 10.752 10.754 10.756 10.758 10.761 10.764 10.768 10.772 10.777 10.782 10.788 10.795 10.803 10.812 10.822 10.834 10.846 10.861 10.877 10.895 10.915 10.937 10.961 10.987 11.016 11.047 11.080 11.114 11.151 11.188 11.226 11.263 11.300 11.335 11.368 11.397 11.421 + -2.701 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.354 10.355 10.355 10.355 10.355 10.355 10.355 10.355 10.355 10.355 10.355 10.356 10.356 10.356 10.356 10.357 10.357 10.357 10.358 10.358 10.359 10.360 10.360 10.361 10.362 10.363 10.365 10.366 10.368 10.370 10.373 10.375 10.378 10.382 10.386 10.391 10.396 10.402 10.409 10.417 10.426 10.436 10.447 10.460 10.475 10.491 10.510 10.531 10.554 10.580 10.608 10.640 10.675 10.712 10.753 10.797 10.843 10.892 10.943 10.995 11.048 11.101 11.152 11.202 11.247 11.289 + -2.635 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.888 9.889 9.889 9.889 9.889 9.889 9.889 9.889 9.889 9.889 9.889 9.890 9.890 9.890 9.890 9.891 9.891 9.891 9.892 9.893 9.893 9.894 9.895 9.896 9.897 9.898 9.900 9.901 9.903 9.905 9.908 9.911 9.914 9.918 9.922 9.927 9.933 9.940 9.947 9.956 9.965 9.976 9.989 10.003 10.019 10.038 10.058 10.082 10.108 10.137 10.170 10.206 10.246 10.290 10.338 10.390 10.447 10.507 10.570 10.637 10.706 10.777 10.849 10.919 10.988 11.054 11.115 + -2.568 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.343 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.344 9.345 9.345 9.345 9.345 9.346 9.346 9.346 9.347 9.347 9.348 9.348 9.349 9.350 9.351 9.352 9.353 9.355 9.357 9.359 9.361 9.363 9.366 9.370 9.373 9.378 9.383 9.389 9.395 9.403 9.412 9.422 9.433 9.446 9.461 9.478 9.497 9.518 9.543 9.571 9.602 9.637 9.677 9.721 9.770 9.824 9.884 9.949 10.019 10.095 10.175 10.260 10.349 10.441 10.534 10.627 10.719 10.809 10.894 + -2.502 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.781 8.782 8.782 8.782 8.782 8.782 8.782 8.783 8.783 8.783 8.784 8.784 8.785 8.785 8.786 8.787 8.787 8.788 8.790 8.791 8.792 8.794 8.796 8.798 8.801 8.804 8.808 8.812 8.816 8.822 8.828 8.835 8.843 8.852 8.863 8.875 8.890 8.906 8.925 8.947 8.972 9.001 9.034 9.071 9.114 9.163 9.218 9.280 9.349 9.425 9.509 9.600 9.699 9.804 9.914 10.029 10.147 10.266 10.385 10.502 10.615 + -2.436 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.402 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.403 8.404 8.404 8.404 8.404 8.404 8.405 8.405 8.406 8.406 8.407 8.407 8.408 8.409 8.410 8.411 8.412 8.414 8.416 8.418 8.420 8.423 8.426 8.430 8.435 8.441 8.447 8.455 8.465 8.476 8.490 8.507 8.528 8.553 8.583 8.620 8.664 8.717 8.779 8.853 8.939 9.036 9.146 9.266 9.396 9.534 9.679 9.827 9.976 10.125 10.270 + -2.369 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.280 8.279 8.279 8.279 8.279 8.279 8.279 8.279 8.279 8.279 8.279 8.279 8.279 8.280 8.280 8.281 8.283 8.285 8.288 8.292 8.298 8.307 8.319 8.336 8.360 8.394 8.442 8.506 8.591 8.696 8.823 8.969 9.131 9.304 9.485 9.667 9.848 + -2.303 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.231 8.230 8.230 8.230 8.230 8.230 8.230 8.230 8.229 8.229 8.229 8.229 8.228 8.228 8.227 8.227 8.226 8.225 8.225 8.224 8.223 8.222 8.220 8.219 8.217 8.216 8.215 8.213 8.212 8.210 8.210 8.214 8.224 8.245 8.282 8.343 8.436 8.565 8.729 8.920 9.127 9.343 + -2.237 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.172 8.171 8.171 8.171 8.171 8.171 8.170 8.170 8.170 8.169 8.168 8.168 8.167 8.166 8.164 8.161 8.157 8.152 8.146 8.141 8.137 8.135 8.135 8.147 8.182 8.252 8.374 8.549 8.767 + -2.170 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.089 8.090 8.090 8.090 8.090 8.090 8.090 8.090 8.090 8.091 8.091 8.091 8.091 8.092 8.092 8.093 8.093 8.094 8.094 8.095 8.096 8.096 8.097 8.098 8.099 8.100 8.101 8.101 8.098 8.095 8.092 8.088 8.084 8.078 8.068 8.058 8.050 8.046 8.059 8.111 8.228 + -2.104 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.976 7.977 7.977 7.977 7.977 7.977 7.977 7.977 7.977 7.977 7.978 7.978 7.978 7.978 7.979 7.979 7.980 7.980 7.981 7.981 7.982 7.983 7.984 7.985 7.986 7.988 7.989 7.991 7.993 7.996 7.998 8.001 8.004 8.008 8.011 8.014 8.014 8.015 8.015 8.015 8.015 8.014 8.008 7.999 7.989 7.973 7.955 7.937 7.931 + -2.038 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.830 7.831 7.831 7.831 7.831 7.831 7.831 7.831 7.831 7.832 7.832 7.832 7.832 7.833 7.833 7.833 7.834 7.834 7.835 7.836 7.837 7.838 7.839 7.840 7.841 7.843 7.845 7.847 7.849 7.852 7.855 7.858 7.862 7.867 7.872 7.877 7.883 7.889 7.895 7.900 7.904 7.908 7.913 7.919 7.923 7.922 7.919 7.915 7.903 7.885 7.858 7.821 + -1.972 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.648 7.649 7.649 7.649 7.649 7.649 7.649 7.649 7.650 7.650 7.650 7.651 7.651 7.651 7.652 7.652 7.653 7.654 7.655 7.656 7.657 7.658 7.660 7.661 7.663 7.666 7.668 7.671 7.674 7.678 7.682 7.687 7.693 7.699 7.706 7.714 7.722 7.731 7.741 7.749 7.757 7.767 7.777 7.788 7.798 7.804 7.807 7.810 7.804 7.793 7.769 7.731 + -1.905 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.428 7.429 7.429 7.429 7.429 7.429 7.429 7.429 7.430 7.430 7.430 7.430 7.431 7.431 7.432 7.432 7.433 7.433 7.434 7.435 7.436 7.438 7.439 7.441 7.442 7.445 7.447 7.450 7.453 7.457 7.461 7.466 7.471 7.477 7.484 7.492 7.500 7.510 7.521 7.533 7.546 7.558 7.571 7.585 7.600 7.617 7.634 7.646 7.657 7.668 7.670 7.665 7.648 7.612 + -1.839 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.171 7.172 7.172 7.172 7.172 7.172 7.172 7.173 7.173 7.173 7.174 7.174 7.175 7.175 7.176 7.177 7.178 7.179 7.180 7.181 7.183 7.185 7.187 7.190 7.192 7.195 7.199 7.203 7.208 7.213 7.220 7.227 7.234 7.243 7.254 7.265 7.278 7.293 7.308 7.324 7.340 7.359 7.380 7.402 7.426 7.445 7.464 7.482 7.492 7.495 7.485 7.453 + -1.773 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.881 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.882 6.883 6.883 6.883 6.883 6.884 6.884 6.884 6.885 6.885 6.886 6.886 6.887 6.888 6.889 6.890 6.891 6.892 6.894 6.896 6.898 6.901 6.904 6.907 6.911 6.915 6.920 6.926 6.932 6.939 6.948 6.957 6.968 6.981 6.995 7.010 7.028 7.046 7.065 7.087 7.112 7.139 7.168 7.194 7.220 7.247 7.266 7.277 7.274 7.247 + -1.706 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.575 6.576 6.576 6.576 6.576 6.576 6.576 6.577 6.577 6.577 6.578 6.578 6.579 6.579 6.580 6.581 6.581 6.582 6.583 6.585 6.586 6.588 6.590 6.592 6.595 6.598 6.601 6.605 6.610 6.615 6.621 6.628 6.636 6.645 6.655 6.667 6.681 6.696 6.713 6.732 6.752 6.775 6.802 6.831 6.864 6.894 6.926 6.959 6.986 7.006 7.010 6.988 + -1.640 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.280 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.281 6.282 6.282 6.282 6.282 6.282 6.283 6.283 6.283 6.284 6.284 6.285 6.285 6.286 6.287 6.288 6.289 6.290 6.292 6.294 6.296 6.298 6.301 6.304 6.307 6.311 6.316 6.321 6.327 6.335 6.343 6.352 6.364 6.376 6.390 6.406 6.423 6.444 6.468 6.495 6.526 6.556 6.590 6.626 6.657 6.682 6.693 6.676 + -1.574 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.038 6.039 6.039 6.039 6.039 6.039 6.040 6.040 6.040 6.041 6.041 6.042 6.043 6.043 6.044 6.045 6.047 6.048 6.049 6.051 6.053 6.056 6.059 6.062 6.066 6.070 6.075 6.081 6.088 6.097 6.106 6.115 6.126 6.139 6.155 6.174 6.197 6.219 6.246 6.275 6.303 6.327 6.338 6.323 + -1.507 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.864 5.865 5.865 5.865 5.865 5.865 5.865 5.866 5.866 5.866 5.866 5.867 5.867 5.868 5.869 5.869 5.870 5.871 5.872 5.874 5.875 5.877 5.879 5.882 5.885 5.889 5.893 5.897 5.900 5.905 5.910 5.917 5.926 5.937 5.946 5.958 5.972 5.985 5.996 5.997 5.976 + -1.441 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.741 5.742 5.742 5.742 5.742 5.742 5.742 5.742 5.743 5.743 5.743 5.744 5.744 5.745 5.746 5.746 5.748 5.749 5.750 5.752 5.754 5.754 5.755 5.756 5.758 5.761 5.764 5.764 5.764 5.765 5.763 5.757 5.742 5.710 + -1.375 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.642 5.643 5.643 5.643 5.642 5.641 5.640 5.640 5.640 5.640 5.637 5.633 5.628 5.619 5.605 5.581 5.541 + -1.308 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.553 5.552 5.552 5.552 5.552 5.552 5.552 5.552 5.552 5.552 5.552 5.552 5.552 5.552 5.552 5.551 5.551 5.551 5.551 5.551 5.550 5.550 5.550 5.549 5.549 5.549 5.548 5.548 5.547 5.546 5.545 5.543 5.541 5.539 5.536 5.535 5.532 5.527 5.521 5.515 5.503 5.486 5.459 5.417 + -1.242 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.461 5.460 5.460 5.460 5.460 5.460 5.459 5.459 5.459 5.458 5.458 5.457 5.457 5.456 5.455 5.454 5.452 5.450 5.446 5.443 5.439 5.436 5.431 5.424 5.416 5.407 5.393 5.374 5.345 5.302 + -1.176 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.360 5.359 5.358 5.357 5.355 5.351 5.348 5.344 5.340 5.334 5.325 5.315 5.303 5.287 5.265 5.235 5.190 + -1.109 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.240 5.241 5.241 5.241 5.241 5.241 5.242 5.242 5.242 5.243 5.243 5.244 5.244 5.245 5.245 5.246 5.247 5.247 5.248 5.249 5.250 5.250 5.248 5.246 5.244 5.241 5.238 5.233 5.225 5.214 5.201 5.183 5.160 5.127 5.082 + -1.043 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.098 5.099 5.099 5.099 5.099 5.099 5.099 5.099 5.099 5.099 5.099 5.100 5.100 5.100 5.101 5.101 5.101 5.102 5.102 5.103 5.104 5.105 5.105 5.106 5.107 5.109 5.110 5.112 5.114 5.116 5.118 5.120 5.122 5.124 5.125 5.125 5.125 5.124 5.124 5.122 5.115 5.107 5.095 5.078 5.055 5.022 4.976 + -0.977 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.940 4.941 4.941 4.941 4.941 4.941 4.942 4.942 4.942 4.943 4.943 4.944 4.944 4.945 4.946 4.947 4.948 4.949 4.950 4.951 4.953 4.955 4.957 4.959 4.961 4.965 4.968 4.972 4.976 4.980 4.982 4.984 4.987 4.989 4.992 4.994 4.991 4.986 4.978 4.964 4.943 4.912 4.867 + -0.910 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.775 4.776 4.776 4.776 4.776 4.776 4.776 4.776 4.776 4.776 4.777 4.777 4.777 4.778 4.778 4.779 4.779 4.779 4.780 4.781 4.782 4.783 4.783 4.785 4.787 4.789 4.791 4.793 4.796 4.798 4.802 4.807 4.812 4.816 4.821 4.824 4.829 4.833 4.838 4.843 4.849 4.850 4.849 4.846 4.836 4.819 4.791 4.747 + -0.844 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.624 4.625 4.625 4.625 4.625 4.625 4.625 4.625 4.625 4.625 4.625 4.625 4.625 4.625 4.626 4.626 4.626 4.627 4.627 4.628 4.629 4.629 4.630 4.631 4.632 4.632 4.633 4.634 4.634 4.635 4.638 4.637 4.640 4.643 4.646 4.651 4.655 4.659 4.663 4.668 4.672 4.678 4.685 4.692 4.695 4.697 4.697 4.691 4.679 4.653 4.612 + -0.778 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.505 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.506 4.507 4.507 4.507 4.507 4.507 4.507 4.507 4.507 4.507 4.506 4.506 4.505 4.505 4.505 4.505 4.505 4.505 4.504 4.505 4.505 4.505 4.506 4.506 4.506 4.508 4.508 4.510 4.513 4.513 4.513 4.513 4.515 4.518 4.522 4.527 4.532 4.534 4.537 4.539 4.534 4.524 4.500 4.459 + -0.711 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.432 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.431 4.430 4.429 4.429 4.428 4.426 4.424 4.425 4.425 4.422 4.419 4.419 4.417 4.416 4.415 4.410 4.408 4.406 4.404 4.401 4.399 4.394 4.394 4.389 4.386 4.387 4.384 4.383 4.382 4.373 4.362 4.338 4.295 + -0.645 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.393 4.392 4.392 4.392 4.392 4.391 4.390 4.390 4.389 4.388 4.387 4.386 4.386 4.386 4.387 4.387 4.386 4.385 4.382 4.378 4.377 4.376 4.375 4.368 4.368 4.367 4.363 4.359 4.356 4.351 4.346 4.341 4.333 4.322 4.313 4.304 4.294 4.285 4.278 4.268 4.255 4.243 4.227 4.209 4.182 4.133 + -0.579 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.370 4.369 4.369 4.369 4.369 4.369 4.369 4.369 4.368 4.368 4.368 4.368 4.368 4.368 4.368 4.368 4.368 4.368 4.368 4.367 4.366 4.364 4.361 4.359 4.359 4.360 4.358 4.355 4.349 4.351 4.349 4.339 4.341 4.333 4.329 4.319 4.313 4.303 4.294 4.281 4.272 4.257 4.244 4.226 4.209 4.189 4.167 4.146 4.118 4.087 4.051 3.997 + -0.512 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.356 4.355 4.355 4.355 4.355 4.355 4.355 4.355 4.355 4.354 4.354 4.354 4.354 4.353 4.353 4.353 4.352 4.352 4.352 4.351 4.351 4.351 4.351 4.350 4.350 4.349 4.349 4.348 4.346 4.345 4.343 4.342 4.339 4.334 4.325 4.325 4.322 4.315 4.312 4.308 4.297 4.293 4.287 4.273 4.261 4.247 4.227 4.212 4.189 4.172 4.139 4.117 4.083 4.046 4.008 3.955 3.894 + -0.446 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.321 4.320 4.320 4.320 4.320 4.320 4.320 4.320 4.320 4.320 4.319 4.319 4.319 4.318 4.317 4.317 4.316 4.314 4.313 4.312 4.312 4.311 4.310 4.310 4.310 4.310 4.310 4.308 4.303 4.298 4.296 4.293 4.285 4.281 4.275 4.266 4.258 4.247 4.236 4.222 4.207 4.192 4.169 4.141 4.112 4.084 4.038 3.997 3.952 3.888 3.827 + -0.380 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.272 4.271 4.271 4.270 4.270 4.269 4.269 4.268 4.266 4.265 4.264 4.263 4.261 4.259 4.255 4.251 4.252 4.250 4.247 4.244 4.235 4.233 4.227 4.222 4.212 4.202 4.188 4.176 4.157 4.140 4.119 4.091 4.056 4.019 3.968 3.919 3.857 3.781 + -0.314 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.212 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.211 4.210 4.210 4.209 4.208 4.207 4.205 4.204 4.204 4.204 4.201 4.200 4.200 4.195 4.192 4.192 4.190 4.188 4.184 4.178 4.170 4.166 4.161 4.152 4.144 4.132 4.116 4.100 4.083 4.055 4.027 3.990 3.941 3.892 3.829 3.752 + -0.247 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.127 4.126 4.126 4.126 4.126 4.126 4.126 4.126 4.126 4.126 4.126 4.126 4.126 4.127 4.127 4.127 4.127 4.127 4.126 4.126 4.125 4.125 4.126 4.126 4.125 4.123 4.122 4.123 4.123 4.121 4.119 4.118 4.116 4.116 4.116 4.112 4.111 4.108 4.102 4.094 4.089 4.082 4.072 4.063 4.051 4.038 4.022 3.998 3.972 3.941 3.903 3.857 3.794 3.724 + -0.181 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.021 4.020 4.020 4.020 4.020 4.020 4.018 4.016 4.015 4.016 4.016 4.015 4.013 4.010 4.006 4.005 4.002 4.000 3.997 3.992 3.988 3.984 3.974 3.962 3.954 3.941 3.926 3.906 3.882 3.851 3.817 3.775 3.722 3.652 + -0.115 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.890 3.889 3.889 3.889 3.889 3.889 3.888 3.888 3.887 3.887 3.888 3.888 3.889 3.890 3.890 3.889 3.887 3.886 3.886 3.884 3.883 3.885 3.884 3.884 3.878 3.880 3.880 3.876 3.874 3.870 3.865 3.861 3.856 3.848 3.839 3.829 3.819 3.800 3.781 3.764 3.732 3.697 3.656 3.600 3.535 + -0.048 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.749 3.748 3.748 3.748 3.748 3.748 3.748 3.748 3.748 3.748 3.748 3.748 3.748 3.748 3.748 3.747 3.747 3.747 3.747 3.747 3.747 3.747 3.748 3.748 3.749 3.749 3.749 3.748 3.747 3.747 3.748 3.747 3.745 3.744 3.744 3.742 3.741 3.734 3.739 3.738 3.738 3.733 3.737 3.733 3.730 3.729 3.724 3.721 3.717 3.704 3.694 3.684 3.674 3.656 3.636 3.604 3.574 3.529 3.473 3.408 + 0.018 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.592 3.593 3.593 3.593 3.593 3.593 3.593 3.593 3.592 3.592 3.592 3.593 3.592 3.592 3.592 3.595 3.595 3.593 3.595 3.593 3.590 3.584 3.580 3.574 3.565 3.568 3.561 3.555 3.545 3.539 3.518 3.496 3.472 3.442 3.399 3.343 3.278 + 0.084 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.416 3.417 3.417 3.417 3.417 3.418 3.418 3.418 3.419 3.419 3.420 3.420 3.421 3.422 3.422 3.423 3.424 3.424 3.425 3.425 3.425 3.425 3.424 3.423 3.422 3.420 3.419 3.416 3.417 3.416 3.415 3.415 3.415 3.415 3.419 3.416 3.417 3.416 3.414 3.406 3.413 3.415 3.412 3.409 3.401 3.394 3.380 3.360 3.340 3.325 3.290 3.251 3.198 3.134 + 0.151 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.233 3.234 3.234 3.234 3.234 3.234 3.234 3.234 3.234 3.234 3.234 3.234 3.234 3.235 3.235 3.235 3.235 3.235 3.236 3.236 3.236 3.237 3.237 3.237 3.238 3.238 3.239 3.238 3.238 3.236 3.235 3.233 3.231 3.231 3.233 3.235 3.237 3.236 3.235 3.230 3.227 3.226 3.228 3.236 3.232 3.241 3.241 3.240 3.235 3.242 3.239 3.239 3.238 3.239 3.229 3.222 3.212 3.192 3.176 3.156 3.126 3.095 3.041 2.979 + 0.217 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.033 3.032 3.032 3.032 3.032 3.032 3.032 3.032 3.032 3.032 3.032 3.031 3.032 3.032 3.032 3.032 3.032 3.032 3.032 3.032 3.031 3.028 3.026 3.026 3.029 3.033 3.034 3.035 3.037 3.041 3.042 3.047 3.052 3.049 3.056 3.051 3.050 3.035 3.063 3.063 3.055 3.061 3.048 3.036 3.033 3.027 2.997 2.977 2.943 2.889 2.827 + 0.283 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.818 2.819 2.819 2.819 2.819 2.819 2.819 2.819 2.819 2.820 2.820 2.820 2.820 2.821 2.821 2.821 2.822 2.822 2.823 2.823 2.824 2.825 2.826 2.827 2.827 2.826 2.825 2.824 2.823 2.822 2.820 2.818 2.818 2.821 2.825 2.828 2.832 2.829 2.825 2.827 2.838 2.843 2.847 2.852 2.854 2.857 2.847 2.841 2.850 2.858 2.866 2.871 2.863 2.847 2.873 2.851 2.868 2.861 2.838 2.818 2.782 2.742 2.687 + 0.350 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.611 2.612 2.612 2.612 2.612 2.612 2.612 2.612 2.612 2.612 2.612 2.613 2.613 2.613 2.613 2.614 2.614 2.614 2.616 2.617 2.618 2.619 2.620 2.622 2.622 2.624 2.623 2.623 2.619 2.615 2.622 2.628 2.633 2.622 2.633 2.626 2.631 2.649 2.659 2.658 2.656 2.669 2.662 2.659 2.653 2.649 2.635 2.616 2.577 2.502 + 0.416 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.383 2.384 2.384 2.384 2.384 2.384 2.384 2.384 2.384 2.384 2.385 2.385 2.385 2.386 2.386 2.387 2.387 2.388 2.388 2.389 2.390 2.392 2.392 2.391 2.390 2.388 2.385 2.386 2.388 2.389 2.390 2.395 2.400 2.406 2.413 2.418 2.425 2.432 2.430 2.432 2.429 2.427 2.432 2.444 2.445 2.449 2.458 2.451 2.454 2.446 2.451 2.449 2.422 2.387 2.336 + 0.482 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.218 2.219 2.219 2.219 2.219 2.219 2.219 2.219 2.219 2.218 2.218 2.217 2.217 2.216 2.216 2.216 2.215 2.214 2.214 2.215 2.220 2.223 2.221 2.221 2.224 2.224 2.228 2.219 2.224 2.225 2.217 2.230 2.232 2.240 2.247 2.244 2.247 2.256 2.244 2.253 2.246 2.241 2.229 2.224 2.216 2.173 2.135 + 0.549 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.058 2.059 2.059 2.059 2.059 2.059 2.059 2.060 2.060 2.061 2.062 2.063 2.064 2.065 2.066 2.066 2.067 2.067 2.066 2.063 2.061 2.059 2.055 2.050 2.046 2.048 2.056 2.068 2.072 2.074 2.081 2.090 2.086 2.086 2.098 2.098 2.094 2.097 2.090 2.096 2.103 2.102 2.113 2.113 2.109 2.089 2.074 2.067 2.025 1.976 + 0.615 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.903 1.902 1.902 1.902 1.902 1.902 1.902 1.902 1.902 1.901 1.901 1.901 1.901 1.900 1.900 1.900 1.899 1.899 1.898 1.898 1.898 1.898 1.899 1.899 1.900 1.902 1.903 1.905 1.907 1.909 1.911 1.915 1.918 1.918 1.918 1.917 1.913 1.910 1.909 1.915 1.922 1.923 1.937 1.938 1.937 1.935 1.923 1.926 1.945 1.964 1.965 1.976 1.980 1.972 1.973 1.992 1.990 1.977 1.978 1.960 1.950 1.914 1.873 + 0.681 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.720 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.719 1.720 1.720 1.721 1.721 1.723 1.725 1.727 1.730 1.733 1.737 1.741 1.748 1.757 1.767 1.777 1.789 1.798 1.796 1.788 1.783 1.795 1.795 1.798 1.816 1.816 1.836 1.814 1.817 1.853 1.847 1.870 1.868 1.855 1.874 1.840 1.844 1.811 1.767 + 0.748 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.541 1.542 1.542 1.542 1.542 1.542 1.542 1.543 1.543 1.543 1.543 1.544 1.544 1.545 1.545 1.546 1.547 1.548 1.548 1.546 1.540 1.534 1.527 1.520 1.512 1.501 1.490 1.498 1.515 1.536 1.553 1.561 1.568 1.579 1.589 1.590 1.587 1.596 1.637 1.673 1.655 1.664 1.666 1.673 1.704 1.723 1.691 1.723 1.744 1.745 1.726 1.732 1.705 1.665 + 0.814 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.196 1.197 1.197 1.197 1.197 1.197 1.197 1.197 1.198 1.198 1.198 1.198 1.199 1.199 1.200 1.200 1.201 1.201 1.202 1.203 1.204 1.205 1.207 1.208 1.210 1.212 1.215 1.217 1.221 1.224 1.228 1.233 1.238 1.246 1.255 1.266 1.280 1.295 1.308 1.312 1.315 1.320 1.338 1.355 1.363 1.372 1.344 1.305 1.369 1.392 1.433 1.454 1.478 1.519 1.515 1.518 1.541 1.561 1.588 1.598 1.568 1.586 1.567 1.537 + 0.880 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.828 0.829 0.829 0.829 0.829 0.829 0.829 0.830 0.830 0.830 0.830 0.831 0.831 0.832 0.833 0.833 0.834 0.835 0.836 0.838 0.839 0.841 0.843 0.845 0.847 0.850 0.853 0.858 0.866 0.874 0.882 0.889 0.898 0.905 0.908 0.914 0.929 0.949 0.963 0.990 1.028 1.090 1.140 1.158 1.144 1.112 1.168 1.234 1.274 1.339 1.403 1.385 1.427 1.395 1.417 1.408 1.383 + 0.947 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.526 0.527 0.527 0.527 0.527 0.527 0.528 0.528 0.528 0.529 0.529 0.530 0.530 0.531 0.532 0.532 0.533 0.534 0.534 0.535 0.536 0.538 0.539 0.541 0.544 0.545 0.546 0.546 0.546 0.545 0.545 0.550 0.558 0.552 0.530 0.517 0.502 0.548 0.596 0.625 0.647 0.682 0.728 0.772 0.783 0.855 0.967 0.957 1.018 1.100 1.163 1.207 1.213 1.216 1.195 + 1.013 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.120 -0.119 -0.119 -0.119 -0.119 -0.119 -0.118 -0.118 -0.118 -0.117 -0.117 -0.116 -0.116 -0.115 -0.114 -0.114 -0.113 -0.111 -0.110 -0.109 -0.107 -0.105 -0.103 -0.100 -0.097 -0.093 -0.088 -0.083 -0.077 -0.070 -0.063 -0.038 -0.025 -0.011 0.005 0.026 0.048 0.071 0.097 0.124 0.153 0.174 0.191 0.208 0.236 0.283 0.329 0.383 0.418 0.361 0.445 0.522 0.636 0.753 0.799 0.820 0.947 0.956 + 1.079 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.917 -0.916 -0.916 -0.916 -0.916 -0.916 -0.916 -0.916 -0.916 -0.915 -0.915 -0.915 -0.914 -0.914 -0.914 -0.913 -0.913 -0.912 -0.911 -0.910 -0.909 -0.908 -0.907 -0.906 -0.904 -0.902 -0.900 -0.897 -0.895 -0.891 -0.888 -0.883 -0.877 -0.871 -0.864 -0.855 -0.846 -0.836 -0.824 -0.812 -0.798 -0.782 -0.765 -0.744 -0.722 -0.697 -0.670 -0.642 -0.611 -0.571 -0.530 -0.489 -0.444 -0.397 -0.349 -0.302 -0.259 -0.212 -0.123 -0.072 -0.026 0.027 0.173 0.311 0.423 + 1.146 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 + 1.212 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 + 1.278 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 + 1.344 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 + 1.411 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 + 1.477 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 -40.000 Index: old/WuerzburgSoft/Thomas/mphys/energyloss.C =================================================================== --- old/WuerzburgSoft/Thomas/mphys/energyloss.C (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/energyloss.C (revision 15424) @@ -0,0 +1,530 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "mbase/MParContainer.h" +#include "mphys/MPhoton.h" +#include "mphys/MElectron.h" +#include "mphys/MPairProduction.h" +#include "mhist/MBinning.h" +#include "mhist/MH.h" + +// 2.96 +// 2.87 +// 2.73 + +// ================================================================ +Double_t DiSum(Double_t *x, Double_t *k=NULL) +{ + Double_t t = x[0]; + + Double_t disum = t; + Double_t add = 0; + + Double_t eps = fabs(t*1e-1); + + Int_t n = 2; + Double_t pow = t*t; // t^2 + + do + { + add = pow/n/n; + + pow *= t; // pow = t^n + n++; + + disum += add; + + } while (fabs(add)>eps); + + return disum; +} + +Double_t F(Double_t *x, Double_t *k=NULL) +{ + Double_t o = x[0]; + Double_t s = -2.*o; + +// if (o<1e-10) +// return 2.125; //-3./8.; + + return -o/4. + (9./4. + 1./o + o/2.) * log(1.+2.*o) + 1./8.*(1.+2.*o) + MElectron::Li2(&s); //- 3./8. +} + +void rkck(Double_t y[], Double_t dydx[], int n, Double_t x, Double_t h, Double_t yout[], + Double_t yerr[], void (*derivs)(Double_t, Double_t[], Double_t[])) +{ + /* + * --------------------------------------------------------- + * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method + * --------------------------------------------------------- + */ + const Double_t a2 = 0.2; + const Double_t a3 = 0.3; + const Double_t a4 = 0.6; + const Double_t a5 = 1.0; + const Double_t a6 = 0.875; + const Double_t b21 = 0.2; + const Double_t b31 = 3.0/40.0; + const Double_t b32 = 9.0/40.0; + const Double_t b41 = 0.3; + const Double_t b42 = -0.9; + const Double_t b43 = 1.2; + const Double_t b51 = -11.0/54.0; + const Double_t b52 = 2.5; + const Double_t b53 = -70.0/27.0; + const Double_t b54 = 35.0/27.0; + const Double_t b61 = 1631.0/55296.0; + const Double_t b62 = 175.0/512.0; + const Double_t b63 = 575.0/13824.0; + const Double_t b64 = 44275.0/110592.0; + const Double_t b65 = 253.0/4096.0; + const Double_t c1 = 37.0/378.0; + const Double_t c3 = 250.0/621.0; + const Double_t c4 = 125.0/594.0; + const Double_t c6 = 512.0/1771.0; + const Double_t dc5 = -277.00/14336.0; + + const Double_t dc1 = c1-2825.0/27648.0; + const Double_t dc3 = c3-18575.0/48384.0; + const Double_t dc4 = c4-13525.0/55296.0; + const Double_t dc6 = c6-0.25; + + Double_t ak2[n]; + Double_t ak3[n]; + Double_t ak4[n]; + Double_t ak5[n]; + Double_t ak6[n]; + Double_t ytemp[n]; + + for (int i=0; ib ? a : b; +} + +void SolvEq(Double_t y[], Double_t dydx[], int n, Double_t *x, Double_t htry, Double_t eps, + Double_t yscal[], Double_t *hdid, Double_t *hnext, + void (*derivs)(Double_t, Double_t[], Double_t[])) // rkqs +{ + /* + * --------------------------------------------------------- + * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method + * --------------------------------------------------------- + */ + const Double_t SAFETY = 0.9; + const Double_t PGROW = -0.2; + const Double_t PSHRNK = -0.25; + const Double_t ERRCON = 1.89e-4; + + Double_t yerr[n]; + Double_t ytemp[n]; + + Double_t h = htry; + Double_t errmax; + while (1) + { + rkck(y, dydx, n, *x, h, ytemp, yerr, derivs); + + errmax=0.0; + + for (int i=0; i= 0.0 ? FMAX(htemp,0.1*h) : FMIN(htemp,0.1*h)); + + Double_t xnew= (*x) + h; + + if (xnew != *x) + continue; + + cout << "stepsize underflow in rkqs" << endl; + break; + } + + *hnext = errmax>ERRCON ? SAFETY*h*pow(errmax,PGROW) : 5.0*h; + + *x += (*hdid=h); + + for (int i=0; i=0 ? MParticle::ZofR(&R) : 0; + + Double_t e = 1.602176462e-19; // [C] + Double_t T = 2.96*(z+1); // [K] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + Double_t kBT = kB*T; // [GeV] + Double_t alpha = 1./137; // [|] + Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + Double_t E0 = 511e-6; // [GeV] + + + Double_t k = TMath::Pi() * alpha * kBT; + + Double_t ln = 4.*kBT/E0/E0; + + return -1./3/h* k*k * (log(ln*E)-1.9805); +} + +void dEdt(Double_t t, Double_t E[], Double_t dedt[]) +{ + dedt[0] = dEdt(E[0], t); + //cout << t << "\t" << E[0] << "\t" << dedt[0] << endl; +} + +void DrawDevelopmentHiLim(Double_t E0, Double_t z, Option_t *opt="") +{ + Double_t t = 0; + Double_t E[1] = { E0 }; + Double_t yscal[1] = { 1 }; + + Double_t dedt[1]; + + Double_t eps;// = 1e5; + Double_t step = 5e6; + + Double_t hdid; + Double_t hnext; + + cout << "Start: " << dedt[0] << endl; + + Int_t n = 15000; + Double_t tres[n]; + Double_t Eres[n]; + int i; + for (i=0; i1.5e14 || E[0]<5e6) + break; +// cout << tres[i] << "\t" << Eres[i] << "\t(" << step << ")" << endl; + } + + cout << i << endl; + + TGraph grp(iDraw(opt); + +} + +Double_t EnergyLossRateLoLim(Double_t *x, Double_t *k=NULL) +{ + Double_t t = x[0]; + Double_t E = k[0]; + Double_t t0 = k[1]; + + Double_t c = 299792458; // [m/s] + Double_t ly = 3600.*24.*365.*c; // [m/ly] + Double_t pc = 1./3.258; // [pc/ly] + + Double_t r = t * c / ly * pc / 1000; // [kpc] + Double_t R = MParticle::RofZ(&k[2]) - r; + Double_t z = k[2]>=0 ? MParticle::ZofR(&R) : 0; + + Double_t T = 2.96*(z+1); // [K] + // Double_t alpha = 1./137; // [1] + Double_t sigma = 6.653e-29; // [m^2] + Double_t E0 = 511e-6; // [GeV] + Double_t e = 1.602176462e-19; // [C] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + Double_t hc = h*c; // [GeVm] + Double_t pi = TMath::Pi(); // [1] + Double_t khc = pi*kB/hc; // [1 / K m] + Double_t a = 8./15 * pi * khc*khc*khc*khc * hc; // [Gev / K^4 / m^3] + Double_t konst = 4./3 * sigma * a * T*T*T*T * c / (E0* E0); // [1 / GeV s] + + Double_t ret = 1./(konst*(t-t0) + 1./E); + + return ret; +} + +void DrawDevelopmentLoLim(Double_t t0, Double_t E0, Double_t z=-1, Option_t *opt="") +{ + // 8.7 + + Double_t val[] = { E0, t0, z }; + Double_t t = 1.5e14; + while (EnergyLossRateLoLim(&t, val)<1e4) + t -= 0.01e14; + + TF1 *f1=new TF1("LoLim", EnergyLossRateLoLim, t, 1.5e14, 2); + f1->SetParameter(0, E0); + f1->SetParameter(1, t0); + + f1->Draw(opt); + f1->SetBit(kCanDelete); +} + +// +// (3) Energy loss rate of electrons and 'high energy particle' +// +Double_t DrawDevelopment(Double_t E, Double_t z, Option_t *opt="", TH1 *hist=NULL) +{ + Double_t ly = 3600.*24.*365.; // [s/ly] + Double_t pc = 1./3.258; // [pc/ly] + + TGraph *graph = new TGraph; + graph->SetPoint(0, 0, E); + graph->SetMaximum(E*3); // *MENU* + + cout << "------ " << endl; + + static TRandom rand; + + Double_t x = 0; + for (int i=1; i<10; i++) + { + Double_t l = rand.Exp(MElectron::InteractionLength(&E, &z)); + + if (z>=0) + { + Double_t r = MParticle::RofZ(&z) - l; + + cout << " " << i << ". R=" << MParticle::RofZ(&z) << " l=" << l << " z=" << MParticle::ZofR(&r) << endl; + + z = r>=0 ? MParticle::ZofR(&r) : 0; + + if (z==0) + cout << "z<0" << endl; + } + + x += l; + + Double_t t = x/pc*ly*1000; + graph->SetPoint(i*2-1, t, E); + + Double_t e1 = MElectron::GetEnergyLoss(E, z<0?0:z); + + E -= e1; + + if (hist) + hist->Fill(e1); + + cout << " Ep=" << e1 << flush; + + graph->SetPoint(i*2, t, E); + } + + graph->SetMinimum(E/3); // *MENU* + graph->Draw(opt); + graph->SetBit(kCanDelete); + + //if (E<31500) + cout << "t=" << x*ly/pc*1000 << "\tE=" << E << "\tz=" << z << endl; + + return E; +} + +void EnergyLossRate() +{ + if (gPad) + gPad->Clear(); + + Double_t E = 1.5e9; // [GeV] + Double_t z = 0.03; + + MBinning bins; + bins.SetEdgesLog(18, 0.1, 1e9); + + TH1D hist; + hist.SetName("Phot"); + hist.SetTitle("Photons from inverse Compton"); + MH::SetBinning(&hist, &bins); + + cout << "Working..." << flush; + + for (int i=0; i<50; i++) + { + cout << i << "." << flush; + DrawDevelopment(E, z, i?"L":"AL", &hist); + } + + //DrawDevelopmentLoLim(2e14, 1.64e2, "Lsame"); // seen + DrawDevelopmentLoLim(1.78e14, 280, z, "Lsame"); // seen + DrawDevelopmentHiLim(E, z, "L"); + gPad->SetLogy(); + + new TCanvas("Photons", "Photons created in inverse Compton"); + hist.DrawCopy(); + gPad->SetLogx(); + + cout << "...done." << endl; +} + +void DrawRZ() +{ + new TCanvas("RZ", "r and z"); + + TF1 f1("ZofR", MParticle::ZofR, 0, 7.1e6, 0); + TF1 f2("RofZ", MParticle::RofZ, 0, 5, 0); + + gPad->Divide(2,2); + + gPad->GetVirtCanvas()->cd(1); + TH1 *h = f1.DrawCopy()->GetHistogram(); + + h->SetTitle("z(r)"); + h->SetXTitle("r [kpc]"); + h->SetYTitle("z"); + + gPad->Modified(); + gPad->Update(); + + gPad->GetVirtCanvas()->cd(2); + h = f2.DrawCopy()->GetHistogram(); + + h->SetTitle("r(z)"); + h->SetXTitle("z"); + h->SetYTitle("r [kpc]"); + + gPad->Modified(); + gPad->Update(); +} + +// ------------------------------------------------------------------- + +Double_t func(Double_t *x, Double_t *k) +{ + return MPhoton::Int2(x, k)*1e68; +} + +void energyloss() +{ +/* Double_t E0 = 511e-6; // [GeV] + + Double_t Eg = 1e4; //3.6e4; + Double_t z = 5; + + Double_t val[2] = { Eg, z }; + + Double_t lolim = E0*E0/Eg; + Double_t inf = Eg<1e6 ? 3e-11*z : 3e-12*z; + + cout << Eg << " " << z << " " << lolim << " " << inf << endl; + + TF1 f("int2", func, lolim, inf, 2); + + Double_t int2 = f.Integral(lolim, inf, val); //[GeV^3 m^2] + + f.SetParameter(0, Eg); + f.SetParameter(1, z); + + cout << int2 << endl; + + new TCanvas("ILPhoton", "Mean Interaction Length Photon"); + + gPad->SetLogx(); +// gPad->SetLogy(); + gPad->SetGrid(); + + f.SetLineWidth(1); + f.DrawCopy(); + + gPad->Modified(); + gPad->Update(); + + // return; + */ MPhoton p; + p.DrawInteractionLength(); + + return; +// EnergyLossRate(); + + DrawRZ(); + + return; + +} Index: old/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- old/WuerzburgSoft/Thomas/mphys/phys.C (revision 15424) +++ old/WuerzburgSoft/Thomas/mphys/phys.C (revision 15424) @@ -0,0 +1,23 @@ +#include "MCascade.h" + +void phys() +{ + MCascade cascade; + + cascade.SetSourceZ(0.1); // Readshift of source + cascade.SetB(0/*1e-6*/); // [G] mean magnetic field + cascade.SetBradius(50); // [Mpc] + cascade.SetRuntime(8*60); // [min] maximum time to run the simulation + cascade.SetEnergyBins(18, 1e2, 1e5); // [GeV] + cascade.SetMaxInvCompton(512); // [#] maximum number of inv. Compton (-1 means infinite) + cascade.SetRatioInvCompton(0); // [%] allowed Emis of electron (0 means disabled) + cascade.SetSpectralIndex(-1); // -1 means a E^2 plot + cascade.SetDisplayIndex(1); // 1 means a E^-2 spectrum + + // + // Run the simulation: filename, eventdisply (on/off) + // + //cascade.Run("cascade_0.03_18_1e2_1e5_B1e-6_50Mpc_256_1.root", kTRUE); + cascade.Run("cascade_0.1_18_1e2_1e5_B0_512_02.root", kFALSE); +} +