Fork me on GitHub

Changeset 59 in svn


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

remove bug

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Delphes.cpp

    r55 r59  
    207207      while( (particle = (TRootGenParticle*) itGen.Next()) )
    208208        {
    209           genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
    210          
    211209          int pid = abs(particle->PID);
    212           float eta = fabs(genMomentum.Eta());
    213 
    214           //subarray of the quarks (i.e. not final particles, i.e status not equal to 1)
    215           // in the tracker (i.e. for those we know the tracks)
    216           // with enough p_T
    217           //// This subarray is needed for the B-jet algorithm
     210         //// This subarray is needed for the B-jet algorithm
    218211          // optimization for speed : put first PID condition, then ETA condition, then either pt or status
    219212          if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ?
    220               eta < DET->MAX_TRACKER &&
     213              fabs(particle->Eta) < DET->MAX_TRACKER &&
    221214              particle->Status != 1 &&
    222215              particle->PT > DET->PT_QUARKS_MIN ) {
    223216            NFCentralQ.Add(particle);
    224217          }
    225          
    226          
     218                 
    227219          // keeps only final particles, visible by the central detector, including the fiducial volume
    228220          // the ordering of conditions have been optimised for speed : put first the STATUS condition
     
    230222          //
    231223          if( (particle->Status == 1)    &&
    232               ( 
    233                (pid == pMU  && eta < DET->MAX_MU) ||
    234                (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD)
    235                 )
    236               ) {
    237           //  TRACP->Propagation(particle,genMomentum);
    238             eta = fabs(genMomentum.Eta());
     224              ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
     225              (fabs(particle->Eta) < DET->MAX_CALO_FWD)
     226            )
     227            {
     228            genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
     229            TRACP->Propagation(particle,genMomentum);
     230            float eta=fabs(genMomentum.Eta());
    239231            switch(pid) {
    240232             
     
    264256              break;
    265257            } // switch (pid)
    266            
     258
    267259            // all final particles but muons and neutrinos     
    268260            // for calorimetric towers and mission PT
     
    275267                // back of the input_particles vector
    276268                input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
    277                 //genMomentumCalo.SetPtEtaPhiE(CaloTower.Et(),CaloTower.iEta(),CaloTower.iPhi(),CaloTower.fourVector.E);
     269               
    278270                genMomentumCalo.SetPxPyPzE(CaloTower.fourVector.px,CaloTower.fourVector.py,CaloTower.fourVector.pz,CaloTower.fourVector.E);
    279271                elementCalo = (TRootCalo*) branchCalo->NewEntry();
     
    301293                TrackCentral.push_back(genMomentum);
    302294              }
    303            
     295
    304296          } // switch
    305297         
     
    340332      elementEtmis->Px = (-PTmis).Px();
    341333      elementEtmis->Py = (-PTmis).Py();
     334
    342335      //*****************************
     336      treeWriter->Fill();
    343337     
    344338      sorted_jets=JETRUN->RunJets(input_particles);
     
    346340      JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral);
    347341
    348 
    349       treeWriter->Fill();
    350 
    351342      // Add here the trigger
    352343      // Should test all the trigger table on the event, based on reconstructed objects
    353      
    354  //     Trigger *TRIG = new Trigger();
    355  //     TRIG->TriggerReader("data/trigger.dat");
    356344     
    357345    } // Loop over all events
  • 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.