Ignore:
Timestamp:
08/02/02 10:50:32 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/WuerzburgSoft/Thomas/mphys
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/WuerzburgSoft/Thomas/mphys/Changelog

    r1475 r1476  
    44   * MCascade.[h,cc]:
    55     - doesn't overwrite existing files anymore
     6     - removed redefinition of default argument
     7
     8   * MFit.cc:
     9     - removed redefinition of default argument
     10
     11   * analp.C:
     12     - added smearing
    613
    714
  • trunk/WuerzburgSoft/Thomas/mphys/MCascade.cc

    r1475 r1476  
    11#include "MCascade.h"
    22
     3#include <math.h>     // fabs, for alpha
    34#include <iostream.h>
    45
     
    215216        return kFALSE;
    216217    }
    217     return kFALSE;
    218218}
    219219
  • trunk/WuerzburgSoft/Thomas/mphys/MFit.cc

    r1449 r1476  
    9696}
    9797
    98 void MFit::SetHist(TH1 *h, Bool_t candelete=kFALSE)
     98void MFit::SetHist(TH1 *h, Bool_t candelete)
    9999{
    100100    DeleteH();
     
    106106}
    107107
    108 void MFit::Fit(Bool_t log=kFALSE)
     108void MFit::Fit(Bool_t log)
    109109{
    110110    fgBinMin = fHist->FindBin(fFunc->GetXmin()); // excluded
     
    128128}
    129129
    130 void MFit::Print(Option_t *opt="") const
     130void MFit::Print(Option_t *opt) const
    131131{
    132132    for (int i=0; i<fFunc->GetNpar(); i++)
     
    159159}
    160160
    161 TF1 *MFit::DrawCopy(Option_t *o="") const
     161TF1 *MFit::DrawCopy(Option_t *o) const
    162162{
    163163    return fFunc->DrawCopy(o);
  • trunk/WuerzburgSoft/Thomas/mphys/analp.C

    r1449 r1476  
    2020#include "MFit.h"
    2121#include "MPhoton.h"
     22
     23TH1D Smear(const TH1 &h, Double_t sigma)
     24{
     25    TH1D ret;
     26
     27    MH::SetBinning((TH1*)&ret, (TH1*)&h);
     28
     29    Int_t n = h.GetNbinsX();
     30
     31    for (int i=1; i<=n; i++)
     32    {
     33        Double_t xi = sqrt(h.GetBinLowEdge(i+1)*h.GetBinLowEdge(i));
     34        for (int j=1; j<=n; j++)
     35        {
     36            Double_t xj = sqrt(h.GetBinLowEdge(j+1)*h.GetBinLowEdge(j));
     37
     38            Double_t dx = log10(xj/xi)/sigma;
     39
     40            Double_t gaus = exp(-dx*dx/2)/(sigma*sqrt(TMath::Pi()*2));
     41
     42            gaus *= h.GetBinContent(i);
     43
     44            ret.Fill(xj, gaus/sqrt(n));
     45        }
     46    }
     47    return ret;
     48}
    2249
    2350void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1, Bool_t zero=kTRUE)
     
    110137    // -------------------------------------------------------------------
    111138
     139    TString dir = "~/data/";
     140
    112141    TChain chain("Photons");
    113     chain.Add("cascade_0.1_18_1e2_1e5_B0_256_01.root");
    114     chain.Add("cascade_0.1_18_1e2_1e5_B0_256_02.root");
    115     chain.Add("cascade_0.1_18_1e2_1e5_B0_256_03.root");
    116     chain.Add("cascade_0.1_18_1e2_1e5_B0_256_04.root");
    117     chain.Add("cascade_0.1_18_1e2_1e5_B0_512_01.root");
    118     chain.Add("cascade_0.1_18_1e2_1e5_B0_512_03.root");
    119     chain.Add("cascade_0.1_18_1e2_1e5_B0_512_02.root");
    120142    /*
    121      chain.Add("cascade_0.03_18_1e2_1e5_B0_256_3.root");
    122      chain.Add("cascade_0.03_18_1e2_1e5_B0_256_4.root");
    123      chain.Add("cascade_0.03_18_1e2_1e5_B0_256_5.root");
    124      chain.Add("cascade_0.03_18_1e2_1e5_B0_256_6.root");
    125      chain.Add("cascade_0.03_18_1e2_1e5_B0_256_7.root");
    126      chain.Add("cascade_0.03_18_1e2_1e5_B0_256_8.root");
    127      chain.Add("cascade_0.03_18_1e2_1e5_B0_256_9.root");
    128 
    129      chain.Add("cascade_0.03_18_1e2_1e5_B0_512_01.root");
    130      chain.Add("cascade_0.03_18_1e2_1e5_B0_512_02.root");
     143     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_01.root");
     144     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_02.root");
     145     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_03.root");
     146     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_04.root");
     147
     148     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_01.root");
     149     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_02.root");
     150     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_03.root");
     151     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_05.root");
    131152     */
     153
     154    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B0_512_01.root");
     155    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_10Mpc_512_01.root");
     156    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_50Mpc_512_01.root");
     157    chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_100Mpc_512_01.root");
     158
     159    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-6_10Mpc_512_01.root");
     160    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-6_50Mpc_512_01.root");
     161    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-7_100Mpc_512_01.root");
     162
     163    /*
     164     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_01.root");
     165     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_02.root");
     166     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_03.root");
     167     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_04.root");
     168     */
     169
     170    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-6_10Mpc_512_01.root");
     171    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-6_50Mpc_512_01.root");
     172    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-7_100Mpc_512_01.root");
    132173
    133174    Int_t nbinse = 18;      // number of bins in your histogram
     
    273314        {
    274315            h3.Fill(Ep, pow(Ep,plot) * weight);
    275             r3.Fill(phi, 0, weight);
    276             a3.Fill(psi, 0, weight);
     316            r3.Fill(phi, 0.0, weight);
     317            a3.Fill(psi, 0.0, weight);
    277318            isprim = kFALSE;
    278319            continue;
     
    456497    // ----------------------------------------------------------------------
    457498
    458   //  delete gROOT->FindObject("Analysis Photons");
     499    //  delete gROOT->FindObject("Analysis Photons");
    459500    c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
    460501    c->Divide(1,2);
     
    463504    gPad->SetLogx();
    464505    gPad->SetLogy();
     506    h3.SetTitle(Form(" E^{%.1f}  Spectra @ z=%.1e (R=%.0fkpc) ", alpha, z, R));
     507    h3.SetXTitle("E\\gamma [GeV]");
     508    h3.SetYTitle(Form("E^{%.1f} * Counts", plot));
     509    h3.GetXaxis()->SetLabelOffset(-0.015);
     510    h3.GetXaxis()->SetTitleOffset(1.1);
     511    h3.SetFillStyle(0);
     512    h3.SetName("PrimPhotons");
     513    h3.SetMarkerStyle(kFullCircle);
     514    h3.SetMarkerColor(kRed);
     515    h3.SetMarkerSize(0.8);
     516    h3.SetLineColor(kRed);
     517    h3.DrawCopy("C"); //E1
     518
     519    TH1D hs = Smear(h, 0.2);
     520    hs.SetLineColor(kGreen);
     521    hs.SetFillStyle(0);
     522    hs.DrawCopy("Csame");
     523
    465524    h.SetFillStyle(0);
    466525    h.SetName("Photons");
    467     h.SetTitle(Form(" E^{%.1f}  Spectra ", alpha));
    468     h.SetXTitle("E\\gamma [GeV]");
    469     h.SetYTitle(Form("E^{%.1f} * Counts", plot));
    470     h.GetXaxis()->SetLabelOffset(-0.015);
    471     h.GetXaxis()->SetTitleOffset(1.1);
    472     h.SetMarkerStyle(kPlus);
     526    h.SetMarkerStyle(kFullCircle);
     527    h.SetMarkerSize(0.8);
     528    h.DrawCopy("Csame E4");
    473529    prim.SetFillStyle(0);
    474530    prim.SetMarkerStyle(kPlus);
    475531    prim.SetMarkerColor(kRed);
    476532    prim.SetLineColor(kRed);
    477     TH1 &g1=*h.DrawCopy("P E4");
    478     g2=prim.DrawCopy("Psame E0");
    479     g2->Draw("Csame");
    480     g1.Draw("Csame");
     533    prim.SetMarkerSize(0.8);
     534    prim.DrawCopy("Csame");
    481535    h2.SetFillStyle(0);
    482536    h2.SetName("SecPhotons");
    483     h2.SetMarkerStyle(kMultiply);
     537    h2.SetMarkerStyle(kFullCircle);
    484538    h2.SetMarkerColor(kBlue);
     539    h2.SetMarkerSize(0.8);
    485540    h2.SetLineColor(kBlue);
    486     TH1 &g3=*h2.DrawCopy("Psame E2");
    487     g3.Draw("Csame");
    488     h3.SetFillStyle(0);
    489     h3.SetName("PrimPhotons");
    490     h3.SetMarkerStyle(kMultiply);
    491     h3.SetMarkerColor(kGreen);
    492     h3.SetLineColor(kGreen);
    493     TH1 &g4=*h3.DrawCopy("Psame E1");
    494     g4.Draw("Csame");
     541    h2.DrawCopy("Csame E4"); //E2
    495542
    496543    draw1h1426(prim.GetBinContent(2),plot);
     
    507554    div.SetXTitle("E\\gamma [GeV]");
    508555    div.SetYTitle("Ratio");
     556    div.SetMarkerStyle(kPlus);
    509557    TH1 *gHist = div.DrawCopy("E4");
    510558
     
    526574        }
    527575    }
    528     fit.SetParameter(0, "t", 750,  10, 1e5);
    529     fit.SetParameter(1, "m", 0.5, 0.1,  10);
     576    fit.SetParameter(0, "t", 750,   10, 1e5);
     577    fit.SetParameter(1, "m", 0.5,  0.1,  10);
    530578    fit.FitLog(gHist);
    531579    fit.Print();
    532580    fit.DrawCopy("same");
    533581
    534     cout << "Cutoff:  " << setprecision(2) << fit[0]/1e3 << "TeV +- ";
    535     cout << fit.GetParError(0) << endl;
     582    cout << "Cutoff:  " << setprecision(2) << fit[0]/1e3 << " +- ";
     583    cout << fit.GetParError(0)/1e3 << " TeV" << endl;
    536584
    537585     // ----------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.