Changeset 59 in svn

Timestamp:

Nov 27, 2008, 8:49:44 PM (16 years ago)

Author:

severine ovyn

Message:

remove bug

Location:

trunk

Files:

• Delphes.cpp (modified) (7 diffs)
• src/BFieldProp.cc (modified) (3 diffs)

trunk/Delphes.cpp

@@ -207 +207 @@
       while( (particle = (TRootGenParticle*) itGen.Next()) )
         {
-          genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
-          
           int pid = abs(particle->PID);
-          float eta = fabs(genMomentum.Eta());
-
-          //subarray of the quarks (i.e. not final particles, i.e status not equal to 1) 
-          // in the tracker (i.e. for those we know the tracks)
-          // with enough p_T
-          //// This subarray is needed for the B-jet algorithm
+         //// This subarray is needed for the B-jet algorithm
           // optimization for speed : put first PID condition, then ETA condition, then either pt or status
           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 ?
-              eta < DET->MAX_TRACKER &&
+              fabs(particle->Eta) < DET->MAX_TRACKER &&
               particle->Status != 1 &&
               particle->PT > DET->PT_QUARKS_MIN ) {
             NFCentralQ.Add(particle);
           }
-          
-          
+                  
           // keeps only final particles, visible by the central detector, including the fiducial volume
           // the ordering of conditions have been optimised for speed : put first the STATUS condition
@@ -230 +222 @@
           // 
           if( (particle->Status == 1)    &&
-              ( 
-               (pid == pMU  && eta < DET->MAX_MU) ||
-               (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) 
-                )
-              ) { 
-          //  TRACP->Propagation(particle,genMomentum);
-            eta = fabs(genMomentum.Eta());
+              ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
+              (fabs(particle->Eta) < DET->MAX_CALO_FWD)
+            ) 
+            { 
+            genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
+            TRACP->Propagation(particle,genMomentum);
+            float eta=fabs(genMomentum.Eta());
             switch(pid) {
               
@@ -264 +256 @@
               break;
             } // switch (pid)
-            
+
             // all final particles but muons and neutrinos      
             // for calorimetric towers and mission PT
@@ -275 +267 @@
                 // back of the input_particles vector
                 input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
-                //genMomentumCalo.SetPtEtaPhiE(CaloTower.Et(),CaloTower.iEta(),CaloTower.iPhi(),CaloTower.fourVector.E);
+                
                 genMomentumCalo.SetPxPyPzE(CaloTower.fourVector.px,CaloTower.fourVector.py,CaloTower.fourVector.pz,CaloTower.fourVector.E);
                 elementCalo = (TRootCalo*) branchCalo->NewEntry();
@@ -301 +293 @@
                 TrackCentral.push_back(genMomentum);
               }
-            
+
           } // switch
           
@@ -340 +332 @@
       elementEtmis->Px = (-PTmis).Px();
       elementEtmis->Py = (-PTmis).Py();
+
       //*****************************
+      treeWriter->Fill();
       
       sorted_jets=JETRUN->RunJets(input_particles);
@@ -346 +340 @@
       JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral);
 
-
-      treeWriter->Fill();
-
       // Add here the trigger
       // Should test all the trigger table on the event, based on reconstructed objects
-      
- //     Trigger *TRIG = new Trigger();
- //     TRIG->TriggerReader("data/trigger.dat");
      
     } // Loop over all events

trunk/src/BFieldProp.cc

@@ -45 +45 @@
   
   //out of trackibg coverage?
-  if(sqrt(pow(Xvertex1,2)+pow(Yvertex1,2)) > TRACKING_RADIUS){return;}
+  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > TRACKING_RADIUS){return;}
   if(fabs(Zvertex1) > TRACKING_LENGTH){return;}
   
@@ -55 +55 @@
   double py = Py / 0.003;
   double pz = Pz / 0.003;
