Changeset 405 in svn for trunk/src/VeryForward.cc

Timestamp:

May 20, 2009, 10:39:01 AM (15 years ago)

Author:

Xavier Rouby

Message:

bug removed in Hector implementation for RP220/FP420

File:

• trunk/src/VeryForward.cc (modified) (13 diffs)

trunk/src/VeryForward.cc

@@ -31 +31 @@
 
 #include "VeryForward.h"
+#include "PdgParticle.h"
 #include "H_RomanPot.h"
-#include "PdgParticle.h"
+
 #include <iostream>
+#include <fstream>
 #include<cmath>
 
 using namespace std;
 
+/* Notes on the correct initialisation for Hector
+ * -- these notes apply to the LHC beamlines
+ *
+ *  beam1 : forward direction is for increasing 's' values
+ *    beamline1 = new H_BeamLine(1,...);
+ *    beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);
+ *  beam2 : forward direction is for decreasing 's' values
+ *    beamline2 = new H_BeamLine(-1,...);
+ *    beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);
+ *
+ *  relative_energy should be false -- kickers_on should be 1
+ *
+ */
 
 //------------------------------------------------------------------------------
 VeryForward::VeryForward() : 
    DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)), 
-   beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
-   relative_energy(true), // should always be true
-   kickers_on(1)         // should always be 1
+   beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(-1,d_max)),
+   rel_energy(true), // should always be true
+   kickers(1)         // should always be 1
    {
    init(); //Initialisation of Hector
@@ -55 +70 @@
    const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s);
    beamline1 = new H_BeamLine(1,d_max);
-   beamline2 = new H_BeamLine(1,d_max);
+   beamline2 = new H_BeamLine(-1,d_max);
    init(); //Initialisation of Hector
-   relative_energy = true; // should always be true
-   kickers_on = 1;         // should always be 1
+   rel_energy = true; // should always be true
+   kickers = 1;         // should always be 1
 }
 
 VeryForward::VeryForward(const RESOLution * DetDatacard) :
         DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
-        beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
-        relative_energy(true), // should always be true
-        kickers_on(1)          // should always be 1
+        beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(-1,d_max)),
+        rel_energy(true), // should always be true
+        kickers(1)          // should always be 1
    {
    init();  //Initialisation of Hector
@@ -73 +88 @@
    DET(new RESOLution(*(vf.DET))),   d_max(vf.d_max),
    beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))),
-   relative_energy(vf.relative_energy),
-   kickers_on(vf.kickers_on) {
+   rel_energy(vf.rel_energy),
+   kickers(vf.kickers) {
 }
 
@@ -83 +98 @@
    beamline1 = new H_BeamLine(*(vf.beamline1));
    beamline2 = new H_BeamLine(*(vf.beamline2));
-   relative_energy =vf.relative_energy;
-   kickers_on = vf.kickers_on;
+   rel_energy =vf.rel_energy;
+   kickers = vf.kickers;
    return *this;
 }
@@ -93 +108 @@
   static unsigned int counter;
   counter =0;
-  relative_energy = true; // should always be true
-  kickers_on = 1;         // should always be 1
+  extern bool relative_energy;
+  extern int kickers_on;
+  relative_energy = rel_energy; // should always be true
+  kickers_on = kickers;         // should always be 1
   beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);                                
-  beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x);
-  H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*1E6); // RP 220m, 2mm, beam 1  
-  H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*1E6); // RP 420m, 4mm, beam 1  
+  beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
+  H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*1E6); 
+  // RP 220m, 2mm, beam 1  
+  H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*1E6); 
+  // RP 420m, 4mm, beam 1  
+  //rp220_1->printProperties();
+  //rp420_1->printProperties();
   beamline1->add(rp220_1); 
   beamline1->add(rp420_1);
   
   beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);                              
