Merge pull request #28 from brownd1978/tweak

Tweak

Data/Mu2e2Tracker.dat

@@ -1,7 +1,7 @@
 #
 #  Parameters describing the Mu2e-ii tracker
 # Volume information: rmin rmax zpos zhalf
-365.0 700.0 0.0 1500.0
+365.0 700.0 0.0 1450.0
 # Cell information: orientation (0=azimuthal, 1=axial), Ncells, cell radius, cell length, wall thickness, wire radius, minimum path
 0 20736 2.5 800.0 0.008 0.0125 100
 # material information: wall, gas, and wire

Data/Mu2eTracker.dat

@@ -1,7 +1,7 @@
 #
 #  Parameters describing the Mu2e tracker
 # Volume information: rmin rmax zpos zhalf
-365.0 700.0 0.0 1500.0
+365.0 700.0 0.0 1450.0
 # Cell information: orientation (0=azimuthal, 1=axial), Ncells, cell radius, cell length, wall thickness, wire radius, minimum path
 0 20736 2.5 800.0 0.015 0.0125 100
 # material information: wall, gas, and wire

Detector/CylindricalShell.hh

@@ -62,7 +62,8 @@ namespace TrackToy {
           double tx = ttest - tstep*fabs(dr/(dr-olddr));
           // compute the crossing time range
           auto vel = pktraj.velocity(tx);
-          double dt = rhalf_/vel.Rho();
+          double vr = vel.Rho();
+          double dt = vr > 1e-8 ?  rhalf()/vr : 2.0*sqrt(2.0*radius()*rhalf())/vel.R();
           trange = TimeRange(tx-dt,tx+dt);
           break;
         }

Detector/HollowCylinder.hh

