Changeset 22dc7fd in git for modules

Timestamp:

Dec 21, 2013, 3:00:11 PM (11 years ago)

Author:

mselvaggi <mselvaggi@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

24d005f

Parents:

af88093

Message:

timing implemented

Location:

modules

Files:

• Calorimeter.cc (modified) (10 diffs)
• Calorimeter.h (modified) (1 diff)
• FastJetFinder.cc (modified) (3 diffs)
• ModulesLinkDef.h (modified) (2 diffs)
• PileUpMerger.cc (modified) (8 diffs)
• PileUpMerger.h (modified) (2 diffs)
• TrackPileUpSubtractor.cc (modified) (4 diffs)
• TrackPileUpSubtractor.h (modified) (1 diff)
• TreeWriter.cc (modified) (24 diffs)

modules/Calorimeter.cc

@@ -330 +330 @@
       fTrackHCalEnergy = 0.0;
 
+      fTowerECalTime = 0.0;
+      fTowerHCalTime = 0.0;
+
+      fTrackECalTime = 0.0;
+      fTrackHCalTime = 0.0;
+
+      fTowerECalWeightTime = 0.0; 
+      fTowerHCalWeightTime = 0.0;
+     
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -343 +352 @@
       track = static_cast<Candidate*>(fTrackInputArray->At(number));
       momentum = track->Momentum;
-
+      position = track->Position;
+
+      
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -349 +360 @@
       fTrackECalEnergy += ecalEnergy;
       fTrackHCalEnergy += hcalEnergy;
+      
+      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+       
+      fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
+      fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
 
       fTowerTrackArray->Add(track);
@@ -360 +377 @@
     particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     momentum = particle->Momentum;
+    position = particle->Position;
 
     // fill current tower
@@ -367 +385 @@
     fTowerECalEnergy += ecalEnergy;
     fTowerHCalEnergy += hcalEnergy;
+
+    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+
+    fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
+    fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
+    
 
     fTower->AddCandidate(particle);
@@ -383 +408 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalSigma, hcalSigma;
+  Double_t ecalTime, hcalTime, time;
 
   if(!fTower) return;
@@ -392 +418 @@
 
   ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
+  ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
 
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
@@ -399 +426 @@
 
   hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
+  hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
 
   energy = ecalEnergy + hcalEnergy;
+  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy)); 
 
 //  eta = fTowerEta;
@@ -410 +439 @@
   pt = energy / TMath::CosH(eta);
 
-  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
+ // fTower->Position.SetXYZT(-time, 0.0, 0.0, time);
+  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
@@ -420 +450 @@
   fTower->Edges[3] = fTowerEdges[3];
 
+
   // fill calorimeter towers and photon candidates
   if(energy > 0.0)

modules/Calorimeter.h

@@ -45 +45 @@
   Double_t fTowerECalEnergy, fTowerHCalEnergy;
   Double_t fTrackECalEnergy, fTrackHCalEnergy;
+  
+  Double_t fTowerECalTime, fTowerHCalTime;
+  Double_t fTrackECalTime, fTrackHCalTime;
+   
+  Double_t fTowerECalWeightTime, fTowerHCalWeightTime;
+  Double_t fTrackECalWeightTime, fTrackHCalWeightTime;
+  
   Int_t fTowerTrackHits, fTowerPhotonHits;
 

modules/FastJetFinder.cc

@@ -193 +193 @@
   TLorentzVector momentum;
   Double_t deta, dphi, detaMax, dphiMax;
+  Double_t time, weightTime, avTime;
   Int_t number;
   Double_t rho = 0;
@@ -259 +260 @@
     candidate = factory->NewCandidate();
 
+    time=0;
+    weightTime=0;
+
     inputList.clear();
     inputList = sequence->constituents(*itOutputList);
@@ -269 +273 @@
       if(deta > detaMax) detaMax = deta;
       if(dphi > dphiMax) dphiMax = dphi;
