Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1370) @@ -1,3 +1,29 @@ -*-*- END -*-*- + 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 Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1370) @@ -34,7 +34,11 @@ #include #include +#include +#include + +#include + #include #include -#include #include "MPhoton.h" @@ -362,2 +366,131 @@ DrawInteractionLength(fZ); } + +Bool_t MElectron::SetNewPositionB(Double_t B) +{ + static TRandom rand(0); + Double_t x = rand.Exp(GetInteractionLength()); + + // ----------------------------- + + 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); + + 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 + // + 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)); + + // ----------------------------- + + 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 + // ----------------------------- + TVectorD p(3); + + Double_t rho = TMath::Pi()/2; // [2pi] + 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; + } + + 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; + s -= p; + + fR = sqrt(s(0)*s(0)+s(1)*s(1)); + fPhi = atan2(s(1), s(0)); + fZ = ZofR(&s(2)); + + // ----------------------------- + + TVectorD w(3); + w(0) = beta_p/beta; + w(1) = beta_o/beta*cos(rho); + w(2) = beta_o/beta*sin(rho); + + w *= N2; + w *= M2; + + fPsi = atan2(w(1), w(0)); + fTheta = asin(sqrt(w(0)*w(0)+w(1)*w(1))); + + // ----------------------------- + + return fZ<0 ? kFALSE : kTRUE; +} Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1370) @@ -11,5 +11,5 @@ { public: - MElectron(Double_t e=0, Double_t z=0) : MParticle(MParticle::kEGamma) + MElectron(Double_t e=0, Double_t z=0, Bool_t type=kFALSE) : MParticle(type?MParticle::kEPositron:MParticle::kEElectron) { fEnergy = e; @@ -39,5 +39,8 @@ Double_t GetEnergyLoss(Double_t *ep) const; + // ---------------------------------------------------------------- + MPhoton *DoInvCompton(Double_t theta); + Bool_t SetNewPositionB(Double_t B); // ---------------------------------------------------------------- Index: /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1370) @@ -119,6 +119,6 @@ // cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush; - MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ()); - MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ()); + 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 Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1370) @@ -107,4 +107,11 @@ } +void MParticle::InitRandom() +{ + TRandom rnd(0); + fPhi = rnd.Uniform(TMath::Pi()*2); + fPsi = rnd.Uniform(TMath::Pi()*2); +} + void MParticle::SetNewDirection(Double_t theta, Double_t phi) { Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1370) @@ -34,4 +34,6 @@ MParticle(ParticleType_t t=kENone, const char *name=NULL, const char *title=NULL); + void InitRandom(); + static Double_t ZofR(Double_t *x, Double_t *k=NULL); static Double_t RofZ(Double_t *x, Double_t *k=NULL); @@ -66,4 +68,6 @@ Bool_t SetNewPosition(); + ParticleType_t GetPType() const { return fPType; } + Double_t GetEnergy() const { return fEnergy; } Index: /trunk/WuerzburgSoft/Thomas/mphys/anale.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/anale.C (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/anale.C (revision 1370) @@ -4,9 +4,13 @@ { TChain chain("Electrons"); - //chain.Add("cascade.root"); - +// 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"); + */ @@ -74,21 +78,15 @@ gPad->SetLogx(); gPad->SetLogy(); - h.GetXaxis()->SetRangeUser(1e4, 1e9); + h.GetXaxis()->SetRangeUser(5e3, 1e9); h.GetXaxis()->SetLabelOffset(-0.015); h.GetXaxis()->SetTitleOffset(1.1); - h.GetXaxis()->SetRangeUser(1e4, 1e9); h.GetYaxis()->SetTitleOffset(1.1); +// h.SetMinimum(1e1); +// h.SetMaximum(1e9); h.DrawCopy(); TH1 *p=h.ProfileX(); + p->SetDirectory(NULL); p->SetBit(kCanDelete); - p->SetTitle("Profile e^{-} / \\gamma "); p->SetLineColor(kRed); - p->GetXaxis()->SetLabelOffset(-0.015); - p->GetXaxis()->SetTitleOffset(1.1); - p->SetXTitle("E_{e} [GeV]"); - p->SetYTitle("E\\gamma [GeV]"); - p->SetMinimum(1e3); - p->SetMaximum(1e9); - p->GetXaxis()->SetRangeUser(1e3, 1e9); p->Draw("same"); line.DrawLine(4.2, 4.2, 8.8, 8.8); @@ -102,4 +100,5 @@ r.DrawCopy(); p=r.ProfileX(); + p->SetDirectory(NULL); p->SetBit(kCanDelete); p->SetLineColor(kRed); @@ -109,4 +108,5 @@ gPad->SetLogx(); p=h.ProjectionX(); + p->SetDirectory(NULL); p->SetBit(kCanDelete); p->SetTitle("e^{-} / \\gamma Distribution"); @@ -118,4 +118,5 @@ p->Draw(); p=h.ProjectionY(); + p->SetDirectory(NULL); p->SetBit(kCanDelete); p->SetTitle("Projection Y"); @@ -126,4 +127,5 @@ gPad->SetLogx(); p=r.ProjectionY(); + p->SetDirectory(NULL); p->SetBit(kCanDelete); p->SetTitle("Ratio E\\gamma / E_{e}"); Index: /trunk/WuerzburgSoft/Thomas/mphys/analp.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/analp.C (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/analp.C (revision 1370) @@ -1,54 +1,122 @@ //Double_t kRad2Deg = 180./TMath::Pi(); + +#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, Double_t &min, Double_t &max) +{ + TString str("MPhoton.MParticle."); + + MPhoton *p = new MPhoton; + + chain->SetBranchAddress("MPhoton.", &p); + + min = FLT_MAX; + max = -FLT_MAX; + + Int_t n = chain->GetEntries(); + for (int i=0; iGetEntry(i); + + TLeaf *leaf = chain->GetLeaf(str+name); + if (!leaf) + return 0; + Double_t val = leaf->GetValue(); + + if (p->IsPrimary()) + continue; + + if (val < min) + min = val; + if (val > max) + max = val; + } +} void analp() { TChain chain("Photons"); - chain.Add("cascade_500kpc_0*.root"); - chain.Add("cascade_500kpc_2*.root"); -// chain.Add("cascade_100kpc_0*.root"); -// chain.Add("cascade_100kpc_1*.root"); + // chain.Add("cascade_100kpc_a_1e6.root"); + // chain.Add("cascade_100kpc_b_1e6.root"); + // chain.Add("cascade_500kpc_0*.root"); + chain.Add("cascade_500kpc_21d_0*.root"); + // chain.Add("cascade_100kpc_0*.root"); + // chain.Add("cascade_100kpc_14*.root"); + // chain.Add("cascade_100kpc_0*.root"); + // chain.Add("cascade_500kpc_21d_*.root"); + // chain.Add("cascade_100kpc_4*.root"); + // chain.Add("cascade_500kpc_21d_B1e-18_*.root"); + + + Int_t n = chain.GetEntries(); + + cout << "Found " << n << " entries." << endl; + + if (n==0) + return; + + MBinning binsx; + binsx.SetEdgesLog(21, 1e4, 1e11); + + MBinning binspolx; + MBinning binspola; + MBinning binspolr; + binspolx.SetEdges(16, -180, 180); + + const Double_t ls = 299792458; // [m/s] + const Double_t ly = 3600.*24.*365.*ls; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + + Double_t max; + Double_t min; + Int_t nbins=20; + + TTreePlayer &play = *chain.GetPlayer(); + + GetRange(&chain, "fTheta", min, max); + Double_t conv = kRad2Deg*60; + min *= conv; + max *= conv; + cout << "fTheta: " << min << " < " << max << " [']" << endl; + play.FindGoodLimits(10, nbins, min, max, kFALSE); + cout << " " << nbins << ": " << min << " < " << max << endl; + binspola.SetEdges(nbins, min, max); + + GetRange(&chain, "fR", min, max); + max/=4; + cout << "fR: " << min << " < " << max << " [kpc]" << endl; + play.FindGoodLimits(10, nbins, min, max, kFALSE); + cout << " " << nbins << ": " << min << " < " << max << endl; + binspolr.SetEdges(nbins, min, max); + + TH1D h; + h.SetName("Photons"); + MH::SetBinning(&h, &binsx); + + TH1D prim; + MH::SetBinning(&prim, &binsx); + + TH2D r; + TH2D a; + MH::SetBinning(&a, &binspolx, &binspola); + MH::SetBinning(&r, &binspolx, &binspolr); MPhoton *p = new MPhoton; chain.SetBranchAddress("MPhoton.", &p); - Int_t n = chain.GetEntries(); - - cout << "Found " << n << " entries." << endl; - - if (n==0) - return; - - MBinning binsx; - binsx.SetEdgesLog(14, 1e4, 1e11); - - MBinning binspolx; - MBinning binspoly1; - MBinning binspoly2; - binspolx.SetEdges(16, -180, 180); - - const Double_t ls = 299792458; // [m/s] - const Double_t ly = 3600.*24.*365.*ls; // [m/ly] - const Double_t pc = 1./3.258; // [pc/ly] - - binspoly1.SetEdges(20, 0, 2e-6*3600); - binspoly2.SetEdges(20, 0, 1e-5*1e3/pc*365); - - TH1D h; - h.SetName("Photons"); - h.SetTitle(" E^{2} Spectra "); - h.SetXTitle("E\\gamma [GeV]"); - h.SetYTitle("E^{2} * Counts"); - MH::SetBinning(&h, &binsx); - - TH1D prim; - MH::SetBinning(&prim, &binsx); - - TH2D r; - TH2D a; - MH::SetBinning(&a, &binspolx, &binspoly1); - MH::SetBinning(&r, &binspolx, &binspoly2); - - Double_t weight = 0; - Double_t alpha = -2; + Double_t weight = 0; + Double_t alpha = -2; + Double_t plot = 2; + Double_t E; for (int i=0; iGetPhi()*kRad2Deg, p->GetR()*1e3/pc*365, weight); - a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg*3600, weight); + //if (Ep==E) + // continue; + + h.Fill(Ep, pow(Ep,plot) * weight); + Double_t v = p->GetPhi()*kRad2Deg; + r.Fill(fmod(v+180, 360)-180, p->GetR(), weight); + v = p->GetPsi()*kRad2Deg; + a.Fill(fmod(v+180, 360)-180, p->GetTheta()*kRad2Deg*60, 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; + gStyle->SetOptStat(10); + + TLine line; // delete gROOT->FindObject("Analysis Arrival"); @@ -85,4 +168,7 @@ gPad->SetTheta(70); gPad->SetPhi(90); + h.SetTitle(Form(" E^{%.1f} Spectra ", alpha)); + h.SetXTitle("E\\gamma [GeV]"); + h.SetYTitle(Form("E^{%.1f} * Counts", plot)); r.SetTitle(" Distance from Observer "); r.GetXaxis()->SetLabelOffset(-0.015); @@ -90,5 +176,5 @@ r.GetXaxis()->SetRangeUser(1e4, 1e9); r.SetXTitle("\\Phi [\\circ]"); - r.SetYTitle("R [light days]"); + r.SetYTitle("R [kpc]"); r.DrawCopy("surf1pol"); @@ -103,10 +189,12 @@ a.GetXaxis()->SetRangeUser(1e4, 1e9); a.SetXTitle("\\Psi [\\circ]"); - a.SetYTitle("\\Theta [\"]"); + a.SetYTitle("\\Theta [']"); a.DrawCopy("surf1pol"); c->cd(3); gPad->SetLogy(); - TH1 *g=r.ProjectionY("p1"); + g=r.ProjectionY(); + g->SetMinimum(GetMin(g)/2); + g->SetDirectory(NULL); g->SetBit(kCanDelete); g->SetTitle(" Hit Observers Plain "); @@ -114,5 +202,5 @@ g->GetXaxis()->SetTitleOffset(1.1); g->GetYaxis()->SetTitleOffset(1.3); - g->SetXTitle("R [light days]"); + g->SetXTitle("R [kpc]"); g->SetYTitle("Counts"); g->Draw(); @@ -120,5 +208,7 @@ c->cd(4); gPad->SetLogy(); - g=a.ProjectionY("p2"); + g=a.ProjectionY(); + g->SetMinimum(GetMin(g)/2); + g->SetDirectory(NULL); g->SetBit(kCanDelete); g->SetTitle("Hit Angle"); @@ -126,5 +216,5 @@ g->GetXaxis()->SetTitleOffset(1.1); g->GetYaxis()->SetTitleOffset(1.3); - g->SetXTitle("\\Phi [\"]"); + g->SetXTitle("\\Phi [']"); g->SetYTitle("Counts"); g->Draw(); @@ -157,9 +247,12 @@ MH::SetBinning(&div, &binsx); div.Divide(&h, &prim); - div.SetTitle("Spectrum / Primara Spectrum"); + div.SetTitle("Spectrum / Primary Spectrum"); div.GetXaxis()->SetLabelOffset(-0.015); div.GetXaxis()->SetTitleOffset(1.1); div.SetXTitle("E\\gamma [GeV]"); div.SetYTitle("Ratio"); + line.SetLineColor(kBlue); + line.SetLineWidth(1); div.DrawCopy(); + line.DrawLine(4, 1, 11, 1); } Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1369) +++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1370) @@ -76,5 +76,5 @@ Double_t w = log10(hi/lo)/n; - Double_t E = lo*pow(10, rnd.Uniform(w) + w*bin); + Double_t E = lo*pow(10, rnd.Uniform(w) + w*(n-bin-1)); //cout << endl << w << " " << n << " " << w*bin << " " << E << endl; @@ -118,6 +118,10 @@ void DoIt() { - Double_t R = 100; //MParticle::RofZ(&startz); // [kpc] + Double_t R = 500; //MParticle::RofZ(&startz); // [kpc] Double_t startz = MParticle::ZofR(&R); + + const char *filename = "cascade_500kpc_21d_B1e-18_03.root"; + + const Double_t B = 1e-18; cout << "R = " << R << "kpc" << endl; @@ -127,5 +131,5 @@ MPairProduction pair; - Double_t runtime = 15*60; // [s] + Double_t runtime = 1; // [s] Double_t lo = 1e4; @@ -208,9 +212,9 @@ MElectron *e4file = NULL; - TFile file("cascade.root", "RECREATE", "Intergalactic cascade", 9); - TTree *T = new TTree("Photons", "Photons from Cascade"); + TFile file(filename, "RECREATE", "Intergalactic cascade", 9); + TTree *T1 = new TTree("Photons", "Photons from Cascade"); TTree *T2 = new TTree("Electrons", "Electrons in the Cascade"); - TBranch *B =T ->Branch("MPhoton.", "MPhoton", &p4file, 32000); - TBranch *B2=T2->Branch("MElectron.", "MElectron", &e4file, 32000); + TBranch *B1 = T1->Branch("MPhoton.", "MPhoton", &p4file, 32000); + TBranch *B2 = T2->Branch("MElectron.", "MElectron", &e4file, 32000); // ------------------------------ @@ -234,8 +238,9 @@ MPhoton *gamma=new MPhoton(E, startz); gamma->SetIsPrimary(); + gamma->InitRandom(); listg.Add(gamma); - B->SetAddress(&gamma); - T->Fill(); + B1->SetAddress(&gamma); + T1->Fill(); cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush; @@ -248,5 +253,5 @@ TIter NextP(&listg); MPhoton *p = NULL; - B->SetAddress(&p); + B1->SetAddress(&p); while ((p=(MPhoton*)NextP())) { @@ -264,5 +269,5 @@ p->SetIsPrimary(kFALSE); - T->Fill(); + T1->Fill(); delete listg.Remove(p); @@ -327,5 +332,5 @@ while(1) { - if (!e->SetNewPosition()) + if (!e->SetNewPositionB(B)) { cout << "!" << flush; @@ -339,5 +344,5 @@ T2->Fill(); - if (fabs(p->GetTheta()*3437)<1 && // < 1min + if (fabs(p->GetTheta()*3437)<60 && // < 60min p->GetEnergy()>lo) { @@ -354,6 +359,6 @@ } - if (fabs(e->GetTheta()*3437)<1 && // < 1min - e->GetEnergy()>lo) + if (fabs(e->GetTheta()*3437)<60 && // < 60min + e->GetEnergy()>5*lo) continue;