Fork me on GitHub

Changeset 1345 in svn for trunk/modules


Ignore:
Timestamp:
Dec 21, 2013, 3:00:11 PM (11 years ago)
Author:
Michele Selvaggi
Message:

timing implemented

Location:
trunk/modules
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/modules/Calorimeter.cc

    r1280 r1345  
    330330      fTrackHCalEnergy = 0.0;
    331331
     332      fTowerECalTime = 0.0;
     333      fTowerHCalTime = 0.0;
     334
     335      fTrackECalTime = 0.0;
     336      fTrackHCalTime = 0.0;
     337
     338      fTowerECalWeightTime = 0.0;
     339      fTowerHCalWeightTime = 0.0;
     340     
    332341      fTowerTrackHits = 0;
    333342      fTowerPhotonHits = 0;
     
    343352      track = static_cast<Candidate*>(fTrackInputArray->At(number));
    344353      momentum = track->Momentum;
    345 
     354      position = track->Position;
     355
     356     
    346357      ecalEnergy = momentum.E() * fTrackECalFractions[number];
    347358      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
     
    349360      fTrackECalEnergy += ecalEnergy;
    350361      fTrackHCalEnergy += hcalEnergy;
     362     
     363      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
     364      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
     365       
     366      fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
     367      fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
    351368
    352369      fTowerTrackArray->Add(track);
     
    360377    particle = static_cast<Candidate*>(fParticleInputArray->At(number));
    361378    momentum = particle->Momentum;
     379    position = particle->Position;
    362380
    363381    // fill current tower
     
    367385    fTowerECalEnergy += ecalEnergy;
    368386    fTowerHCalEnergy += hcalEnergy;
     387
     388    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
     389    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
     390
     391    fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
     392    fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
     393   
    369394
    370395    fTower->AddCandidate(particle);
     
    383408  Double_t ecalEnergy, hcalEnergy;
    384409  Double_t ecalSigma, hcalSigma;
     410  Double_t ecalTime, hcalTime, time;
    385411
    386412  if(!fTower) return;
     
    392418
    393419  ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
     420  ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
    394421
    395422  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
     
    399426
    400427  hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
     428  hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
    401429
    402430  energy = ecalEnergy + hcalEnergy;
     431  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
    403432
    404433//  eta = fTowerEta;
     
    410439  pt = energy / TMath::CosH(eta);
    411440
    412   fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
     441 // fTower->Position.SetXYZT(-time, 0.0, 0.0, time);
     442  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
    413443  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    414444  fTower->Eem = ecalEnergy;
     
    420450  fTower->Edges[3] = fTowerEdges[3];
    421451
     452
    422453  // fill calorimeter towers and photon candidates
    423454  if(energy > 0.0)
  • trunk/modules/Calorimeter.h

    r1273 r1345  
    4545  Double_t fTowerECalEnergy, fTowerHCalEnergy;
    4646  Double_t fTrackECalEnergy, fTrackHCalEnergy;
     47 
     48  Double_t fTowerECalTime, fTowerHCalTime;
     49  Double_t fTrackECalTime, fTrackHCalTime;
     50   
     51  Double_t fTowerECalWeightTime, fTowerHCalWeightTime;
     52  Double_t fTrackECalWeightTime, fTrackHCalWeightTime;
     53 
    4754  Int_t fTowerTrackHits, fTowerPhotonHits;
    4855
  • trunk/modules/FastJetFinder.cc

    r1335 r1345  
    193193  TLorentzVector momentum;
    194194  Double_t deta, dphi, detaMax, dphiMax;
     195  Double_t time, weightTime, avTime;
    195196  Int_t number;
    196197  Double_t rho = 0;
     
    259260    candidate = factory->NewCandidate();
    260261
     262    time=0;
     263    weightTime=0;
     264
    261265    inputList.clear();
    262266    inputList = sequence->constituents(*itOutputList);
     
    269273      if(deta > detaMax) detaMax = deta;
    270274      if(dphi > dphiMax) dphiMax = dphi;
    271 
     275     
     276      time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
     277      weightTime += TMath::Sqrt(constituent->Momentum.E());
     278   
    272279      candidate->AddCandidate(constituent);
    273280    }
     281   
     282    avTime = time/weightTime;
    274283
    275284    candidate->Momentum = momentum;
     285    candidate->Position.SetT(avTime);
    276286    candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
    277287
  • trunk/modules/ModulesLinkDef.h

    r1177 r1345  
    1919#include "modules/EnergySmearing.h"
    2020#include "modules/MomentumSmearing.h"
     21#include "modules/TimeSmearing.h"
    2122#include "modules/Calorimeter.h"
    2223#include "modules/Isolation.h"
     
    5051#pragma link C++ class EnergySmearing+;
    5152#pragma link C++ class MomentumSmearing+;
     53#pragma link C++ class TimeSmearing+;
    5254#pragma link C++ class Calorimeter+;
    5355#pragma link C++ class Isolation+;
  • trunk/modules/PileUpMerger.cc

    r1323 r1345  
    1616#include "classes/DelphesClasses.h"
    1717#include "classes/DelphesFactory.h"
    18 #include "classes/DelphesFormula.h"
     18#include "classes/DelphesTF2.h"
    1919#include "classes/DelphesPileUpReader.h"
    2020
     
    4141
    4242PileUpMerger::PileUpMerger() :
    43   fReader(0), fItInputArray(0)
    44 {
    45 }
     43  fFunction(0), fReader(0), fItInputArray(0)
     44{
     45  fFunction = new DelphesTF2;
     46}
     47
    4648
    4749//------------------------------------------------------------------------------
     
    4951PileUpMerger::~PileUpMerger()
    5052{
     53  delete fFunction;
    5154}
    5255
     
    6063
    6164  fMeanPileUp  = GetDouble("MeanPileUp", 10);
    62   fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
    63 
     65 
     66  fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
     67  fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
     68
     69  // read vertex smearing formula
     70       
     71  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
     72  fFunction->SetRange(-fZVertexSpread,-fTVertexSpread,fZVertexSpread,fTVertexSpread);
     73 
    6474  fileName = GetString("PileUpFile", "MinBias.pileup");
    6575  fReader = new DelphesPileUpReader(fileName);
     
    90100  Float_t x, y, z, t;
    91101  Float_t px, py, pz, e;
    92   Double_t dz, dphi;
     102  Double_t dz, dphi, dt;
    93103  Int_t numberOfEvents, event;
    94104  Long64_t allEntries, entry;
    95   Candidate *candidate;
     105  Candidate *candidate, *vertexcandidate;
    96106  DelphesFactory *factory;
     107 
     108  const Double_t c_light = 2.99792458E8;
    97109
    98110  fItInputArray->Reset();
     111 
     112    // --- Deal with Primary vertex first  ------
     113 
     114  fFunction->GetRandom2(dz,dt);
     115 
     116  dt *= c_light*1.0E3; //necessary in order to make t in mm/c
     117  dz *= 1.0E3; //necessary in order to make z in mm
     118 
    99119  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    100120  {
     121    candidate->Position.SetXYZT(x, y, z+dz, t+dt);
    101122    fParticleOutputArray->Add(candidate);
    102123  }
    103124
    104125  factory = GetFactory();
    105 
     126 
     127  vertexcandidate = factory->NewCandidate();
     128  vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
     129  fVertexOutputArray->Add(vertexcandidate);
     130 
     131 
     132  // --- Then with pile-up vertices  ------
     133 
    106134  switch(fPileUpDistribution)
    107135  {
     
    118146
    119147  allEntries = fReader->GetEntries();
    120 
     148 
    121149  for(event = 0; event < numberOfEvents; ++event)
    122150  {
     
    128156
    129157    fReader->ReadEntry(entry);
    130 
    131     dz = gRandom->Gaus(0.0, fZVertexSpread);
     158 
     159   // --- Pile-up vertex smearing
     160   
     161    fFunction->GetRandom2(dz,dt);
     162 
     163    dt *= c_light*1.0E3; //necessary in order to make t in mm/c
     164    dz *= 1.0E3; //necessary in order to make z in mm
     165 
    132166    dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    133167
    134     candidate = factory->NewCandidate();
    135     candidate->Position.SetXYZT(0.0, 0.0, dz, 0.0);
    136     fVertexOutputArray->Add(candidate);
     168    vertexcandidate = factory->NewCandidate();
     169    vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
     170    vertexcandidate->IsPU = 1;
     171   
     172    fVertexOutputArray->Add(vertexcandidate);
    137173
    138174    while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
     
    149185
    150186      candidate->IsPU = 1;
    151 
     187   
    152188      candidate->Momentum.SetPxPyPzE(px, py, pz, e);
    153189      candidate->Momentum.RotateZ(dphi);
    154190
    155       candidate->Position.SetXYZT(x, y, z + dz, t);
     191      candidate->Position.SetXYZT(x, y, z+dz, t+dt);
    156192      candidate->Position.RotateZ(dphi);
    157 
     193     
    158194      fParticleOutputArray->Add(candidate);
    159195    }
  • trunk/modules/PileUpMerger.h

    r1323 r1345  
    1919class TObjArray;
    2020class DelphesPileUpReader;
     21class DelphesTF2;
    2122
    2223class PileUpMerger: public DelphesModule
     
    3536  Int_t fPileUpDistribution;
    3637  Double_t fMeanPileUp;
     38 
    3739  Double_t fZVertexSpread;
     40  Double_t fTVertexSpread;
    3841
    3942  DelphesPileUpReader *fReader;
    40 
     43 
     44  DelphesTF2 *fFunction; //!
     45 
    4146  TIterator *fItInputArray; //!
    4247
  • trunk/modules/TrackPileUpSubtractor.cc

    r1314 r1345  
    5353void TrackPileUpSubtractor::Init()
    5454{
     55// import input array
     56
     57  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
     58  fItVertexInputArray = fVertexInputArray->MakeIterator();
     59 
    5560  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
    5661
     
    8590    if(iterator) delete iterator;
    8691  }
     92
     93  if(fItVertexInputArray) delete fItVertexInputArray;
    8794}
    8895
     
    95102  TIterator *iterator;
    96103  TObjArray *array;
    97   Double_t z;
     104  Double_t z, zvtx;
     105
     106 
     107  // find z position of primary vertex
     108 
     109  fItVertexInputArray->Reset();
     110  while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
     111  {
     112    if(candidate->IsPU == 0)
     113    {
     114    zvtx = candidate->Position.Z();
     115    break;
     116    }
     117  }
     118
    98119
    99120  // loop over all input arrays
     
    112133      // apply pile-up subtraction
    113134      // assume perfect pile-up subtraction for tracks outside fZVertexResolution
    114       if(candidate->IsPU && TMath::Abs(z) > fZVertexResolution) continue;
     135      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
    115136
    116137      array->Add(candidate);
  • trunk/modules/TreeWriter.cc

    r1325 r1345  
    9595      continue;
    9696    }
    97 
     97   
    9898    itClassMap = fClassMap.find(branchClass);
    9999    if(itClassMap == fClassMap.end())
     
    115115void TreeWriter::Finish()
    116116{
    117 
    118117}
    119118
     
    161160  GenParticle *entry = 0;
    162161  Double_t pt, signPz, cosTheta, eta, rapidity;
    163 
     162 
     163  const Double_t c_light = 2.99792458E8;
     164 
    164165  // loop over all particles
    165166  iterator.Reset();
     
    208209    entry->Y = position.Y();
    209210    entry->Z = position.Z();
    210     entry->T = position.T();
     211    entry->T = position.T()*1.0E-3/c_light;
    211212  }
    212213}
     
    219220  Candidate *candidate = 0;
    220221  Vertex *entry = 0;
     222 
     223  const Double_t c_light = 2.99792458E8;
    221224
    222225  // loop over all vertices
     
    231234    entry->Y = position.Y();
    232235    entry->Z = position.Z();
    233     entry->T = position.T();
     236    entry->T = position.T()*1.0E-3/c_light;
    234237  }
    235238}
     
    244247  Track *entry = 0;
    245248  Double_t pt, signz, cosTheta, eta, rapidity;
    246 
     249  const Double_t c_light = 2.99792458E8;
     250 
    247251  // loop over all tracks
    248252  iterator.Reset();
     
    271275    entry->YOuter = position.Y();
    272276    entry->ZOuter = position.Z();
     277    entry->TOuter = position.T()*1.0E-3/c_light;
    273278
    274279    const TLorentzVector &momentum = candidate->Momentum;
     
    290295    entry->Y = initialPosition.Y();
    291296    entry->Z = initialPosition.Z();
     297    entry->T = initialPosition.T()*1.0E-3/c_light;
    292298
    293299    entry->Particle = particle;
     
    303309  Tower *entry = 0;
    304310  Double_t pt, signPz, cosTheta, eta, rapidity;
    305 
     311  const Double_t c_light = 2.99792458E8;
     312 
    306313  // loop over all towers
    307314  iterator.Reset();
     
    309316  {
    310317    const TLorentzVector &momentum = candidate->Momentum;
    311 
     318    const TLorentzVector &position = candidate->Position;
     319   
    312320    pt = momentum.Pt();
    313321    cosTheta = TMath::Abs(momentum.CosTheta());
     
    331339    entry->Edges[2] = candidate->Edges[2];
    332340    entry->Edges[3] = candidate->Edges[3];
    333 
     341   
     342    entry->T = position.T()*1.0E-3/c_light;
     343   
    334344    FillParticles(candidate, &entry->Particles);
    335345  }
     
    344354  Photon *entry = 0;
    345355  Double_t pt, signPz, cosTheta, eta, rapidity;
    346 
     356  const Double_t c_light = 2.99792458E8;
     357 
    347358  array->Sort();
    348359
     
    353364    TIter it1(candidate->GetCandidates());
    354365    const TLorentzVector &momentum = candidate->Momentum;
     366    const TLorentzVector &position = candidate->Position;
     367   
    355368
    356369    pt = momentum.Pt();
     
    366379    entry->PT = pt;
    367380    entry->E = momentum.E();
    368 
     381   
     382    entry->T = position.T()*1.0E-3/c_light;
     383   
    369384    entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9;
    370385
     
    381396  Electron *entry = 0;
    382397  Double_t pt, signPz, cosTheta, eta, rapidity;
    383 
     398  const Double_t c_light = 2.99792458E8;
     399 
    384400  array->Sort();
    385401
     
    389405  {
    390406    const TLorentzVector &momentum = candidate->Momentum;
    391 
     407    const TLorentzVector &position = candidate->Position;
     408   
    392409    pt = momentum.Pt();
    393410    cosTheta = TMath::Abs(momentum.CosTheta());
     
    401418    entry->Phi = momentum.Phi();
    402419    entry->PT = pt;
    403 
     420   
     421    entry->T = position.T()*1.0E-3/c_light;
     422   
    404423    entry->Charge = candidate->Charge;
    405424
     
    418437  Muon *entry = 0;
    419438  Double_t pt, signPz, cosTheta, eta, rapidity;
    420 
     439 
     440  const Double_t c_light = 2.99792458E8;
     441 
    421442  array->Sort();
    422443
     
    426447  {
    427448    const TLorentzVector &momentum = candidate->Momentum;
     449    const TLorentzVector &position = candidate->Position;
     450   
    428451
    429452    pt = momentum.Pt();
     
    442465    entry->PT = pt;
    443466
     467    entry->T = position.T()*1.0E-3/c_light;
     468   
     469   // cout<<entry->PT<<","<<entry->T<<endl;
     470
    444471    entry->Charge = candidate->Charge;
    445472
     
    457484  Double_t pt, signPz, cosTheta, eta, rapidity;
    458485  Double_t ecalEnergy, hcalEnergy;
    459 
     486  const Double_t c_light = 2.99792458E8;
     487 
    460488  array->Sort();
    461489
     
    465493  {
    466494    TIter itConstituents(candidate->GetCandidates());
    467     const TLorentzVector &momentum = candidate->Momentum;
    468 
     495   
     496    const TLorentzVector &momentum = candidate->Momentum;
     497    const TLorentzVector &position = candidate->Position;
     498   
    469499    pt = momentum.Pt();
    470500    cosTheta = TMath::Abs(momentum.CosTheta());
     
    479509    entry->PT = pt;
    480510
     511    entry->T = position.T()*1.0E-3/c_light;
     512   
    481513    entry->Mass = momentum.M();
    482514
Note: See TracChangeset for help on using the changeset viewer.