-
+      
+      time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
+      weightTime += TMath::Sqrt(constituent->Momentum.E());
+    
       candidate->AddCandidate(constituent);
     }
+   
+    avTime = time/weightTime;
 
     candidate->Momentum = momentum;
+    candidate->Position.SetT(avTime);
     candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
 

modules/ModulesLinkDef.h

@@ -19 +19 @@
 #include "modules/EnergySmearing.h"
 #include "modules/MomentumSmearing.h"
+#include "modules/TimeSmearing.h"
 #include "modules/Calorimeter.h"
 #include "modules/Isolation.h"
@@ -50 +51 @@
 #pragma link C++ class EnergySmearing+;
 #pragma link C++ class MomentumSmearing+;
+#pragma link C++ class TimeSmearing+;
 #pragma link C++ class Calorimeter+;
 #pragma link C++ class Isolation+;

modules/PileUpMerger.cc

@@ -16 +16 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesFormula.h"
+#include "classes/DelphesTF2.h"
 #include "classes/DelphesPileUpReader.h"
 
@@ -41 +41 @@
 
 PileUpMerger::PileUpMerger() :
-  fReader(0), fItInputArray(0)
-{
-}
+  fFunction(0), fReader(0), fItInputArray(0)
+{
+  fFunction = new DelphesTF2;
+}
+
 
 //------------------------------------------------------------------------------
@@ -49 +51 @@
 PileUpMerger::~PileUpMerger()
 {
+  delete fFunction;
 }
 
@@ -60 +63 @@
 
   fMeanPileUp  = GetDouble("MeanPileUp", 10);
-  fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
-
+  
+  fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
+  fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
+
+  // read vertex smearing formula
+        
+  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
+  fFunction->SetRange(-fZVertexSpread,-fTVertexSpread,fZVertexSpread,fTVertexSpread);
+ 
   fileName = GetString("PileUpFile", "MinBias.pileup");
   fReader = new DelphesPileUpReader(fileName);
@@ -90 +100 @@
   Float_t x, y, z, t;
   Float_t px, py, pz, e;
-  Double_t dz, dphi;
+  Double_t dz, dphi, dt;
   Int_t numberOfEvents, event;
   Long64_t allEntries, entry;
-  Candidate *candidate;
+  Candidate *candidate, *vertexcandidate;
   DelphesFactory *factory;
+  
+  const Double_t c_light = 2.99792458E8;
 
   fItInputArray->Reset();
+  
+    // --- Deal with Primary vertex first  ------
+  
+  fFunction->GetRandom2(dz,dt);
+ 
+  dt *= c_light*1.0E3; //necessary in order to make t in mm/c
+  dz *= 1.0E3; //necessary in order to make z in mm
+  
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
+    candidate->Position.SetXYZT(x, y, z+dz, t+dt);
     fParticleOutputArray->Add(candidate);
   }
 
   factory = GetFactory();
