← Back to branch summary

~jan.greis/maus/1811

« back to all changes in this revision

Viewing changes to src/common_cpp/Recon/Kalman/MAUSSciFiPropagators.cc

  • Committer: Adam Dobbs
  • Date: 2016-01-13 11:31:12 UTC
  • mfrom: (659.2.18 release-candidate)
  • Revision ID: phuccj@gmail.com-20160113113112-bhbguupn50eyvd0z
Tags: MAUS-v1.4, MAUS-v1.4.0
MAUS-v1.4.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Recon/Kalman/MAUSSciFiPropagators.cc
252
252
//  double old_y      = old_vec(2, 0);
253
253
    double old_py     = old_vec(3, 0);
254
254
    double old_kappa  = old_vec(4, 0);
 
255
    double old_pz     = fabs(1.0 / old_kappa);
255
256
    double charge = old_kappa/fabs(old_kappa);
 
257
    double old_momentum = sqrt(old_px*old_px + old_py*old_py + old_pz*old_pz);
256
258
 
257
259
    _Bz = _geometry_helper->GetFieldValue((start.GetId() > 0 ? 1 : 0));
258
260
 
265
267
 
266
268
    if (start.GetId() < 0) u = - u;
267
269
 
 
270
    double reduction_factor = 1.0;
 
271
 
 
272
    if (_subtract_eloss) {
 
273
      SciFiMaterialsList materials;
 
274
      _geometry_helper->FillMaterialsList(start.GetId(), end.GetId(), materials);
 
275
 
 
276
      // Reduce/increase momentum vector accordingly.
 
277
      double e_loss_sign = 1.0;
 
278
      if (end.GetId() > 0) {
 
279
        e_loss_sign = -1.0;
 
280
      }
 
281
 
 
282
      double delta_p = 0.0;
 
283
      double momentum = old_momentum;
 
284
      int n_steps = materials.size();
 
285
      for (int i = 0; i < n_steps; i++) {  // In mm => times by 0.1;
 
286
//        delta_p = _geometry_helper->BetheBlochStoppingPower(momentum,
 
287
//                                               materials.at(i).first)*materials.at(i).second*0.1;
 
288
        double SP = _geometry_helper->BetheBlochStoppingPower(momentum, materials.at(i).first);
 
289
        delta_p = SP*materials.at(i).second*0.1;
 
290
        momentum += e_loss_sign*delta_p;
 
291
      }
 
292
 
 
293
      reduction_factor = momentum/old_momentum;
 
294
    }
 
295
 
268
296
    TMatrixD new_prop(5, 5);
269
297
 
270
298
    new_prop(0, 0) = 1.;
271
 
    new_prop(0, 1) = sine/u;
 
299
    new_prop(0, 1) = (sine/u);
272
300
    new_prop(0, 2) = 0.;
273
 
    new_prop(0, 3) = (cosine-1.)/u;
 
301
    new_prop(0, 3) = ((cosine-1.)/u);
274
302
 
275
303
    new_prop(1, 0) = 0.;
276
 
    new_prop(1, 1) = cosine;
 
304
    new_prop(1, 1) = cosine*reduction_factor;
277
305
    new_prop(1, 2) = 0.;
278
 
    new_prop(1, 3) = -sine;
 
306
    new_prop(1, 3) = -sine*reduction_factor;
279
307
 
280
308
    new_prop(2, 0) = 0.;
281
 
    new_prop(2, 1) = (1.-cosine)/u;
 
309
    new_prop(2, 1) = ((1.-cosine)/u);
282
310
    new_prop(2, 2) = 1.;
283
 
    new_prop(2, 3) = sine/u;
 
311
    new_prop(2, 3) = (sine/u);
284
312
 
285
313
    new_prop(3, 0) = 0.;
286
 
    new_prop(3, 1) = sine;
 
314
    new_prop(3, 1) = sine*reduction_factor;
287
315
    new_prop(3, 2) = 0.;
288
 
    new_prop(3, 3) = cosine;
 
316
    new_prop(3, 3) = cosine*reduction_factor;
289
317
 
290
318
    new_prop(4, 0) = 0.;
291
319
    new_prop(4, 1) = 0.;
292
320
    new_prop(4, 2) = 0.;
293
321
    new_prop(4, 3) = 0.;
294
 
    new_prop(4, 4) = 1.;
 
322
    new_prop(4, 4) = 1./reduction_factor;
295
323
 
296
324
    if (_correct_Pz) {
297
325
      if (start.GetId() < 0) {
298
326
        new_prop(0, 4) = delta_z*(old_px*cosine - old_py*sine);
299
 
        new_prop(1, 4) = -u*delta_z*(old_px*sine + old_py*cosine);
 
327
        new_prop(1, 4) = -u*delta_z*(old_px*sine + old_py*cosine)*reduction_factor;
300
328
        new_prop(2, 4) = delta_z*(old_px*sine + old_py*cosine);
301
 
        new_prop(3, 4) = u*delta_z*(old_px*cosine - old_py*sine);
 
329
        new_prop(3, 4) = u*delta_z*(old_px*cosine - old_py*sine)*reduction_factor;
302
330
      } else {
303
331
        new_prop(0, 4) = delta_z*(old_px*cosine - old_py*sine);
304
 
        new_prop(1, 4) = -u*delta_z*(old_px*sine + old_py*cosine);
 
332
        new_prop(1, 4) = -u*delta_z*(old_px*sine + old_py*cosine)*reduction_factor;
305
333
        new_prop(2, 4) = delta_z*(old_px*sine + old_py*cosine);
306
 
        new_prop(3, 4) = u*delta_z*(old_px*cosine - old_py*sine);
 
334
        new_prop(3, 4) = u*delta_z*(old_px*cosine - old_py*sine)*reduction_factor;
307
335
      }
308
336
    } else {
309
337
      new_prop(0, 4) = 0.0;
311
339
      new_prop(2, 4) = 0.0;
312
340
      new_prop(3, 4) = 0.0;
313
341
    }
 
342
 
314
343
//    new_prop.Print();
315
344
 
316
345
    return new_prop;