Changeset 1370 for trunk/WuerzburgSoft/Thomas/mphys
- Timestamp:
- 06/19/02 12:09:25 (22 years ago)
- Location:
- trunk/WuerzburgSoft/Thomas/mphys
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/WuerzburgSoft/Thomas/mphys/Changelog
r1369 r1370 1 1 -*-*- 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 2 28 2002/06/17: Thomas Bretz 3 29 -
trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc
r1369 r1370 34 34 #include <TF1.h> 35 35 #include <TH1.h> 36 #include <TMatrixD.h> 37 #include <TVectorD.h> 38 39 #include <TRandom.h> 40 36 41 #include <TPad.h> 37 42 #include <TCanvas.h> 38 #include <TRandom.h>39 43 40 44 #include "MPhoton.h" … … 362 366 DrawInteractionLength(fZ); 363 367 } 368 369 Bool_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 11 11 { 12 12 public: 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) 14 14 { 15 15 fEnergy = e; … … 39 39 Double_t GetEnergyLoss(Double_t *ep) const; 40 40 41 // ---------------------------------------------------------------- 42 41 43 MPhoton *DoInvCompton(Double_t theta); 44 Bool_t SetNewPositionB(Double_t B); 42 45 43 46 // ---------------------------------------------------------------- -
trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc
r1367 r1370 119 119 // cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush; 120 120 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); 123 123 p0 = *gamma; // Set Position and direction 124 124 p1 = *gamma; // Set Position and direction -
trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc
r1369 r1370 107 107 } 108 108 109 void MParticle::InitRandom() 110 { 111 TRandom rnd(0); 112 fPhi = rnd.Uniform(TMath::Pi()*2); 113 fPsi = rnd.Uniform(TMath::Pi()*2); 114 } 115 109 116 void MParticle::SetNewDirection(Double_t theta, Double_t phi) 110 117 { -
trunk/WuerzburgSoft/Thomas/mphys/MParticle.h
r1364 r1370 34 34 MParticle(ParticleType_t t=kENone, const char *name=NULL, const char *title=NULL); 35 35 36 void InitRandom(); 37 36 38 static Double_t ZofR(Double_t *x, Double_t *k=NULL); 37 39 static Double_t RofZ(Double_t *x, Double_t *k=NULL); … … 66 68 Bool_t SetNewPosition(); 67 69 70 ParticleType_t GetPType() const { return fPType; } 71 68 72 Double_t GetEnergy() const { return fEnergy; } 69 73 -
trunk/WuerzburgSoft/Thomas/mphys/anale.C
r1364 r1370 4 4 { 5 5 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 /* 8 11 chain.Add("cascade_100kpc_01.root"); 9 12 chain.Add("cascade_100kpc_02.root"); 10 13 chain.Add("cascade_100kpc_03.root"); 14 */ 11 15 12 16 … … 74 78 gPad->SetLogx(); 75 79 gPad->SetLogy(); 76 h.GetXaxis()->SetRangeUser( 1e4, 1e9);80 h.GetXaxis()->SetRangeUser(5e3, 1e9); 77 81 h.GetXaxis()->SetLabelOffset(-0.015); 78 82 h.GetXaxis()->SetTitleOffset(1.1); 79 h.GetXaxis()->SetRangeUser(1e4, 1e9);80 83 h.GetYaxis()->SetTitleOffset(1.1); 84 // h.SetMinimum(1e1); 85 // h.SetMaximum(1e9); 81 86 h.DrawCopy(); 82 87 TH1 *p=h.ProfileX(); 88 p->SetDirectory(NULL); 83 89 p->SetBit(kCanDelete); 84 p->SetTitle("Profile e^{-} / \\gamma ");85 90 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);93 91 p->Draw("same"); 94 92 line.DrawLine(4.2, 4.2, 8.8, 8.8); … … 102 100 r.DrawCopy(); 103 101 p=r.ProfileX(); 102 p->SetDirectory(NULL); 104 103 p->SetBit(kCanDelete); 105 104 p->SetLineColor(kRed); … … 109 108 gPad->SetLogx(); 110 109 p=h.ProjectionX(); 110 p->SetDirectory(NULL); 111 111 p->SetBit(kCanDelete); 112 112 p->SetTitle("e^{-} / \\gamma Distribution"); … … 118 118 p->Draw(); 119 119 p=h.ProjectionY(); 120 p->SetDirectory(NULL); 120 121 p->SetBit(kCanDelete); 121 122 p->SetTitle("Projection Y"); … … 126 127 gPad->SetLogx(); 127 128 p=r.ProjectionY(); 129 p->SetDirectory(NULL); 128 130 p->SetBit(kCanDelete); 129 131 p->SetTitle("Ratio E\\gamma / E_{e}"); -
trunk/WuerzburgSoft/Thomas/mphys/analp.C
r1364 r1370 1 1 //Double_t kRad2Deg = 180./TMath::Pi(); 2 3 #include <float.h> 4 Double_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 15 void 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 } 2 45 3 46 void analp() 4 47 { 5 48 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); 10 113 11 114 MPhoton *p = new MPhoton; 12 115 chain.SetBranchAddress("MPhoton.", &p); 13 116 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; 53 121 for (int i=0; i<n; i++) 54 122 { … … 63 131 { 64 132 weight = pow(Ep, alpha); 65 prim.Fill(Ep, Ep*Ep * weight); 133 prim.Fill(Ep, pow(Ep,plot) * weight); 134 E = Ep; 66 135 continue; 67 136 } 68 137 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); 72 146 } 73 147 74 148 delete p; 75 149 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 76 157 gStyle->SetOptStat(10); 158 159 TLine line; 77 160 78 161 // delete gROOT->FindObject("Analysis Arrival"); … … 85 168 gPad->SetTheta(70); 86 169 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)); 87 173 r.SetTitle(" Distance from Observer "); 88 174 r.GetXaxis()->SetLabelOffset(-0.015); … … 90 176 r.GetXaxis()->SetRangeUser(1e4, 1e9); 91 177 r.SetXTitle("\\Phi [\\circ]"); 92 r.SetYTitle("R [ light days]");178 r.SetYTitle("R [kpc]"); 93 179 r.DrawCopy("surf1pol"); 94 180 … … 103 189 a.GetXaxis()->SetRangeUser(1e4, 1e9); 104 190 a.SetXTitle("\\Psi [\\circ]"); 105 a.SetYTitle("\\Theta [ \"]");191 a.SetYTitle("\\Theta [']"); 106 192 a.DrawCopy("surf1pol"); 107 193 108 194 c->cd(3); 109 195 gPad->SetLogy(); 110 TH1 *g=r.ProjectionY("p1"); 196 g=r.ProjectionY(); 197 g->SetMinimum(GetMin(g)/2); 198 g->SetDirectory(NULL); 111 199 g->SetBit(kCanDelete); 112 200 g->SetTitle(" Hit Observers Plain "); … … 114 202 g->GetXaxis()->SetTitleOffset(1.1); 115 203 g->GetYaxis()->SetTitleOffset(1.3); 116 g->SetXTitle("R [ light days]");204 g->SetXTitle("R [kpc]"); 117 205 g->SetYTitle("Counts"); 118 206 g->Draw(); … … 120 208 c->cd(4); 121 209 gPad->SetLogy(); 122 g=a.ProjectionY("p2"); 210 g=a.ProjectionY(); 211 g->SetMinimum(GetMin(g)/2); 212 g->SetDirectory(NULL); 123 213 g->SetBit(kCanDelete); 124 214 g->SetTitle("Hit Angle"); … … 126 216 g->GetXaxis()->SetTitleOffset(1.1); 127 217 g->GetYaxis()->SetTitleOffset(1.3); 128 g->SetXTitle("\\Phi [ \"]");218 g->SetXTitle("\\Phi [']"); 129 219 g->SetYTitle("Counts"); 130 220 g->Draw(); … … 157 247 MH::SetBinning(&div, &binsx); 158 248 div.Divide(&h, &prim); 159 div.SetTitle("Spectrum / Primar aSpectrum");249 div.SetTitle("Spectrum / Primary Spectrum"); 160 250 div.GetXaxis()->SetLabelOffset(-0.015); 161 251 div.GetXaxis()->SetTitleOffset(1.1); 162 252 div.SetXTitle("E\\gamma [GeV]"); 163 253 div.SetYTitle("Ratio"); 254 line.SetLineColor(kBlue); 255 line.SetLineWidth(1); 164 256 div.DrawCopy(); 257 line.DrawLine(4, 1, 11, 1); 165 258 } -
trunk/WuerzburgSoft/Thomas/mphys/phys.C
r1367 r1370 76 76 Double_t w = log10(hi/lo)/n; 77 77 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)); 79 79 80 80 //cout << endl << w << " " << n << " " << w*bin << " " << E << endl; … … 118 118 void DoIt() 119 119 { 120 Double_t R = 100; //MParticle::RofZ(&startz); // [kpc]120 Double_t R = 500; //MParticle::RofZ(&startz); // [kpc] 121 121 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; 122 126 123 127 cout << "R = " << R << "kpc" << endl; … … 127 131 MPairProduction pair; 128 132 129 Double_t runtime = 1 5*60; // [s]133 Double_t runtime = 1; // [s] 130 134 131 135 Double_t lo = 1e4; … … 208 212 MElectron *e4file = NULL; 209 213 210 TFile file( "cascade.root", "RECREATE", "Intergalactic cascade", 9);211 TTree *T 214 TFile file(filename, "RECREATE", "Intergalactic cascade", 9); 215 TTree *T1 = new TTree("Photons", "Photons from Cascade"); 212 216 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); 215 219 216 220 // ------------------------------ … … 234 238 MPhoton *gamma=new MPhoton(E, startz); 235 239 gamma->SetIsPrimary(); 240 gamma->InitRandom(); 236 241 listg.Add(gamma); 237 242 238 B ->SetAddress(&gamma);239 T ->Fill();243 B1->SetAddress(&gamma); 244 T1->Fill(); 240 245 241 246 cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush; … … 248 253 TIter NextP(&listg); 249 254 MPhoton *p = NULL; 250 B ->SetAddress(&p);255 B1->SetAddress(&p); 251 256 while ((p=(MPhoton*)NextP())) 252 257 { … … 264 269 265 270 p->SetIsPrimary(kFALSE); 266 T ->Fill();271 T1->Fill(); 267 272 268 273 delete listg.Remove(p); … … 327 332 while(1) 328 333 { 329 if (!e->SetNewPosition ())334 if (!e->SetNewPositionB(B)) 330 335 { 331 336 cout << "!" << flush; … … 339 344 T2->Fill(); 340 345 341 if (fabs(p->GetTheta()*3437)< 1 && // < 1min346 if (fabs(p->GetTheta()*3437)<60 && // < 60min 342 347 p->GetEnergy()>lo) 343 348 { … … 354 359 } 355 360 356 if (fabs(e->GetTheta()*3437)< 1 && // < 1min357 e->GetEnergy()> lo)361 if (fabs(e->GetTheta()*3437)<60 && // < 60min 362 e->GetEnergy()>5*lo) 358 363 continue; 359 364
Note:
See TracChangeset
for help on using the changeset viewer.