-
+  
+  vertexcandidate = factory->NewCandidate();
+  vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
+  fVertexOutputArray->Add(vertexcandidate);
+ 
+  
+  // --- Then with pile-up vertices  ------
+ 
   switch(fPileUpDistribution)
   {
@@ -118 +146 @@
 
   allEntries = fReader->GetEntries();
-
+ 
   for(event = 0; event < numberOfEvents; ++event)
   {
@@ -128 +156 @@
 
     fReader->ReadEntry(entry);
-
-    dz = gRandom->Gaus(0.0, fZVertexSpread);
+ 
+   // --- Pile-up vertex smearing 
+    
+    fFunction->GetRandom2(dz,dt);
+ 
+    dt *= c_light*1.0E3; //necessary in order to make t in mm/c
+    dz *= 1.0E3; //necessary in order to make z in mm
+ 
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    candidate = factory->NewCandidate();
-    candidate->Position.SetXYZT(0.0, 0.0, dz, 0.0);
-    fVertexOutputArray->Add(candidate);
+    vertexcandidate = factory->NewCandidate();
+    vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
+    vertexcandidate->IsPU = 1;
+    
+    fVertexOutputArray->Add(vertexcandidate);
 
     while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
@@ -149 +185 @@
 
       candidate->IsPU = 1;
-
+    
       candidate->Momentum.SetPxPyPzE(px, py, pz, e);
       candidate->Momentum.RotateZ(dphi);
 
-      candidate->Position.SetXYZT(x, y, z + dz, t);
+      candidate->Position.SetXYZT(x, y, z+dz, t+dt);
       candidate->Position.RotateZ(dphi);
-
+      
       fParticleOutputArray->Add(candidate);
     }

modules/PileUpMerger.h

@@ -19 +19 @@
 class TObjArray;
 class DelphesPileUpReader;
+class DelphesTF2;
 
 class PileUpMerger: public DelphesModule
@@ -35 +36 @@
   Int_t fPileUpDistribution;
   Double_t fMeanPileUp;
+  
   Double_t fZVertexSpread;
+  Double_t fTVertexSpread;
 
   DelphesPileUpReader *fReader;
-
+  
+  DelphesTF2 *fFunction; //!
+  
   TIterator *fItInputArray; //!
 

modules/TrackPileUpSubtractor.cc

@@ -53 +53 @@
 void TrackPileUpSubtractor::Init()
 {
+// import input array
+
+  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
+  fItVertexInputArray = fVertexInputArray->MakeIterator();
+  
   fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
 
@@ -85 +90 @@
     if(iterator) delete iterator;
   }
+
+  if(fItVertexInputArray) delete fItVertexInputArray;
 }
 
@@ -95 +102 @@
   TIterator *iterator;
   TObjArray *array;
-  Double_t z;
+  Double_t z, zvtx;
+
+  
+  // find z position of primary vertex
+  
+  fItVertexInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
+  {
+    if(candidate->IsPU == 0)
+    {
+    zvtx = candidate->Position.Z();
+    break;
+    }
+  }
+
 
   // loop over all input arrays
@@ -112 +133 @@
       // apply pile-up subtraction
       // assume perfect pile-up subtraction for tracks outside fZVertexResolution
-      if(candidate->IsPU && TMath::Abs(z) > fZVertexResolution) continue;
+      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
 
       array->Add(candidate);

modules/TrackPileUpSubtractor.h

@@ -39 +39 @@
 
   ClassDef(TrackPileUpSubtractor, 1)
+
+  TIterator *fItVertexInputArray; //!
+
+  const TObjArray *fVertexInputArray; //!
+
 };
 

modules/TreeWriter.cc

@@ -95 +95 @@
       continue;
     }
-
+    
     itClassMap = fClassMap.find(branchClass);
     if(itClassMap == fClassMap.end())
@@ -115 +115 @@
 void TreeWriter::Finish()
 {
-
 }
 
@@ -161 +160 @@
   GenParticle *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  
+  const Double_t c_light = 2.99792458E8;
+  
   // loop over all particles
   iterator.Reset();
@@ -208 +209 @@
     entry->Y = position.Y();
     entry->Z = position.Z();
-    entry->T = position.T();
+    entry->T = position.T()*1.0E-3/c_light;
   }
 }
@@ -219 +220 @@
   Candidate *candidate = 0;
   Vertex *entry = 0;
+  
+  const Double_t c_light = 2.99792458E8;
 
   // loop over all vertices
@@ -231 +234 @@
     entry->Y = position.Y();
     entry->Z = position.Z();
-    entry->T = position.T();
+    entry->T = position.T()*1.0E-3/c_light;
   }
 }
@@ -244 +247 @@
   Track *entry = 0;
   Double_t pt, signz, cosTheta, eta, rapidity;
