Fork me on GitHub

Changeset 59 in svn for trunk/src/BFieldProp.cc


Ignore:
Timestamp:
Nov 27, 2008, 8:49:44 PM (16 years ago)
Author:
severine ovyn
Message:

remove bug

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/BFieldProp.cc

    r53 r59  
    4545 
    4646  //out of trackibg coverage?
    47   if(sqrt(pow(Xvertex1,2)+pow(Yvertex1,2)) > TRACKING_RADIUS){return;}
     47  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > TRACKING_RADIUS){return;}
    4848  if(fabs(Zvertex1) > TRACKING_LENGTH){return;}
    4949 
    … …  
    5555  double py = Py / 0.003;
    5656  double pz = Pz / 0.003;
    57   double pt = sqrt(pow(px,2)+pow(py,2));
     57  double pt = sqrt(px*px+py*py);
     58  double p  = sqrt(px*px+py*py+pz*pz);
     59 
     60  if(q!=0){
     61     double e  = Part->E  / 0.003;
     62//     double M  = sqrt(e*e - (px*px + py*py + pz*pz) );
     63     double M  = Part->M;
     64/*if(fabs(Part->PID)==11)
     65{
     66     cout<<"genMomentum.M() "<<genMomentum.M()<<"  "<<Part->M<<endl;
     67     cout<<"e*e - (px*px + py*py + pz*pz) "<<e*e - (px*px + py*py + pz*pz)<<endl;
     68}*/
     69     double vx = px/M;
     70     double vy = py/M;
     71     double vz = pz/M;
    5872
    59   double p  = sqrt(pow(px,2)+pow(py,2)+pow(pz,2));
    60  
    61   double m = sqrt(pow(Part->E,2) -  pow(Part->Px,2)+ pow(Part->Py,2)+ pow(Part->Pz,2));
    62  
    63   if(q!=0)
    64     {
    65       double e  = Part->E  / 0.003;
    66       double M  = sqrt(e*e - (px*px + py*py + pz*pz) );
    67       double vx = px/M;
    68       double vy = py/M;
    69       double vz = pz/M;
    70      
    71       /*double Bx = BField_[0];
    72         double By = BField_[1];
    73         double Bz = BField_[2];
    74       */
    75      
    76       double Bx = 0;
    77       double By = 0;
    78       double Bz = 3.8;
    79      
    80       double ax =  (q/M)*(Bz*vy - By*vz);
    81       double ay =  (q/M)*(Bx*vz - Bz*vx);
    82       double az =  (q/M)*(By*vx - Bx*vy);
    83      
    84       double xold=Xvertex1;       
    85       double yold=Yvertex1;       
    86       double zold=Zvertex1;       
    87      
    88       double dt = 1/p;
    89            
    90       double x = xold;
    91       double y = yold;       
    92       double z = zold; 
    93      
    94       double VTold = sqrt(pow(vx,2)+pow(vy,2));
    95      
    96       int k = 0;
    97      
    98       while(k < MAXITERATION)
    99         {
    100           k++;
    101          
    102           //
    103           vx += ax*dt;
    104           vy += ay*dt;
    105           vz += az*dt;
    106          
    107           double VTratio = VTold/sqrt(pow(vx,2)+pow(vy,2));
    108           vx *= VTratio;
    109           vy *= VTratio;
    110          
    111           ax = (q/M)*(Bz*vy - By*vz);
    112           ay = (q/M)*(Bx*vz - Bz*vx);
    113           az = (q/M)*(By*vx - Bx*vy);
    114          
    115           x  += vx*dt;
    116           y  += vy*dt;
    117           z  += vz*dt;
    118          
    119           if( (x*x+y*y) > (TRACKING_RADIUS*TRACKING_RADIUS) )
    120             {
    121               x /= (x*x+y*y)/ (TRACKING_RADIUS*TRACKING_RADIUS);
    122               y /= (x*x+y*y)/ (TRACKING_RADIUS*TRACKING_RADIUS);
    123               break;
    124             }
    125           if( fabs(z) > TRACKING_LENGTH)break;
    126           if( fabs((pow(x,2)+pow(y,2)+pow(z,2)) - (pow(xold,2)+pow(yold,2)+pow(zold,2))) < MINSEGLENGTH )continue;
    127          
    128           xold = x;
    129           yold = y;
    130           zold = z;
    131         }       
     73     double Bx = 0;
     74     double By = 0;
     75     double Bz = 3.8;
     76
     77     double ax =  (q/M)*(Bz*vy - By*vz);
     78     double ay =  (q/M)*(Bx*vz - Bz*vx);
     79     double az =  (q/M)*(By*vx - Bx*vy);
     80/*
     81cout<<"M "<<M<<" et le pid "<<Part->PID<<endl;
     82cout<<"ax "<<ax<<" ay "<<ay<<" az "<<az<<endl;
     83*/
     84     double dt = 1/p;
     85     if(pt<266 && vz < 0.0012)       dt = fabs(0.001/vz);
     86 
     87     double xold=Xvertex1;         double x=xold;
     88     double yold=Yvertex1;         double y=yold;
     89     double zold=Zvertex1;         double z=zold;
     90
     91     double VTold = sqrt(vx*vx+vy*vy);
     92 
     93     int k = 0;
     94     
     95     while(k < MAXITERATION){
     96        k++;
     97
     98        vx += ax*dt;
     99        vy += ay*dt;
     100        vz += az*dt;
     101
     102        double VTratio = VTold/sqrt(vx*vx+vy*vy);
     103        vx *= VTratio;
     104        vy *= VTratio;
     105
     106        ax = (q/M)*(Bz*vy - By*vz);
     107        ay = (q/M)*(Bx*vz - Bz*vx);
     108        az = (q/M)*(By*vx - Bx*vy);
     109
     110        x  += vx*dt;
     111        y  += vy*dt;
     112        z  += vz*dt;
     113
     114       if( (x*x+y*y) > TRACKING_RADIUS*TRACKING_RADIUS ){ x /= (x*x+y*y)/(TRACKING_RADIUS*TRACKING_RADIUS); y /= (x*x+y*y)/(TRACKING_RADIUS*TRACKING_RADIUS); break;}
     115       if( fabs(z)>TRACKING_LENGTH)break;
     116
     117       xold = x;
     118       yold = y;
     119       zold = z;
     120     }
     121 
    132122      if(x!=0 && y!=0 && z!=0)
    133123        {
    … …  
    139129          genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
    140130        }
    141       cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
     131     // cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
    142132    }
    143133}
Note: See TracChangeset for help on using the changeset viewer.