Changeset 1370 for trunk/WuerzburgSoft/Thomas/mphys

Timestamp:

06/19/02 12:09:25 (23 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

• Changelog (modified) (1 diff)
• MElectron.cc (modified) (2 diffs)
• MElectron.h (modified) (2 diffs)
• MPairProduction.cc (modified) (1 diff)
• MParticle.cc (modified) (1 diff)
• MParticle.h (modified) (2 diffs)
• anale.C (modified) (6 diffs)
• analp.C (modified) (9 diffs)
• phys.C (modified) (10 diffs)

trunk/WuerzburgSoft/Thomas/mphys/Changelog

@@ -1 +1 @@
                                                                   -*-*- 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
 

trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

@@ -34 +34 @@
 #include <TF1.h>
 #include <TH1.h>
+#include <TMatrixD.h>
+#include <TVectorD.h>
+
+#include <TRandom.h>
+
 #include <TPad.h>
 #include <TCanvas.h>
-#include <TRandom.h>
 
 #include "MPhoton.h"
@@ -362 +366 @@
     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;
+}

trunk/WuerzburgSoft/Thomas/mphys/MElectron.h

@@ -11 +11 @@
 {
 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 +39 @@
     Double_t GetEnergyLoss(Double_t *ep) const;
 
+    // ----------------------------------------------------------------
+
     MPhoton *DoInvCompton(Double_t theta);
+    Bool_t SetNewPositionB(Double_t B);
 
     // ----------------------------------------------------------------

trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc

@@ -119 +119 @@
     // 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

trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

@@ -107 +107 @@
 }
 
+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)
 {

trunk/WuerzburgSoft/Thomas/mphys/MParticle.h

@@ -34 +34 @@
     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 +68 @@
     Bool_t SetNewPosition();
 
+    ParticleType_t GetPType() const { return fPType; }
+
     Double_t GetEnergy() const { return fEnergy; }
 

trunk/WuerzburgSoft/Thomas/mphys/anale.C

@@ -4 +4 @@
 {
     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 +78 @@
     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 +100 @@
     r.DrawCopy();
     p=r.ProfileX();
+    p->SetDirectory(NULL);
     p->SetBit(kCanDelete);
     p->SetLineColor(kRed);
@@ -109 +108 @@
     gPad->SetLogx();
     p=h.ProjectionX();
+    p->SetDirectory(NULL);
     p->SetBit(kCanDelete);
     p->SetTitle("e^{-} / \\gamma Distribution");
@@ -118 +118 @@
     p->Draw();
     p=h.ProjectionY();
+    p->SetDirectory(NULL);
     p->SetBit(kCanDelete);
     p->SetTitle("Projection Y");
@@ -126 +127 @@
     gPad->SetLogx();
     p=r.ProjectionY();
+    p->SetDirectory(NULL);
     p->SetBit(kCanDelete);
     p->SetTitle("Ratio E\\gamma / E_{e}");

trunk/WuerzburgSoft/Thomas/mphys/analp.C

@@ -1 +1 @@
 //Double_t kRad2Deg = 180./TMath::Pi();
+
+#include <float.h>
+Double_t GetMin(TH1 *h)
+{
+    Double_t min = FLT_MAX;
+    for (int i=1; i<=h->GetXaxis()->GetNbins(); i++)
+    {
+        if (h->GetBinContent(i)<min && h->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; i<n; i++)
+    {
+        chain->GetEntry(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; i<n; i++)
     {
@@ -63 +131 @@
         {
             weight = pow(Ep, alpha);
-            prim.Fill(Ep, Ep*Ep * weight);
+            prim.Fill(Ep, pow(Ep,plot) * weight);
+            E = Ep;
             continue;
         }
 
-        h.Fill(Ep, Ep*Ep * weight);
-        r.Fill(p->GetPhi()*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 +168 @@
     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 +176 @@
     r.GetXaxis()->SetRangeUser(1e4, 1e9);
     r.SetXTitle("\\Phi [\\circ]");
-    r.SetYTitle("R [light days]");
+    r.SetYTitle("R [kpc]");
     r.DrawCopy("surf1pol");
 
@@ -103 +189 @@
     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 +202 @@
     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 +208 @@
     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 +216 @@
     g->GetXaxis()->SetTitleOffset(1.1);
     g->GetYaxis()->SetTitleOffset(1.3);
-    g->SetXTitle("\\Phi [\"]");
+    g->SetXTitle("\\Phi [']");
     g->SetYTitle("Counts");
     g->Draw();
@@ -157 +247 @@
     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);
 }

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -76 +76 @@
     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 +118 @@
 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 +131 @@
     MPairProduction pair;
 
-    Double_t runtime = 15*60; // [s]
+    Double_t runtime = 1; // [s]
 
     Double_t lo = 1e4;
@@ -208 +212 @@
     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 +238 @@
         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 +253 @@
             TIter NextP(&listg);
             MPhoton *p = NULL;
-            B->SetAddress(&p);
+            B1->SetAddress(&p);
             while ((p=(MPhoton*)NextP()))
             {
@@ -264 +269 @@
 
                     p->SetIsPrimary(kFALSE);
-                    T->Fill();
+                    T1->Fill();
 
                     delete listg.Remove(p);
@@ -327 +332 @@
                 while(1)
                 {
-                    if (!e->SetNewPosition())
+                    if (!e->SetNewPositionB(B))
                     {
                         cout << "!" << flush;
@@ -339 +344 @@
                     T2->Fill();
 
-                    if (fabs(p->GetTheta()*3437)<1 &&  // < 1min
+                    if (fabs(p->GetTheta()*3437)<60 &&  // < 60min
                         p->GetEnergy()>lo)
                     {
@@ -354 +359 @@
                     }
 
-                    if (fabs(e->GetTheta()*3437)<1 &&  // < 1min
-                        e->GetEnergy()>lo)
+                    if (fabs(e->GetTheta()*3437)<60 &&  // < 60min
+                        e->GetEnergy()>5*lo)
                         continue;