Fork me on GitHub

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


Ignore:
Files:
4 deleted
15 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r4acf2fd r46d3442  
    412412        modules/ImpactParameterSmearing.h \
    413413        modules/TimeSmearing.h \
    414         modules/TimeOfFlight.h \
    415414        modules/SimpleCalorimeter.h \
    416415        modules/DenseTrackFilter.h \
    … …  
    449448        modules/DecayFilter.h \
    450449        modules/ParticleDensity.h \
    451         modules/TruthVertexFinder.h \
    452450        modules/ExampleModule.h
    453451tmp/modules/ModulesDict$(PcmSuf): \
    … …  
    10141012        classes/DelphesFactory.h \
    10151013        classes/DelphesFormula.h
    1016 tmp/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
    10251014tmp/modules/TimeSmearing.$(ObjSuf): \
    10261015        modules/TimeSmearing.$(SrcSuf) \
    … …  
    10821071        external/ExRootAnalysis/ExRootResult.h \
    10831072        external/ExRootAnalysis/ExRootTreeBranch.h
    1084 tmp/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
    10941073tmp/modules/UniqueObjectFinder.$(ObjSuf): \
    10951074        modules/UniqueObjectFinder.$(SrcSuf) \
    … …  
    12361215        tmp/modules/TaggingParticlesSkimmer.$(ObjSuf) \
    12371216        tmp/modules/TauTagging.$(ObjSuf) \
    1238         tmp/modules/TimeOfFlight.$(ObjSuf) \
    12391217        tmp/modules/TimeSmearing.$(ObjSuf) \
    12401218        tmp/modules/TrackCountingBTagging.$(ObjSuf) \
    … …  
    12441222        tmp/modules/TrackSmearing.$(ObjSuf) \
    12451223        tmp/modules/TreeWriter.$(ObjSuf) \
    1246         tmp/modules/TruthVertexFinder.$(ObjSuf) \
    12471224        tmp/modules/UniqueObjectFinder.$(ObjSuf) \
    12481225        tmp/modules/VertexFinder.$(ObjSuf) \
    … …  
    20902067        @touch $@
    20912068
    2092 modules/TimeOfFlight.h: \
    2093         classes/DelphesModule.h
    2094         @touch $@
    2095 
    20962069external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh: \
    20972070        external/fastjet/ClusterSequence.hh \
    … …  
    22122185        @touch $@
    22132186
    2214 modules/TruthVertexFinder.h: \
    2215         classes/DelphesModule.h
    2216         @touch $@
    2217 
    22182187classes/DelphesSTDHEPReader.h: \
    22192188        classes/DelphesXDRReader.h
    … …  
    23312300        @touch $@
    23322301
     2302modules/TimeSmearing.h: \
     2303        classes/DelphesModule.h
     2304        @touch $@
     2305
    23332306modules/TreeWriter.h: \
    2334         classes/DelphesModule.h
    2335         @touch $@
    2336 
    2337 modules/TimeSmearing.h: \
    23382307        classes/DelphesModule.h
    23392308        @touch $@

  • cards/delphes_card_IDEA.tcl

    r4acf2fd r46d3442  
    1111set B 2.0
    1212
    13 ## Drift chamber coordinates
    14 set DCHZMIN -2.125
    15 set DCHZMAX 2.125
    16 set DCHRMIN 0.345
    17 set DCHRMAX 2.02
    18 
    19 
    2013#######################################
    2114# Order of execution of various modules
    … …  
    2316
    2417set ExecutionPath {
    25 
    26   TruthVertexFinder
    2718  ParticlePropagator
    2819
    … …  
    3425  TrackSmearing
    3526  ClusterCounting
    36   TimeSmearing
    37   TimeOfFlight
    3827
    3928  TrackMerger
    … …  
    7261  ScalarHT
    7362  TreeWriter
    74 }
    75 
    76 #################################
    77 # Truth Vertex Finder
    78 #################################
    79 
    80 module 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
    8763}
    8864
    … …  
    125101    }
    126102}
     103
     104#       (pt <= 0.1)                                     * (0.00) +
     105#       (abs(eta) <= 3.0)               * (pt > 0.1)    * (1.00) +
     106#       (abs(eta) > 3)                                  * (0.00)
    127107
    128108
    … …  
    232212      2        VTXDSK      0.141  0.3    0.92     0.00028   0.0937     2        0          1.5708           7e-006        7e-006         1
    233213
    234       1 DCHCANI $DCHZMIN $DCHZMAX $DCHRMIN 0.0002 0.237223 0 0 0 0 0 0
     214      1 DCHCANI -2.125 2.125 0.345 0.0002 0.237223 0 0 0 0 0 0
    235215      1 DCH -2 2 0.36 0.0147748 1400 1 0.0203738 0 0.0001 0 1
    236216      1 DCH -2 2 0.374775 0.0147748 1400 1 -0.0212097 0 0.0001 0 1
    … …  
    345325      1 DCH -2 2 1.98523 0.0147748 1400 1 0.111898 0 0.0001 0 1
    346326      1 DCH -2 2 2 0.0147748 1400 1 -0.112723 0 0.0001 0 1
    347       1 DCHCANO $DCHZMIN $DCHZMAX $DCHRMAX $DCHRMAX 0.02 1.667 0 0 0 0 0 0
     327      1 DCHCANO -2.125 2.125 2.02 0.02 1.667 0 0 0 0 0 0
    348328      1 BSILWRP -2.35 2.35 2.04 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    349329      1 BSILWRP -2.35 2.35 2.06 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    350330      1 MAG -2.5 2.5 2.25 0.05 0.0658 0 0 0 0 0 0
    351331      1 BPRESH -2.55 2.55 2.45 0.02 1 2 0 1.5708 7e-005 0.01 1
    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
     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
    354334      2 FSILWRP 0.354 2.02 -2.32 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    355335      2 FSILWRP 0.35 2.02 -2.3 0.00047 0.0937 2 0 1.5708 7e-006 9e-005 1
    … …  
    375355  set Bz $B
    376356
    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:
     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
    384363  # 0:  Helium 90% - Isobutane 10%
    385364  # 1:  Helium 100%
    … …  
    392371
    393372
    394 ########################################
    395 #   Time Smearing MIP
    396 ########################################
    397 
    398 module 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 
    412 module 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 
    425373
    426374##############
    … …  
    430378module Merger TrackMerger {
    431379# add InputArray InputArray
    432   add InputArray TimeOfFlight/tracks
     380  add InputArray ClusterCounting/tracks
    433381  set OutputArray tracks
    434382}
    … …  
    831779
    832780module UniqueObjectFinder UniqueObjectFinder {
    833 
    834   set UseUniqueID true
    835781# earlier arrays take precedence over later ones
    836782# add InputArray InputArray OutputArray
    … …  
    842788
    843789
     790
    844791##################
    845792# ROOT tree writer

  • classes/DelphesClasses.cc

    r4acf2fd r46d3442  
    231231  Xd(0), Yd(0), Zd(0),
    232232  Nclusters(0.0),
    233   dNdx(0.0),
    234233  TrackResolution(0),
    235234  NCharged(0),
    … …  
    409408  object.Zd = Zd;
    410409  object.Nclusters = Nclusters;
    411   object.dNdx = dNdx;
    412410  object.TrackResolution = TrackResolution;
    413411  object.NCharged = NCharged;
    … …  
    543541  Zd = 0.0;
    544542  Nclusters = 0.0;
    545   dNdx = 0.0;
    546543  TrackResolution = 0.0;
    547544  NCharged = 0;

  • classes/DelphesClasses.h

    r4acf2fd r46d3442  
    462462  Float_t DZ; // track longitudinal impact parameter
    463463  Float_t Nclusters; // Number of ionization clusters
    464   Float_t dNdx; // Number of ionization clusters
    465464
    466465  Float_t ErrorP; // track momentum error
    … …  
    568567  Float_t DZ; // track longitudinal impact parameter
    569568  Float_t Nclusters; // Number of ionization clusters
    570   Float_t dNdx; // Number of ionization clusters
    571569
    572570  Float_t ErrorP; // track momentum error
    … …  
    699697
    700698  Float_t Nclusters; // Number of ionization clusters
    701   Float_t dNdx; // Number of ionization clusters per unit length
    702699
    703700  Float_t Xd;

  • examples/Example5.C

    r4acf2fd r46d3442  
    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<<", 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;
     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;
    5151    }
    5252  }

  • examples/Pythia8/ee_zh.cmnd

    r4acf2fd r46d3442  
    11
    2 ! number of events to generate
    3 Main:numberOfEvents = 1000         ! number of events to generate
     2Main:numberOfEvents = 10         ! number of events to generate
    43
    54Beams:idA = 11                   ! first beam, e- = -11
    … …  
    87
    98! Vertex smearing :
     9Beams: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ï¿Œ
    1015
    11 Beams:allowVertexSpread = on
    12 Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
    13 Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
    14 Beams:sigmaVertexZ = 0.64      !  0.64 mm
    15 Beams:sigmaTime    = 0.64      !  0.64 mm
    1616
    1717! Higgsstrahlung process

  • modules/ClusterCounting.cc

    r4acf2fd r46d3442  
    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
    3034#include "modules/ClusterCounting.h"
     35
    3136#include "classes/DelphesClasses.h"
    3237#include "TrackCovariance/TrkUtil.h"
    … …  
    103108{
    104109  Candidate *candidate, *mother, *particle;
    105   Double_t mass, trackLength, Ncl;
     110  Double_t mass, Ncl;
    106111
    107112  fItInputArray->Reset();
    … …  
    119124    mass = candidateMomentum.M();
    120125
    121     trackLength = fTrackUtil->TrkLen(Par);
    122 
    123126    mother    = candidate;
    124127    candidate = static_cast<Candidate*>(candidate->Clone());
    125128
    126     Ncl = 0.;
    127     if (fTrackUtil->IonClusters(Ncl, mass, Par))
     129    Ncl = -999;
     130    // computation of Nclusters is not supported for electrons
     131    if (TMath::Abs(particle->PID) == 11)
    128132    {
    129133      candidate->Nclusters = Ncl;
    130       candidate->dNdx = (trackLength > 0.) ? Ncl/trackLength : -1;
    131134    }
     135    else if (fTrackUtil->IonClusters(Ncl, mass, Par))
     136    {
     137      candidate->Nclusters = Ncl;
     138    }
     139    //cout<<candidate->PID<<", "<<mass<<", "<<candidate->Nclusters<<endl;
    132140
    133141    candidate->AddCandidate(mother);

  • modules/DualReadoutCalorimeter.cc

    r4acf2fd r46d3442  
    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       fEFlowPhotonOutputArray->Add(tower);
    626     }
     625    }
     626
    627627    // if hadronic fraction > 0, use HCAL resolution
    628628    else
    … …  
    631631      tower->Ehad = neutralEnergy;
    632632      tower->PID = 130;
    633       fEFlowNeutralHadronOutputArray->Add(tower);
    634     }
     633    }
     634
     635    fEFlowPhotonOutputArray->Add(tower);
     636
    635637
    636638    //clone tracks

  • modules/Efficiency.cc

    r4acf2fd r46d3442  
    100100  Double_t pt, eta, phi, e;
    101101
    102 
    103102  fItInputArray->Reset();
    104103  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    … …  
    116115    pt = candidateMomentum.Pt();
    117116    e = candidateMomentum.E();
    118 
     117   
    119118    // apply an efficency formula
    120119    if(gRandom->Uniform() > fFormula->Eval(pt, eta, phi, e, candidate)) continue;

  • modules/ModulesLinkDef.h

    r4acf2fd r46d3442  
    4040#include "modules/ImpactParameterSmearing.h"
    4141#include "modules/TimeSmearing.h"
    42 #include "modules/TimeOfFlight.h"
    4342#include "modules/SimpleCalorimeter.h"
    4443#include "modules/DenseTrackFilter.h"
    … …  
    7776#include "modules/DecayFilter.h"
    7877#include "modules/ParticleDensity.h"
    79 #include "modules/TruthVertexFinder.h"
    8078#include "modules/ExampleModule.h"
    8179
    … …  
    10098#pragma link C++ class ImpactParameterSmearing+;
    10199#pragma link C++ class TimeSmearing+;
    102 #pragma link C++ class TimeOfFlight+;
    103100#pragma link C++ class SimpleCalorimeter+;
    104101#pragma link C++ class DenseTrackFilter+;
    … …  
    137134#pragma link C++ class DecayFilter+;
    138135#pragma link C++ class ParticleDensity+;
    139 #pragma link C++ class TruthVertexFinder+;
    140136#pragma link C++ class ExampleModule+;
    141137

  • modules/ParticlePropagator.cc

    r4acf2fd r46d3442  
    125125  TLorentzVector particlePosition, particleMomentum, beamSpotPosition;
    126126  Double_t px, py, pz, pt, pt2, e, q;
    127   Double_t x, y, z, t, r;
     127  Double_t x, y, z, t, r, phi;
    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 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;
     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;
    135136  Double_t bsx, bsy, bsz;
    136   Double_t rxp, rdp, t_R;
    137   Double_t td, pio, phid, sign_pz, vz;
    138 
     137  Double_t s0, s1, sd;
     138 
    139139  const Double_t c_light = 2.99792458E8;
    140140
    … …  
    161161    particlePosition = particle->Position;
    162162    particleMomentum = particle->Momentum;
    163 
    164     // Constants
    165 
    166163    x = particlePosition.X() * 1.0E-3;
    167164    y = particlePosition.Y() * 1.0E-3;
    … …  
    208205    else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
    209206    {
    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));
     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));
    227236
    228237      mother = candidate;
    229       candidate = static_cast<Candidate*>(candidate->Clone());
     238      candidate = static_cast<Candidate *>(candidate->Clone());
    230239
    231240      candidate->InitialPosition = particlePosition;
    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;
     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;
    234243
    235244      candidate->Momentum = particleMomentum;
    … …  
    237246
    238247      fOutputArray->Add(candidate);
    239 
    240248      if(TMath::Abs(q) > 1.0E-9)
    241249      {
    … …  
    266274      //     helix radius r = p_{T0} / (omega gamma m)
    267275
    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]
     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]
    271279
    272280      phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
    273281
    274282      // 2. helix axis coordinates
    275       x_c = x + r*TMath::Sin(phi_0);
    276       y_c = y - r*TMath::Cos(phi_0);
     283      x_c = x + r * TMath::Sin(phi_0);
     284      y_c = y - r * TMath::Cos(phi_0);
    277285      r_c = TMath::Hypot(x_c, y_c);
    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());
     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());
    303318
    304319      // calculate additional track parameters (correct for beamspot position)
    305       d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
    306       dz = zd - bsz;
    307       ctgTheta  = 1.0 / TMath::Tan (particleMomentum.Theta());
     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());
    308325
    309326      // 3. time evaluation t = TMath::Min(t_r, t_z)
    310327      //    t_r : time to exit from the sides
    311328      //    t_z : time to exit from the front or the back
    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       {
     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;
     333      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
    320339        t = t_z;
    321340      }
    322341      else
    323342      {
    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 
     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);
    328364        t = TMath::Min(t_r, t_z);
    329365      }
    330366
    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);
     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);
    335372
    336373      // compute path length for an helix
    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);
     374
     375      alpha = pz * 1.0E9 / c_light / gammam;
     376      l = t * TMath::Sqrt(alpha * alpha + r * r * omega * omega);
    340377
    341378      if(r_t > 0.0)
    342379      {
     380
    343381        // store these variables before cloning
    344382        if(particle == candidate)
    … …  
    346384          particle->D0 = d0 * 1.0E3;
    347385          particle->DZ = dz * 1.0E3;
    348           particle->P = particleMomentum.P();
     386          particle->P = p;
    349387          particle->PT = pt;
    350388          particle->CtgTheta = ctgTheta;
    351           particle->Phi = particleMomentum.Phi();
     389          particle->Phi = phip;
    352390        }
    353391

  • modules/TimeSmearing.cc

    r4acf2fd r46d3442  
    1919/** \class TimeSmearing
    2020 *
    21  *  Performs time smearing.
     21 *  Performs transverse momentum resolution smearing.
    2222 *
    23  *  \author M. Selvaggi - CERN
     23 *  \author P. Demin - UCL, Louvain-la-Neuve
    2424 *
    2525 */
    … …  
    4949
    5050using namespace std;
     51
    5152//------------------------------------------------------------------------------
    5253
    5354TimeSmearing::TimeSmearing() :
    54   fItTrackInputArray(0), fResolutionFormula(0)
     55  fItInputArray(0)
    5556{
    56         fResolutionFormula = new DelphesFormula;
    5757}
    5858
    … …  
    6161TimeSmearing::~TimeSmearing()
    6262{
    63         if(fResolutionFormula) delete fResolutionFormula;
    6463}
    6564
    … …  
    7069  // read resolution formula
    7170
    72   // read time resolution formula in seconds
    73   fResolutionFormula->Compile(GetString("TimeResolution", "30e-12"));
     71  fTimeResolution = GetDouble("TimeResolution", 1.0E-10);
     72  // import input array
    7473
    75   // import track input array
    76   fTrackInputArray = ImportArray(GetString("TrackInputArray", "MuonMomentumSmearing/muons"));
    77   fItTrackInputArray = fTrackInputArray->MakeIterator();
    78 
     74  fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
     75  fItInputArray = fInputArray->MakeIterator();
    7976
    8077  // create output array
    81   fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
     78
     79  fOutputArray = ExportArray(GetString("OutputArray", "muons"));
    8280}
    8381
    … …  
    8684void TimeSmearing::Finish()
    8785{
    88   if(fItTrackInputArray) delete fItTrackInputArray;
     86  if(fItInputArray) delete fItInputArray;
    8987}
    9088
    … …  
    9492{
    9593  Candidate *candidate, *mother;
    96   Double_t tf_smeared, tf;
    97   Double_t eta, energy;
    98   Double_t timeResolution;
    99 
     94  Double_t ti, tf_smeared, tf;
    10095  const Double_t c_light = 2.99792458E8;
    10196
    102   fItTrackInputArray->Reset();
    103   while((candidate = static_cast<Candidate *>(fItTrackInputArray->Next())))
     97  fItInputArray->Reset();
     98  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    10499  {
    105     // converting to meters
     100    const TLorentzVector &candidateInitialPosition = candidate->InitialPosition;
    106101    const TLorentzVector &candidateFinalPosition = candidate->Position;
    107     const TLorentzVector &candidateMomentum = candidate->Momentum;
    108102
     103    ti = candidateInitialPosition.T() * 1.0E-3 / c_light;
    109104    tf = candidateFinalPosition.T() * 1.0E-3 / c_light;
    110105
    111     eta = candidateMomentum.Eta();
    112     energy = candidateMomentum.E();
    113 
    114106    // apply smearing formula
    115     timeResolution = fResolutionFormula->Eval(0.0, eta, 0.0, energy);
    116     tf_smeared = gRandom->Gaus(tf, timeResolution);
     107    tf_smeared = gRandom->Gaus(tf, fTimeResolution);
     108    ti = ti + tf_smeared - tf;
     109    tf = tf_smeared;
    117110
    118111    mother = candidate;
    119112    candidate = static_cast<Candidate *>(candidate->Clone());
     113    candidate->InitialPosition.SetT(ti * 1.0E3 * c_light);
     114    candidate->Position.SetT(tf * 1.0E3 * c_light);
    120115
    121     candidate->Position.SetT(tf_smeared * 1.0E3 * c_light);
    122     candidate->ErrorT = timeResolution * 1.0E3 * c_light;
     116    candidate->ErrorT = fTimeResolution * 1.0E3 * c_light;
    123117
    124118    candidate->AddCandidate(mother);
     119
    125120    fOutputArray->Add(candidate);
    126121  }
    127122}
     123
     124//------------------------------------------------------------------------------

  • modules/TimeSmearing.h

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

  • modules/TrackCovariance.cc

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

  • modules/TreeWriter.cc

    r4acf2fd r46d3442  
    356356    entry->DZ = candidate->DZ;
    357357    entry->Nclusters = candidate->Nclusters;
    358     entry->dNdx = candidate->dNdx;
    359358
    360359    entry->ErrorP = candidate->ErrorP;
    … …  
    405404    entry->Mass = m;
    406405
    407     //particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    408     //const TLorentzVector &initialPosition = particle->Position;
    409     const TLorentzVector &initialPosition = candidate->InitialPosition;
     406    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
     407    const TLorentzVector &initialPosition = particle->Position;
    410408
    411409    entry->X = initialPosition.X();
    … …  
    413411    entry->Z = initialPosition.Z();
    414412    entry->T = initialPosition.T() * 1.0E-3 / c_light;
    415     entry->ErrorT =candidate-> ErrorT * 1.0E-3 / c_light;
    416413
    417414    entry->Particle = particle;
    … …  
    512509    entry->DZ = candidate->DZ;
    513510    entry->Nclusters = candidate->Nclusters;
    514     entry->dNdx = candidate->dNdx;
    515511
    516512    entry->ErrorP = candidate->ErrorP;
    … …  
    561557    entry->Mass = m;
    562558
    563     //particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    564     //const TLorentzVector &initialPosition = particle->Position;
    565     const TLorentzVector &initialPosition = candidate->InitialPosition;
     559    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
     560    const TLorentzVector &initialPosition = particle->Position;
    566561
    567562    entry->X = initialPosition.X();
    … …  
    569564    entry->Z = initialPosition.Z();
    570565    entry->T = initialPosition.T() * 1.0E-3 / c_light;
    571     entry->ErrorT = candidate-> ErrorT * 1.0E-3 / c_light;
    572566
    573567    entry->VertexIndex = candidate->ClusterIndex;
    … …  
    580574    entry->Edges[3] = candidate->Edges[3];
    581575
    582     //entry->T = position.T() * 1.0E-3 / c_light;
     576    entry->T = position.T() * 1.0E-3 / c_light;
    583577    entry->NTimeHits = candidate->NTimeHits;
    584578
Note: See TracChangeset for help on using the changeset viewer.