Ignore:
Timestamp:
06/19/02 12:09:25 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/WuerzburgSoft/Thomas/mphys
Files:
9 edited

Legend:

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

    r1369 r1370  
    11                                                                  -*-*- END -*-*-
     2 2002/06/17: Thomas Bretz
     3
     4   * MElectron.[h,cc]:
     5     - Added SetNewPositionB
     6     - changed the constructor to support positrons
     7
     8   * MPairProduction.cc:
     9     - use new constructor of MElectron
     10
     11   * MParticle.[h,cc]:
     12     - added InitRandom
     13     - added GetType
     14
     15   * anale.C:
     16     - added SetDirectory(NULL)
     17
     18   * analp.C:
     19     - changed to an automatic scaling
     20
     21   * phys.C:
     22     - changed limit to 60min
     23     - added lower limit for electrons 5*lo
     24     - Changed to use a B-field
     25
     26
     27
    228 2002/06/17: Thomas Bretz
    329
  • trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

    r1369 r1370  
    3434#include <TF1.h>
    3535#include <TH1.h>
     36#include <TMatrixD.h>
     37#include <TVectorD.h>
     38
     39#include <TRandom.h>
     40
    3641#include <TPad.h>
    3742#include <TCanvas.h>
    38 #include <TRandom.h>
    3943
    4044#include "MPhoton.h"
     
    362366    DrawInteractionLength(fZ);
    363367}
     368
     369Bool_t MElectron::SetNewPositionB(Double_t B)
     370{
     371    static TRandom rand(0);
     372    Double_t x = rand.Exp(GetInteractionLength());
     373
     374    // -----------------------------
     375
     376    const Double_t E0 = 511e-6;            //[GeV]
     377
     378    Double_t B_theta = rand.Uniform(TMath::Pi());
     379    Double_t B_phi   = rand.Uniform(TMath::Pi()*2);
     380
     381    TMatrixD M(3,3);
     382
     383    M(0, 0) =  sin(B_theta)*cos(B_phi);
     384    M(1, 0) =  cos(B_theta)*cos(B_phi);
     385    M(2, 0) = -sin(B_phi);
     386
     387    M(0, 1) =  sin(B_theta)*sin(B_phi);
     388    M(1, 1) =  cos(B_theta)*sin(B_phi);
     389    M(2, 1) =  cos(B_phi);
     390
     391    M(0, 2) =  cos(B_theta);
     392    M(1, 2) = -sin(B_theta);
     393    M(2, 2) =  0;
     394
     395    const Double_t beta = sqrt(1-E0/fEnergy*E0/fEnergy);
     396
     397    //
     398    // rotate vector of velocity into a system in
     399    // which the x-axis is identical with the B-Field
     400    //
     401    TVectorD v(3);
     402    v(0) = beta*sin(fTheta)*cos(fPsi);
     403    v(1) = beta*sin(fTheta)*sin(fPsi);
     404    v(2) = beta*cos(fTheta);
     405    v *= M;
     406
     407    //
     408    // Now rotate the system so, that the velocity vector
     409    // lays in the y-z-plain
     410    //
     411    Double_t chi = atan(-v(2)/v(1));
     412
     413    // -----------------------------
     414
     415    TMatrixD N(3,3);
     416    N(0, 0) =  1;
     417    N(1, 0) =  0;
     418    N(2, 0) =  0;
     419
     420    N(0, 1) =  0;
     421    N(1, 1) = -cos(chi);
     422    N(2, 1) = -sin(chi);
     423
     424    N(0, 2) =  0;
     425    N(1, 2) =  sin(chi);
     426    N(2, 2) = -cos(chi);
     427    v *= N;
     428
     429    Double_t beta_p = v(0);         // beta parallel
     430    Double_t beta_o = v(1);         // beta orthogonal
     431                                    // v(2) = 0
     432    // -----------------------------
     433    TVectorD p(3);
     434
     435    Double_t rho = TMath::Pi()/2;    // [2pi]
     436    if (B>0)
     437    {
     438        const Double_t c  = 299792458;               // [m/s]
     439        const Double_t ly = 3600.*24.*365.*c;        // [m/ly]
     440        const Double_t pc = 1./3.258;                // [pc/ly]
     441
     442        Double_t r = (fEnergy*1e9)*beta_o/(c*B)*(pc/ly/1e3);  // [kpc]
     443
     444        rho = beta_o/beta*x/r;                                // [2pi]
     445
     446        // -----------------------------
     447
     448        if (fabs(rho*3437)>5*60)  // > 5*60min
     449        {
     450            cout << "r" << flush;
     451            return kFALSE;
     452        }
     453
     454        p(0) = beta_p/beta*x;
     455        p(1) = GetPType()==kEElectron?r*sin(rho):-r*sin(rho);
     456        p(2) = r*(1-cos(rho));
     457    }
     458    else
     459    {
     460        p(0) = beta_p/beta*x;
     461        p(1) = beta_o/beta*x;
     462        p(2) = 0;
     463    }
     464
     465    TVectorD s(3);
     466    s(0) = fR*cos(fPhi);
     467    s(1) = fR*sin(fPhi);
     468    s(2) = RofZ(&fZ);
     469
     470    TMatrixD N2(TMatrixD::kTransposed, N);
     471    TMatrixD M2(TMatrixD::kTransposed, M);
     472    p *= N2;
     473    p *= M2;
     474    s -= p;
     475
     476    fR   = sqrt(s(0)*s(0)+s(1)*s(1));
     477    fPhi = atan2(s(1), s(0));
     478    fZ   = ZofR(&s(2));
     479
     480    // -----------------------------
     481
     482    TVectorD w(3);
     483    w(0) = beta_p/beta;
     484    w(1) = beta_o/beta*cos(rho);
     485    w(2) = beta_o/beta*sin(rho);
     486
     487    w *= N2;
     488    w *= M2;
     489
     490    fPsi   = atan2(w(1), w(0));
     491    fTheta = asin(sqrt(w(0)*w(0)+w(1)*w(1)));
     492
     493    // -----------------------------
     494
     495    return fZ<0 ? kFALSE : kTRUE;
     496}
  • trunk/WuerzburgSoft/Thomas/mphys/MElectron.h

    r1364 r1370  
    1111{
    1212public:
    13     MElectron(Double_t e=0, Double_t z=0) : MParticle(MParticle::kEGamma)
     13    MElectron(Double_t e=0, Double_t z=0, Bool_t type=kFALSE) : MParticle(type?MParticle::kEPositron:MParticle::kEElectron)
    1414    {
    1515        fEnergy = e;
     
    3939    Double_t GetEnergyLoss(Double_t *ep) const;
    4040
     41    // ----------------------------------------------------------------
     42
    4143    MPhoton *DoInvCompton(Double_t theta);
     44    Bool_t SetNewPositionB(Double_t B);
    4245
    4346    // ----------------------------------------------------------------
  • trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc

    r1367 r1370  
    119119    // cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush;
    120120
    121     MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ());
    122     MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ());
     121    MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ(), kTRUE);
     122    MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ(), kFALSE);
    123123    p0 = *gamma; // Set Position and direction
    124124    p1 = *gamma; // Set Position and direction
  • trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

    r1369 r1370  
    107107}
    108108
     109void MParticle::InitRandom()
     110{
     111    TRandom rnd(0);
     112    fPhi = rnd.Uniform(TMath::Pi()*2);
     113    fPsi = rnd.Uniform(TMath::Pi()*2);
     114}
     115
    109116void MParticle::SetNewDirection(Double_t theta, Double_t phi)
    110117{
  • trunk/WuerzburgSoft/Thomas/mphys/MParticle.h

    r1364 r1370  
    3434    MParticle(ParticleType_t t=kENone, const char *name=NULL, const char *title=NULL);
    3535
     36    void InitRandom();
     37
    3638    static Double_t ZofR(Double_t *x, Double_t *k=NULL);
    3739    static Double_t RofZ(Double_t *x, Double_t *k=NULL);
     
    6668    Bool_t SetNewPosition();
    6769
     70    ParticleType_t GetPType() const { return fPType; }
     71
    6872    Double_t GetEnergy() const { return fEnergy; }
    6973
  • trunk/WuerzburgSoft/Thomas/mphys/anale.C

    r1364 r1370  
    44{
    55    TChain chain("Electrons");
    6     //chain.Add("cascade.root");
    7 
     6//    chain.Add("cascade_100kpc_0*.root");
     7//    chain.Add("cascade_100kpc_14*.root");
     8//    chain.Add("cascade_100kpc_4*.root");
     9    chain.Add("cascade_100kpc_21d.root");
     10/*
    811    chain.Add("cascade_100kpc_01.root");
    912    chain.Add("cascade_100kpc_02.root");
    1013    chain.Add("cascade_100kpc_03.root");
     14    */
    1115
    1216
     
    7478    gPad->SetLogx();
    7579    gPad->SetLogy();
    76     h.GetXaxis()->SetRangeUser(1e4, 1e9);
     80    h.GetXaxis()->SetRangeUser(5e3, 1e9);
    7781    h.GetXaxis()->SetLabelOffset(-0.015);
    7882    h.GetXaxis()->SetTitleOffset(1.1);
    79     h.GetXaxis()->SetRangeUser(1e4, 1e9);
    8083    h.GetYaxis()->SetTitleOffset(1.1);
     84//    h.SetMinimum(1e1);
     85//    h.SetMaximum(1e9);
    8186    h.DrawCopy();
    8287    TH1 *p=h.ProfileX();
     88    p->SetDirectory(NULL);
    8389    p->SetBit(kCanDelete);
    84     p->SetTitle("Profile e^{-} / \\gamma ");
    8590    p->SetLineColor(kRed);
    86     p->GetXaxis()->SetLabelOffset(-0.015);
    87     p->GetXaxis()->SetTitleOffset(1.1);
    88     p->SetXTitle("E_{e} [GeV]");
    89     p->SetYTitle("E\\gamma [GeV]");
    90     p->SetMinimum(1e3);
    91     p->SetMaximum(1e9);
    92     p->GetXaxis()->SetRangeUser(1e3, 1e9);
    9391    p->Draw("same");
    9492    line.DrawLine(4.2, 4.2, 8.8, 8.8);
     
    102100    r.DrawCopy();
    103101    p=r.ProfileX();
     102    p->SetDirectory(NULL);
    104103    p->SetBit(kCanDelete);
    105104    p->SetLineColor(kRed);
     
    109108    gPad->SetLogx();
    110109    p=h.ProjectionX();
     110    p->SetDirectory(NULL);
    111111    p->SetBit(kCanDelete);
    112112    p->SetTitle("e^{-} / \\gamma Distribution");
     
    118118    p->Draw();
    119119    p=h.ProjectionY();
     120    p->SetDirectory(NULL);
    120121    p->SetBit(kCanDelete);
    121122    p->SetTitle("Projection Y");
     
    126127    gPad->SetLogx();
    127128    p=r.ProjectionY();
     129    p->SetDirectory(NULL);
    128130    p->SetBit(kCanDelete);
    129131    p->SetTitle("Ratio E\\gamma / E_{e}");
  • trunk/WuerzburgSoft/Thomas/mphys/analp.C

    r1364 r1370  
    11//Double_t kRad2Deg = 180./TMath::Pi();
     2
     3#include <float.h>
     4Double_t GetMin(TH1 *h)
     5{
     6    Double_t min = FLT_MAX;
     7    for (int i=1; i<=h->GetXaxis()->GetNbins(); i++)
     8    {
     9        if (h->GetBinContent(i)<min && h->GetBinContent(i)!=0)
     10            min = h->GetBinContent(i);
     11    }
     12    return min;
     13}
     14
     15void GetRange(TChain *chain, const char *name, Double_t &min, Double_t &max)
     16{
     17    TString str("MPhoton.MParticle.");
     18
     19    MPhoton *p = new MPhoton;
     20
     21    chain->SetBranchAddress("MPhoton.", &p);
     22
     23    min =  FLT_MAX;
     24    max = -FLT_MAX;
     25
     26    Int_t n = chain->GetEntries();
     27    for (int i=0; i<n; i++)
     28    {
     29        chain->GetEntry(i);
     30
     31        TLeaf *leaf = chain->GetLeaf(str+name);
     32        if (!leaf)
     33            return 0;
     34        Double_t val = leaf->GetValue();
     35
     36        if (p->IsPrimary())
     37            continue;
     38
     39        if (val < min)
     40            min = val;
     41        if (val > max)
     42            max = val;
     43    }
     44}
    245
    346void analp()
    447{
    548    TChain chain("Photons");
    6     chain.Add("cascade_500kpc_0*.root");
    7     chain.Add("cascade_500kpc_2*.root");
    8 //    chain.Add("cascade_100kpc_0*.root");
    9 //    chain.Add("cascade_100kpc_1*.root");
     49    // chain.Add("cascade_100kpc_a_1e6.root");
     50    // chain.Add("cascade_100kpc_b_1e6.root");
     51    // chain.Add("cascade_500kpc_0*.root");
     52    chain.Add("cascade_500kpc_21d_0*.root");
     53    // chain.Add("cascade_100kpc_0*.root");
     54    // chain.Add("cascade_100kpc_14*.root");
     55    // chain.Add("cascade_100kpc_0*.root");
     56    // chain.Add("cascade_500kpc_21d_*.root");
     57    // chain.Add("cascade_100kpc_4*.root");
     58    // chain.Add("cascade_500kpc_21d_B1e-18_*.root");
     59
     60
     61    Int_t n = chain.GetEntries();
     62
     63    cout << "Found " << n << " entries." << endl;
     64
     65    if (n==0)
     66        return;
     67
     68    MBinning binsx;
     69    binsx.SetEdgesLog(21, 1e4, 1e11);
     70
     71    MBinning binspolx;
     72    MBinning binspola;
     73    MBinning binspolr;
     74    binspolx.SetEdges(16, -180, 180);
     75
     76    const Double_t ls = 299792458;         // [m/s]
     77    const Double_t ly = 3600.*24.*365.*ls; // [m/ly]
     78    const Double_t pc = 1./3.258;          // [pc/ly]
     79
     80    Double_t max;
     81    Double_t min;
     82    Int_t nbins=20;
     83
     84    TTreePlayer &play = *chain.GetPlayer();
     85
     86    GetRange(&chain, "fTheta", min, max);
     87    Double_t conv = kRad2Deg*60;
     88    min *= conv;
     89    max *= conv;
     90    cout << "fTheta: " << min  << " < " << max  << " [']" << endl;
     91    play.FindGoodLimits(10, nbins, min, max, kFALSE);
     92    cout << "    " << nbins << ": " << min << " < " << max << endl;
     93    binspola.SetEdges(nbins, min, max);
     94
     95    GetRange(&chain, "fR", min, max);
     96    max/=4;
     97    cout << "fR:     " << min  << " < " << max  << " [kpc]" << endl;
     98    play.FindGoodLimits(10, nbins, min, max, kFALSE);
     99    cout << "    " << nbins << ": " << min << " < " << max << endl;
     100    binspolr.SetEdges(nbins, min, max);
     101
     102    TH1D h;
     103    h.SetName("Photons");
     104    MH::SetBinning(&h, &binsx);
     105
     106    TH1D prim;
     107    MH::SetBinning(&prim, &binsx);
     108
     109    TH2D r;
     110    TH2D a;
     111    MH::SetBinning(&a, &binspolx, &binspola);
     112    MH::SetBinning(&r, &binspolx, &binspolr);
    10113
    11114    MPhoton *p = new MPhoton;
    12115    chain.SetBranchAddress("MPhoton.", &p);
    13116
    14     Int_t n = chain.GetEntries();
    15 
    16     cout << "Found " << n << " entries." << endl;
    17 
    18     if (n==0)
    19         return;
    20 
    21     MBinning binsx;
    22     binsx.SetEdgesLog(14, 1e4, 1e11);
    23 
    24     MBinning binspolx;
    25     MBinning binspoly1;
    26     MBinning binspoly2;
    27     binspolx.SetEdges(16, -180, 180);
    28 
    29     const Double_t ls = 299792458;        // [m/s]
    30     const Double_t ly = 3600.*24.*365.*ls; // [m/ly]
    31     const Double_t pc = 1./3.258;         // [pc/ly]
    32 
    33     binspoly1.SetEdges(20, 0, 2e-6*3600);
    34     binspoly2.SetEdges(20, 0, 1e-5*1e3/pc*365);
    35 
    36     TH1D h;
    37     h.SetName("Photons");
    38     h.SetTitle(" E^{2} Spectra ");
    39     h.SetXTitle("E\\gamma [GeV]");
    40     h.SetYTitle("E^{2} * Counts");
    41     MH::SetBinning(&h, &binsx);
    42 
    43     TH1D prim;
    44     MH::SetBinning(&prim, &binsx);
    45 
    46     TH2D r;
    47     TH2D a;
    48     MH::SetBinning(&a, &binspolx, &binspoly1);
    49     MH::SetBinning(&r, &binspolx, &binspoly2);
    50 
    51     Double_t weight = 0;
    52     Double_t alpha = -2;
     117    Double_t weight =  0;
     118    Double_t alpha  = -2;
     119    Double_t plot   =  2;
     120    Double_t E;
    53121    for (int i=0; i<n; i++)
    54122    {
     
    63131        {
    64132            weight = pow(Ep, alpha);
    65             prim.Fill(Ep, Ep*Ep * weight);
     133            prim.Fill(Ep, pow(Ep,plot) * weight);
     134            E = Ep;
    66135            continue;
    67136        }
    68137
    69         h.Fill(Ep, Ep*Ep * weight);
    70         r.Fill(p->GetPhi()*kRad2Deg, p->GetR()*1e3/pc*365, weight);
    71         a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg*3600, weight);
     138        //if (Ep==E)
     139        //    continue;
     140
     141        h.Fill(Ep, pow(Ep,plot) * weight);
     142        Double_t v = p->GetPhi()*kRad2Deg;
     143        r.Fill(fmod(v+180, 360)-180, p->GetR(), weight);
     144        v = p->GetPsi()*kRad2Deg;
     145        a.Fill(fmod(v+180, 360)-180, p->GetTheta()*kRad2Deg*60, weight);
    72146    }
    73147
    74148    delete p;
    75149
     150    TH1 *g=r.ProjectionX();
     151    cout << "Mean fPhi: " << g->GetMean() << " +- " << g->GetRMS() << endl;
     152    delete g;
     153    g=a.ProjectionX();
     154    cout << "Mean fPsi: " << g->GetMean() << " +- " << g->GetRMS() <<  endl;
     155    delete g;
     156
    76157    gStyle->SetOptStat(10);
     158
     159    TLine line;
    77160
    78161//    delete gROOT->FindObject("Analysis Arrival");
     
    85168    gPad->SetTheta(70);
    86169    gPad->SetPhi(90);
     170    h.SetTitle(Form(" E^{%.1f}  Spectra ", alpha));
     171    h.SetXTitle("E\\gamma [GeV]");
     172    h.SetYTitle(Form("E^{%.1f} * Counts", plot));
    87173    r.SetTitle(" Distance from Observer ");
    88174    r.GetXaxis()->SetLabelOffset(-0.015);
     
    90176    r.GetXaxis()->SetRangeUser(1e4, 1e9);
    91177    r.SetXTitle("\\Phi [\\circ]");
    92     r.SetYTitle("R [light days]");
     178    r.SetYTitle("R [kpc]");
    93179    r.DrawCopy("surf1pol");
    94180
     
    103189    a.GetXaxis()->SetRangeUser(1e4, 1e9);
    104190    a.SetXTitle("\\Psi [\\circ]");
    105     a.SetYTitle("\\Theta [\"]");
     191    a.SetYTitle("\\Theta [']");
    106192    a.DrawCopy("surf1pol");
    107193
    108194    c->cd(3);
    109195    gPad->SetLogy();
    110     TH1 *g=r.ProjectionY("p1");
     196    g=r.ProjectionY();
     197    g->SetMinimum(GetMin(g)/2);
     198    g->SetDirectory(NULL);
    111199    g->SetBit(kCanDelete);
    112200    g->SetTitle(" Hit Observers Plain ");
     
    114202    g->GetXaxis()->SetTitleOffset(1.1);
    115203    g->GetYaxis()->SetTitleOffset(1.3);
    116     g->SetXTitle("R [light days]");
     204    g->SetXTitle("R [kpc]");
    117205    g->SetYTitle("Counts");
    118206    g->Draw();
     
    120208    c->cd(4);
    121209    gPad->SetLogy();
    122     g=a.ProjectionY("p2");
     210    g=a.ProjectionY();
     211    g->SetMinimum(GetMin(g)/2);
     212    g->SetDirectory(NULL);
    123213    g->SetBit(kCanDelete);
    124214    g->SetTitle("Hit Angle");
     
    126216    g->GetXaxis()->SetTitleOffset(1.1);
    127217    g->GetYaxis()->SetTitleOffset(1.3);
    128     g->SetXTitle("\\Phi [\"]");
     218    g->SetXTitle("\\Phi [']");
    129219    g->SetYTitle("Counts");
    130220    g->Draw();
     
    157247    MH::SetBinning(&div, &binsx);
    158248    div.Divide(&h, &prim);
    159     div.SetTitle("Spectrum / Primara Spectrum");
     249    div.SetTitle("Spectrum / Primary Spectrum");
    160250    div.GetXaxis()->SetLabelOffset(-0.015);
    161251    div.GetXaxis()->SetTitleOffset(1.1);
    162252    div.SetXTitle("E\\gamma [GeV]");
    163253    div.SetYTitle("Ratio");
     254    line.SetLineColor(kBlue);
     255    line.SetLineWidth(1);
    164256    div.DrawCopy();
     257    line.DrawLine(4, 1, 11, 1);
    165258}
  • trunk/WuerzburgSoft/Thomas/mphys/phys.C

    r1367 r1370  
    7676    Double_t w = log10(hi/lo)/n;
    7777
    78     Double_t E = lo*pow(10, rnd.Uniform(w) + w*bin);
     78    Double_t E = lo*pow(10, rnd.Uniform(w) + w*(n-bin-1));
    7979
    8080    //cout << endl << w << " " << n << " " << w*bin << " " << E << endl;
     
    118118void DoIt()
    119119{
    120     Double_t R      = 100; //MParticle::RofZ(&startz); // [kpc]
     120    Double_t R      = 500; //MParticle::RofZ(&startz); // [kpc]
    121121    Double_t startz = MParticle::ZofR(&R);
     122
     123    const char *filename = "cascade_500kpc_21d_B1e-18_03.root";
     124
     125    const Double_t B = 1e-18;
    122126
    123127    cout << "R = " << R << "kpc" << endl;
     
    127131    MPairProduction pair;
    128132
    129     Double_t runtime = 15*60; // [s]
     133    Double_t runtime = 1; // [s]
    130134
    131135    Double_t lo = 1e4;
     
    208212    MElectron *e4file  = NULL;
    209213
    210     TFile file("cascade.root", "RECREATE", "Intergalactic cascade", 9);
    211     TTree *T  = new TTree("Photons",   "Photons from Cascade");
     214    TFile file(filename, "RECREATE", "Intergalactic cascade", 9);
     215    TTree *T1 = new TTree("Photons",   "Photons from Cascade");
    212216    TTree *T2 = new TTree("Electrons", "Electrons in the Cascade");
    213     TBranch *B =T ->Branch("MPhoton.",   "MPhoton",   &p4file, 32000);
    214     TBranch *B2=T2->Branch("MElectron.", "MElectron", &e4file, 32000);
     217    TBranch *B1 = T1->Branch("MPhoton.",   "MPhoton",   &p4file, 32000);
     218    TBranch *B2 = T2->Branch("MElectron.", "MElectron", &e4file, 32000);
    215219
    216220    // ------------------------------
     
    234238        MPhoton *gamma=new MPhoton(E, startz);
    235239        gamma->SetIsPrimary();
     240        gamma->InitRandom();
    236241        listg.Add(gamma);
    237242
    238         B->SetAddress(&gamma);
    239         T->Fill();
     243        B1->SetAddress(&gamma);
     244        T1->Fill();
    240245
    241246        cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush;
     
    248253            TIter NextP(&listg);
    249254            MPhoton *p = NULL;
    250             B->SetAddress(&p);
     255            B1->SetAddress(&p);
    251256            while ((p=(MPhoton*)NextP()))
    252257            {
     
    264269
    265270                    p->SetIsPrimary(kFALSE);
    266                     T->Fill();
     271                    T1->Fill();
    267272
    268273                    delete listg.Remove(p);
     
    327332                while(1)
    328333                {
    329                     if (!e->SetNewPosition())
     334                    if (!e->SetNewPositionB(B))
    330335                    {
    331336                        cout << "!" << flush;
     
    339344                    T2->Fill();
    340345
    341                     if (fabs(p->GetTheta()*3437)<1 &&  // < 1min
     346                    if (fabs(p->GetTheta()*3437)<60 &&  // < 60min
    342347                        p->GetEnergy()>lo)
    343348                    {
     
    354359                    }
    355360
    356                     if (fabs(e->GetTheta()*3437)<1 &&  // < 1min
    357                         e->GetEnergy()>lo)
     361                    if (fabs(e->GetTheta()*3437)<60 &&  // < 60min
     362                        e->GetEnergy()>5*lo)
    358363                        continue;
    359364
Note: See TracChangeset for help on using the changeset viewer.