Fork me on GitHub

Changeset 374 in svn for trunk/src


Ignore:
Timestamp:
May 10, 2009, 8:23:32 PM (15 years ago)
Author:
Xavier Rouby
Message:

added Resolution terms for energy and timing for ZDC

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SmearUtil.cc

    r350 r374  
    6262  ELG_Nfwd          = 0.0;               // N term for FCAL
    6363  ELG_Cfwd          = 0.107;             // C term for FCAL
     64  ELG_Szdc          = 0.70;              // S term for ZDC
     65  ELG_Nzdc          = 0.0;               // N term for ZDC
     66  ELG_Czdc          = 0.08;              // C term for ZDC
    6467
    6568  // Energy resolution for hadrons in ecal/hcal/hf
    … …  
    7174  HAD_Nhf           = 0.;                // N term for FCAL
    7275  HAD_Chf           = 0.13;              // C term for FCAL
     76  HAD_Szdc          = 1.38;              // S term for ZDC
     77  HAD_Nzdc          = 0.;                // N term for ZDC
     78  HAD_Czdc          = 0.13;              // C term for ZDC
    7379
    7480  // Muon smearing
    7581  MU_SmearPt        = 0.01;
     82
     83  // time resolution
     84  ZDC_T_resolution   = 0;       // resolution for time measurement [s]
     85  RP220_T_resolution = 0;
     86  RP420_T_resolution = 0;
    7687
    7788  // Tracking efficiencies
    … …  
    95106
    96107
    97   // Thresholds for reconstructed objetcs
     108  // Thresholds for reconstructed objetcs (GeV)
    98109  PTCUT_elec      = 10.0;
    99110  PTCUT_muon      = 10.0;
    … …  
    101112  PTCUT_gamma     = 10.0;
    102113  PTCUT_taujet    = 10.0;
     114
     115  ZDC_gamma_E     = 20; // GeV
     116  ZDC_n_E         = 50; // GeV
    103117
    104118  // Isolation
    … …  
    199213  ELG_Sfwd          = DET.ELG_Sfwd;
    200214  ELG_Nfwd          = DET.ELG_Nfwd;
    201 
    202   // Energy resolution for hadrons in ecal/hcal/hf
     215  ELG_Czdc          = DET.ELG_Czdc;
     216  ELG_Szdc          = DET.ELG_Szdc;
     217  ELG_Nzdc          = DET.ELG_Nzdc;
     218
     219  // Energy resolution for hadrons in ecal/hcal/hf/zdc
    203220  HAD_Shcal         = DET.HAD_Shcal;
    204221  HAD_Nhcal         = DET.HAD_Nhcal;
    … …  
    207224  HAD_Nhf           = DET.HAD_Nhf;
    208225  HAD_Chf           = DET.HAD_Chf;
     226  HAD_Szdc          = DET.HAD_Szdc;
     227  HAD_Nzdc          = DET.HAD_Nzdc;
     228  HAD_Czdc          = DET.HAD_Czdc;
     229
     230  // time resolution
     231  ZDC_T_resolution   = DET.ZDC_T_resolution;       // resolution for time measurement [s]
     232  RP220_T_resolution = DET.RP220_T_resolution;
     233  RP420_T_resolution = DET.RP420_T_resolution;
    209234
    210235  // Muon smearing
    … …  
    229254  PTCUT_gamma     = DET.PTCUT_gamma;
    230255  PTCUT_taujet    = DET.PTCUT_taujet;
     256
     257  ZDC_gamma_E     = DET.ZDC_gamma_E;
     258  ZDC_n_E         = DET.ZDC_n_E;
    231259
    232260  // Isolation
    … …  
    321349  ELG_Sfwd          = DET.ELG_Sfwd;
    322350  ELG_Nfwd          = DET.ELG_Nfwd;
     351  ELG_Czdc          = DET.ELG_Czdc;
     352  ELG_Szdc          = DET.ELG_Szdc;
     353  ELG_Nzdc          = DET.ELG_Nzdc;
    323354
    324355  // Energy resolution for hadrons in ecal/hcal/hf
    … …  
    329360  HAD_Nhf           = DET.HAD_Nhf;
    330361  HAD_Chf           = DET.HAD_Chf;
     362  HAD_Szdc          = DET.HAD_Szdc;
     363  HAD_Nzdc          = DET.HAD_Nzdc;
     364  HAD_Czdc          = DET.HAD_Czdc;
     365
     366  // time resolution
     367  ZDC_T_resolution   = DET.ZDC_T_resolution;       // resolution for time measurement [s]
     368  RP220_T_resolution = DET.RP220_T_resolution;
     369  RP420_T_resolution = DET.RP420_T_resolution;
    331370
    332371  // Muon smearing
    … …  
    351390  PTCUT_gamma     = DET.PTCUT_gamma;
    352391  PTCUT_taujet    = DET.PTCUT_taujet;
     392
     393  ZDC_gamma_E     = DET.ZDC_gamma_E;
     394  ZDC_n_E         = DET.ZDC_n_E;
    353395
    354396  // Isolation
    … …  
    475517    else if(strstr(temp_string.c_str(),"ELG_Cfwd"))         {curstring >> varname >> value; ELG_Cfwd          = value;}
    476518    else if(strstr(temp_string.c_str(),"ELG_Nfwd"))         {curstring >> varname >> value; ELG_Nfwd          = value;}
     519    else if(strstr(temp_string.c_str(),"ELG_Szdc"))         {curstring >> varname >> value; ELG_Szdc          = value;}
     520    else if(strstr(temp_string.c_str(),"ELG_Czdc"))         {curstring >> varname >> value; ELG_Czdc          = value;}
     521    else if(strstr(temp_string.c_str(),"ELG_Nzdc"))         {curstring >> varname >> value; ELG_Nzdc          = value;}
     522
    477523    else if(strstr(temp_string.c_str(),"HAD_Shcal"))        {curstring >> varname >> value; HAD_Shcal         = value;}
    478524    else if(strstr(temp_string.c_str(),"HAD_Nhcal"))        {curstring >> varname >> value; HAD_Nhcal         = value;}
    … …  
    481527    else if(strstr(temp_string.c_str(),"HAD_Nhf"))          {curstring >> varname >> value; HAD_Nhf           = value;}
    482528    else if(strstr(temp_string.c_str(),"HAD_Chf"))          {curstring >> varname >> value; HAD_Chf           = value;}
     529    else if(strstr(temp_string.c_str(),"HAD_Szdc"))         {curstring >> varname >> value; HAD_Szdc          = value;}
     530    else if(strstr(temp_string.c_str(),"HAD_Nzdc"))         {curstring >> varname >> value; HAD_Nzdc          = value;}
     531    else if(strstr(temp_string.c_str(),"HAD_Czdc"))         {curstring >> varname >> value; HAD_Czdc          = value;}
     532    else if(strstr(temp_string.c_str(),"ZDC_T_resolution")) {curstring >> varname >> value; ZDC_T_resolution  = value;}
     533    else if(strstr(temp_string.c_str(),"RP220_T_resolution")) {curstring >> varname >> value; RP220_T_resolution  = value;}
     534    else if(strstr(temp_string.c_str(),"RP420_T_resolution")) {curstring >> varname >> value; RP420_T_resolution  = value;}
    483535    else if(strstr(temp_string.c_str(),"MU_SmearPt"))       {curstring >> varname >> value; MU_SmearPt        = value;}
    484536   
    … …  
    503555    else if(strstr(temp_string.c_str(),"PTCUT_gamma"))      {curstring >> varname >> value; PTCUT_gamma       = value;}
    504556    else if(strstr(temp_string.c_str(),"PTCUT_taujet"))     {curstring >> varname >> value; PTCUT_taujet      = value;}
     557    else if(strstr(temp_string.c_str(),"ZDC_gamma_E"))      {curstring >> varname >> value; ZDC_gamma_E      = value;}
     558    else if(strstr(temp_string.c_str(),"ZDC_n_E"))          {curstring >> varname >> value; ZDC_n_E      = value;}
    505559
    506560    else if(strstr(temp_string.c_str(),"ISOL_PT"))          {curstring >> varname >> value; ISOL_PT           = value;}
    … …  
    694748  f_out << left << setw(30) <<"* C term for FCAL: "<<""
    695749        << left << setw(30) <<ELG_Cfwd       <<""<< right << setw(10)<<"*"<<"\n";
     750  f_out << left << setw(30) <<"* S term for ZDC: "<<""
     751        << left << setw(30) <<ELG_Szdc       <<""<< right << setw(10)<<"*"<<"\n";
     752  f_out << left << setw(30) <<"* N term for ZDC: "<<""
     753        << left << setw(30) <<ELG_Nzdc       <<""<< right << setw(10)<<"*"<<"\n";
     754  f_out << left << setw(30) <<"* C term for ZDC: "<<""
     755        << left << setw(30) <<ELG_Czdc       <<""<< right << setw(10)<<"*"<<"\n";
     756
    696757  f_out<<"*                                                                    *"<<"\n";
    697758  f_out<<"#*****************************                                       *"<<"\n";
    … …  
    711772  f_out << left << setw(30) <<"* C term for FCAL: "<<""
    712773        << left << setw(30) <<HAD_Chf         <<""<< right << setw(10)<<"*"<<"\n";
     774  f_out << left << setw(30) <<"* S term for ZDC: "<<""
     775        << left << setw(30) <<HAD_Szdc        <<""<< right << setw(10)<<"*"<<"\n";
     776  f_out << left << setw(30) <<"* N term for ZDC: "<<""
     777        << left << setw(30) <<HAD_Nzdc        <<""<< right << setw(10)<<"*"<<"\n";
     778  f_out << left << setw(30) <<"* C term for ZDC: "<<""
     779        << left << setw(30) <<HAD_Czdc        <<""<< right << setw(10)<<"*"<<"\n";
     780
     781  f_out<<"*                                                                    *"<<"\n";
     782  f_out<<"#*************************                                           *"<<"\n";
     783  f_out<<"# Time smearing parameters                                           *"<<"\n";
     784  f_out<<"#*************************                                           *"<<"\n";
     785  f_out<<"*                                                                    *"<<"\n";
     786  f_out << left << setw(55) <<"* Time resolution for ZDC        :              "<<""
     787        << left << setw(5) <<ZDC_T_resolution       <<""<< right << setw(10)<<"*"<<"\n";
     788  f_out << left << setw(55) <<"* Time resolution for RP220      :              "<<""
     789        << left << setw(5) <<RP220_T_resolution     <<""<< right << setw(10)<<"*"<<"\n";
     790  f_out << left << setw(55) <<"* Time resolution for RP420      :              "<<""
     791        << left << setw(5) <<RP420_T_resolution     <<""<< right << setw(10)<<"*"<<"\n";
     792  f_out<<"*                                                                    *"<<"\n";
     793
    713794  f_out<<"*                                                                    *"<<"\n";
    714795  f_out<<"#*************************                                           *"<<"\n";
    … …  
    786867  f_out << left << setw(40) <<"* Minimum pT for photons: "<<""
    787868        << left << setw(20) <<PTCUT_gamma        <<""<< right << setw(10)<<"*"<<"\n";
     869  f_out << left << setw(40) <<"* Minimum E for photons in ZDC: "<<""
     870        << left << setw(20) <<ZDC_gamma_E        <<""<< right << setw(10)<<"*"<<"\n";
     871  f_out << left << setw(40) <<"* Minimum E for neutrons in ZDC: "<<""
     872        << left << setw(20) <<ZDC_n_E        <<""<< right << setw(10)<<"*"<<"\n";
     873
    788874  f_out<<"*                                                                    *"<<"\n";
    789875  f_out<<"#*******************                                                 *"<<"\n";

  • trunk/src/VeryForward.cc

    r361 r374  
    44**         |  Delphes, a framework for the fast simulation  |         **
    55**         |       of a  generic collider experiment        |         **
    6 **          \----------------------------------------------/          **
     6**          \-------------  arXiv:0903.2225v1  ------------/          **
    77**                                                                    **
    88**                                                                    **
    … …  
    3939
    4040//------------------------------------------------------------------------------
    41 
    42 VeryForward::VeryForward() {
    43    DET = new RESOLution();
    44    beamline1 = new H_BeamLine(1,500.);
    45    beamline2 = new H_BeamLine(1,500.);
    46    init();
    47   //Initialisation of Hector
    48   relative_energy = true; // should always be true
    49   kickers_on = 1;         // should always be 1
    50 
    51 }
    52 
    53 VeryForward::VeryForward(const string& DetDatacard) {
    54    DET = new RESOLution();
     41VeryForward::VeryForward() :
     42   DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
     43   beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
     44   relative_energy(true), // should always be true
     45   kickers_on(1)         // should always be 1
     46   {
     47   init(); //Initialisation of Hector
     48}
     49
     50VeryForward::VeryForward(const string& DetDatacard) :
     51           DET(new RESOLution())
     52   {
    5553   DET->ReadDataCard(DetDatacard);
    56    beamline1 = new H_BeamLine(1,500.);
    57    beamline2 = new H_BeamLine(1,500.);
    58    init();
    59   //Initialisation of Hector
    60   relative_energy = true; // should always be true
    61   kickers_on = 1;         // should always be 1
    62 
    63 }
    64 
    65 VeryForward::VeryForward(const RESOLution * DetDatacard) {
    66    DET = new RESOLution(*DetDatacard);
    67    beamline2 = new H_BeamLine(1,500.);
    68    beamline1 = new H_BeamLine(1,500.);
    69 
    70    init();
    71   //Initialisation of Hector
    72   relative_energy = true; // should always be true
    73   kickers_on = 1;         // should always be 1
    74 
    75 }
    76 
    77 VeryForward::VeryForward(const VeryForward& vf) {
    78    DET = new RESOLution(*(vf.DET));
    79    beamline1 = new H_BeamLine(*(vf.beamline1));
    80    beamline2 = new H_BeamLine(*(vf.beamline2));
     54   const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s);
     55   beamline1 = new H_BeamLine(1,d_max);
     56   beamline2 = new H_BeamLine(1,d_max);
     57   init(); //Initialisation of Hector
     58   relative_energy = true; // should always be true
     59   kickers_on = 1;         // should always be 1
     60}
     61
     62VeryForward::VeryForward(const RESOLution * DetDatacard) :
     63        DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
     64        beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
     65        relative_energy(true), // should always be true
     66        kickers_on(1)          // should always be 1
     67   {
     68   init();  //Initialisation of Hector
     69}
     70
     71VeryForward::VeryForward(const VeryForward& vf) :
     72   DET(new RESOLution(*(vf.DET))),   d_max(vf.d_max),
     73   beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))),
     74   relative_energy(vf.relative_energy),
     75   kickers_on(vf.kickers_on) {
    8176}
    8277
    … …  
    8479   if (this==&vf) return *this;
    8580   DET = new RESOLution(*(vf.DET));
     81   d_max = vf.d_max;
    8682   beamline1 = new H_BeamLine(*(vf.beamline1));
    8783   beamline2 = new H_BeamLine(*(vf.beamline2));
     84   relative_energy =vf.relative_energy;
     85   kickers_on = vf.kickers_on;
    8886   return *this;
    8987}
    … …  
    9492  relative_energy = true; // should always be true
    9593  kickers_on = 1;         // should always be 1
    96   // user should provide : (1) optics file for each beamline, and IPname,
    97   // and offset data (s,x) for optical elements
    9894  beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);                               
    9995  beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x);
    … …  
    114110
    115111
    116 void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC,TRootGenParticle *particle)
     112void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle)
    117113{
    118   int pid=abs(particle->PID);
    119114  TRootZdcHits *elementZdc;
    120   TLorentzVector genMomentum;
     115  float energy = particle->E;
     116  //TLorentzVector genMomentum;
     117
    121118  // Zero degree calorimeter, for forward neutrons and photons
    122   if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && fabs(particle->Eta) > DET->VFD_min_zdc ) {
    123     genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
    124     // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
    125     // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
    126     // !!!!!!!!! pas de smearing ! on garde trop d'info !
     119  if (particle->Status ==1 && ( (particle->PID==pN     && energy>DET->ZDC_n_E) ||
     120                                (particle->PID==pGAMMA && energy>DET->ZDC_gamma_E) )
     121                           && fabs(particle->Eta) > DET->VFD_min_zdc ) {
    127122    elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
    128     elementZdc->Set(genMomentum);
    129    
    130     // time of flight t is t = T + d/[ cos(theta) v ]
    131     //double tx = acos(particle->Px/particle->Pz);
    132     //double ty = acos(particle->Py/particle->Pz);
    133     //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
    134     //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
    135     double flight_distance = DET->VFD_s_zdc - particle->Z*(1E-6);
    136     // assumes also that the emission angle is so small that 1/(cos theta) = 1
    137     elementZdc->T = particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
    138 //cout << "ZDC: T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
     123 
     124    // 1) energy smearing
     125    float energyS = -1.;
     126    if (particle->PID == pGAMMA)
     127        energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->ELG_Nzdc,2) +
     128                                                   pow(DET->ELG_Czdc*particle->E,2) +
     129                                                   pow(DET->ELG_Szdc*sqrt(particle->E),2) ));
     130    else // smearing with hadronic resolution
     131        energyS =  gRandom->Gaus(particle->E, sqrt( pow(DET->HAD_Nzdc,2) +
     132                                                    pow(DET->HAD_Czdc*particle->E,2) +
     133                                                    pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
     134    elementZdc->E = energyS;
     135 
     136
     137    // 2) time of flight t is t = T + d/[ cos(theta) v ]
     138    float cos_theta = 1;      //very good approximation, if eta_zdc >3
     139    if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation
     140       double tx = atan(particle->Px/particle->Pz);
     141       double ty = atan(particle->Py/particle->Pz);
     142       double theta = sqrt( pow(tx,2) + pow(ty,2) );
     143       //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
     144       // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
     145       // eta = 2.6 <-> cos(theta) = 0.99
     146       // eta = 3.0 <-> cos(theta) = 0.995
     147       cos_theta = cos(theta);
     148    }
     149    // units from StdHEP : Z [mm] T[mm/c]
     150    // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
     151    double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ;
     152    double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
     153    double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution);
     154    elementZdc->T = timeS;
     155
     156    // 3) side: which ZDC has been hit?
    139157    elementZdc->side = sign(particle->Eta);
    140    
    141    
     158
     159    // 4) object nature : e.m. (photon) or had (neutron) ?
     160    elementZdc->hadronic_hit = (bool) (particle->PID==pN);
    142161  }
    143162 
    … …  
    148167 
    149168  TRootRomanPotHits* elementRP220;
    150   TRootRomanPotHits* elementFP420;
     169  TRootForwardTaggerHits* elementFP420;
    151170 
    152171  TLorentzVector genMomentum;
    … …  
    177196          elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
    178197          elementRP220->S = p1.getS();          // [m]
    179           // in first approximation only ! this number is always lower than the real distance-of-flight
    180           double flight_distance = p1.getS() - particle->Z*(1E-6);
    181           elementRP220->T = particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
    182 //cout << "T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
    183           elementRP220->E = p1.getE();          // not yet implemented
    184           elementRP220->q2 = -1;                // not yet implemented
    185           elementRP220->side = sign(particle->Eta);
     198
     199            /* time of flight t is t = T + d/[ cos(theta) v ]
     200            // nb: here we assume a straight path to the detector, which is not the case!
     201            // this time estimate is always underestimated (while exact for the ZDC case)
     202            float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
     203            if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
     204               double tx = atan(particle->Px/particle->Pz);
     205               double ty = atan(particle->Py/particle->Pz);
     206               double theta = sqrt( pow(tx,2) + pow(ty,2) );
     207               //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
     208               // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
     209               // eta = 2.6 <-> cos(theta) = 0.99
     210               // eta = 3.0 <-> cos(theta) = 0.995
     211               cos_theta = cos(theta);
     212            }
     213            // units from StdHEP : Z [mm] T[mm/c]
     214            // units from Delphes : p1.getS [m]   speed_of_light [m/s]
     215            //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
     216            //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
     217            */
     218            elementRP220->E = p1.getE();        // not yet implemented
     219            elementRP220->q2 = -1;                // not yet implemented
     220            elementRP220->side = sign(particle->Eta);
    186221         
    187222        } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
    188223          p1.propagate(DET->RP_420_s);
    189           elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
     224          elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
    190225          elementFP420->X  = (1E-6)*p1.getX();  // [m]
    191226          elementFP420->Y  = (1E-6)*p1.getY();  // [m]
    … …  
    193228          elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
    194229          elementFP420->S = p1.getS();          // [m]
    195           // in first approximation only ! this number is always lower than the real distance-of-flight
    196           double flight_distance = p1.getS() - particle->Z*(1E-6);
    197 //cout << "T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
    198           elementFP420->T = particle->T + flight_distance/speed_of_light;
     230
     231            // time of flight t is t = T + d/[ cos(theta) v ]
     232            // nb: here we assume a straight path to the detector, which is not the case!
     233            // this time estimate is always underestimated (while exact for the ZDC case)
     234            float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
     235            if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
     236               double tx = atan(particle->Px/particle->Pz);
     237               double ty = atan(particle->Py/particle->Pz);
     238               double theta = sqrt( pow(tx,2) + pow(ty,2) );
     239               //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
     240               // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
     241               // eta = 2.6 <-> cos(theta) = 0.99
     242               // eta = 3.0 <-> cos(theta) = 0.995
     243               cos_theta = cos(theta);
     244            }
     245            // units from StdHEP : Z [mm] T[mm/c]
     246            // units from Delphes : p1.getS [m]   speed_of_light [m/s]
     247            double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
     248            elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
    199249          elementFP420->E = p1.getE();          // not yet implemented
    200250          elementFP420->q2 = -1;                // not yet implemented
Note: See TracChangeset for help on using the changeset viewer.