Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1373) +++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1374) @@ -1,3 +1,16 @@ -*-*- END -*-*- + 2002/06/20: 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 Index: /trunk/WuerzburgSoft/Thomas/mphys/analp.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/analp.C (revision 1373) +++ /trunk/WuerzburgSoft/Thomas/mphys/analp.C (revision 1374) @@ -2,4 +2,5 @@ #include + Double_t GetMin(TH1 *h) { @@ -13,5 +14,5 @@ } -void GetRange(TChain *chain, const char *name, Double_t &min, Double_t &max) +void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1) { TString str("MPhoton.MParticle."); @@ -32,4 +33,5 @@ if (!leaf) return 0; + Double_t val = leaf->GetValue(); @@ -42,4 +44,14 @@ max = val; } + + min *= conv; + max *= conv; + + chain.GetPlayer(); + TTreePlayer &play = *chain.GetPlayer(); + + Int_t n; + play.FindGoodLimits(nbins, n, min, max, kFALSE); + nbins = n; } @@ -47,25 +59,22 @@ { TChain chain("Photons"); - // 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"); - + chain.Add("cascade_500kpc_*.root"); + + MBinning binsx; + binsx.SetEdgesLog(21, 1e4, 1e11); + + Double_t alpha = -1.8; + Double_t plot = 1.8; + + Int_t nbins = 20; + + // ------------------------------ Int_t n = chain.GetEntries(); - cout << "Found " << n << " entries." << endl; + cout << endl << "Found " << n << " entries." << endl; if (n==0) return; - - MBinning binsx; - binsx.SetEdgesLog(21, 1e4, 1e11); MBinning binspolx; @@ -74,35 +83,16 @@ 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; + Double_t min, max; + GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60); 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; + cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl; + + GetRange(&chain, "fR", nbins, min, max); binspolr.SetEdges(nbins, min, max); + cout << "fR, " << nbins << ": " << min << " ... " << max << " [kpc]" << endl; TH1D h; - h.SetName("Photons"); + TH1D prim; MH::SetBinning(&h, &binsx); - - TH1D prim; MH::SetBinning(&prim, &binsx); @@ -115,8 +105,16 @@ chain.SetBranchAddress("MPhoton.", &p); - Double_t weight = 0; - Double_t alpha = -2; - Double_t plot = 2; - Double_t E; + 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; + + Double_t weight = 0; + for (int i=0; iGetEnergy(); - if (i==0) - weight = pow(Ep, alpha); - if (p->IsPrimary()) { weight = pow(Ep, alpha); - prim.Fill(Ep, pow(Ep,plot) * weight); - E = Ep; + prim.Fill(Ep, pow(Ep, plot+alpha)); continue; } - //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); + + Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180; + Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180; + + r.Fill(phi, p->GetR(), weight); + a.Fill(psi, p->GetTheta()*kRad2Deg*60, weight); } @@ -163,12 +156,7 @@ c->Divide(2,2); + gStyle->SetPalette(1, 0); + c->cd(1); - gPad->SetLogy(); - gPad->SetLogz(); - 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); @@ -177,11 +165,32 @@ r.SetXTitle("\\Phi [\\circ]"); r.SetYTitle("R [kpc]"); - r.DrawCopy("surf1pol"); + TPad &pad = *gPad; + pad.Divide(2,2); + pad.cd(1); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(0); + gPad->SetPhi(90); + 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); - gPad->SetLogy(); - gPad->SetLogz(); - gPad->SetTheta(70); - gPad->SetPhi(90); a.SetTitle(" Hit Angle for Observer "); a.GetXaxis()->SetLabelOffset(-0.015); @@ -190,5 +199,30 @@ a.SetXTitle("\\Psi [\\circ]"); a.SetYTitle("\\Theta [']"); - a.DrawCopy("surf1pol"); + TPad &pad = *gPad; + pad.Divide(2,2); + pad.cd(1); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(0); + gPad->SetPhi(90); + g = a.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(3); @@ -227,18 +261,20 @@ gPad->SetLogx(); gPad->SetLogy(); - prim.SetFillStyle(0); h.SetFillStyle(0); + h.SetName("Photons"); + h.SetTitle(Form(" E^{%.1f} Spectra ", alpha)); + h.SetXTitle("E\\gamma [GeV]"); + h.SetYTitle(Form("E^{%.1f} * Counts", plot)); h.GetXaxis()->SetLabelOffset(-0.015); h.GetXaxis()->SetTitleOffset(1.1); -// h.GetXaxis()->SetRangeUser(1e4, 1e9); + h.SetMarkerStyle(kMultiply); + prim.SetFillStyle(0); prim.SetMarkerStyle(kPlus); - h.SetMarkerStyle(kMultiply); prim.SetMarkerColor(kRed); prim.SetLineColor(kRed); - - h.DrawCopy("P"); - prim.DrawCopy("Psame"); - prim.DrawCopy("Csame"); - h.DrawCopy("Csame"); + TH1 &g1=*h.DrawCopy("P"); + TH1 &g2=*prim.DrawCopy("Psame"); + g2.Draw("Csame"); + g1.Draw("Csame"); c->cd(2); Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1373) +++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1374) @@ -10,6 +10,7 @@ #include #include +#include +#include #include -#include #include #include @@ -121,5 +122,5 @@ Double_t startz = MParticle::ZofR(&R); - const char *filename = "data/cascade_delme.root"; + const char *filename = "cascade_delme.root"; const Double_t B = 0; @@ -128,12 +129,12 @@ cout << "Z = " << startz << endl; - static TRandom rand(0); MPairProduction pair; - Double_t runtime = 60; // [s] + Double_t runtime = 5*60; // [s] Double_t lo = 1e4; Double_t hi = 1e11; Double_t alpha = -2; + Double_t plot = 2; Int_t nbins = 21; //4*log10(hi/lo); @@ -142,6 +143,6 @@ TH1D histsrc; - hist.SetMinimum(lo*lo * pow(lo, alpha)/10); - histsrc.SetMinimum(lo*lo * pow(lo, alpha)/10); + hist.SetMinimum(pow(lo, plot+alpha)/10); + histsrc.SetMinimum(pow(lo, plot+alpha)/10); TList listg; @@ -152,25 +153,4 @@ gStyle->SetOptStat(10); - - TH2D position; - TH2D angle; - TH1D histpos; - TH1D hista; - - MBinning binspolx; - MBinning binspoly1; - MBinning binspoly2; - binspolx.SetEdges(16, -180, 180); - binspoly1.SetEdges(20, 0, 5e-6); - binspoly2.SetEdges(20, 0, 1e-5); - MH::SetBinning(&angle, &binspolx, &binspoly1); - MH::SetBinning(&position, &binspolx, &binspoly2); - MH::SetBinning(&hista, &binspoly1); - MH::SetBinning(&histpos, &binspoly2); - - hista.SetTitle("Particle Angle"); - angle.SetTitle("Particle Angle"); - histpos.SetTitle("Particle Position"); - position.SetTitle("Particle Position"); // @@ -190,5 +170,5 @@ MH::SetBinning(&histsrc, &bins); - MH::MakeDefCanvas(); + TCanvas *c=MH::MakeDefCanvas(); gPad->SetGrid(); @@ -196,32 +176,31 @@ gPad->SetLogy(); - histsrc.SetName("Spectrum"); - histsrc.SetXTitle("E [GeV]"); - histsrc.SetYTitle("E^{2}\\dotCounts"); - histsrc.GetXaxis()->SetLabelOffset(-0.015); - histsrc.GetXaxis()->SetTitleOffset(1.1); - - histsrc.Draw("P"); - hist.Draw("Psame"); - histsrc.Draw("Csame"); - hist.Draw("Csame"); + 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"); // ------------------------------ - MPhoton *p4file = NULL; - MElectron *e4file = NULL; - + void *ptr = NULL; 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 *B1 = T1->Branch("MPhoton.", "MPhoton", &p4file, 32000); - TBranch *B2 = T2->Branch("MElectron.", "MElectron", &e4file, 32000); + TTree *T1 = new TTree ("Photons", "Photons from Cascade"); + TTree *T2 = new TTree ("Electrons", "Electrons in the Cascade"); + TBranch *B1 = T1->Branch("MPhoton.", "MPhoton", &ptr); + TBranch *B2 = T2->Branch("MElectron.", "MElectron", &ptr); // ------------------------------ + + TTimer timer("gSystem->ProcessEvents();", 250, kFALSE); + timer.TurnOn(); TStopwatch clock; clock.Start(); - - Int_t timer = 0; Int_t n=0; @@ -234,5 +213,5 @@ Double_t E = GetEnergy(nbins, lo, hi); Double_t weight = pow(E, alpha); - histsrc.Fill(E, E*E * weight); + histsrc.Fill(E, pow(E, plot) * weight); MPhoton *gamma=new MPhoton(E, startz); @@ -248,5 +227,4 @@ Double_t Esum=0; Double_t Emis=0; - while (1) { @@ -265,9 +243,5 @@ cout << "!" << flush; - hist.Fill(Eg, Eg*Eg*weight); - position.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight); - angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight); - histpos.Fill(p->GetR(), weight); - hista.Fill(p->GetTheta()*kRad2Deg, weight); + hist.Fill(Eg, pow(Eg, plot)*weight); p->SetIsPrimary(kFALSE); @@ -286,11 +260,5 @@ phot.SetParameter(0, Eg); phot.SetParameter(1, p->GetZ()); - /* - if (phot.Integral(-log10(Eg)-8, -10.5, (Double_t*)NULL, 1e-2)==0) - { - cout << "z" << flush; - continue; - } - */ + Double_t Ep = pow(10, phot.GetRandom()); if (Ep==pow(10,0)) @@ -327,5 +295,5 @@ while ((e=(MElectron*)Next())) { - e->SetIsPrimary(); + e->SetIsPrimary(kTRUE); T2->Fill(); e->SetIsPrimary(kFALSE); @@ -340,12 +308,9 @@ { cout << "!" << flush; - Esum += e->GetEnergy(); - Emis += e->GetEnergy(); break; } // WRONG! - Double_t theta = 0; //rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4; - MPhoton *p = e->DoInvCompton(theta); + MPhoton *p = e->DoInvCompton(0); T2->Fill(); @@ -371,6 +336,4 @@ { cout << "T" << flush; - Esum += e->GetEnergy(); - Emis += e->GetEnergy(); break; } @@ -379,9 +342,9 @@ { cout << "E" << flush; - Esum += e->GetEnergy(); - Emis += e->GetEnergy(); break; } } + Esum += e->GetEnergy(); + Emis += e->GetEnergy(); } liste.Delete(); @@ -389,13 +352,10 @@ cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl; - Int_t now = (int)(TStopwatch::GetRealTime()- starttime)/5; - if (now!=timer || n<10) - { - histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz, - (int)runtime/60, (int)runtime%60, hist.GetEntries())); - gPad->Modified(); - gPad->Update(); - timer = now; - } + 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); } @@ -405,4 +365,6 @@ clock.Print(); + timer.Stop(); + file.Write(); @@ -414,40 +376,10 @@ gPad->Clear(); - hist.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n)); - gPad->Modified(); - gPad->Update(); - - hist.DrawCopy("P"); - histsrc.DrawCopy("Psame"); - hist.DrawCopy("Csame"); - histsrc.DrawCopy("Csame"); - - gPad->SetGrid(); - gPad->SetLogx(); - gPad->SetLogy(); - - MH::MakeDefCanvas(); - position.SetXTitle("Pos: \\Phi [\\circ]"); - position.SetYTitle("Pos: R [kpc]"); - position.DrawCopy("surf1pol"); - //gPad->SetLogy(); - - MH::MakeDefCanvas(); - angle.SetXTitle("Angle: \\Psi [\\circ]"); - angle.SetYTitle("Angle: \\Theta [\\circ]"); - angle.DrawCopy("surf1pol"); - //gPad->SetLogy(); - - MH::MakeDefCanvas(); - histpos.SetXTitle("R [kpc]"); - histpos.SetYTitle("Counts"); - histpos.DrawCopy(); - gPad->SetLogx(); - - MH::MakeDefCanvas(); - hista.SetXTitle("\\Theta [\\circ]"); - hista.SetYTitle("Counts"); - hista.DrawCopy(); - gPad->SetLogx(); + 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"); }