-
+  const Double_t c_light = 2.99792458E8;
+  
   // loop over all tracks
   iterator.Reset();
@@ -271 +275 @@
     entry->YOuter = position.Y();
     entry->ZOuter = position.Z();
+    entry->TOuter = position.T()*1.0E-3/c_light;
 
     const TLorentzVector &momentum = candidate->Momentum;
@@ -290 +295 @@
     entry->Y = initialPosition.Y();
     entry->Z = initialPosition.Z();
+    entry->T = initialPosition.T()*1.0E-3/c_light;
 
     entry->Particle = particle;
@@ -303 +309 @@
   Tower *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  const Double_t c_light = 2.99792458E8;
+  
   // loop over all towers
   iterator.Reset();
@@ -309 +316 @@
   {
     const TLorentzVector &momentum = candidate->Momentum;
-
+    const TLorentzVector &position = candidate->Position;
+    
     pt = momentum.Pt();
     cosTheta = TMath::Abs(momentum.CosTheta());
@@ -331 +339 @@
     entry->Edges[2] = candidate->Edges[2];
     entry->Edges[3] = candidate->Edges[3];
-
+    
+    entry->T = position.T()*1.0E-3/c_light;
+    
     FillParticles(candidate, &entry->Particles);
   }
@@ -344 +354 @@
   Photon *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  const Double_t c_light = 2.99792458E8;
+  
   array->Sort();
 
@@ -353 +364 @@
     TIter it1(candidate->GetCandidates());
     const TLorentzVector &momentum = candidate->Momentum;
+    const TLorentzVector &position = candidate->Position;
+    
 
     pt = momentum.Pt();
@@ -366 +379 @@
     entry->PT = pt;
     entry->E = momentum.E();
-
+    
+    entry->T = position.T()*1.0E-3/c_light;
+   
     entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9;
 
@@ -381 +396 @@
   Electron *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  const Double_t c_light = 2.99792458E8;
+  
   array->Sort();
 
@@ -389 +405 @@
   {
     const TLorentzVector &momentum = candidate->Momentum;
-
+    const TLorentzVector &position = candidate->Position;
+    
     pt = momentum.Pt();
     cosTheta = TMath::Abs(momentum.CosTheta());
@@ -401 +418 @@
     entry->Phi = momentum.Phi();
     entry->PT = pt;
-
+    
+    entry->T = position.T()*1.0E-3/c_light;
+    
     entry->Charge = candidate->Charge;
 
@@ -418 +437 @@
   Muon *entry = 0;
   Double_t pt, signPz, cosTheta, eta, rapidity;
-
+  
+  const Double_t c_light = 2.99792458E8;
+  
   array->Sort();
 
@@ -426 +447 @@
   {
     const TLorentzVector &momentum = candidate->Momentum;
+    const TLorentzVector &position = candidate->Position;
+    
 
     pt = momentum.Pt();
@@ -442 +465 @@
     entry->PT = pt;
 
+    entry->T = position.T()*1.0E-3/c_light;
+    
+   // cout<<entry->PT<<","<<entry->T<<endl;
+
     entry->Charge = candidate->Charge;
 
@@ -457 +484 @@
   Double_t pt, signPz, cosTheta, eta, rapidity;
   Double_t ecalEnergy, hcalEnergy;
-
+  const Double_t c_light = 2.99792458E8;
+  
   array->Sort();
 
@@ -465 +493 @@
   {
     TIter itConstituents(candidate->GetCandidates());
-    const TLorentzVector &momentum = candidate->Momentum;
-
+    
+    const TLorentzVector &momentum = candidate->Momentum;
+    const TLorentzVector &position = candidate->Position;
+   
     pt = momentum.Pt();
     cosTheta = TMath::Abs(momentum.CosTheta());
@@ -479 +509 @@
     entry->PT = pt;
 
+    entry->T = position.T()*1.0E-3/c_light;
+    
     entry->Mass = momentum.M();