Fork me on GitHub

Changeset 193 in svn


Ignore:
Timestamp:
Jan 26, 2009, 3:43:42 PM (16 years ago)
Author:
Xavier Rouby
Message:

new algorithm inserted but not yet active

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • TabularUnified trunk/interface/BFieldProp.h

    r100 r193  
    1212 *             */
    1313
    14 #include <vector>
    1514#include "TLorentzVector.h"
    16 
    1715#include "interface/SmearUtil.h"
    18 
    1916#include "Utilities/ExRootAnalysis/interface/BlockClasses.h"
    2017#include "Utilities/ExRootAnalysis/interface/TSimpleArray.h"
     
    2623 public:
    2724  // Constructor
    28     TrackPropagation(string DetDatacard);
     25    TrackPropagation(const string DetDatacard);
    2926 
    3027    void Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum);
    3128
    32     int MAXITERATION;
    33     int MINSEGLENGTH;
    3429
    3530 private:
     31   unsigned int MAXITERATION;
     32   RESOLution *DET;
    3633
    37    RESOLution *DET; 
     34
     35   /// radial/longitudinal extensions of magnetic field volume
     36   double R_max, z_max;
     37   /// magnetic field components
     38   double B_x, B_y, B_z;
     39   /// particle charge
     40   double q;
     41   /// initial coordinate
     42   double phi_0;
     43   /// relativistic gamma \times m
     44   double gammam;
     45   /// giration frequency and radius
     46   double omega, r;
     47   /// center of the helix in the transverse plane
     48   double x_c, y_c, R_c, Phi_c;
     49   // variable for an intermediate computing
     50   double rr;
     51   /// times for exiting the magnetic field volume
     52   double t, t_z, t_T;
     53   /// coordinates of the exit point
     54   double x_t, y_t, z_t, R_t, Phi_t, Theta_t, Eta_t;
     55
     56   unsigned int loop_overflow_counter;
    3857};
    3958 
  • TabularUnified trunk/src/BFieldProp.cc

    r189 r193  
    1010
    1111#include "interface/BFieldProp.h"
    12 #include "TRandom.h"
    13 
    14 #include <iostream>
    15 #include <sstream>
    16 #include <fstream>
    17 #include <iomanip>
    18 
    1912#include<cmath>
    20 
    21 
     13#include "TMath.h"
    2214using namespace std;
    2315
     
    2517//------------------------------------------------------------------------------
    2618
    27 TrackPropagation::TrackPropagation(string DetDatacard) {
    28 
     19TrackPropagation::TrackPropagation(const string DetDatacard):
     20 MAXITERATION(10000), q(-9999.), phi_0(-9999.), gammam(-9999.), omega(-9999.), r(-9999.),
     21 x_c(-9999.), y_c(-9999.), R_c(-9999.), Phi_c(-9999.),
     22 rr(-9999.), t(-9999.), t_z(-9999.), t_T(-9999.),
     23 x_t(-9999.), y_t(-9999.), z_t(-9999.),
     24 R_t(-9999.), Phi_t(-9999.), Theta_t(-9999.), Eta_t(-9999.) {
     25
     26 // if(DetDatacard="") { DET = new RESOLution(); }
     27 // else DET = new RESOLution(DetDatacard);
    2928 DET = new RESOLution();
    3029 DET->ReadDataCard(DetDatacard);
    31  MAXITERATION = 10000;
    32  MINSEGLENGTH = 70;
    33 
     30
     31 // magnetic field parameters
     32   R_max = DET->TRACK_radius;
     33   z_max = DET->TRACK_length/2.;
     34   B_x = DET->TRACK_bfield_x;
     35   B_y = DET->TRACK_bfield_y;
     36   B_z = DET->TRACK_bfield_z;
     37
     38   loop_overflow_counter=0;
    3439}
    3540
    36 void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum)
    37 {
    38 
    39   double q  = Charge(Part->PID);
    40   if(q==0)return;
    41  
     41void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &momentum) {
     42
     43  q  = Charge(Part->PID);
     44  if(q==0) return;
     45
     46  if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
     47  if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
     48
     49  if (0 && B_x== 0 && B_y== 0) { // faster if only B_z
     50    if (B_z==0) return; // nothing to do
     51
     52    // initial conditions:
     53      // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
     54      // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z;
     55
     56    // 1.  initial transverse momentum p_{T0} : Part->PT
     57    //     initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
     58    //     relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
     59    //     giration frequency \omega = q/(gamma m) B_z
     60    //     helix radius r = p_T0 / (omega gamma m)
     61      phi_0 = -atan2(Part->Px,Part->Py);
     62      gammam = Part->E; // here c==1
     63      //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl;
     64      omega = q * B_z /gammam;
     65      r = Part->PT / (omega * gammam);
     66
     67    // 2.  Helix parameters : center coordinates in transverse plane
     68    //     x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0)
     69    //     R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c}
     70      x_c = Part->X - r*cos(phi_0);
     71      y_c = Part->Y - r*sin(phi_0);
     72      R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) );
     73      Phi_c = atan2(y_c,x_c);
     74
     75    // 3. time evaluation t = min(t_T, t_z)
     76    //    t_T : time to exit from the sides
     77    //    t_z : time to exit from the front or the back
     78      rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
     79      t_T=0;
     80      t_z = gammam / Part->Pz * (-Part->Z + DET-> TRACK_length/2. * TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ;
     81      if ( fabs(R_c - r) > R_max || R_c + r  < R_max ) t = t_z;
     82      else {
     83         t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega;
     84         t = min(t_T,t_z);
     85      }
     86
     87    // 4. position in terms of x(t), y(t), z(t)
     88    //    x(t) = x_c + r cos (omega t + phi_0)
     89    //    y(t) = y_c + r sin (omega t + phi_0)
     90    //    z(t) = z_0 + (p_Z0/gammam) t
     91      x_t = x_c + r * cos(omega * t + phi_0);
     92      y_t = y_c + r * sin(omega * t + phi_0);
     93      z_t = Part->Z + Part->Pz / gammam * t;
     94
     95    // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
     96    //    R(t) = sqrt(x(t)² + y(t)²)
     97    //    Phi(t) = atan(y(t)/x(t))
     98    //    Theta(t) = atan(R(t)/z(t))
     99    //    Eta(t) = -ln tan (Theta(t)/2)
     100      R_t   = sqrt( pow(x_t,2.) + pow(y_t,2.)  );
     101      double R_t2  = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
     102      if(fabs(R_t - R_t2) > 1e-11)
     103                cout << "ERROR-BUG: R_t != R_t2:  R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl;
     104      Phi_t = atan2( y_t, x_t);
     105      Theta_t = atan2( R_t , z_t);
     106      Eta_t = - log(tan(Theta_t/2.));
     107
     108        double Px_t = - Part->PT * sin(omega*t + phi_0);
     109        double Py_t =   Part->PT * cos(omega*t + phi_0);
     110        double Pz_t =   Part->Pz;
     111        double PT_t =   sqrt(Px_t*Px_t + Py_t*Py_t);
     112        double E=sqrt(Part->M*Part->M + PT_t*PT_t + Pz_t*Pz_t); // cross-check
     113        momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E);
     114        TLorentzVector mom_temp;
     115        mom_temp.SetPtEtaPhiM(PT_t,Eta_t,Phi_t+PI,Part->M);
     116
     117        if (momentum != mom_temp) {
     118//              cout << "momentum et mom_temp different\n";
     119/*              if (momentum.Px() != mom_temp.Px() ) cout << "  px: " << momentum.Px() << " & " << mom_temp.Px() << endl;
     120                if (momentum.Py() != mom_temp.Py() ) cout << "  py: " << momentum.Py() << " & " << mom_temp.Py() << endl;
     121                if (momentum.Pz() != mom_temp.Pz() ) cout << "  pz: " << momentum.Pz() << " & " << mom_temp.Pz() << endl;
     122                if (momentum.Pt() != mom_temp.Pt() ) cout << "  pt: " << momentum.Pt() << " & " << mom_temp.Pt() << endl;*/
     123//              cout << "eta: " << momentum.Eta() << " " << mom_temp.Eta() << endl;
     124//              cout << "the: " << momentum.Theta() << " " << mom_temp.Theta() << endl;
     125//              cout << "phi: " << momentum.Phi() << " " << mom_temp.Phi() << endl;
     126        }
     127        //cout << Part->PT << " or " << momentum.Pt() << " or " << mom_temp.Pt() << endl;
     128
     129        if( fabs(E - gammam) > 1e-4 ) {
     130                cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
     131                cout << "E - momentum.E() = " <<  fabs(E - momentum.E()) <<  " gammam - E " << fabs(gammam -E) << endl; }
     132        if( fabs(PT_t - Part->PT) > 1e-10 ) {
     133                cout << "ERROR-BUG: PT is not conserved in src/BFieldProp.cc. ";
     134                cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
     135                }
     136        if(momentum.Pz() != Pz_t)
     137                cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
     138
     139      //momentum.SetPtEtaPhiE(Part->PT,Eta_t,Phi_t,Part->E);
     140   
     141  } else { // if B_x or B_y are non zero: longer computation
     142
    42143  float Xvertex1 = Part->X;
    43144  float Yvertex1 = Part->Y;
    44145  float Zvertex1 = Part->Z;
    45146 
    46   //out of trackibg coverage?
    47   if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > DET->TRACK_radius){return;}
    48   if(fabs(Zvertex1) > DET->TRACK_length/2.){return;}
     147  //out of tracking coverage?
     148  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;}
     149  if(fabs(Zvertex1) > z_max){return;}
    49150 
    50   float Px = Part->Px;
    51   float Py = Part->Py;
    52   float Pz = Part->Pz;
    53    
    54   double px = Px / 0.003;
    55   double py = Py / 0.003;
    56   double pz = Pz / 0.003;
    57   double pt = sqrt(px*px+py*py);
     151  double px = Part->Px / 0.003;
     152  double py = Part->Py / 0.003;
     153  double pz = Part->Pz / 0.003;
     154  double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
    58155  double p  = sqrt(px*px+py*py+pz*pz);
    59156 
    60   if(q!=0){
    61      double M  = Part->M;
    62      double vx = px/M;
    63      double vy = py/M;
    64      double vz = pz/M;
    65 
    66      double Bx = DET->TRACK_bfield_x;
    67      double By = DET->TRACK_bfield_y;
    68      double Bz = DET->TRACK_bfield_z;
    69 
    70      double ax =  (q/M)*(Bz*vy - By*vz);
    71      double ay =  (q/M)*(Bx*vz - Bz*vx);
    72      double az =  (q/M)*(By*vx - Bx*vy);
    73      double dt = 1/p;
    74      if(pt<266 && vz < 0.0012)       dt = fabs(0.001/vz);
    75  
    76      double xold=Xvertex1;         double x=xold;
    77      double yold=Yvertex1;         double y=yold;
    78      double zold=Zvertex1;         double z=zold;
    79 
    80      double VTold = sqrt(vx*vx+vy*vy);
    81  
    82      int k = 0;
    83      
    84      double Radius=DET->TRACK_radius;
    85      double Length=DET->TRACK_length/2.;
    86 
    87      while(k < MAXITERATION){
     157  double M  = Part->M;
     158  double vx = px/M;
     159  double vy = py/M;
     160  double vz = pz/M;
     161  double qm = q/M;
     162
     163  double ax =  qm*(B_z*vy - B_y*vz);
     164  double ay =  qm*(B_x*vz - B_z*vx);
     165  double az =  qm*(B_y*vx - B_x*vy);
     166  double dt = 1/p;
     167  if(pt<266 && vz < 0.0012)       dt = fabs(0.001/vz); // ?????
     168 
     169  double xold=Xvertex1;         double x=xold;
     170  double yold=Yvertex1;         double y=yold;
     171  double zold=Zvertex1;         double z=zold;
     172
     173  double VTold = pt/M; //=sqrt(vx*vx+vy*vy);
     174 
     175  unsigned int k = 0;
     176  double VTratio=0;
     177  double R_max2 = R_max*R_max;
     178  double r2=0; // will be x*x+y*y
     179
     180  while(k < MAXITERATION){
    88181        k++;
    89182
     
    92185        vz += az*dt;
    93186
    94         double VTratio = VTold/sqrt(vx*vx+vy*vy);
     187        VTratio = VTold/sqrt(vx*vx+vy*vy);
    95188        vx *= VTratio;
    96189        vy *= VTratio;
    97190
    98         ax = (q/M)*(Bz*vy - By*vz);
    99         ay = (q/M)*(Bx*vz - Bz*vx);
    100         az = (q/M)*(By*vx - Bx*vy);
     191        ax = qm*(B_z*vy - B_y*vz);
     192        ay = qm*(B_x*vz - B_z*vx);
     193        az = qm*(B_y*vx - B_x*vy);
    101194
    102195        x  += vx*dt;
    103196        y  += vy*dt;
    104197        z  += vz*dt;
    105 
    106        if( (x*x+y*y) > Radius*Radius ){ x /= (x*x+y*y)/(Radius*Radius); y /= (x*x+y*y)/(Radius*Radius); break;}
    107        if( fabs(z)>Length)break;
     198        r2  = x*x + y*y;
     199
     200       if( r2 > R_max2 ){
     201                x /= r2/R_max2;
     202                y /= r2/R_max2;
     203                break;
     204       }
     205       if( fabs(z)>z_max)break;
    108206
    109207       xold = x;
    110208       yold = y;
    111209       zold = z;
    112      }
    113  
    114       if(x!=0 && y!=0 && z!=0)
    115         {
    116           double eta;
    117           float Theta = atan2(sqrt(x*x+y*y),z);
    118           eta  = -log(tan(Theta/2));
     210  } // while loop
     211
     212  if(k == MAXITERATION) loop_overflow_counter++;
     213  //cout << "too short loop in " << loop_overflow_counter << " cases" << endl;
     214 
     215  if(x!=0 && y!=0 && z!=0) {
     216          float Theta = atan2(sqrt(r2),z);
     217          double eta  = -log(tan(Theta/2.));
    119218          double phi = atan2(y,x);
    120           genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
    121         }
    122     }
     219          momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
     220  }
     221   
     222 } // if b_x or b_y non zero
    123223}
Note: See TracChangeset for help on using the changeset viewer.