Index: /trunk/MagicSoft/Mars/mhist/MH.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MH.cc (revision 1448) +++ /trunk/MagicSoft/Mars/mhist/MH.cc (revision 1449) @@ -52,4 +52,5 @@ #include #include +#include #include @@ -335,2 +336,76 @@ ScaleAxis(*h->GetXaxis()->GetXbins(), fx); } + +void MH::FindGoodLimits(Int_t nbins, Int_t &newbins, Double_t &xmin, Double_t &xmax, Bool_t isInteger) +{ +//*-*-*-*-*-*-*-*-*Find reasonable bin values*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +//*-* ========================== + + Double_t dx = 0.1*(xmax-xmin); + Double_t umin = xmin - dx; + Double_t umax = xmax + dx; + + if (umin < 0 && xmin >= 0) + umin = 0; + + if (umax > 0 && xmax <= 0) + umax = 0; + + Int_t n=0; + Double_t binlow =0; + Double_t binhigh =0; + Double_t binwidth=0; + TGaxis::Optimize(umin, umax, nbins, binlow, binhigh, n, binwidth, ""); + + if (binwidth <= 0 || binwidth > 1.e+39) + { + xmin = -1; + xmax = 1; + } + else + { + xmin = binlow; + xmax = binhigh; + } + + if (isInteger) + { + Int_t ixmin = (Int_t)xmin; + Int_t ixmax = (Int_t)xmax; + Double_t dxmin = (Double_t)ixmin; + Double_t dxmax = (Double_t)ixmax; + + xmin = xmin<0 && xmin!=dxmin ? dxmin - 1 : dxmin; + xmax = xmax>0 && xmax!=dxmax ? dxmax + 1 : dxmax; + + if (xmin>=xmax) + xmax = xmin+1; + + Int_t bw = 1 + (Int_t)((xmax-xmin)/nbins); + + nbins = (Int_t)((xmax-xmin)/bw); + + if (xmin+nbins*bw < xmax) + { + nbins++; + xmax = xmin +nbins*bw; + } + } + + newbins = nbins; +} + +Double_t MH::GetMinimumGT(const TH1 &h, Double_t gt) +{ + Double_t min = FLT_MAX; + + const TAxis &axe = *((TH1&)h).GetXaxis(); + + for (int i=1; i<=axe.GetNbins(); i++) + { + Double_t x = h.GetBinContent(i); + if (gt + +#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 RandomTheta(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("RndTheta", 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, "RECREATE", "Intergalactic cascade", 9); + + 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; + } + + // WRONG! + MPhoton *p = e.DoInvCompton(0); + + 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); + + 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() +{ + TRandom r(0); + delete gRandom; + 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; +} + +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) +{ + 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 (draw) + { + 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 << "Output File... " << flush; + + TFile *file=OpenFile(filename); + + // ------------------------------ + + cout << "Timers... " << flush; + + TTimer timer("gSystem->ProcessEvents();", 333, kFALSE); + if (draw) + 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 (!draw) + 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 (draw) + { + 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: /trunk/WuerzburgSoft/Thomas/mphys/MCascade.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MCascade.h (revision 1449) +++ /trunk/WuerzburgSoft/Thomas/mphys/MCascade.h (revision 1449) @@ -0,0 +1,77 @@ +#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; //! + + 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: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1449) @@ -129,8 +129,28 @@ const Double_t E = k[0]; - Double_t omega = epsilon*E/(E0*E0); + 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(&omega)*n; + return flim*n; } @@ -250,9 +270,13 @@ const Double_t E0 = 511e-6; //[GeV] - const Double_t lolim = -log10(E)/7*4-13; - - TF1 fP("p_e", p_e, lolim, -10.6, 2); - TF1 fQ("G", G_q, 0, 1., 1); - + 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); @@ -288,6 +312,4 @@ MPhoton *MElectron::DoInvCompton(Double_t theta) { - static TRandom rand(0); - const Double_t E0 = 511e-6; //[GeV] @@ -304,18 +326,48 @@ const Double_t f = fEnergy/e; - Double_t t=-1; + Double_t t=-10; Double_t arg; do { - if (t>0) + if (t>-10) cout << "~" << flush; - t = rand.Uniform(TMath::Pi())+TMath::Pi()/2; - Double_t er = epsilon*(gamma-gammabeta*cos(t)); // photon energy rest frame - arg = (f - E0/er - 1)/(sqrt(fEnergy*fEnergy-E0*E0)/e+1); + + // + // Create an angle uniformly distributed in the range of possible + // interaction angles due to the known energies. + // + // The distibution is a theta-function which is incorrect, but + // should be correct in a first order... + // + t = gRandom->Uniform(TMath::Pi())+TMath::Pi()/2; + Double_t er = epsilon*(gamma-gammabeta*cos(t)); // photon energy rest frame + Double_t n = sqrt(fEnergy*fEnergy-E0*E0)/e+1; + arg = (f - E0/er - 1)/n; + + /* ------------------ some hints ------------------------------ + -1 < arg < 1 + -n < f - r/er - 1 < n + 1-f-n < r/er < 1-f+n + r/(1-f+n) < er < r/(1-f-n) + r/(1-f+n) < gamma-gammabeta*cos(t) < r/(1-f-n) + r/(1-f+n)-gamma < -gammabeta*cos(t) < r/(1-f-n)-gamma + gamma-r/(1-f-n) < gammabeta*cos(t) < gamma-r/(1-f+n) + (gamma-r/(1-f-n))/gammabeta < cos(t) < (gamma-r/(1-f+n))/gammabeta + acos((gamma-r/(1-f+n))/gammabeta) < t < acos((gamma-r/(1-f-n))/gammabeta) + + + (gamma+E0/(n*arg+1-f)/epsilon)/gammabeta = cos(t); + 1 = (epsilon/E0*cos(t)*gammabeta-gamma)*(n*arg+1-f); + + arg=1: (gamma+E0/epsilon*(1-f+n))/gammabeta = cos(t); + arg=-1 : (gamma+E0/epsilon*(1-f-n))/gammabeta = cos(t); + + (E0/(1-f-n)/epsilon+gamma)/gammabeta < cos(t) < (E0/(1-f+n)/epsilon+gamma)/gammabeta + */ } while (arg<-1 || arg>1); const Double_t theta1s = acos(arg); - const Double_t thetas = atan(sin(t)/(gamma*cos(t)-gammabeta)); + const Double_t thetas = atan(sin(t)/(gamma*cos(t)-gammabeta)); const Double_t thetastar = thetas-theta1s; @@ -325,8 +377,11 @@ fEnergy -= e; - const Double_t phi = rand.Uniform(TMath::Pi()*2); - + const Double_t phi = gRandom->Uniform(TMath::Pi()*2); + + /* MPhoton &p = *new MPhoton(e, fZ); p = *this; + */ + MPhoton &p = *new MPhoton(*this, e); p.SetNewDirection(theta1, phi); @@ -374,6 +429,5 @@ return SetNewPosition(); - static TRandom rand(0); - Double_t x = rand.Exp(GetInteractionLength()); + Double_t x = gRandom->Exp(GetInteractionLength()); // ----------------------------- @@ -381,8 +435,8 @@ const Double_t E0 = 511e-6; //[GeV] - Double_t B_theta = rand.Uniform(TMath::Pi()); - Double_t B_phi = rand.Uniform(TMath::Pi()*2); - - TMatrixD M(3,3); + 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); @@ -398,5 +452,5 @@ M(2, 2) = 0; - const Double_t beta = sqrt(1-E0/fEnergy*E0/fEnergy); + const Double_t beta = sqrt(1.-E0/fEnergy*E0/fEnergy); // @@ -404,5 +458,5 @@ // which the x-axis is identical with the B-Field // - TVectorD v(3); + static TVectorD v(3); v(0) = beta*sin(fTheta)*cos(fPsi); v(1) = beta*sin(fTheta)*sin(fPsi); @@ -418,5 +472,5 @@ // ----------------------------- - TMatrixD N(3,3); + static TMatrixD N(3,3); N(0, 0) = 1; N(1, 0) = 0; @@ -436,5 +490,5 @@ // v(2) = 0 // ----------------------------- - TVectorD p(3); + static TVectorD p(3); Double_t rho = 0; @@ -469,5 +523,5 @@ } - TVectorD s(3); + static TVectorD s(3); s(0) = fR*cos(fPhi); s(1) = fR*sin(fPhi); @@ -494,5 +548,5 @@ // ----------------------------- - TVectorD w(3); + static TVectorD w(3); w(0) = beta_p/beta; w(1) = beta_o/beta*cos(rho); Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1449) @@ -16,4 +16,7 @@ fZ = z; } + + MElectron(MParticle &p, Bool_t type=kFALSE) : MParticle(p, type?MParticle::kEPositron:MParticle::kEElectron) { } + MElectron(MParticle &p, Double_t e, Bool_t type=kFALSE) : MParticle(p, e, type?MParticle::kEPositron:MParticle::kEElectron) {} void operator=(MParticle &p) { MParticle::operator=(p); } Index: /trunk/WuerzburgSoft/Thomas/mphys/MFit.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MFit.cc (revision 1449) +++ /trunk/WuerzburgSoft/Thomas/mphys/MFit.cc (revision 1449) @@ -0,0 +1,165 @@ +#include "MFit.h" + +#include + +#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=kFALSE) +{ + DeleteH(); + + fHist = h; + + if (candelete) + fHist->SetBit(kIsOwner); +} + +void MFit::Fit(Bool_t log=kFALSE) +{ + 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: /trunk/WuerzburgSoft/Thomas/mphys/MFit.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MFit.h (revision 1449) +++ /trunk/WuerzburgSoft/Thomas/mphys/MFit.h (revision 1449) @@ -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, Double_t max); + + 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: /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1449) @@ -119,11 +119,8 @@ // cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush; - MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ(), kTRUE); - MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ(), kFALSE); - p0 = *gamma; // Set Position and direction - p1 = *gamma; // Set Position and direction + MElectron &p0 = *new MElectron(*gamma, E*(1.-f), kTRUE); + MElectron &p1 = *new MElectron(*gamma, E*(1.+f), kFALSE); - static TRandom rand(0); - Double_t rnd = rand.Uniform(TMath::Pi() * 2); + Double_t rnd = gRandom->Uniform(TMath::Pi() * 2); p0.SetNewDirection(atan(+tan1), rnd); p1.SetNewDirection(atan(-tan2), rnd+TMath::Pi()); Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1449) @@ -9,4 +9,5 @@ #include #include +#include ClassImp(MParticle); @@ -39,5 +40,5 @@ const Double_t r = x[0] /pc*ly*1e3; // [m] - const Double_t R = 1./(1-r*H0/c/2); // [1] + const Double_t R = 1./(1.-r*H0/c/2); // [1] return R*R - 1; @@ -67,5 +68,5 @@ const Double_t pc = 1./3.258; // [pc/ly] - const Double_t R = c/H0 * 2 * (1 - 1./sqrt(z1)); // [m] + const Double_t R = c/H0 * 2 * (1. - 1./sqrt(z1)); // [m] return R * pc/ly/1e3; // [kpc] @@ -85,8 +86,8 @@ 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 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); @@ -123,4 +124,7 @@ } + static TAxis &axez = *hist2->GetXaxis(); + static TAxis &axee = *hist2->GetYaxis(); + // // y = (y1-y0)/(x1-x0) * (x-x0) + y0 @@ -129,12 +133,12 @@ Double_t E = x[0]; - Int_t i = hist2->GetXaxis()->FindFixBin(z); - Int_t j = hist2->GetYaxis()->FindFixBin(E); - - Double_t z1 = hist2->GetXaxis()->GetBinLowEdge(i+1); - Double_t z0 = hist2->GetXaxis()->GetBinLowEdge(i); - - Double_t E1 = hist2->GetYaxis()->GetBinLowEdge(j+1); - Double_t E0 = hist2->GetYaxis()->GetBinLowEdge(j); + 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); @@ -207,31 +211,50 @@ void MParticle::InitRandom() { - static TRandom rnd(0); - fPhi = rnd.Uniform(TMath::Pi()*2); - fPsi = rnd.Uniform(TMath::Pi()*2); + fPhi = gRandom->Uniform(TMath::Pi()*2); + fPsi = gRandom->Uniform(TMath::Pi()*2); } void MParticle::SetNewDirection(Double_t theta, Double_t phi) { - TMatrixD B(3, 3); - - B(0, 0) = cos(fTheta)*cos(fPsi); - B(1, 0) = cos(fTheta)*sin(fPsi); - B(2, 0) = -sin(fTheta); - - B(0, 1) = -sin(fPsi); - B(1, 1) = cos(fPsi); + 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; - B(0, 2) = sin(fTheta)*cos(fPsi); - B(1, 2) = sin(fTheta)*sin(fPsi); - B(2, 2) = cos(fTheta); - - // ------------------------------ - - TVectorD r(3); - - r(0) = sin(theta)*cos(phi); - r(1) = sin(theta)*sin(phi); + // 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); @@ -253,42 +276,32 @@ Bool_t rc=kTRUE; - TVectorD x(3); - - x(0) = sin(fTheta)*cos(fPsi); - x(1) = sin(fTheta)*sin(fPsi); - x(2) = cos(fTheta); - - // ------------------------------ - - const Double_t R = RofZ(&fZ); - - const Double_t dx = R - dr*x(2); - + 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/x(2); // R>0 --> x(2)>0 + dr = R/d.z(); // R>0 --> x(2)>0 rc = kFALSE; } - else - fX += dr*(1.-x(2)); - - x *= dr; - - // ------------------------------ - - TVectorD r(3); - - r(0) = fR*cos(fPhi); - r(1) = fR*sin(fPhi); - r(2) = R; - - // ------------------------------ - - r -= x; + + fX += dr*(1.-d(2)); + d *= dr; + + // ------------------------------ + + TVector3 r(fR*cos(fPhi), fR*sin(fPhi), R); + + r -= d; + + // ------------------------------ if (fR!=0) - fPhi = atan2(r(1), r(0)); - - fR = sqrt(r(0)*r(0)+r(1)*r(1)); + fPhi = atan2(r.y(), r.x()); + fR = sqrt(r.x()*r.x()+r.y()*r.y()); fZ = ZofR(&r(2)); @@ -298,6 +311,5 @@ Bool_t MParticle::SetNewPosition() { - static TRandom rand(0); - Double_t r = rand.Exp(GetInteractionLength()); + Double_t r = gRandom->Exp(GetInteractionLength()); return SetNewPosition(r); } Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1449) @@ -37,4 +37,6 @@ 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(); @@ -89,5 +91,5 @@ Double_t GetPsi() const { return fPsi; } - ClassDef(MParticle, 1) // Container which holds hostograms for the Hillas parameters + ClassDef(MParticle, 2) // Container which holds hostograms for the Hillas parameters }; Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1449) @@ -200,2 +200,7 @@ DrawInteractionLength(fZ); } + +void MPhoton::Fill(TH1 &h, Double_t idx, Double_t w) const +{ + h.Fill(fEnergy, pow(fEnergy, idx)*w); +} Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1449) @@ -5,4 +5,6 @@ #include "MParticle.h" #endif + +class TH1; class MPhoton : public MParticle @@ -17,5 +19,10 @@ } + 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; // ---------------------------------------------------------------- Index: /trunk/WuerzburgSoft/Thomas/mphys/Makefile =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Makefile (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/Makefile (revision 1449) @@ -22,5 +22,5 @@ # connect the include files defined in the config.mk file # -INCLUDES = -I. -I../mbase +INCLUDES = -I. -I../mbase -I../mhist #------------------------------------------------------------------------------ @@ -29,7 +29,9 @@ SRCFILES = MParticle.cc \ + MPairProduction.cc \ MPhoton.cc \ - MElectron.cc \ - MPairProduction.cc + MFit.cc \ + MCascade.cc \ + MElectron.cc SRCS = $(SRCFILES) Index: /trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h (revision 1449) @@ -8,4 +8,6 @@ #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+; Index: /trunk/WuerzburgSoft/Thomas/mphys/analp.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/analp.C (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/analp.C (revision 1449) @@ -2,35 +2,55 @@ #include - -Double_t GetMin(TH1 *h) -{ - Double_t min = FLT_MAX; - for (int i=1; i<=h->GetXaxis()->GetNbins(); i++) - { - if (h->GetBinContent(i)GetBinContent(i)!=0) - min = h->GetBinContent(i); - } - return min; -} - -void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1) +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "MH.h" +#include "MBinning.h" + +#include "MFit.h" +#include "MPhoton.h" + +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 n = chain->GetEntries(); + 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); - TLeaf *leaf = chain->GetLeaf(str+name); - if (!leaf) - return 0; + if (numt != chain->GetTreeNumber()) + { + if (!(leaf = chain->GetLeaf(str))) + return; + + numt = chain->GetTreeNumber(); + } Double_t val = leaf->GetValue(); @@ -39,5 +59,5 @@ continue; - if (val < min) + if (val < min && (zero || val!=0)) min = val; if (val > max) @@ -48,73 +68,88 @@ max *= conv; - chain.GetPlayer(); - TTreePlayer &play = *chain.GetPlayer(); - - Int_t n; - play.FindGoodLimits(nbins, n, min, max, kFALSE); - nbins = n; + Int_t newn=0; + MH::FindGoodLimits(nbins, newn, min, max, kFALSE); + nbins = newn; } void draw1h1426(Double_t E, Double_t plot=1) { - plot=2; - Double_t level = log10(E); - level -= log10(5.73e-0*pow( 0.39e3, plot)); - - cout << "Level: " << level << endl; - - TPolyMarker *m = new TPolyMarker(7); - /* - m->SetPoint(0, 0.28e3, 1.86e-0*pow( 0.28e3, plot)*level); - m->SetPoint(1, 0.39e3, 5.73e-0*pow( 0.39e3, plot)*level); - m->SetPoint(2, 0.80e3, 1.22e-0*pow( 0.80e3, plot)*level); - m->SetPoint(3, 1.55e3, 5.10e-2*pow( 1.55e3, plot)*level); - m->SetPoint(4, 2.82e3, 1.50e-2*pow( 2.82e3, plot)*level); - m->SetPoint(5, 5.33e3, 7.70e-3*pow( 5.33e3, plot)*level); - m->SetPoint(6, 10.00e3, 1.20e-4*pow(10.00e3, plot)*level); - */ - m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot))+level); - m->SetPoint(1, log10( 0.39e3), log10(5.73e-0*pow( 0.39e3, plot))+level); - m->SetPoint(2, log10( 0.80e3), log10(1.22e-0*pow( 0.80e3, plot))+level); - m->SetPoint(3, log10( 1.55e3), log10(5.10e-2*pow( 1.55e3, plot))+level); - m->SetPoint(4, log10( 2.82e3), log10(1.50e-2*pow( 2.82e3, plot))+level); - m->SetPoint(5, log10( 5.33e3), log10(7.70e-3*pow( 5.33e3, plot))+level); - m->SetPoint(6, log10(10.00e3), log10(1.20e-4*pow(10.00e3, plot))+level); + 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->SetMarkerSize(1); + 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() { + // ------------------------------------------------------------------- + TChain chain("Photons"); + chain.Add("cascade_0.1_18_1e2_1e5_B0_256_01.root"); + chain.Add("cascade_0.1_18_1e2_1e5_B0_256_02.root"); + chain.Add("cascade_0.1_18_1e2_1e5_B0_256_03.root"); + chain.Add("cascade_0.1_18_1e2_1e5_B0_256_04.root"); + chain.Add("cascade_0.1_18_1e2_1e5_B0_512_01.root"); + chain.Add("cascade_0.1_18_1e2_1e5_B0_512_03.root"); + chain.Add("cascade_0.1_18_1e2_1e5_B0_512_02.root"); /* - chain.Add("delme3.root"); - chain.Add("delme3B.root"); - chain.Add("delme3C.root"); - chain.Add("delme3D.root"); - chain.Add("delme3E.root"); + chain.Add("cascade_0.03_18_1e2_1e5_B0_256_3.root"); + chain.Add("cascade_0.03_18_1e2_1e5_B0_256_4.root"); + chain.Add("cascade_0.03_18_1e2_1e5_B0_256_5.root"); + chain.Add("cascade_0.03_18_1e2_1e5_B0_256_6.root"); + chain.Add("cascade_0.03_18_1e2_1e5_B0_256_7.root"); + chain.Add("cascade_0.03_18_1e2_1e5_B0_256_8.root"); + chain.Add("cascade_0.03_18_1e2_1e5_B0_256_9.root"); + + chain.Add("cascade_0.03_18_1e2_1e5_B0_512_01.root"); + chain.Add("cascade_0.03_18_1e2_1e5_B0_512_02.root"); */ - chain.Add("delme3H.root"); - - Int_t nbins = 24; - Double_t lo = 1e2; - Double_t hi = 1e6; - - MBinning binsx; - binsx.SetEdgesLog(nbins, lo, hi); - - // ------------------------------ - - Double_t cutoff = 1e12; - - Double_t alpha = -1; - Double_t plot = 1; - - Int_t nbins = 16; - - // ------------------------------ - - Int_t n = chain.GetEntries(); + + 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; @@ -123,36 +158,53 @@ return; - MBinning binspolx; - MBinning binspola; - MBinning binspolr; - binspolx.SetEdges(16, -180, 180); + // ----------------------------------- + + MBinning bins; + + cout << "Determin min/max..." << endl; Double_t min, max; GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60); - binspola.SetEdges(nbins, min, max*1.1); 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); - binspolr.SetEdges(nbins, min, max*1.1); cout << "fR, " << nbins << ": " << min << " ... " << max << " [kpc]" << endl; - TH1D h, h2, h3; - TH1D prim; - MH::SetBinning(&h, &binsx); - MH::SetBinning(&h2, &binsx); - MH::SetBinning(&h3, &binsx); + 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); - - TH2D r, r2, r3; - TH2D a, a2, a3; - MH::SetBinning(&a, &binspolx, &binspola); - MH::SetBinning(&r, &binspolx, &binspolr); - MH::SetBinning(&a2, &binspolx, &binspola); - MH::SetBinning(&r2, &binspolx, &binspolr); - MH::SetBinning(&a3, &binspolx, &binspola); - MH::SetBinning(&r3, &binspolx, &binspolr); +// prim.Sumw2(); + h.Sumw2(); + h2.Sumw2(); + h3.Sumw2(); + + // ----------------------------------- MPhoton *p = new MPhoton; chain.SetBranchAddress("MPhoton.", &p); + chain.SetBranchStatus("*", 1); chain.GetEntry(0); @@ -165,4 +217,26 @@ 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; @@ -193,16 +267,11 @@ continue; - h.Fill(Ep, pow(Ep,plot) * weight); - Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180; Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180; - - r.Fill(phi, isprim?0:p->GetR(), weight); - a.Fill(psi, isprim?0:p->GetTheta()*kRad2Deg*60, weight); if (isprim) { h3.Fill(Ep, pow(Ep,plot) * weight); - r3.Fill(phi, 0, weight); + r3.Fill(phi, 0, weight); a3.Fill(psi, 0, weight); isprim = kFALSE; @@ -211,26 +280,61 @@ h2.Fill(Ep, pow(Ep,plot) * weight); - r2.Fill(phi, p->GetR(), 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; - TH1 *g=r.ProjectionX(); - cout << "Mean fPhi: " << g->GetMean() << " +- " << g->GetRMS() << endl; - delete g; - g=a.ProjectionX(); - cout << "Mean fPsi: " << g->GetMean() << " +- " << g->GetRMS() << endl; - delete g; + 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"); - c = MH::MakeDefCanvas("Analysis Arrival", ""); + TCanvas *c = MH::MakeDefCanvas("Analysis Arrival", ""); c->Divide(2,2); - - gStyle->SetPalette(1, 0); c->cd(1); @@ -241,5 +345,5 @@ r.SetXTitle("\\Phi [\\circ]"); r.SetYTitle("R [kpc]"); - TPad &pad = *gPad; + TVirtualPad &pad = *gPad; pad.Divide(2,2); pad.cd(1); @@ -248,5 +352,5 @@ gPad->SetTheta(0); gPad->SetPhi(90); - g = r.DrawCopy("surf1pol"); + TH1 *g = r.DrawCopy("surf1pol"); pad.cd(2); gPad->SetLogy(); @@ -275,7 +379,7 @@ a.SetXTitle("\\Psi [\\circ]"); a.SetYTitle("\\Theta [']"); - TPad &pad = *gPad; - pad.Divide(2,2); - pad.cd(1); + TVirtualPad &pad2 = *gPad; + pad2.Divide(2,2); + pad2.cd(1); gPad->SetLogy(); gPad->SetLogz(); @@ -283,5 +387,5 @@ gPad->SetPhi(90); g = a.DrawCopy("surf1pol"); - pad.cd(2); + pad2.cd(2); gPad->SetLogy(); gPad->SetLogz(); @@ -289,5 +393,5 @@ gPad->SetPhi(90); g->Draw("surf1pol"); - pad.cd(3); + pad2.cd(3); gPad->SetLogy(); gPad->SetLogz(); @@ -295,5 +399,5 @@ gPad->SetPhi(90); g->Draw("surf1pol"); - pad.cd(4); + pad2.cd(4); gPad->SetLogy(); gPad->SetLogz(); @@ -307,5 +411,5 @@ g->SetDirectory(NULL); TH1 *g2=r2.ProjectionY(); - g->SetMinimum(GetMin(g2)/2); + g->SetMinimum(MH::GetMinimumGT(*g2)/2); g->SetBit(kCanDelete); g->SetTitle(" Hit Observers Plain "); @@ -329,5 +433,5 @@ gPad->SetLogy(); g=a.ProjectionY(); - g->SetMinimum(GetMin(g)/2); + g->SetMinimum(MH::GetMinimumGT(*g)/2); g->SetDirectory(NULL); g->SetBit(kCanDelete); @@ -350,6 +454,8 @@ g->Draw("same"); - // delete gROOT->FindObject("Analysis Photons"); - TCanvas *c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870); + // ---------------------------------------------------------------------- + + // delete gROOT->FindObject("Analysis Photons"); + c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870); c->Divide(1,2); @@ -364,11 +470,11 @@ h.GetXaxis()->SetLabelOffset(-0.015); h.GetXaxis()->SetTitleOffset(1.1); - h.SetMarkerStyle(kMultiply); + h.SetMarkerStyle(kPlus); prim.SetFillStyle(0); prim.SetMarkerStyle(kPlus); prim.SetMarkerColor(kRed); prim.SetLineColor(kRed); - TH1 &g1=*h.DrawCopy("P"); - g2=*prim.DrawCopy("Psame"); + TH1 &g1=*h.DrawCopy("P E4"); + g2=prim.DrawCopy("Psame E0"); g2->Draw("Csame"); g1.Draw("Csame"); @@ -378,5 +484,5 @@ h2.SetMarkerColor(kBlue); h2.SetLineColor(kBlue); - TH1 &g3=*h2.DrawCopy("Psame"); + TH1 &g3=*h2.DrawCopy("Psame E2"); g3.Draw("Csame"); h3.SetFillStyle(0); @@ -385,5 +491,5 @@ h3.SetMarkerColor(kGreen); h3.SetLineColor(kGreen); - TH1 &g4=*h3.DrawCopy("Psame"); + TH1 &g4=*h3.DrawCopy("Psame E1"); g4.Draw("Csame"); @@ -394,13 +500,53 @@ TH1D div; MH::SetBinning(&div, &binsx); + div.Sumw2(); div.Divide(&h, &prim); - div.SetTitle("Spectrum / Primary Spectrum"); + div.SetTitle(" Spectrum / Primary Spectrum "); div.GetXaxis()->SetLabelOffset(-0.015); div.GetXaxis()->SetTitleOffset(1.1); div.SetXTitle("E\\gamma [GeV]"); div.SetYTitle("Ratio"); + 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.SetLineWidth(1); - div.DrawCopy(); 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 << "TeV +- "; + cout << fit.GetParError(0) << 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: /trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1448) +++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1449) @@ -1,506 +1,23 @@ -#include +#include "MCascade.h" -#include +void phys() +{ + MCascade cascade; -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include + 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 -#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 PrimSpect(Double_t *x, Double_t *k) -{ - return pow(pow(10, x[0]), k[0]); + // + // 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); } -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 RandomTheta(Double_t Eg, Double_t Ep) -{ - Double_t E0 = 511e-6; // [GeV] - - Double_t f = Eg/E0*Ep/E0; - - if (f<1) - return 0; - - TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0); - - func.SetNpx(50); - - Double_t sbar = pow(10, func.GetRandom()); - Double_t theta = acos(1.-sbar*2/f); - - return theta; -} - -Double_t GetEnergy(Int_t n, Double_t lo, Double_t hi) -{ - static int bin=0; - static TRandom rnd(0); - - Double_t w = log10(hi/lo)/n; - - Double_t E = lo*pow(10, rnd.Uniform(w) + w*(n-bin-1)); - - //cout << endl << w << " " << n << " " << w*bin << " " << E << endl; - - if (++bin==n) - bin=0; - - return E; - - /* - // TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV - // src.SetParameter(0, 0); - // Double_t E = pow(10, src.GetRandom()); - - // - // 0: 11111111111111|22222222222222 2 2^0 + 1 1 - // 1: 11111111|22222222222|33333333 3 2^1 + 1 2 - // 2: 11111|22222|33333|44444|55555 5 2^2 + 1 7 - // 3: | | | | | | | | 15 - // - - static int run=0; - static int num=1; - - Double_t w = log10(hi/lo)/((1< 5 - // b --> 10 - // c --> 2 - - // delme, delmeB starting integration at lolim - // delme2, delme2B starting integration at -log10(Eg)-8 - // delme3, lolim, 24, 1e2-1e6, 2x inv.Compton - // delme3B, lolim, 24, 1e2-1e6, 4x inv.Compton - // delme3C, lolim, 24, 1e2-1e6, 8x inv.Compton - // delme3D, lolim, 24, 1e2-1e6, 16x inv.Compton - // delme3E, lolim, 24, 1e2-1e6, 32x inv.Compton - - // delme3H, lolim, 24, 1e2-1e6, 256x inv.Compton 8*60 - - Double_t startz = 0.03; //MParticle::ZofR(&R); - Double_t R = MParticle::RofZ(&startz); // [kpc] - - const char *filename = "cascade_0.03_18_1e2_1e5_B0_256_1.root"; - //const char *filename = "test.root"; - - const Double_t B = 0;//1e-6*1e-4; // [T=1e4G] mean magnetic field - - Double_t runtime = 45*60; // [s] maximum time to run the simulation - - Int_t nbins = 18; // number of bins produced in energy spectrum - - Double_t lo = 1e2; // lower limit of displayed spectrum - Double_t hi = 1e5; // upper limit of spectrum (cutoff) - - Int_t counter = 256; // maximum number of inv. Compton (-1 means infinite) - - Double_t alpha = -1; // -1 means a E^2 plot - Double_t plot = 1; // 1 means a E^-2 spectrum - - Double_t bubbler = R-1e3*10; // [kpc] - Double_t bubblez = MParticle::ZofR(&bubbler); - - // -------------------------------------------------------------------- - - cout << "R = " << R << "kpc" << endl; - cout << "Z = " << startz << endl; - - cout << "Setting up: Histograms... " << flush; - - MPairProduction pair; - - TH1D hist; - TH1D histsrc; - - hist.SetMinimum(pow(lo, plot+alpha)/100); - histsrc.SetMinimum(pow(lo, plot+alpha)/100); - - TList listg; - TList liste; - - listg.SetOwner(); - liste.SetOwner(); - - gStyle->SetOptStat(10); - - // - // Don't change the order!!! - // - hist.SetFillStyle(0); - hist.SetMarkerStyle(kPlus); - histsrc.SetFillStyle(0); - histsrc.SetMarkerStyle(kMultiply); - histsrc.SetMarkerColor(kRed); - histsrc.SetLineColor(kRed); - - MBinning bins; - bins.SetEdgesLog(nbins, lo, hi); - - MH::SetBinning(&hist, &bins); - MH::SetBinning(&histsrc, &bins); - - TCanvas *c=MH::MakeDefCanvas(); - - gPad->SetGrid(); - gPad->SetLogx(); - gPad->SetLogy(); - - hist.SetName("Spectrum"); - hist.SetXTitle("E [GeV]"); - hist.SetYTitle(Form("E^{%.1f} Counts", plot)); - hist.GetXaxis()->SetLabelOffset(-0.015); - hist.GetXaxis()->SetTitleOffset(1.1); - - hist.Draw("P"); - histsrc.Draw("Psame"); - histsrc.Draw("same"); - hist.Draw("same"); - - // ------------------------------ - - cout << "Output File... " << flush; - - TFile file(filename, "RECREATE", "Intergalactic cascade", 9); - 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; - TBranch *B1 = T1->Branch("MPhoton.", "MPhoton", &ptr); - MPhoton dummye; - ptr = &dummye; - TBranch *B2 = T2->Branch("MElectron.", "MElectron", &ptr); - - // ------------------------------ - - cout << "Timers... " << flush; - - TTimer timer("gSystem->ProcessEvents();", 250, kFALSE); - timer.TurnOn(); - - TStopwatch clock; - clock.Start(); - - cout << "Done. " << endl; - - Int_t n=0; - Double_t starttime = TStopwatch::GetRealTime(); // s - while (TStopwatch::GetRealTime()SetIsPrimary(); - gamma->SetSrcR(R); - gamma->InitRandom(); - listg.Add(gamma); - - B1->SetAddress(&gamma); - T1->Fill(); - - cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush; - - Double_t Esum=0; - Double_t Emis=0; - while (1) - { - if (listg.GetSize()!=0) - cout << " |P" << flush; - - TIter NextP(&listg); - MPhoton *p = NULL; - B1->SetAddress(&p); - while ((p=(MPhoton*)NextP())) - { - cout << ":" << flush; - Double_t Eg = p->GetEnergy(); - - Double_t E0 = 511e-6; - Double_t z = p->GetZ(); - 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, 9.4-log10(Eg)*2); - - //TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2); - //Double_t E0 = 511e-6; // [GeV] - TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2); - phot.SetParameter(0, Eg); - while (1) - { - if (!p->SetNewPosition()) - { - cout << "!" << flush; - - hist.Fill(Eg, pow(Eg, plot)*weight); - - p->SetIsPrimary(kFALSE); - T1->Fill(); - - Esum += Eg; - - break; - } - - // - // Sample phtoton from background and interaction angle - // - phot.SetParameter(1, p->GetZ()); - Double_t pe = phot.GetRandom(); - if (pe==0) - { - cout << "z" << flush; - continue; - } - - Double_t Ep = pow(10, pe); - Double_t theta = RandomTheta(Eg, Ep); - if (theta==0) - { - cout << "t" << flush; - continue; - } - - if (!pair.Process(p, Ep, theta, &liste)) - { - cout << "0" << flush; - continue; - } - - cout << "." << flush; - break; - } - delete listg.Remove(p); - } - - if (liste.GetSize()==0 && listg.GetSize()==0) - break; - - if (liste.GetSize()!=0) - cout << " |E" << flush; - - TIter Next(&liste); - MElectron *e = NULL; - B2->SetAddress(&e); - while ((e=(MElectron*)Next())) - { - e->SetIsPrimary(kTRUE); - T2->Fill(); - e->SetIsPrimary(kFALSE); - - Double_t Ee = e->GetEnergy(); - - if (EeSetNewPositionB(e->GetZ()>bubblez ? B : 0)) - { - cout << "!" << flush; - break; - } - - // WRONG! - MPhoton *p = e->DoInvCompton(0); - - T2->Fill(); - - cout << "." << flush; - listg.Add(p); - /* - if (fabs(p->GetTheta()*3437)<60 && // < 60min - p->GetEnergy()>lo) - { - cout << "." << flush; - listg.Add(p); - } - else - { - if (p->GetEnergy()GetEnergy()<=lo) - cout << "e" << flush; - else - cout << "t" << flush; // ignored - delete p; - Esum += p->GetEnergy(); - Emis += p->GetEnergy(); - } - */ - if (fabs(e->GetTheta()*3437)>60) // < 60min - { - cout << "T" << flush; - break; - } - - if (e->GetEnergy()GetEnergy(); - Emis += e->GetEnergy(); - } - liste.Delete(); - } - cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl; - - hist.SetTitle(Form("E^{%.1f} z=%f T=%d'%d\" N=%d", alpha, startz, - (int)runtime/60, (int)runtime%60, n)); - - timer.Stop(); - c->Update(); - timer.Start(250); - - } - cout << endl; - - clock.Stop(); - clock.Print(); - - timer.Stop(); - - file.Write(); - - cout << "Wrote: " << filename << endl; - cout << "Created " << n << " gammas from source with E^" << alpha << endl; - cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma." << endl; - - // ------------------------------ - - gPad->Clear(); - - hist.SetTitle(Form("E^{%.1f} z=%f T=%d'%d\" N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n)); - - TH1 &h1 = *hist.DrawCopy("P"); - TH1 &h2 = *histsrc.DrawCopy("Psame"); - h2.Draw("Csame"); - h1.Draw("Csame"); -} - -// ------------------------------------------------------------------- -void phys() -{ - /* - Double_t Eg = 10; - Double_t E0 = 511e-6; - Double_t z = 0.03; - 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, 9.4-log10(Eg)*2); - TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2); - phot.SetParameter(0, Eg); - phot.SetParameter(1, z); - - Double_t val[2] = {Eg,z}; - cout << phot.Integral(log10(lolim), log10(inf), val, 1e-2) << endl; - - cout << lolim << " " << inf << endl; - - phot.DrawCopy(); - return; -*/ - - DoIt(); - - /* - Double_t E = 1e10; - TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2); - phot.SetParameter(0, E); - phot.SetParameter(1, 5); - phot.DrawCopy(); - return; - */ - - /* - Double_t Eg = 1e5; - - Double_t E0 = 511e-6; // [GeV] - Double_t lolim = E0*E0/Eg; - Double_t inf = 4e-12; //E0*E0/Eg * sqrt(Eg); - - // 1e5: 5e-13, 4e-12 - // 1e6: 5e-14, 2e-12 - // 1e8: 5e-16, 1e-12 - // 1e10: 5e-18, 1e-12 - - // 1e5: 2.6e-12, 4e-12 - // 1e6: 2.6e-13, 2e-12 - // 1e8: 2.6e-15, 1e-12 - // 1e10: 1.6e-17, 1e-12 - - TF1 f("int2", MPhoton::Int2, lolim, inf, 2); - f.SetParameter(0, Eg); - f.SetParameter(1, 0); - - MH::MakeDefCanvas(); - gPad->SetLogx(); - f.DrawCopy(); - */ -} -