Fork me on GitHub

Changeset 199 in svn


Ignore:
Timestamp:
Jan 28, 2009, 12:03:03 AM (16 years ago)
Author:
Xavier Rouby
Message:

fast bfield is working

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/interface/BFieldProp.h

    r193 r199  
    5353   /// coordinates of the exit point
    5454   double x_t, y_t, z_t, R_t, Phi_t, Theta_t, Eta_t;
     55   /// energy-momentum of the exit point [Gev] et c=1
     56   double Px_t, Py_t, Pz_t, PT_t, p_t, E_t;
    5557
    5658   unsigned int loop_overflow_counter;
  • trunk/src/BFieldProp.cc

    r193 r199  
    2222 rr(-9999.), t(-9999.), t_z(-9999.), t_T(-9999.),
    2323 x_t(-9999.), y_t(-9999.), z_t(-9999.),
    24  R_t(-9999.), Phi_t(-9999.), Theta_t(-9999.), Eta_t(-9999.) {
     24 R_t(-9999.), Phi_t(-9999.), Theta_t(-9999.), Eta_t(-9999.),
     25 Px_t(-9999), Py_t(-9999), Pz_t(-9999), PT_t(-9999), p_t(-9999), E_t(-9999) {
    2526
    2627 // if(DetDatacard="") { DET = new RESOLution(); }
     
    4748  if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
    4849
    49   if (0 && B_x== 0 && B_y== 0) { // faster if only B_z
     50  if (B_x== 0 && B_y== 0) { // faster if only B_z
    5051    if (B_z==0) return; // nothing to do
    5152
     
    6869    //     x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0)
    6970    //     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      x_c = Part->X - r*cos(phi_0);  /// TEST !!
    7172      y_c = Part->Y - r*sin(phi_0);
    7273      R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) );
     
    7576    // 3. time evaluation t = min(t_T, t_z)
    7677    //    t_T : time to exit from the sides
     78    //    t_T= [ Phi_c - phi_0 + atan( (R_max^2 - (R_c^2 + r^2))/(2rR_c) )  ]/omega
    7779    //    t_z : time to exit from the front or the back
     80    //    t_z = gamma * m /p_z0 \times (-z_0 + z_max * sign(p_z0))
    7881      rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
    7982      t_T=0;
    80       t_z = gammam / Part->Pz * (-Part->Z + DET-> TRACK_length/2. * TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ;
     83      t_z = gammam / Part->Pz * (-Part->Z + z_max* TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ;
    8184      if ( fabs(R_c - r) > R_max || R_c + r  < R_max ) t = t_z;
    8285      else {
     
    99102    //    Eta(t) = -ln tan (Theta(t)/2)
    100103      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)
     104      Phi_t = atan2( y_t, x_t);
     105/*      if(R_t>0) {
     106              Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t));
     107              Eta_t = - log(tan(Theta_t/2.));
     108      } else{
     109                Theta_t=0; Eta_t = 9999;
     110      }
     111*/
     112        Px_t = - Part->PT * sin(omega*t + phi_0);
     113        Py_t =   Part->PT * cos(omega*t + phi_0);
     114        Pz_t =   Part->Pz;
     115        PT_t =   sqrt(Px_t*Px_t + Py_t*Py_t);
     116        p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t);
     117        E_t=sqrt(Part->M*Part->M +p_t);
     118        if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.;
     119        if(p_t>0) Theta_t = acos(Pz_t/p_t);
     120        momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t);
     121
     122// test zone ---
     123/*
     124        cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl;
     125        double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.;
     126        double Eta_t2 = log( (temp_p+Pz_t)/(temp_p-Pz_t) )/2.;
     127        if(0 &&  fabs(Eta_t -Eta_t2)>1e-310) {
     128                cout << "ERROR-BUG: Eta_t != Eta_t2\n";
     129                cout << "Eta_t= " << Eta_t << "\t Eta_t1= " << Eta_t1 << "\t Eta_t2= " << Eta_t2 << endl;
     130        }
     131
     132        double R_t2  = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
     133        if(fabs(R_t - R_t2) > 1e-7)
    103134                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 ) {
     135        if( fabs(E - gammam) > 1e-3 ) {
    130136                cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
    131137                cout << "E - momentum.E() = " <<  fabs(E - momentum.E()) <<  " gammam - E " << fabs(gammam -E) << endl; }
    132138        if( fabs(PT_t - Part->PT) > 1e-10 ) {
    133                 cout << "ERROR-BUG: PT is not conserved in src/BFieldProp.cc. ";
     139                cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. ";
    134140                cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
    135141                }
     
    137143                cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
    138144
    139       //momentum.SetPtEtaPhiE(Part->PT,Eta_t,Phi_t,Part->E);
    140    
     145        double temp_p0=sqrt(Part->PT*Part->PT + Part->Pz*Part->Pz);
     146        if(fabs( (temp_p-temp_p0)*(temp_p+temp_p0) )>1e-10  ) {
     147                cout << "ERROR-BUG: momentum |vec{p}| is not conserved in src/BFieldProp.cc\n";
     148                cout << temp_p << "\t" << temp_p0 << endl;
     149        }
     150
     151   // if x_c == y_c ==0 (set it by hand!), easy cross-check
     152   //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl;
     153*/
     154
    141155  } else { // if B_x or B_y are non zero: longer computation
    142156
     
    153167  double pz = Part->Pz / 0.003;
    154168  double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
    155   double p  = sqrt(px*px+py*py+pz*pz);
     169  double p  = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
    156170 
    157171  double M  = Part->M;
Note: See TracChangeset for help on using the changeset viewer.