-  beamline2->offsetElements(DET->RP_offsetEl_s,+DET->RP_offsetEl_x);
-  H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6);// RP 220m, 2mm, beam 2     
-  H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6);// RP 420m, 4mm, beam 2     
+  beamline2->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
+  H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6);
+  // RP 220m, 2mm, beam 2     
+  H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6);
+  // RP 420m, 4mm, beam 2     
+  //rp220_2->printProperties();
+  //rp420_2->printProperties();
   beamline2->add(rp220_2);
   beamline2->add(rp420_2); 
   // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine
   // do not put explicit delete
+}
+
+
+float VeryForward::time_of_flight(TRootGenParticle *particle, const float detector_s, const float detector_etamin, const float detector_t_resolution) {
+    // time of flight t is t = T + d/[ cos(theta) v ]
+    float cos_theta = 1;     //very good approximation, if detector_etamin >3
+    if (detector_etamin<3) { // if smaller eta -> make the complete calculation
+       double tx = atan(particle->Px/particle->Pz);
+       double ty = atan(particle->Py/particle->Pz);
+       double theta = sqrt( pow(tx,2) + pow(ty,2) );
+       // cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
+       // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
+       // eta = 2.6 <-> cos(theta) = 0.99
+       // eta = 3.0 <-> cos(theta) = 0.995
+       cos_theta = cos(theta);
+    }
+    // units from StdHEP : Z [mm] T[mm/c]
+    // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
+    double flight_distance = (detector_s - particle->Z*(1E-3))/cos_theta ;
+    // assumes highly relativistic particles
+    double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; 
+    double timeS = gRandom->Gaus(flight_time,detector_t_resolution);
+    return timeS;
 }
 
@@ -125 +173 @@
     elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
 
-
-    elementZdc->pid = particle->PID;
-
-     
-    // for compatibility with 'old' version
     TLorentzVector genMomentum;
     genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
-    elementZdc->Set(genMomentum); 
-    // ******************
-    
- 
-    //particle->print();
+    elementZdc->Set(genMomentum);  // initialises the gen-level data
+    elementZdc->pid = particle->PID;
  
     // 1) energy smearing
@@ -149 +189 @@
                                                     pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
     elementZdc->E = energyS;
- 
-
-    // 2) time of flight t is t = T + d/[ cos(theta) v ]
-    float cos_theta = 1;      //very good approximation, if eta_zdc >3
-    if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation
-       double tx = atan(particle->Px/particle->Pz);
-       double ty = atan(particle->Py/particle->Pz);
-       double theta = sqrt( pow(tx,2) + pow(ty,2) );
-       //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
-       // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
-       // eta = 2.6 <-> cos(theta) = 0.99
-       // eta = 3.0 <-> cos(theta) = 0.995
-       cos_theta = cos(theta);
-    }
-    // units from StdHEP : Z [mm] T[mm/c] 
-    // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
-    double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ;
-    double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
-    double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution);
-    elementZdc->T = timeS;
+
+    // 2) time of flight t is t = T + d/[ cos(theta) v ] + detector smearing on time
+    elementZdc->T = time_of_flight(particle, DET->VFD_s_zdc, DET->VFD_min_zdc, DET->ZDC_T_resolution);
 
     // 3) side: which ZDC has been hit? 
@@ -174 +197 @@
 
     // 4) object nature : e.m. (photon) or had (neutron) ?
-    //elementZdc->hadronic_hit = (bool) (particle->PID==pN);
-  }
+    elementZdc->hadronic_hit = (bool) (particle->PID!=pGAMMA);
+  } // if neutrons or photons over E_threshold
+
 }
 