-  double pt = sqrt(pow(px,2)+pow(py,2));
+  double pt = sqrt(px*px+py*py);
+  double p  = sqrt(px*px+py*py+pz*pz);
+ 
+  if(q!=0){
+     double e  = Part->E  / 0.003;
+//     double M  = sqrt(e*e - (px*px + py*py + pz*pz) );
+     double M  = Part->M;
+/*if(fabs(Part->PID)==11)
+{
+     cout<<"genMomentum.M() "<<genMomentum.M()<<"  "<<Part->M<<endl;
+     cout<<"e*e - (px*px + py*py + pz*pz) "<<e*e - (px*px + py*py + pz*pz)<<endl;
+}*/
+     double vx = px/M;
+     double vy = py/M;
+     double vz = pz/M;
 
-  double p  = sqrt(pow(px,2)+pow(py,2)+pow(pz,2));
-  
-  double m = sqrt(pow(Part->E,2) -  pow(Part->Px,2)+ pow(Part->Py,2)+ pow(Part->Pz,2));
-  
-  if(q!=0)
-    {
-      double e  = Part->E  / 0.003;
-      double M  = sqrt(e*e - (px*px + py*py + pz*pz) );
-      double vx = px/M;
-      double vy = py/M;
-      double vz = pz/M;
-      
-      /*double Bx = BField_[0];
-        double By = BField_[1];
-        double Bz = BField_[2];
-      */
-      
-      double Bx = 0;
-      double By = 0;
-      double Bz = 3.8;
-      
-      double ax =  (q/M)*(Bz*vy - By*vz);
-      double ay =  (q/M)*(Bx*vz - Bz*vx);
-      double az =  (q/M)*(By*vx - Bx*vy);
-      
-      double xold=Xvertex1;        
-      double yold=Yvertex1;        
-      double zold=Zvertex1;       
-      
-      double dt = 1/p;
-            
-      double x = xold;
-      double y = yold;        
-      double z = zold;  
-      
-      double VTold = sqrt(pow(vx,2)+pow(vy,2));
-      
-      int k = 0;
-      
-      while(k < MAXITERATION)
-        {
-          k++;
-          
-          //
-          vx += ax*dt;
-          vy += ay*dt;
-          vz += az*dt;
-          
-          double VTratio = VTold/sqrt(pow(vx,2)+pow(vy,2));
-          vx *= VTratio;
-          vy *= VTratio;
-          
-          ax = (q/M)*(Bz*vy - By*vz);
-          ay = (q/M)*(Bx*vz - Bz*vx);
-          az = (q/M)*(By*vx - Bx*vy);
-          
-          x  += vx*dt;
-          y  += vy*dt;
-          z  += vz*dt;
-          
-          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;
-            }
-          if( fabs(z) > TRACKING_LENGTH)break;
-          if( fabs((pow(x,2)+pow(y,2)+pow(z,2)) - (pow(xold,2)+pow(yold,2)+pow(zold,2))) < MINSEGLENGTH )continue;
-          
-          xold = x;
-          yold = y;
-          zold = z;
-        }       
+     double Bx = 0;
+     double By = 0;
+     double Bz = 3.8;
+
+     double ax =  (q/M)*(Bz*vy - By*vz);
+     double ay =  (q/M)*(Bx*vz - Bz*vx);
+     double az =  (q/M)*(By*vx - Bx*vy);
+/*
+cout<<"M "<<M<<" et le pid "<<Part->PID<<endl;
+cout<<"ax "<<ax<<" ay "<<ay<<" az "<<az<<endl;
+*/
+     double dt = 1/p;
+     if(pt<266 && vz < 0.0012)       dt = fabs(0.001/vz);
+ 
+     double xold=Xvertex1;         double x=xold;
+     double yold=Yvertex1;         double y=yold;
+     double zold=Zvertex1;         double z=zold;
+
+     double VTold = sqrt(vx*vx+vy*vy);
+ 
+     int k = 0;
+     
+     while(k < MAXITERATION){
+        k++;
+
+        vx += ax*dt;
+        vy += ay*dt;
+        vz += az*dt;
+
+        double VTratio = VTold/sqrt(vx*vx+vy*vy);
+        vx *= VTratio;
+        vy *= VTratio;
+
+        ax = (q/M)*(Bz*vy - By*vz);
+        ay = (q/M)*(Bx*vz - Bz*vx);
+        az = (q/M)*(By*vx - Bx*vy);
+
+        x  += vx*dt;
+        y  += vy*dt;
+        z  += vz*dt;
+
+       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;}
+       if( fabs(z)>TRACKING_LENGTH)break;
+
+       xold = x;
+       yold = y;
+       zold = z;
+     }
+ 
       if(x!=0 && y!=0 && z!=0)
         {
@@ -139 +129 @@
           genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
         }
-      cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
+     // cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
     }
 }