@@ -55,8 +55,7 @@ namespace TrackToy {
     double ttest = tstart;
     auto pos = pktraj.position3(ttest);
     bool inside = isInside(pos);
-    TimeRange trange;
-    if(inside) trange = TimeRange(ttest,ttest);
+    TimeRange trange(tstart,tstart);
     while(ttest < pktraj.range().end()){
       //      cout << "particle enters at " << pos << endl;
       auto vel = pktraj.velocity(ttest);

Detector/IPA.hh

@@ -32,6 +32,7 @@ namespace TrackToy {
   template<class PKTRAJ> bool IPA::extendTrajectory(KinKal::BFieldMap const& bfield, PKTRAJ& pktraj, TimeRanges& intersections,double tol) const {
     using KinKal::VEC3;
     using KinKal::TimeRange;
+    static unsigned moyalterms_(20); // number of terms in Moyal expansion
     intersections.clear();
     // record the end of the previous extension; this is where new extensions start
     double tstart = pktraj.back().range().begin();
@@ -47,18 +48,21 @@ namespace TrackToy {
       while( (!trange.null()) && trange.end() < pktraj.range().end()) {
         intersections.push_back(trange);
         double energy = pktraj.energy(trange.mid());
+        auto momvec = pktraj.momentum3(trange.mid());
+        double mom = momvec.R();
         double plen = pktraj.speed(trange.mid())*trange.range();
         // Moyal dist. models ionization loss
-        double demean = mat_->energyLoss(energy,plen,pktraj.mass());
-        double derms = mat_->energyLossRMS(energy,plen,pktraj.mass());
+        double demean = mat_->energyLoss(mom,plen,pktraj.mass());
+        double derms = mat_->energyLossRMS(mom,plen,pktraj.mass());
         // model ionization energy loss using a Moyal distribution
-        MoyalDist edist(MoyalDist::MeanRMS(fabs(demean), derms),10);
+//        std::cout << "IPA demean " << demean << " derms " << derms << std::endl;
+        MoyalDist edist(MoyalDist::MeanRMS(fabs(demean), derms),moyalterms_);
         double ionloss = edist.sample(tr_.Uniform(0.0,1.0));
         // add radiative energy loss.  note we have to convert to cm!!!
         double radFrac = mat_->radiationFraction(trange.range())/10;
         BremssLoss bLoss;
         double bremloss = bLoss.sampleSSPGamma(energy,radFrac);
-        // delta energy loss
+        // delta energy loss; note unit change mm->cm!
         DeltaRayLoss dLoss(mat_, energy,plen/10, pktraj.mass());
         double dloss = dLoss.sampleDRL();
         double totloss = ionloss + bremloss + dloss;
@@ -67,8 +71,6 @@ namespace TrackToy {
         energy -= totloss;
         //          std::cout << "old energy " << oldenergy << " new energy " << energy << std::endl;
         // scattering
-        auto momvec = pktraj.momentum3(trange.mid());
-        double mom = momvec.R();
         double scatterRMS = mat_->scatterAngleRMS(mom,plen,pktraj.mass());
         // generate random momentum scatter
         VEC3 phidir = VEC3(momvec.Y(),-momvec.X(),0.0).Unit();

Detector/Target.hh

@@ -45,9 +45,10 @@ namespace TrackToy {
   template<class PKTRAJ> bool Target::extendTrajectory(KinKal::BFieldMap const& bfield, PKTRAJ& pktraj, TimeRanges& intersections,double tol) const {
     using KinKal::VEC3;
     using KinKal::TimeRange;
+    static unsigned moyalterms_(20); // number of terms in Moyal expansion
     intersections.clear();
     static double pfactor = 0.001*density_/mat_->density(); // unit conversion cm->mm, and scale for the effective density
-//    std::cout << "density factor" << pfactor << std::endl;
+    //    std::cout << "density factor" << pfactor << std::endl;
     // extend to the  of the target or exiting the BField (backwards)
     bool retval = extendZ(pktraj,bfield, cyl_.zmax(), tol);
     //    cout << "Z target extend " << ztgt << endl;
@@ -58,40 +59,41 @@ namespace TrackToy {
       TimeRange trange = cyl_.intersect(pktraj,tstart,tstep);
       while( (!trange.null()) && trange.end() < pktraj.range().end()) {
         double speed = pktraj.speed(trange.mid());
-        double plen = trange.range()*speed;
-        if(plen > minpath_){
-          intersections.push_back(trange);
-          double energy = pktraj.energy(trange.mid());
-          // to get physical results, scale the path
-          plen *= pfactor;
-          double demean = mat_->energyLoss(energy,plen,pktraj.mass());
-          double derms = mat_->energyLossRMS(energy,plen,pktraj.mass());
-          MoyalDist edist(MoyalDist::MeanRMS(fabs(demean),derms),10);
-          double ionloss = edist.sample(tr_.Uniform(0.0,1.0));
-          // add radiative energy loss.  note we have to convert to cm!!!
-          double radFrac = mat_->radiationFraction(trange.range())/10;
-          BremssLoss bLoss;
-          double bremloss = bLoss.sampleSSPGamma(energy,radFrac);
-          // delta energy loss
-          DeltaRayLoss dLoss(mat_, energy,plen/10, pktraj.mass());
-          double dloss = dLoss.sampleDRL();
-          double totloss = ionloss + bremloss + dloss;
-          //        std::cout << "Target Ionization eloss = " << ionloss << " Delta eloss " << dloss << " rad eloss "  << bremloss << " tot " << totloss << std::endl;
-          energy -= totloss;
-          // scattering
-          auto momvec = pktraj.momentum3(trange.mid());
-          double mom = momvec.R();
-          double scatterRMS = mat_->scatterAngleRMS(mom,plen,pktraj.mass());
-//          std::cout << "scatterRMS " << scatterRMS << std::endl;
-          // generate random momentum scatter
-          VEC3 phidir = VEC3(momvec.Y(),-momvec.X(),0.0).Unit();
-          double momtan = tan(momvec.Theta());
-          VEC3 thedir = VEC3(momvec.X()/momtan,momvec.Y()/momtan,-momvec.Z()*momtan).Unit();
-          auto dmom = tr_.Gaus(0.0,scatterRMS)*mom*phidir + tr_.Gaus(0.0,scatterRMS)*mom*thedir;
-//          std::cout << "phidot " << momvec.Dot(dmom) << " dmag = " << mom - (momvec + dmom).R() << std::endl;
-          retval = updateEnergy(pktraj,trange.mid(),energy,dmom);
-          if(!retval)break;
-        }
+        double plen = std::max(trange.range()*speed, minpath_);
+        // require a physical minimum
+        intersections.push_back(trange);
+        double energy = pktraj.energy(trange.mid());
+        auto momvec = pktraj.momentum3(trange.mid());
+        double mom = momvec.R();
+        // to get physical results, scale the path
+//        double estarde = electronEnergyLoss(energy-pktraj.mass(),plen);
+        plen *= pfactor;
+        double demean = mat_->energyLoss(mom,plen,pktraj.mass());
+        double derms = mat_->energyLossRMS(mom,plen,pktraj.mass());
+        MoyalDist edist(MoyalDist::MeanRMS(fabs(demean),derms),moyalterms_);
+        double ionloss = edist.sample(tr_.Uniform(0.0,1.0));
+        // add radiative energy loss.  note we have to convert to cm!!!
+        double radFrac = mat_->radiationFraction(trange.range())/10;
+        BremssLoss bLoss;
+        double bremloss = bLoss.sampleSSPGamma(energy,radFrac);
+        // delta energy loss
+        DeltaRayLoss dLoss(mat_, energy,plen/10, pktraj.mass());
+        double dloss = dLoss.sampleDRL();
+        double totloss = ionloss + bremloss + dloss;
+        //        std::cout << "Target Ionization eloss = " << ionloss << " Delta eloss " << dloss << " rad eloss "  << bremloss << " tot " << totloss << std::endl;
+        energy -= totloss;
+//        std::cout << "Target demean " << demean << " derms " << derms << " estar de " << estarde << " totlos " << totloss << std::endl;
+        // scattering
+        double scatterRMS = mat_->scatterAngleRMS(mom,plen,pktraj.mass());
+        //          std::cout << "scatterRMS " << scatterRMS << std::endl;
+        // generate random momentum scatter
+        VEC3 phidir = VEC3(momvec.Y(),-momvec.X(),0.0).Unit();
+        double momtan = tan(momvec.Theta());
+        VEC3 thedir = VEC3(momvec.X()/momtan,momvec.Y()/momtan,-momvec.Z()*momtan).Unit();
+        auto dmom = tr_.Gaus(0.0,scatterRMS)*mom*phidir + tr_.Gaus(0.0,scatterRMS)*mom*thedir;
+        //          std::cout << "phidot " << momvec.Dot(dmom) << " dmag = " << mom - (momvec + dmom).R() << std::endl;
+        retval = updateEnergy(pktraj,trange.mid(),energy,dmom);
+        if(!retval)break;
         trange = cyl_.intersect(pktraj,trange.end(),tstep);
       }
     }

Detector/Tracker.hh

@@ -81,19 +81,19 @@ namespace TrackToy {
       std::vector<std::shared_ptr<KinKal::Hit<KTRAJ>>>& hits,
       std::vector<std::shared_ptr<KinKal::ElementXing<KTRAJ>>>& xings,
       std::vector<KinKal::TimeRange>& tinters, double tol) const {
+    using KinKal::TimeRange;
     double tstart = mctraj.back().range().begin();
     double speed = mctraj.speed(tstart);
     double tstep = cellRadius()/speed;
     // find intersections with tracker
-    cylinder().intersect(mctraj,tinters,tstart,tstep);
-    //    std::cout << "ninters " << tinters.size() << std::endl;
-    for(auto const& tinter : tinters) {
-      double clen = tinter.range()*speed;
-      //
+    TimeRange trange = cylinder().intersect(mctraj,tstart,tstep);
+    while( (!trange.null()) && trange.end() < mctraj.range().end()) {
+      double clen = trange.range()*speed;
       if(clen > minpath_){
+        tinters.push_back(trange);
         unsigned ncells = (unsigned)floor(clen*cellDensity_);
-        double hstep = tinter.range()/(ncells+1);
-        double htime = tinter.begin()+0.5*hstep;
+        double hstep = trange.range()/(ncells+1);
+        double htime = trange.begin()+0.5*hstep;
         for(unsigned icell=0;icell<ncells;++icell){
           // extend the trajectory to this time
           extendTraj(bfield,mctraj,htime,tol);
@@ -107,6 +107,7 @@ namespace TrackToy {
           htime += hstep;
         }
       }
+      trange = cylinder().intersect(mctraj,trange.end(),tstep);
     }
   }
 
@@ -127,7 +128,7 @@ namespace TrackToy {
       static double pi_over_3 = M_PI/3.0;
       // generate azimuthal angle WRT the hit
       double wphi = tr_.Uniform(0.0,pi_over_3);
-// acceptance test, using a triangular inner hole
+      // acceptance test, using a triangular inner hole
       double rmax = 0.5*cyl_.rmax()/cos(wphi);
       if(pos.Rho() < rmax)return false;
       auto rdir = PerpVector(pos,zdir).Unit(); // radial direction
@@ -214,6 +215,7 @@ namespace TrackToy {
 
   template <class KTRAJ> void Tracker::updateTraj(KinKal::BFieldMap const& bfield,
       KinKal::ParticleTrajectory<KTRAJ>& mctraj, const KinKal::ElementXing<KTRAJ>* sxing) const {
+    static unsigned moyalterms_(20); // number of terms in Moyal expansion
     // simulate energy loss and multiple scattering from this xing
     auto txing = sxing->crossingTime();
     auto const& endpiece = mctraj.nearestPiece(txing);
@@ -225,7 +227,7 @@ namespace TrackToy {
     sxing->materialEffects(mctraj,KinKal::TimeDir::forwards, dmom, momvar);
 //   std::cout << "DMOM " << dmom[0] << std::endl;
    // note materialEffects returns the normalized energy change
-    MoyalDist edist(MoyalDist::MeanRMS(fabs(dmom[KinKal::MomBasis::momdir_])*mom,sqrt(momvar[KinKal::MomBasis::momdir_])),10);
+    MoyalDist edist(MoyalDist::MeanRMS(fabs(dmom[KinKal::MomBasis::momdir_])*mom,sqrt(momvar[KinKal::MomBasis::momdir_])),moyalterms_);
     // radiation energy loss model
     BremssLoss bLoss;
     for(int idir=0;idir<=KinKal::MomBasis::phidir_; idir++) {
@@ -234,9 +236,8 @@ namespace TrackToy {
       double ionloss, bremloss, dloss;
       double dm;
       double radFrac = sxing->radiationFraction()/10; // convert to cm
-      // only include wall material for now TODO
+      // only include wall material for now; gas can be added later TODO
       DeltaRayLoss dLoss(&(sxing->matXings()[0].dmat_), mom,sxing->matXings()[0].plen_/10.0, endpiece.mass());
-      dLoss.setCutOffEnergy(0.001); // 1 KeV
       // generate a random effect given this variance and mean.  Note momEffect is scaled to momentum
       switch( mdir ) {
         case KinKal::MomBasis::perpdir_: case KinKal::MomBasis::phidir_ :

General/ELossDistributions.cc

@@ -45,7 +45,7 @@ namespace TrackToy {
     }
 
     double y = sum;
-    double x = _mode - 2 * _sigma * std::log ((std::sqrt(2.) * y));
+    double x = _mode - 2 * _sigma * std::log (M_SQRT2 * y);
     return x;
   }
 

General/ELossDistributions.hh

@@ -26,7 +26,7 @@ namespace TrackToy {
 
       MoyalDist(MeanRMS const& meanrms, int max = 20):_mean(meanrms.Mean_),_rms(meanrms.RMS_),_kmax(max) {
         //Variance of moyal = (pi * sigma)^2 /2
-        _sigma = std::sqrt(2.0) * _rms / M_PI ;
+        _sigma = M_SQRT2 * _rms / M_PI ;
         _mode = _mean - _sigma * MFACTOR;
         setCoeffs(_kmax);
       }

Mains/CMakeLists.txt

@@ -6,8 +6,8 @@ set( MAIN_FILES
   BField_test.cc
   MuStops_test.cc
   Tracks_test.cc
-  ELossDists_unit.cc
-  MoyalDist_unit.cc
+  ELossDists_test.cc
+  MoyalDist_test.cc
 )
 # Generate test targets
 foreach( mainfile ${MAIN_FILES} )

Tests/MuStops.C

@@ -6,7 +6,7 @@
 #include "TCanvas.h"
 #include "TStyle.h"
 #include "TLegend.h"
-void MuStop() {
+void MuStops() {
   TFile *mubeam = TFile::Open("/Users/brownd/data/sim.mu2e.MuBeamCat.MDC2020k.001201_00000000.art");
   TTree* bevents = (TTree*)mubeam->Get("Events");
   TH1F* mumom = new TH1F("mumom","Muon momentum;Momentum (MeV/c)",120,0,120.0);