← Back to branch summary

~chris-rogers/maus/emr_mc_digitization

« back to all changes in this revision

Viewing changes to src/legacy/Optics/TransferMap.cc

  • Committer: Chris Rogers
  • Date: 2014-04-16 11:48:45 UTC
  • mfrom: (707 merge)
  • mto: This revision was merged to the branch mainline in revision 711.
  • Revision ID: chris.rogers@stfc.ac.uk-20140416114845-h3u3q7pdcxkxvovs
Update to trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/legacy/Optics/TransferMap.cc
136
136
Tensor  TransferMap::GetMap(int order)    const
137
137
{
138
138
        if(order > GetOrder() || order - 2 > int(_higherOrderMaps.size()) ) return Tensor(6, order+1, 0);
139
 
        if(order < 3) throw(Squeal(Squeal::recoverable, "Use GetMap for maps of order > 2", "TransferMap::GetOrder(int)"));
 
139
        if(order < 3) throw(MAUS::Exception(MAUS::Exception::recoverable, "Use GetMap for maps of order > 2", "TransferMap::GetOrder(int)"));
140
140
        return *_higherOrderMaps[order-3];
141
141
}
142
142
 
296
296
double TransferMap::PhaseAdvance(int axis) const
297
297
{
298
298
        CLHEP::HepMatrix subMatrix = GetFirstOrderMap().sub(axis*2+1, axis*2+2, axis*2+1, axis*2+2);
299
 
        double    phaseAdvance = acos( (subMatrix[0][0] + subMatrix[1][1])/2. ) ;
300
 
        if(phaseAdvance != phaseAdvance) throw(Squeal(Squeal::recoverable, "Complex phase advance", "TransferMap::PhaseAdvance")); 
301
 
        return    phaseAdvance;
 
299
  // TM given by
 
300
  // cos(\phi) + \alpha sin(\phi), \beta sin(\phi) / p
 
301
  // - \gamma sin(\phi),            cos(\phi) - \alpha sin(\phi)
 
302
  // we can get the phase advance quadrant by considering both sin(phi) and
 
303
  // cos(phi) and using atan2(phi)
 
304
  double n00 = (subMatrix[0][0]-subMatrix[1][1])/2.; // alpha sin(phi)
 
305
  double sin_phi = sqrt(-n00*n00 - subMatrix[0][1]*subMatrix[1][0]); // ((-alpha^2 + beta gamma) sin^2(phi) )^0.5
 
306
  double cos_phi = (subMatrix[0][0] + subMatrix[1][1])/2.; // cos(phi)
 
307
        double phaseAdvance = atan2(sin_phi, cos_phi);
 
308
        if(phaseAdvance != phaseAdvance) throw(MAUS::Exception(MAUS::Exception::recoverable, "Complex phase advance", "TransferMap::PhaseAdvance")); 
 
309
        return phaseAdvance;
302
310
}
303
311
 
304
312
MAUS::PolynomialMap* TransferMap::GetPolynomialMap() {