Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1428) +++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1429) @@ -1,3 +1,40 @@ -*-*- END -*-*- + 2002/06/21: Thomas Bretz + + * MElectron.cc: + - changed upper limit in InteractionLength from 1e-11 to 2e-11 + - changed upper limit of p_e in EnergyLoss from -10.8 to -10.6 + + * MParticle.[h,cc]: + - replaced Planck by tanjas background + + * background.txt: + - added + + * MPhoton.[h,cc]: + - changed lower limit of Int2-Integration for Interaction Length + in case Eg<5e4 to fit better its needs + - changed max/min of drawing interaction length + - changed lower limit of interaction length plot + + * analp.C: + - added plotting 1h1426 data + - changed alpha/plot from -2/2 to -1/1 (which means I do the well + known E^2 plot of an E^-2 spectrum) + - added a cutoff to simulate cutoff spectra + - some small corrections for the layout of the plots + + * phys.C: + - added a maximum number of inv. Compton processes + - don't throw away photons which were reemitted lower than the lower + limit of the plot (so they are still written to the file) + - changed the lower limit of the PhotonSpect simulation function + to a physically more realistic value. No photon which don't have + an corresponding angle are produced anymore. + - stop the electron only if its energy is less than a thousand + of its primary energy + + + 2002/06/21: Thomas Bretz Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1428) +++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1429) @@ -34,8 +34,9 @@ { Double_t Ep = pow(10, x[0]); + Double_t res = MPhoton::Int2(&Ep, k); - return res*1e55; //65/k[0]; - // return MPhoton::Planck(&Ep, &k[1]); + + //return MPhoton::Planck(&Ep, &k[1]); } @@ -119,10 +120,38 @@ void DoIt() { - Double_t R = 100; //MParticle::RofZ(&startz); // [kpc] - Double_t startz = MParticle::ZofR(&R); - - const char *filename = "cascade_delme.root"; - - const Double_t B = 0; + // a --> 5 + // b --> 10 + // c --> 2 + + // delme, delmeB starting integration at lolim + // delme2, delme2B starting integration at -log10(Eg)-8 + // delme3, lolim, 24, 1e2-1e6, 2x inv.Compton + // delme3B, lolim, 24, 1e2-1e6, 4x inv.Compton + // delme3C, lolim, 24, 1e2-1e6, 8x inv.Compton + // delme3D, lolim, 24, 1e2-1e6, 16x inv.Compton + // delme3E, lolim, 24, 1e2-1e6, 32x inv.Compton + + // delme3H, lolim, 24, 1e2-1e6, 256x inv.Compton 8*60 + + Double_t startz = 0.03; //MParticle::ZofR(&R); + Double_t R = MParticle::RofZ(&startz); // [kpc] + + const char *filename = "delme3H.root"; + + const Double_t B = 0; // [T] mean magnetic field + + Double_t runtime = 8*60*60; // [s] maximum time to run the simulation + + Int_t nbins = 24; // number of bins produced in energy spectrum + + Double_t lo = 1e2; // lower limit of displayed spectrum + Double_t hi = 1e6; // upper limit of spectrum (cutoff) + + Int_t counter = 256; // maximum number of inv. Compton (-1 means infinite) + + Double_t alpha = -1; // -1 means a E^2 plot + Double_t plot = 1; // 1 means a E^-2 spectrum + + // -------------------------------------------------------------------- cout << "R = " << R << "kpc" << endl; @@ -131,18 +160,9 @@ MPairProduction pair; - Double_t runtime = 5*60; // [s] - - Double_t lo = 1e4; - Double_t hi = 1e11; - Double_t alpha = -2; - Double_t plot = 2; - - Int_t nbins = 21; //4*log10(hi/lo); - TH1D hist; TH1D histsrc; - hist.SetMinimum(pow(lo, plot+alpha)/10); - histsrc.SetMinimum(pow(lo, plot+alpha)/10); + hist.SetMinimum(pow(lo, plot+alpha)/100); + histsrc.SetMinimum(pow(lo, plot+alpha)/100); TList listg; @@ -230,5 +250,5 @@ { if (listg.GetSize()!=0) - cout << " |P:" << flush; + cout << " |P" << flush; TIter NextP(&listg); @@ -237,49 +257,63 @@ while ((p=(MPhoton*)NextP())) { + cout << ":" << flush; Double_t Eg = p->GetEnergy(); - if (!p->SetNewPosition()) + Double_t E0 = 511e-6; + Double_t z = p->GetZ(); + Double_t lolim = E0*E0/Eg; + Double_t inf = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1)); + if (Eg<5e4) + inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2); + + //TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2); + //Double_t E0 = 511e-6; // [GeV] + TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2); + phot.SetParameter(0, Eg); + while (1) { - cout << "!" << flush; - - hist.Fill(Eg, pow(Eg, plot)*weight); - - p->SetIsPrimary(kFALSE); - T1->Fill(); - - Esum += Eg; - - delete listg.Remove(p); - continue; + if (!p->SetNewPosition()) + { + cout << "!" << flush; + + hist.Fill(Eg, pow(Eg, plot)*weight); + + p->SetIsPrimary(kFALSE); + T1->Fill(); + + Esum += Eg; + + break; + } + + // + // Sample phtoton from background and interaction angle + // + phot.SetParameter(1, p->GetZ()); + + Double_t Ep = pow(10, phot.GetRandom()); + if (Ep==pow(10,0)) + { + cout << "z" << flush; + continue; + } + + Double_t theta = RandomTheta(Eg, Ep); + if (theta==0) + { + cout << "t" << flush; + continue; + } + + if (!pair.Process(p, Ep, theta, &liste)) + { + cout << "0" << flush; + continue; + } + + cout << "." << flush; + break; } - - // - // Sample phtoton from background and interaction angle - // - TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2); - phot.SetParameter(0, Eg); - phot.SetParameter(1, p->GetZ()); - - Double_t Ep = pow(10, phot.GetRandom()); - if (Ep==pow(10,0)) - { - cout << "z" << flush; - continue; - } - Double_t theta = RandomTheta(Eg, Ep); - if (theta==0) - { - cout << "t" << flush; - continue; - } - - if (!pair.Process(p, Ep, theta, &liste)) - { - cout << "0" << flush; - continue; - } - delete listg.Remove(p); - cout << "." << flush; } @@ -299,10 +333,14 @@ e->SetIsPrimary(kFALSE); - if (e->GetEnergy()GetEnergy(); + + if (EeSetNewPositionB(B)) { @@ -316,4 +354,7 @@ T2->Fill(); + cout << "." << flush; + listg.Add(p); + /* if (fabs(p->GetTheta()*3437)<60 && // < 60min p->GetEnergy()>lo) @@ -324,5 +365,5 @@ else { - if (p->GetEnergy()<=lo) + if (p->GetEnergy()GetEnergy()<=lo) cout << "e" << flush; else @@ -332,5 +373,5 @@ Emis += p->GetEnergy(); } - + */ if (fabs(e->GetTheta()*3437)>60) // < 60min { @@ -339,5 +380,5 @@ } - if (e->GetEnergy()GetEnergy()