@@ -181 +205 @@
 void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle) 
 {
+  if(particle->Status != 1) return; // reject particles that are not final ones
+  extern bool relative_energy;
+  relative_energy = rel_energy;
+  extern int kickers_on;
+  kickers_on = kickers;
+
   float charge = particle->Charge, mass = particle->M;
+  //float charge, mass, ctau;
+  //charge = mass = ctau = UNDEFINED;
   if (mass<-999) { // unitialised!
           PdgParticle pdg_part = DET->PDGtable[particle->PID];
-          charge  = pdg_part.charge();
-          mass    = pdg_part.mass();
+          charge  = pdg_part.charge();  // e+
+          mass    = pdg_part.mass();    // GeV
+  //        ctau    = pdg_part.ctau();    // m
+  // cout << "ctau = " << ctau << endl;
   }
-  //if(particle->Charge!=1) return; // only particles with Q=+1 can hope to reach RP200/FP420
- 
+
+
+  if(particle->Charge==0) return; // only particles with Q=+1 can hope to reach RP200/FP420
+  //cout << "particle ("<< particle->PID << "): m = " << mass << " \t Q= " << charge << endl;
   TRootRomanPotHits* elementRP220;
   //TRootForwardTaggerHits* elementFP420;
@@ -197 +233 @@
   
   // to go faster, why not rejecting particles already going into the ZDC?
-  if( particle->PID == pP) 
-  if( (particle->Status == 1) &&  (fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd) )
+  
+  // K_L^0 has a ctau of ~16m ; only pi+ and p+ can beat it ; 
+  // so if RP/FP too far away, the particle must be a proton or a mu+
+  if( std::min(DET->RP_420_s,DET->RP_220_s) > 17 && (particle->PID != pP && particle->PID != 13)) return;
+
+  if( fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd )
     {
-      //cout << "VeryForward :: M = " << mass << "\t Q = " << charge << "\t\t " << particle->PID << endl;
       H_BeamParticle p1(mass,charge); 
-      p1.smearAng();   p1.smearPos(); // vertex smearing
-        cout << "x = " << p1.getX() + DET->RP_cross_x 
-             << " y= " << p1.getY() + DET->RP_cross_y 
-             << " tx= " << p1.getTX() - kickers_on*DET->RP_cross_ang 
-             << " ty=" << p1.getTY() << endl;
-      p1.setPosition(p1.getX()+DET->RP_cross_x,p1.getY()+DET->RP_cross_y,p1.getTX()-kickers_on*DET->RP_cross_ang,p1.getTY(),0);
-      //p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
+      double tx = 1E6*atan(particle->Px/particle->Pz); // in microrad
+      double ty = 1E6*atan(particle->Py/particle->Pz); // in microrad
+      //p1.smearAng();   p1.smearPos(); // vertex smearing do not put it here !!!
+
+      /* cout << "x = "  << particle->X << " + " << p1.getX()   << " + " <<  DET->RP_cross_x 
+             << " y= "  << particle->Y << " + " << p1.getY()   << " + " << DET->RP_cross_y 
+             << " tx= " << tx << " + " << p1.getTX() << " - " << kickers_on*DET->RP_cross_ang_x 
+             << " ty= " << ty << " + " << p1.getTY()  << " - " << kickers_on*DET->RP_cross_ang_y  
+             << " z= "  << particle->Z << endl;*/
+
+      // here below, p1.getX(), p1.getY(), p1.getTX() and p1.getTY() =0 unless some smearing is done
+      // all in micrometers or microradians
+      p1.setPosition((1E3)*particle->X + p1.getX() + DET->RP_cross_x, 
+                     (1E3)*particle->Y + p1.getY() + DET->RP_cross_y, 
+                     tx + p1.getTX()- kickers_on*DET->RP_cross_ang_x,
+                     ty + p1.getTY()- kickers_on*DET->RP_cross_ang_y,
+                     (1E3)*particle->Z);
       p1.setE(particle->E);
        
@@ -218 +267 @@
       
       if(p1.stopped(beamline)) {
+        // roman pots at 220 m
         if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
-          static unsigned int counter;
-          counter++;
-          if (counter==1) {
-                //p1.getPath(0,"p1path.txt");
-                cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge 
-                     << "\t" << mass << "=" << particle->M 
-                     << "\t E=" << particle->E 
-                     << endl;           
-          }
+
+        /*static unsigned int counter;
+        counter++;
+        if (counter==1) {
+           p1.getPath(0,"p1path.txt");
+           cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge 
+                << "\t" << mass << "=" << particle->M << "\t E=" << particle->E << endl;
+          }*/
+
           p1.propagate(DET->RP_220_s);
           elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
+         
+          // detector measurements
           elementRP220->X  = (1E-6)*p1.getX();  // [m]
           elementRP220->Y  = (1E-6)*p1.getY();  // [m]
@@ -235 +287 @@
           elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
           elementRP220->S = p1.getS();          // [m]
-
-            /* time of flight t is t = T + d/[ cos(theta) v ]
-            // nb: here we assume a straight path to the detector, which is not the case!
-            // this time estimate is always underestimated (while exact for the ZDC case)
-            float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
-            if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
-               double tx = atan(particle->Px/particle->Pz);
-               double ty = atan(particle->Py/particle->Pz);
-               double theta = sqrt( pow(tx,2) + pow(ty,2) );
-               //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
-               // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
-               // eta = 2.6 <-> cos(theta) = 0.99
-               // eta = 3.0 <-> cos(theta) = 0.995
-               cos_theta = cos(theta);
-            }
-            // units from StdHEP : Z [mm] T[mm/c]
-            // units from Delphes : p1.getS [m]   speed_of_light [m/s]
-            //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
-            //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
-            */
-            elementRP220->E = p1.getE();        // not yet implemented
-            elementRP220->q2 = -1;                // not yet implemented
-            elementRP220->side = sign(particle->Eta);
-
-
-                elementRP220->pid = particle->PID;
+          elementRP220->T = time_of_flight(particle, DET->RP_220_s, DET->CEN_max_calo_fwd, DET->RP220_T_resolution);
+          elementRP220->side = sign(particle->Eta);
           
-        } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
+          // reconstructed data  
+          float sE = p1.getE(); // apply the smearing here!!!
+          elementRP220->E = sE;           // not yet implemented
+          elementRP220->q2 = UNDEFINED;   // not yet implemented
+
+          // generator level data
+          elementRP220->pid = particle->PID;
+          elementRP220->Set(genMomentum);
+        } // if RP220
+
+        // proton taggers at 420 m
+        else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
+
           p1.propagate(DET->RP_420_s);
+          elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
           //elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
-          elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
+
+          // detector measurements
           elementFP420->X  = (1E-6)*p1.getX();  // [m]
           elementFP420->Y  = (1E-6)*p1.getY();  // [m]
@@ -271 +313 @@
           elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
           elementFP420->S = p1.getS();          // [m]
-
-            // time of flight t is t = T + d/[ cos(theta) v ]
-            // nb: here we assume a straight path to the detector, which is not the case!
-            // this time estimate is always underestimated (while exact for the ZDC case)
-            float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
-            if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
-               double tx = atan(particle->Px/particle->Pz);
-               double ty = atan(particle->Py/particle->Pz);
-               double theta = sqrt( pow(tx,2) + pow(ty,2) );
-               //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
-               // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
-               // eta = 2.6 <-> cos(theta) = 0.99
-               // eta = 3.0 <-> cos(theta) = 0.995
-               cos_theta = cos(theta);
-            }
-            // units from StdHEP : Z [mm] T[mm/c]
-            // units from Delphes : p1.getS [m]   speed_of_light [m/s]
-            double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
-            elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
+          elementFP420->T = time_of_flight(particle, DET->RP_420_s, DET->CEN_max_calo_fwd, DET->RP420_T_resolution);
+          elementFP420->side = sign(particle->Eta);
+
+          // reconstructed data
           elementFP420->E = p1.getE();          // not yet implemented
-          elementFP420->q2 = -1;                // not yet implemented
-          elementFP420->side = sign(particle->Eta);
-
-
-                elementFP420->pid = particle->PID;
-        }
-
-      }
-      //  if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
-      //     cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
+          elementFP420->q2 = UNDEFINED;         // not yet implemented
+
+          // generator level data
+          elementFP420->pid = particle->PID;
+          elementFP420->Set(genMomentum);
+        }  // if FP420
+
+      } // if stopped
     } // if forward proton
-}
-
-    // Forward particles in CASTOR ?
-    //    if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
-    //                              && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
-    //
-    //
-    //    } // CASTOR
-    // */
-
+
+}