Fork me on GitHub

Changes in / [46d3442:4acf2fd] in git


Ignore:
Files:
4 added
15 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r46d3442 r4acf2fd  
    412412        modules/ImpactParameterSmearing.h \
    413413        modules/TimeSmearing.h \
     414        modules/TimeOfFlight.h \
    414415        modules/SimpleCalorimeter.h \
    415416        modules/DenseTrackFilter.h \
    … …  
    448449        modules/DecayFilter.h \
    449450        modules/ParticleDensity.h \
     451        modules/TruthVertexFinder.h \
    450452        modules/ExampleModule.h
    451453tmp/modules/ModulesDict$(PcmSuf): \
    … …  
    10121014        classes/DelphesFactory.h \
    10131015        classes/DelphesFormula.h
     1016tmp/modules/TimeOfFlight.$(ObjSuf): \
     1017        modules/TimeOfFlight.$(SrcSuf) \
     1018        modules/TimeOfFlight.h \
     1019        classes/DelphesClasses.h \
     1020        classes/DelphesFactory.h \
     1021        classes/DelphesFormula.h \
     1022        external/ExRootAnalysis/ExRootClassifier.h \
     1023        external/ExRootAnalysis/ExRootFilter.h \
     1024        external/ExRootAnalysis/ExRootResult.h
    10141025tmp/modules/TimeSmearing.$(ObjSuf): \
    10151026        modules/TimeSmearing.$(SrcSuf) \
    … …  
    10711082        external/ExRootAnalysis/ExRootResult.h \
    10721083        external/ExRootAnalysis/ExRootTreeBranch.h
     1084tmp/modules/TruthVertexFinder.$(ObjSuf): \
     1085        modules/TruthVertexFinder.$(SrcSuf) \
     1086        modules/TruthVertexFinder.h \
     1087        classes/DelphesClasses.h \
     1088        classes/DelphesFactory.h \
     1089        classes/DelphesPileUpReader.h \
     1090        classes/DelphesTF2.h \
     1091        external/ExRootAnalysis/ExRootClassifier.h \
     1092        external/ExRootAnalysis/ExRootFilter.h \
     1093        external/ExRootAnalysis/ExRootResult.h
    10731094tmp/modules/UniqueObjectFinder.$(ObjSuf): \
    10741095        modules/UniqueObjectFinder.$(SrcSuf) \
    … …  
    12151236        tmp/modules/TaggingParticlesSkimmer.$(ObjSuf) \
    12161237        tmp/modules/TauTagging.$(ObjSuf) \
     1238        tmp/modules/TimeOfFlight.$(ObjSuf) \
    12171239        tmp/modules/TimeSmearing.$(ObjSuf) \
    12181240        tmp/modules/TrackCountingBTagging.$(ObjSuf) \
    … …  
    12221244        tmp/modules/TrackSmearing.$(ObjSuf) \
    12231245        tmp/modules/TreeWriter.$(ObjSuf) \
     1246        tmp/modules/TruthVertexFinder.$(ObjSuf) \
    12241247        tmp/modules/UniqueObjectFinder.$(ObjSuf) \
    12251248        tmp/modules/VertexFinder.$(ObjSuf) \
    … …  
    20672090        @touch $@
    20682091
     2092modules/TimeOfFlight.h: \
     2093        classes/DelphesModule.h
     2094        @touch $@
     2095
    20692096external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh: \
    20702097        external/fastjet/ClusterSequence.hh \
    … …  
    21852212        @touch $@
    21862213
     2214modules/TruthVertexFinder.h: \
     2215        classes/DelphesModule.h
     2216        @touch $@
     2217
    21872218classes/DelphesSTDHEPReader.h: \
    21882219        classes/DelphesXDRReader.h
    … …  
    23002331        @touch $@
    23012332
     2333modules/TreeWriter.h: \
     2334        classes/DelphesModule.h
     2335        @touch $@
     2336
    23022337modules/TimeSmearing.h: \
    2303         classes/DelphesModule.h
    2304         @touch $@
    2305 
    2306 modules/TreeWriter.h: \
    23072338        classes/DelphesModule.h
    23082339        @touch $@

  • cards/delphes_card_IDEA.tcl

    r46d3442 r4acf2fd  
    1111set B 2.0
    1212
     13## Drift chamber coordinates
     14set DCHZMIN -2.125
     15set DCHZMAX 2.125
     16set DCHRMIN 0.345
     17set DCHRMAX 2.02
     18
     19
    1320#######################################
    1421# Order of execution of various modules
    … …  
    1623
    1724set ExecutionPath {
     25
     26  TruthVertexFinder
    1827  ParticlePropagator
    1928
    … …  
    2534  TrackSmearing
    2635  ClusterCounting
     36  TimeSmearing
     37  TimeOfFlight
    2738
    2839  TrackMerger
    … …  
    6172  ScalarHT
    6273  TreeWriter
     74}
     75
     76#################################
     77# Truth Vertex Finder
     78#################################
     79
     80module TruthVertexFinder TruthVertexFinder {
     81
     82  ## below this distance two vertices are assumed to be the same
     83  set Resolution 1E-06
     84
     85  set InputArray Delphes/stableParticles
     86  set VertexOutputArray vertices
    6387}
    6488
    … …  
    101125    }
    102126}
    103 
    104 #       (pt <= 0.1)                                     * (0.00) +
    105 #       (abs(eta) <= 3.0)               * (pt > 0.1)    * (1.00) +
    106 #       (abs(eta) > 3)                                  * (0.00)
    107127
    108128
    … …  
    212232      2        VTXDSK      0.141  0.3    0.92     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    213233
    214       1 DCHCANI -2.125 2.125 0.345 0.0002 0.237223 0 0 0 0 0 0
     234      1 DCHCANI $DCHZMIN $DCHZMAX $DCHRMIN 0.0002 0.237223 0 0 0 0 0 0
    215235      1 DCH -2 2 0.36 0.0147748 1400 1 0.0203738 0 0.0001 0 1
    216236      1 DCH -2 2 0.374775 0.0147748 1400 1 -0.0212097 0 0.0001 0 1
    … …  
    325345      1 DCH -2 2 1.98523 0.0147748 1400 1 0.111898 0 0.0001 0 1
    326346      1 DCH -2 2 2 0.0147748 1400 1 -0.112723 0 0.0001 0 1
    327       1 DCHCANO -2.125 2.125 2.02 0.02 1.667 0 0 0 0 0 0
     347      1 DCHCANO $DCHZMIN $DCHZMAX $DCHRMAX $DCHRMAX 0.02 1.667 0 0 0 0 0 0
    328348      1 BSILWRP -2.35 2.35 2.04 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    329349      1 BSILWRP -2.35 2.35 2.06 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    330350      1 MAG -2.5 2.5 2.25 0.05 0.0658 0 0 0 0 0 0
    331351      1 BPRESH -2.55 2.55 2.45 0.02 1 2 0 1.5708 7e-005 0.01 1
    332       2 DCHWALL 0.345 2.02 2.125 0.25 5.55 0 0 0 0 0 0
    333       2 DCHWALL 0.345 2.02 -2.125 0.25 5.55 0 0 0 0 0 0
     352      2 DCHWALL $DCHRMIN $DCHRMAX $DCHZMAX 0.25 5.55 0 0 0 0 0 0
     353      2 DCHWALL $DCHRMIN $DCHRMAX $DCHZMIN 0.25 5.55 0 0 0 0 0 0
    334354      2 FSILWRP 0.354 2.02 -2.32 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    335355      2 FSILWRP 0.35 2.02 -2.3 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    … …  
    355375  set Bz $B
    356376
    357   set Rmin 0.35
    358   set Rmax 2.0
    359   set Zmin -2.0
    360   set Zmax 2.0
    361 
    362   # gas mix option: 0
     377  ## check that these are consistent with DCHCANI/DCHNANO parameters in TrackCovariance module
     378  set Rmin $DCHRMIN
     379  set Rmax $DCHRMAX
     380  set Zmin $DCHZMIN
     381  set Zmax $DCHZMAX
     382
     383  # gas mix option:
    363384  # 0:  Helium 90% - Isobutane 10%
    364385  # 1:  Helium 100%
    … …  
    371392
    372393
     394########################################
     395#   Time Smearing MIP
     396########################################
     397
     398module TimeSmearing TimeSmearing {
     399  set TrackInputArray ClusterCounting/tracks
     400  set OutputArray tracks
     401
     402  # assume constant 30 ps resolution for now
     403  set TimeResolution {
     404                       (abs(eta) > 0.0 && abs(eta) <= 3.0)* 30E-12
     405                     }
     406}
     407
     408########################################
     409#   Time Of Flight Measurement
     410########################################
     411
     412module TimeOfFlight TimeOfFlight {
     413  set TrackInputArray TimeSmearing/tracks
     414  set VertexInputArray TruthVertexFinder/vertices
     415
     416  set OutputArray tracks
     417
     418  # 0: assume vertex time tV from MC Truth (ideal case)
     419  # 1: assume vertex time tV=0
     420  # 2: calculate vertex time as vertex TOF, assuming tPV=0
     421  set VertexTimeMode 2
     422
     423}
     424
    373425
    374426##############
    … …  
    378430module Merger TrackMerger {
    379431# add InputArray InputArray
    380   add InputArray ClusterCounting/tracks
     432  add InputArray TimeOfFlight/tracks
    381433  set OutputArray tracks
    382434}
    … …  
    779831
    780832module UniqueObjectFinder UniqueObjectFinder {
     833
     834  set UseUniqueID true
    781835# earlier arrays take precedence over later ones
    782836# add InputArray InputArray OutputArray
    … …  
    788842
    789843
    790 
    791844##################
    792845# ROOT tree writer

  • classes/DelphesClasses.cc

    r46d3442 r4acf2fd  
    231231  Xd(0), Yd(0), Zd(0),
    232232  Nclusters(0.0),
     233  dNdx(0.0),
    233234  TrackResolution(0),
    234235  NCharged(0),
    … …  
    408409  object.Zd = Zd;
    409410  object.Nclusters = Nclusters;
     411  object.dNdx = dNdx;
    410412  object.TrackResolution = TrackResolution;
    411413  object.NCharged = NCharged;
    … …  
    541543  Zd = 0.0;
    542544  Nclusters = 0.0;
     545  dNdx = 0.0;
    543546  TrackResolution = 0.0;
    544547  NCharged = 0;

  • classes/DelphesClasses.h

    r46d3442 r4acf2fd  
    462462  Float_t DZ; // track longitudinal impact parameter
    463463  Float_t Nclusters; // Number of ionization clusters
     464  Float_t dNdx; // Number of ionization clusters
    464465
    465466  Float_t ErrorP; // track momentum error
    … …  
    567568  Float_t DZ; // track longitudinal impact parameter
    568569  Float_t Nclusters; // Number of ionization clusters
     570  Float_t dNdx; // Number of ionization clusters
    569571
    570572  Float_t ErrorP; // track momentum error
    … …  
    697699
    698700  Float_t Nclusters; // Number of ionization clusters
     701  Float_t dNdx; // Number of ionization clusters per unit length
    699702
    700703  Float_t Xd;

  • examples/Example5.C

    r46d3442 r4acf2fd  
    3737    // Load selected branches with data from specified event
    3838    treeReader->ReadEntry(entry);
    39    
     39
    4040    if(entry>10) break;
    41    
     41
    4242    cout<<"" <<endl;
    4343    cout<<"--------- New Event ---------" <<endl;
    4444    cout<<"" <<endl;
    45  
     45
    4646    // loop over all input particles in the event
    4747    for(Int_t i=0; i < branchParticle->GetEntriesFast(); i++)
    48     {   
    49      GenParticle *gen = (GenParticle*) branchParticle->At(i);     
    50      cout<<"N: "<<i<<", St: "<<gen->Status<<", PID: "<<gen->PID<<", E: "<<gen->E<<", Px: "<<gen->Px<<", Py: "<<gen->Py<<", Pz: "<<gen->Pz<<", M: "<<gen->Mass<<", M1: "<<gen->M1<<", M2: "<<gen->M2<<", D1: "<<gen->D1<<", D2: "<<gen->D2<<endl;
     48    {
     49     GenParticle *gen = (GenParticle*) branchParticle->At(i);
     50     cout<<"N: "<<i<<", St: "<<gen->Status<<", PID: "<<gen->PID<<", E: "<<gen->E<<", Px: "<<gen->Px<<", Py: "<<gen->Py<<", Pz: "<<gen->Pz<<", M: "<<gen->Mass<<", T: "<<gen->T<<", X: "<<gen->X<<", Y: "<<gen->Y<<", Z: "<<gen->Z<<", M1: "<<gen->M1<<", M2: "<<gen->M2<<", D1: "<<gen->D1<<", D2: "<<gen->D2<<endl;
    5151    }
    5252  }

  • examples/Pythia8/ee_zh.cmnd

    r46d3442 r4acf2fd  
    11
    2 Main:numberOfEvents = 10         ! number of events to generate
     2! number of events to generate
     3Main:numberOfEvents = 1000         ! number of events to generate
    34
    45Beams:idA = 11                   ! first beam, e- = -11
    … …  
    78
    89! Vertex smearing :
    9 Beams:allowVertexSpread = off
    10 !Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
    11 !Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
    12 !Beams:sigmaVertexZ = 20.64      !  0.64 mm
    13 !ï¿ŒBeams:sigmaVertexT = 0.64      !  0.64 mm
    14 ï¿Œ
    1510
     11Beams:allowVertexSpread = on
     12Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
     13Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
     14Beams:sigmaVertexZ = 0.64      !  0.64 mm
     15Beams:sigmaTime    = 0.64      !  0.64 mm
    1616
    1717! Higgsstrahlung process

  • modules/ClusterCounting.cc

    r46d3442 r4acf2fd  
    1   /*
     1   /*
    22 *  Delphes: a framework for fast simulation of a generic collider experiment
    33 *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), Belgium
    … …  
    2828 */
    2929
    30 //FIXME add reference to Bedeschi-code
    31 //FIXME make sure about units of P, X
    32 //FIXME fix pt > 200 GeV issue and angle > 6.41
    33 
    3430#include "modules/ClusterCounting.h"
    35 
    3631#include "classes/DelphesClasses.h"
    3732#include "TrackCovariance/TrkUtil.h"
    … …  
    108103{
    109104  Candidate *candidate, *mother, *particle;
    110   Double_t mass, Ncl;
     105  Double_t mass, trackLength, Ncl;
    111106
    112107  fItInputArray->Reset();
    … …  
    124119    mass = candidateMomentum.M();
    125120
     121    trackLength = fTrackUtil->TrkLen(Par);
     122
    126123    mother    = candidate;
    127124    candidate = static_cast<Candidate*>(candidate->Clone());
    128125
    129     Ncl = -999;
    130     // computation of Nclusters is not supported for electrons
    131     if (TMath::Abs(particle->PID) == 11)
     126    Ncl = 0.;
     127    if (fTrackUtil->IonClusters(Ncl, mass, Par))
    132128    {
    133129      candidate->Nclusters = Ncl;
     130      candidate->dNdx = (trackLength > 0.) ? Ncl/trackLength : -1;
    134131    }
    135     else if (fTrackUtil->IonClusters(Ncl, mass, Par))
    136     {
    137       candidate->Nclusters = Ncl;
    138     }
    139     //cout<<candidate->PID<<", "<<mass<<", "<<candidate->Nclusters<<endl;
    140132
    141133    candidate->AddCandidate(mother);

  • modules/DualReadoutCalorimeter.cc

    r46d3442 r4acf2fd  
    613613    // create new photon tower
    614614    tower = static_cast<Candidate*>(fTower->Clone());
    615     pt =  neutralEnergy / TMath::CosH(eta);
     615    pt = neutralEnergy / TMath::CosH(eta);
    616616    //cout<<"Creating tower with Pt, Eta, Phi, Energy: "<<pt<<","<<eta<<","<<phi<<","<<neutralEnergy<<endl;
    617617    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
    … …  
    623623      tower->Ehad = 0.0;
    624624      tower->PID = 22;
    625     }
    626 
     625      fEFlowPhotonOutputArray->Add(tower);
     626    }
    627627    // if hadronic fraction > 0, use HCAL resolution
    628628    else
    … …  
    631631      tower->Ehad = neutralEnergy;
    632632      tower->PID = 130;
    633     }
    634 
    635     fEFlowPhotonOutputArray->Add(tower);
    636 
     633      fEFlowNeutralHadronOutputArray->Add(tower);
     634    }
    637635
    638636    //clone tracks

  • modules/Efficiency.cc

    r46d3442 r4acf2fd  
    100100  Double_t pt, eta, phi, e;
    101101
     102
    102103  fItInputArray->Reset();
    103104  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    … …  
    115116    pt = candidateMomentum.Pt();
    116117    e = candidateMomentum.E();
    117    
     118
    118119    // apply an efficency formula
    119120    if(gRandom->Uniform() > fFormula->Eval(pt, eta, phi, e, candidate)) continue;

  • modules/ModulesLinkDef.h

    r46d3442 r4acf2fd  
    4040#include "modules/ImpactParameterSmearing.h"
    4141#include "modules/TimeSmearing.h"
     42#include "modules/TimeOfFlight.h"
    4243#include "modules/SimpleCalorimeter.h"
    4344#include "modules/DenseTrackFilter.h"
    … …  
    7677#include "modules/DecayFilter.h"
    7778#include "modules/ParticleDensity.h"
     79#include "modules/TruthVertexFinder.h"
    7880#include "modules/ExampleModule.h"
    7981
    … …  
    98100#pragma link C++ class ImpactParameterSmearing+;
    99101#pragma link C++ class TimeSmearing+;
     102#pragma link C++ class TimeOfFlight+;
    100103#pragma link C++ class SimpleCalorimeter+;
    101104#pragma link C++ class DenseTrackFilter+;
    … …  
    134137#pragma link C++ class DecayFilter+;
    135138#pragma link C++ class ParticleDensity+;
     139#pragma link C++ class TruthVertexFinder+;
    136140#pragma link C++ class ExampleModule+;
    137141

  • modules/ParticlePropagator.cc

    r46d3442 r4acf2fd  
    125125  TLorentzVector particlePosition, particleMomentum, beamSpotPosition;
    126126  Double_t px, py, pz, pt, pt2, e, q;
    127   Double_t x, y, z, t, r, phi;
     127  Double_t x, y, z, t, r;
    128128  Double_t x_c, y_c, r_c, phi_c, phi_0;
    129129  Double_t x_t, y_t, z_t, r_t;
    130   Double_t t1, t2, t3, t4, t5, t6;
    131   Double_t t_z, t_r, t_ra, t_rb;
    132   Double_t tmp, discr, discr2;
    133   Double_t delta, gammam, omega, asinrho;
    134   Double_t rcu, rc2, xd, yd, zd;
    135   Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
     130  Double_t t_z, t_r;
     131  Double_t discr;
     132  Double_t gammam, omega;
     133  Double_t xd, yd, zd;
     134  Double_t l, d0, dz, ctgTheta, alpha;
    136135  Double_t bsx, bsy, bsz;
    137   Double_t s0, s1, sd;
    138  
     136  Double_t rxp, rdp, t_R;
     137  Double_t td, pio, phid, sign_pz, vz;
     138
    139139  const Double_t c_light = 2.99792458E8;
    140140
    … …  
    161161    particlePosition = particle->Position;
    162162    particleMomentum = particle->Momentum;
     163
     164    // Constants
     165
    163166    x = particlePosition.X() * 1.0E-3;
    164167    y = particlePosition.Y() * 1.0E-3;
    … …  
    205208    else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
    206209    {
    207       // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
    208       tmp = px * y - py * x;
    209       discr2 = pt2 * fRadius2 - tmp * tmp;
    210 
    211       if(discr2 < 0.0)
    212       {
    213         // no solutions
    214         continue;
    215       }
    216 
    217       tmp = px * x + py * y;
    218       discr = TMath::Sqrt(discr2);
    219       t1 = (-tmp + discr) / pt2;
    220       t2 = (-tmp - discr) / pt2;
    221       t = (t1 < 0.0) ? t2 : t1;
    222 
    223       z_t = z + pz * t;
    224       if(TMath::Abs(z_t) > fHalfLength)
    225       {
    226         t3 = (+fHalfLength - z) / pz;
    227         t4 = (-fHalfLength - z) / pz;
    228         t = (t3 < 0.0) ? t4 : t3;
    229       }
    230 
    231       x_t = x + px * t;
    232       y_t = y + py * t;
    233       z_t = z + pz * t;
    234 
    235       l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z));
     210
     211      rxp = x*py - y*px;
     212      rdp = x*px + y*py;
     213
     214      discr = fRadius*fRadius*pt*pt - rxp*rxp;
     215
     216      t_R = e * (sqrt(discr) - rdp) / (c_light * pt * pt);
     217      t_z = e * (TMath::Sign(fHalfLengthMax, pz) - z) / ( c_light * pz);
     218
     219      t = TMath::Min(t_R, t_z);
     220
     221      x_t = x + px*t*c_light/e;
     222      y_t = y + py*t*c_light/e;
     223      z_t = z + pz*t*c_light/e;
     224      r_t = TMath::Hypot(x_t, y_t);
     225
     226      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
    236227
    237228      mother = candidate;
    238       candidate = static_cast<Candidate *>(candidate->Clone());
     229      candidate = static_cast<Candidate*>(candidate->Clone());
    239230
    240231      candidate->InitialPosition = particlePosition;
    241       candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3);
    242       candidate->L = l * 1.0E3;
     232      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*c_light*1.0E3);
     233      candidate->L = l*1.0E3;
    243234
    244235      candidate->Momentum = particleMomentum;
    … …  
    246237
    247238      fOutputArray->Add(candidate);
     239
    248240      if(TMath::Abs(q) > 1.0E-9)
    249241      {
    … …  
    274266      //     helix radius r = p_{T0} / (omega gamma m)
    275267
    276       gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
    277       omega = q * fBz / (gammam); // omega is here in [89875518/s]
    278       r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
     268      gammam = e*1.0E9 / (c_light*c_light);      // gammam in [eV/c^2]
     269      omega = q * fBz / (gammam);                // omega is here in [89875518/s]
     270      r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
    279271
    280272      phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
    281273
    282274      // 2. helix axis coordinates
    283       x_c = x + r * TMath::Sin(phi_0);
    284       y_c = y - r * TMath::Cos(phi_0);
     275      x_c = x + r*TMath::Sin(phi_0);
     276      y_c = y - r*TMath::Cos(phi_0);
    285277      r_c = TMath::Hypot(x_c, y_c);
    286       phi_c = TMath::ATan2(y_c, x_c);
    287       phi = phi_c;
    288       if(x_c < 0.0) phi += TMath::Pi();
    289 
    290       rcu = TMath::Abs(r);
    291       rc2 = r_c * r_c;
    292 
    293       // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
    294       xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c;
    295       xd = (rc2 > 0.0) ? xd / rc2 : -999;
    296       yd = y_c * (-rcu * r_c + rc2);
    297       yd = (rc2 > 0.0) ? yd / rc2 : -999;
    298       zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt;
    299 
    300       // proper calculation of the DCAz coordinate
    301       // s0: track circle parameter at the track origin
    302       // s1: track circle parameter at the closest approach to beam pipe
    303       // sd: s1-s0 signed angular difference
    304       s0 = atan2(y - y_c, x - x_c);
    305       s1 = atan2(yd - y_c, xd - x_c);
    306       sd = atan2(sin(s1 - s0), cos(s1 - s0));
    307       zd = z - r * pz / pt * sd;
    308      
    309       // use perigee momentum rather than original particle
    310       // momentum, since the orignal particle momentum isn't known
    311 
    312       px = TMath::Sign(1.0, r) * pt * (-y_c / r_c);
    313       py = TMath::Sign(1.0, r) * pt * (x_c / r_c);
    314       etap = particleMomentum.Eta();
    315       phip = TMath::ATan2(py, px);
    316 
    317       particleMomentum.SetPtEtaPhiE(pt, etap, phip, particleMomentum.E());
     278      phi_c = TMath::ATan(y_c/x_c);
     279      if(x_c < 0.0) phi_c -= TMath::Sign(1., phi_c)*TMath::Pi();
     280
     281      //Find the time of closest approach
     282      td = (phi_0 - TMath::ATan(-x_c/y_c))/omega;
     283
     284      //Remove all the modulo pi that might have come from the atan
     285      pio = fabs(TMath::Pi()/omega);
     286      while(fabs(td) > 0.5*pio)
     287      {
     288        td -= TMath::Sign(1., td)*pio;
     289      }
     290
     291      //Compute the coordinate of closed approach to z axis
     292      //if wants wtr beamline need to be changedto re-center with a traslation of the z axis
     293      phid = phi_0 - omega*td;
     294      xd = x_c - r*TMath::Sin(phid);
     295      yd = y_c + r*TMath::Cos(phid);
     296      zd = z + c_light*(pz/e)*td;
     297
     298      //Compute momentum at closest approach (perigee??)
     299      px = pt*TMath::Cos(phid);
     300      py = pt*TMath::Sin(phid);
     301
     302      particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E());
    318303
    319304      // calculate additional track parameters (correct for beamspot position)
    320 
    321       d0 = ((x - bsx) * py - (y - bsy) * px) / pt;
    322       dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
    323       p = particleMomentum.P();
    324       ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
     305      d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
     306      dz = zd - bsz;
     307      ctgTheta  = 1.0 / TMath::Tan (particleMomentum.Theta());
    325308
    326309      // 3. time evaluation t = TMath::Min(t_r, t_z)
    327310      //    t_r : time to exit from the sides
    328311      //    t_z : time to exit from the front or the back
    329       t_r = 0.0; // in [ns]
    330       int sign_pz = (pz > 0.0) ? 1 : -1;
    331       if(pz == 0.0)
    332         t_z = 1.0E99;
     312      t = 0;
     313      t_z = 0;
     314      sign_pz = (pz > 0.0) ? 1 : -1;
     315      if(pz == 0.0) t_z = 1.0E99;
     316      else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
     317
     318      if(r_c + TMath::Abs(r)  < fRadius)   // helix does not cross the cylinder sides
     319      {
     320        t = t_z;
     321      }
    333322      else
    334         t_z = gammam / (pz * 1.0E9 / c_light) * (-z + fHalfLength * sign_pz);
    335 
    336       if(r_c + TMath::Abs(r) < fRadius)
    337       {
    338         // helix does not cross the cylinder sides
    339         t = t_z;
    340       }
    341       else
    342       {
    343         asinrho = TMath::ASin((fRadius * fRadius - r_c * r_c - r * r) / (2 * TMath::Abs(r) * r_c));
    344         delta = phi_0 - phi;
    345         if(delta < -TMath::Pi()) delta += 2 * TMath::Pi();
    346         if(delta > TMath::Pi()) delta -= 2 * TMath::Pi();
    347         t1 = (delta + asinrho) / omega;
    348         t2 = (delta + TMath::Pi() - asinrho) / omega;
    349         t3 = (delta + TMath::Pi() + asinrho) / omega;
    350         t4 = (delta - asinrho) / omega;
    351         t5 = (delta - TMath::Pi() - asinrho) / omega;
    352         t6 = (delta - TMath::Pi() + asinrho) / omega;
    353 
    354         if(t1 < 0.0) t1 = 1.0E99;
    355         if(t2 < 0.0) t2 = 1.0E99;
    356         if(t3 < 0.0) t3 = 1.0E99;
    357         if(t4 < 0.0) t4 = 1.0E99;
    358         if(t5 < 0.0) t5 = 1.0E99;
    359         if(t6 < 0.0) t6 = 1.0E99;
    360 
    361         t_ra = TMath::Min(t1, TMath::Min(t2, t3));
    362         t_rb = TMath::Min(t4, TMath::Min(t5, t6));
    363         t_r = TMath::Min(t_ra, t_rb);
     323      {
     324        alpha = -(fRadius*fRadius - r*r - r_c*r_c)/(2*fabs(r)*r_c);
     325        alpha = fabs(TMath::ACos(alpha));
     326        t_r = td + alpha/fabs(omega);
     327
    364328        t = TMath::Min(t_r, t_z);
    365329      }
    366330
    367       // 4. position in terms of x(t), y(t), z(t)
    368       x_t = x_c + r * TMath::Sin(omega * t - phi_0);
    369       y_t = y_c + r * TMath::Cos(omega * t - phi_0);
    370       z_t = z + pz * 1.0E9 / c_light / gammam * t;
    371       r_t = TMath::Hypot(x_t, y_t);
     331      x_t = x_c - r*TMath::Sin(phi_0 - omega*t);
     332      y_t = y_c + r*TMath::Cos(phi_0 - omega*t);
     333      z_t = z + c_light*t*pz/e;
     334      r_t =  TMath::Hypot(x_t, y_t);
    372335
    373336      // compute path length for an helix
    374 
    375       alpha = pz * 1.0E9 / c_light / gammam;
    376       l = t * TMath::Sqrt(alpha * alpha + r * r * omega * omega);
     337      vz = pz*1.0E9 / c_light / gammam;
     338      //lenght of the path from production to tracker
     339      l = t * TMath::Sqrt(vz*vz + r*r*omega*omega);
    377340
    378341      if(r_t > 0.0)
    379342      {
    380 
    381343        // store these variables before cloning
    382344        if(particle == candidate)
    … …  
    384346          particle->D0 = d0 * 1.0E3;
    385347          particle->DZ = dz * 1.0E3;
    386           particle->P = p;
     348          particle->P = particleMomentum.P();
    387349          particle->PT = pt;
    388350          particle->CtgTheta = ctgTheta;
    389           particle->Phi = phip;
     351          particle->Phi = particleMomentum.Phi();
    390352        }
    391353

  • modules/TimeSmearing.cc

    r46d3442 r4acf2fd  
    1919/** \class TimeSmearing
    2020 *
    21  *  Performs transverse momentum resolution smearing.
     21 *  Performs time smearing.
    2222 *
    23  *  \author P. Demin - UCL, Louvain-la-Neuve
     23 *  \author M. Selvaggi - CERN
    2424 *
    2525 */
    … …  
    4949
    5050using namespace std;
    51 
    5251//------------------------------------------------------------------------------
    5352
    5453TimeSmearing::TimeSmearing() :
    55   fItInputArray(0)
     54  fItTrackInputArray(0), fResolutionFormula(0)
    5655{
     56        fResolutionFormula = new DelphesFormula;
    5757}
    5858
    … …  
    6161TimeSmearing::~TimeSmearing()
    6262{
     63        if(fResolutionFormula) delete fResolutionFormula;
    6364}
    6465
    … …  
    6970  // read resolution formula
    7071
    71   fTimeResolution = GetDouble("TimeResolution", 1.0E-10);
    72   // import input array
     72  // read time resolution formula in seconds
     73  fResolutionFormula->Compile(GetString("TimeResolution", "30e-12"));
    7374
    74   fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
    75   fItInputArray = fInputArray->MakeIterator();
     75  // import track input array
     76  fTrackInputArray = ImportArray(GetString("TrackInputArray", "MuonMomentumSmearing/muons"));
     77  fItTrackInputArray = fTrackInputArray->MakeIterator();
     78
    7679
    7780  // create output array
    78 
    79   fOutputArray = ExportArray(GetString("OutputArray", "muons"));
     81  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
    8082}
    8183
    … …  
    8486void TimeSmearing::Finish()
    8587{
    86   if(fItInputArray) delete fItInputArray;
     88  if(fItTrackInputArray) delete fItTrackInputArray;
    8789}
    8890
    … …  
    9294{
    9395  Candidate *candidate, *mother;
    94   Double_t ti, tf_smeared, tf;
     96  Double_t tf_smeared, tf;
     97  Double_t eta, energy;
     98  Double_t timeResolution;
     99
    95100  const Double_t c_light = 2.99792458E8;
    96101
    97   fItInputArray->Reset();
    98   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
     102  fItTrackInputArray->Reset();
     103  while((candidate = static_cast<Candidate *>(fItTrackInputArray->Next())))
    99104  {
    100     const TLorentzVector &candidateInitialPosition = candidate->InitialPosition;
     105    // converting to meters
    101106    const TLorentzVector &candidateFinalPosition = candidate->Position;
     107    const TLorentzVector &candidateMomentum = candidate->Momentum;
    102108
    103     ti = candidateInitialPosition.T() * 1.0E-3 / c_light;
    104109    tf = candidateFinalPosition.T() * 1.0E-3 / c_light;
    105110
     111    eta = candidateMomentum.Eta();
     112    energy = candidateMomentum.E();
     113
    106114    // apply smearing formula
    107     tf_smeared = gRandom->Gaus(tf, fTimeResolution);
    108     ti = ti + tf_smeared - tf;
    109     tf = tf_smeared;
     115    timeResolution = fResolutionFormula->Eval(0.0, eta, 0.0, energy);
     116    tf_smeared = gRandom->Gaus(tf, timeResolution);
    110117
    111118    mother = candidate;
    112119    candidate = static_cast<Candidate *>(candidate->Clone());
    113     candidate->InitialPosition.SetT(ti * 1.0E3 * c_light);
    114     candidate->Position.SetT(tf * 1.0E3 * c_light);
    115120
    116     candidate->ErrorT = fTimeResolution * 1.0E3 * c_light;
     121    candidate->Position.SetT(tf_smeared * 1.0E3 * c_light);
     122    candidate->ErrorT = timeResolution * 1.0E3 * c_light;
    117123
    118124    candidate->AddCandidate(mother);
    119 
    120125    fOutputArray->Add(candidate);
    121126  }
    122127}
    123 
    124 //------------------------------------------------------------------------------

  • modules/TimeSmearing.h

    r46d3442 r4acf2fd  
    2222/** \class TimeSmearing
    2323 *
    24  *  Performs transverse time smearing.
     24 *  Performs time smearing.
    2525 *
    26  *  \author Michele Selvaggi - UCL, Louvain-la-Neuve
     26 *  \author Michele Selvaggi - CERN
    2727 *
    2828 */
    … …  
    3232class TIterator;
    3333class TObjArray;
     34class DelphesFormula;
    3435
    3536class TimeSmearing: public DelphesModule
    … …  
    4445
    4546private:
    46   Double_t fTimeResolution;
    4747
    48   TIterator *fItInputArray; //!
     48  DelphesFormula *fResolutionFormula;
     49  Int_t fVertexTimeMode;
    4950
    50   const TObjArray *fInputArray; //!
     51  TIterator *fItTrackInputArray; //!
     52
     53  const TObjArray *fTrackInputArray; //!
    5154
    5255  TObjArray *fOutputArray; //!

  • modules/TrackCovariance.cc

    r46d3442 r4acf2fd  
    104104  Double_t dd0, ddz, dphi, dct, dp, dpt, dC;
    105105
    106 
    107106  fItInputArray->Reset();
    108107  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))

  • modules/TreeWriter.cc

    r46d3442 r4acf2fd  
    356356    entry->DZ = candidate->DZ;
    357357    entry->Nclusters = candidate->Nclusters;
     358    entry->dNdx = candidate->dNdx;
    358359
    359360    entry->ErrorP = candidate->ErrorP;
    … …  
    404405    entry->Mass = m;
    405406
    406     particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    407     const TLorentzVector &initialPosition = particle->Position;
     407    //particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
     408    //const TLorentzVector &initialPosition = particle->Position;
     409    const TLorentzVector &initialPosition = candidate->InitialPosition;
    408410
    409411    entry->X = initialPosition.X();
    … …  
    411413    entry->Z = initialPosition.Z();
    412414    entry->T = initialPosition.T() * 1.0E-3 / c_light;
     415    entry->ErrorT =candidate-> ErrorT * 1.0E-3 / c_light;
    413416
    414417    entry->Particle = particle;
    … …  
    509512    entry->DZ = candidate->DZ;
    510513    entry->Nclusters = candidate->Nclusters;
     514    entry->dNdx = candidate->dNdx;
    511515
    512516    entry->ErrorP = candidate->ErrorP;
    … …  
    557561    entry->Mass = m;
    558562
    559     particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    560     const TLorentzVector &initialPosition = particle->Position;
     563    //particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
     564    //const TLorentzVector &initialPosition = particle->Position;
     565    const TLorentzVector &initialPosition = candidate->InitialPosition;
    561566
    562567    entry->X = initialPosition.X();
    … …  
    564569    entry->Z = initialPosition.Z();
    565570    entry->T = initialPosition.T() * 1.0E-3 / c_light;
     571    entry->ErrorT = candidate-> ErrorT * 1.0E-3 / c_light;
    566572
    567573    entry->VertexIndex = candidate->ClusterIndex;
    … …  
    574580    entry->Edges[3] = candidate->Edges[3];
    575581
    576     entry->T = position.T() * 1.0E-3 / c_light;
     582    //entry->T = position.T() * 1.0E-3 / c_light;
    577583    entry->NTimeHits = candidate->NTimeHits;
    578584
Note: See TracChangeset for help on using the changeset viewer.