Changeset 1428 for trunk/WuerzburgSoft


Ignore:
Timestamp:
07/24/02 09:17:52 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/WuerzburgSoft/Thomas/mphys
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

    r1371 r1428  
    174174
    175175    const Double_t from   = 1e-17;
    176     const Double_t to     = 1e-11;
     176    const Double_t to     = 2e-11;
    177177
    178178    Double_t val[3] = { E[0], z };                        // E[GeV]
    … …  
    251251    const Double_t lolim = -log10(E)/7*4-13;
    252252
    253     TF1 fP("p_e", p_e, lolim, -10.8, 2);
     253    TF1 fP("p_e", p_e, lolim, -10.6, 2);
    254254    TF1 fQ("G",   G_q, 0, 1., 1);
    255255

  • trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

    r1371 r1428  
    7272}
    7373
    74 Double_t MParticle::Planck(Double_t *x, Double_t *k)
     74#include <fstream.h>
     75#include <TH2.h>
     76#include "../mhist/MBinning.h"
     77#include "../mhist/MH.h"
     78
     79TH2D *hist2;
     80
     81Double_t MParticle::Planck(Double_t *x, Double_t *k=NULL)
     82{
     83    static Bool_t isloaded = kFALSE;
     84
     85    if (!isloaded)
     86    {
     87        Double_t c   = 299792458;             // [m/s]
     88        Double_t e   = 1.602176462e-19;       // [C]
     89        Double_t h   = 1e-9/e*6.62606876e-34; // [GeVs]
     90        Double_t hc  = h*c;                   // [GeVm]
     91        Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
     92
     93        ifstream fin("background.txt");
     94
     95        hist2 = new TH2D;
     96
     97        MBinning binsz;
     98        MBinning binse;
     99        binsz.SetEdgesLog(100, 1e-6,  1);    // --> 101 Edges / 100 bins
     100        binse.SetEdgesLog(100, 7e-15, 3e-8); // --> 101 Edges / 100 bins
     101
     102        MH::SetBinning(hist2, &binsz, &binse);
     103
     104        for (int y=0; y<101; y++)
     105        {
     106            Double_t val;
     107            fin >> val;
     108            for (int x=0; x<101; x++)
     109            {
     110                fin >> val;
     111
     112                val += 9;
     113
     114                Double_t z = binsz.GetEdges()[x];
     115                Double_t f = z+1;
     116
     117                Double_t newval = pow(10, val)/konst;
     118                hist2->SetBinContent(x, y, newval*f*f*f);
     119
     120            }
     121        }
     122        isloaded = kTRUE;
     123    }
     124
     125    //
     126    // y = (y1-y0)/(x1-x0) * (x-x0) + y0
     127    //
     128    Double_t z = k ? k[0] : 0;
     129    Double_t E = x[0];
     130
     131    Int_t i = hist2->GetXaxis()->FindFixBin(z);
     132    Int_t j = hist2->GetYaxis()->FindFixBin(E);
     133
     134    Double_t z1 = hist2->GetXaxis()->GetBinLowEdge(i+1);
     135    Double_t z0 = hist2->GetXaxis()->GetBinLowEdge(i);
     136
     137    Double_t E1 = hist2->GetYaxis()->GetBinLowEdge(j+1);
     138    Double_t E0 = hist2->GetYaxis()->GetBinLowEdge(j);
     139
     140    Double_t n00 = hist2->GetBinContent(i,   j);
     141    Double_t n01 = hist2->GetBinContent(i+1, j);
     142    Double_t n10 = hist2->GetBinContent(i,   j+1);
     143    Double_t n11 = hist2->GetBinContent(i+1, j+1);
     144
     145    Double_t dz = (z-z0)/(z1-z0);
     146    Double_t dE = (E-E0)/(E1-E0);
     147
     148    Double_t n0 = dz*(n01-n00)+n00;
     149    Double_t n1 = dz*(n11-n10)+n10;
     150
     151    Double_t n = dE*(n1-n0)+n0;
     152
     153    return n;
     154    /*
     155     //
     156     // TANJA2
     157     //
     158     Double_t E1 = x[0];
     159     Double_t E2 = x[0]/8;
     160     return (MParticle::Planck0(&E1, k)+MParticle::Planck0(&E2, k)/40e3)*0.7/0.4;
     161     */
     162    /*
     163     //
     164     // TANJA
     165     //
     166     Double_t E1 = x[0];
     167     Double_t E2 = x[0]/8;
     168     return Planck0(&E1, k)+Planck0(&E2, k)/5e3;
     169     */
     170}
     171
     172Double_t MParticle::Planck0(Double_t *x, Double_t *k)
    75173{
    76174    //

  • trunk/WuerzburgSoft/Thomas/mphys/MParticle.h

    r1370 r1428  
    5050    // ----------------------------------------------------------------
    5151
     52    static Double_t Planck0(Double_t *x, Double_t *k=NULL);
    5253    static Double_t Planck(Double_t *x, Double_t *k=NULL);
    5354

  • trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc

    r1367 r1428  
    126126
    127127    Double_t lolim = E0*E0/Eg;
    128     Double_t inf   = Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1);
     128    Double_t inf   = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1));
     129
     130    if (Eg<5e4)
     131        //inf = 3e-11*(z+1)*pow(10, 4.7*0.5-log10(Eg)*0.5);
     132        inf = 3e-11*(z+1)*pow(10, 4.7-log10(Eg));
     133        //inf = 3e-11*(z+1)*pow(10, 7.0-log10(Eg)*1.5);
     134        //inf = 3e-11*(z+1)*pow(10, 8.2-log10(Eg)*1.75);
     135        //inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);
     136
    129137    Double_t int2  = f.Integral(lolim, inf, val, 1e-2); //[GeV^3 m^2]
    130138
    … …  
    158166}
    159167
    160 void MPhoton::DrawInteractionLength(Double_t z, Double_t from, Double_t to)
     168void MPhoton::DrawInteractionLength(Double_t z, Double_t from, Double_t to, Option_t *opt)
    161169{
    162170    if (!gPad)
    … …  
    171179    gPad->SetLogy();
    172180    gPad->SetGrid();
    173     f1.SetMaximum(1e5);
     181    f1.SetMinimum(1);
     182    f1.SetMaximum(1e9);
    174183    f1.SetLineWidth(1);
    175184
    176     TH1 &h=*f1.DrawCopy()->GetHistogram();
     185    TH1 &h=*f1.DrawCopy(opt)->GetHistogram();
    177186
    178187    h.SetTitle("Mean Interaction Length (Photon)");

  • trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h

    r1364 r1428  
    3131    // ----------------------------------------------------------------
    3232
    33     static void DrawInteractionLength(Double_t z, Double_t from=1e4, Double_t to=1e11);
     33    static void DrawInteractionLength(Double_t z, Double_t from=2e2, Double_t to=1e11, Option_t*opt="");
    3434    void DrawInteractionLength() const;
    3535

  • trunk/WuerzburgSoft/Thomas/mphys/analp.C

    r1374 r1428  
    5656}
    5757
     58void draw1h1426(Double_t E, Double_t plot=1)
     59{
     60    plot=2;
     61    Double_t level = log10(E);
     62    level -= log10(5.73e-0*pow( 0.39e3, plot));
     63
     64    cout << "Level: " << level << endl;
     65
     66    TPolyMarker *m = new TPolyMarker(7);
     67    /*
     68     m->SetPoint(0,  0.28e3, 1.86e-0*pow( 0.28e3, plot)*level);
     69     m->SetPoint(1,  0.39e3, 5.73e-0*pow( 0.39e3, plot)*level);
     70     m->SetPoint(2,  0.80e3, 1.22e-0*pow( 0.80e3, plot)*level);
     71     m->SetPoint(3,  1.55e3, 5.10e-2*pow( 1.55e3, plot)*level);
     72     m->SetPoint(4,  2.82e3, 1.50e-2*pow( 2.82e3, plot)*level);
     73     m->SetPoint(5,  5.33e3, 7.70e-3*pow( 5.33e3, plot)*level);
     74     m->SetPoint(6, 10.00e3, 1.20e-4*pow(10.00e3, plot)*level);
     75     */
     76    m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot))+level);
     77    m->SetPoint(1, log10( 0.39e3), log10(5.73e-0*pow( 0.39e3, plot))+level);
     78    m->SetPoint(2, log10( 0.80e3), log10(1.22e-0*pow( 0.80e3, plot))+level);
     79    m->SetPoint(3, log10( 1.55e3), log10(5.10e-2*pow( 1.55e3, plot))+level);
     80    m->SetPoint(4, log10( 2.82e3), log10(1.50e-2*pow( 2.82e3, plot))+level);
     81    m->SetPoint(5, log10( 5.33e3), log10(7.70e-3*pow( 5.33e3, plot))+level);
     82    m->SetPoint(6, log10(10.00e3), log10(1.20e-4*pow(10.00e3, plot))+level);
     83    m->SetMarkerStyle(kOpenStar);
     84//    m->SetMarkerSize(1);
     85    m->Draw();
     86}
     87
    5888void analp()
    5989{
    6090    TChain chain("Photons");
    61     chain.Add("cascade_500kpc_*.root");
     91    /*
     92     chain.Add("delme3.root");
     93     chain.Add("delme3B.root");
     94     chain.Add("delme3C.root");
     95     chain.Add("delme3D.root");
     96     chain.Add("delme3E.root");
     97     */
     98    chain.Add("delme3H.root");
     99
     100    Int_t nbins =  24;
     101    Double_t lo = 1e2;
     102    Double_t hi = 1e6;
    62103
    63104    MBinning binsx;
    64     binsx.SetEdgesLog(21, 1e4, 1e11);
    65 
    66     Double_t alpha = -1.8;
    67     Double_t plot  =  1.8;
    68 
    69     Int_t nbins = 20;
     105    binsx.SetEdgesLog(nbins, lo, hi);
     106
     107    // ------------------------------
     108
     109    Double_t cutoff = 1e12;
     110
     111    Double_t alpha = -1;
     112    Double_t plot  =  1;
     113
     114    Int_t nbins = 16;
    70115
    71116    // ------------------------------
    … …  
    85130    Double_t min, max;
    86131    GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60);
    87     binspola.SetEdges(nbins, min, max);
     132    binspola.SetEdges(nbins, min, max*1.1);
    88133    cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl;
    89134
    90135    GetRange(&chain, "fR", nbins, min, max);
    91     binspolr.SetEdges(nbins, min, max);
     136    binspolr.SetEdges(nbins, min, max*1.1);
    92137    cout << "fR,     " << nbins << ": " << min << " ... " << max << " [kpc]" << endl;
    93138
    94     TH1D h;
     139    TH1D h, h2, h3;
    95140    TH1D prim;
    96141    MH::SetBinning(&h, &binsx);
     142    MH::SetBinning(&h2, &binsx);
     143    MH::SetBinning(&h3, &binsx);
    97144    MH::SetBinning(&prim, &binsx);
    98145
    99     TH2D r;
    100     TH2D a;
     146    TH2D r, r2, r3;
     147    TH2D a, a2, a3;
    101148    MH::SetBinning(&a, &binspolx, &binspola);
    102149    MH::SetBinning(&r, &binspolx, &binspolr);
     150    MH::SetBinning(&a2, &binspolx, &binspola);
     151    MH::SetBinning(&r2, &binspolx, &binspolr);
     152    MH::SetBinning(&a3, &binspolx, &binspola);
     153    MH::SetBinning(&r3, &binspolx, &binspolr);
    103154
    104155    MPhoton *p = new MPhoton;
    … …  
    117168    Double_t weight = 0;
    118169
     170    Bool_t isprim = kFALSE;
     171    Bool_t ignore = kFALSE;
    119172    for (int i=0; i<n; i++)
    120173    {
    … …  
    125178        if (p->IsPrimary())
    126179        {
     180            if (Ep>cutoff)
     181            {
     182                ignore = kTRUE;
     183                continue;
     184            }
     185            ignore = kFALSE;
    127186            weight = pow(Ep, alpha);
    128187            prim.Fill(Ep, pow(Ep, plot+alpha));
     188            isprim = kTRUE;
    129189            continue;
    130190        }
    131191
     192        if (ignore)
     193            continue;
     194
    132195        h.Fill(Ep, pow(Ep,plot) * weight);
    133196
    … …  
    135198        Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180;
    136199
    137         r.Fill(phi, p->GetR(),                 weight);
    138         a.Fill(psi, p->GetTheta()*kRad2Deg*60, weight);
     200        r.Fill(phi, isprim?0:p->GetR(),                 weight);
     201        a.Fill(psi, isprim?0:p->GetTheta()*kRad2Deg*60, weight);
     202
     203        if (isprim)
     204        {
     205            h3.Fill(Ep, pow(Ep,plot) * weight);
     206            r3.Fill(phi, 0,                 weight);
     207            a3.Fill(psi, 0, weight);
     208            isprim = kFALSE;
     209            continue;
     210        }
     211
     212        h2.Fill(Ep, pow(Ep,plot) * weight);
     213        r2.Fill(phi, p->GetR(),                 weight);
     214        a2.Fill(psi, p->GetTheta()*kRad2Deg*60, weight);
    139215    }
    140216
    … …  
    229305    gPad->SetLogy();
    230306    g=r.ProjectionY();
    231     g->SetMinimum(GetMin(g)/2);
    232307    g->SetDirectory(NULL);
     308    TH1 *g2=r2.ProjectionY();
     309    g->SetMinimum(GetMin(g2)/2);
    233310    g->SetBit(kCanDelete);
    234311    g->SetTitle(" Hit Observers Plain ");
    … …  
    239316    g->SetYTitle("Counts");
    240317    g->Draw();
     318    g2->SetDirectory(NULL);
     319    g2->SetBit(kCanDelete);
     320    g2->SetLineColor(kBlue);
     321    g2->Draw("same");
     322    g=r3.ProjectionY();
     323    g->SetDirectory(NULL);
     324    g->SetBit(kCanDelete);
     325    g->SetLineColor(kGreen);
     326    g->Draw("same");
    241327
    242328    c->cd(4);
    … …  
    253339    g->SetYTitle("Counts");
    254340    g->Draw();
     341    g=a2.ProjectionY();
     342    g->SetDirectory(NULL);
     343    g->SetBit(kCanDelete);
     344    g->SetLineColor(kBlue);
     345    g->Draw("same");
     346    g=a3.ProjectionY();
     347    g->SetDirectory(NULL);
     348    g->SetBit(kCanDelete);
     349    g->SetLineColor(kGreen);
     350    g->Draw("same");
    255351
    256352  //  delete gROOT->FindObject("Analysis Photons");
    … …  
    274370    prim.SetLineColor(kRed);
    275371    TH1 &g1=*h.DrawCopy("P");
    276     TH1 &g2=*prim.DrawCopy("Psame");
    277     g2.Draw("Csame");
     372    g2=*prim.DrawCopy("Psame");
     373    g2->Draw("Csame");
    278374    g1.Draw("Csame");
    279  
     375    h2.SetFillStyle(0);
     376    h2.SetName("SecPhotons");
     377    h2.SetMarkerStyle(kMultiply);
     378    h2.SetMarkerColor(kBlue);
     379    h2.SetLineColor(kBlue);
     380    TH1 &g3=*h2.DrawCopy("Psame");
     381    g3.Draw("Csame");
     382    h3.SetFillStyle(0);
     383    h3.SetName("PrimPhotons");
     384    h3.SetMarkerStyle(kMultiply);
     385    h3.SetMarkerColor(kGreen);
     386    h3.SetLineColor(kGreen);
     387    TH1 &g4=*h3.DrawCopy("Psame");
     388    g4.Draw("Csame");
     389
     390    draw1h1426(prim.GetBinContent(2),plot);
     391
    280392    c->cd(2);
    281393    gPad->SetLogx();
    … …  
    291403    line.SetLineWidth(1);
    292404    div.DrawCopy();
    293     line.DrawLine(4, 1, 11, 1);
     405    line.DrawLine(log10(lo), 1, log10(hi), 1);
    294406}
Note: See TracChangeset for help on using the changeset viewer.