Changeset d612dec in git for modules

Timestamp:

Dec 9, 2021, 7:52:15 AM (3 years ago)

Author:

christinaw97 <christina.wang@…>

Children:

29b722a

Parents:

a5af1df (diff), 0c0c9af (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'master' of github.com:Christinaw97/delphes into HEAD

Location:

modules

Files:

• AngularSmearing.cc (modified) (2 diffs)
• Calorimeter.cc (modified) (7 diffs)
• Calorimeter.h (modified) (1 diff)
• ClusterCounting.cc (added)
• ClusterCounting.h (added)
• DecayFilter.cc (added)
• DecayFilter.h (added)
• Delphes.cc (modified) (2 diffs)
• DenseTrackFilter.cc (modified) (2 diffs)
• DualReadoutCalorimeter.cc (modified) (22 diffs)
• DualReadoutCalorimeter.h (modified) (2 diffs)
• Efficiency.cc (modified) (3 diffs)
• Efficiency.h (modified) (1 diff)
• EnergySmearing.cc (modified) (3 diffs)
• FastJetFinder.cc (modified) (7 diffs)
• FastJetFinder.h (modified) (1 diff)
• IdentificationMap.cc (modified) (1 diff)
• JetFakeParticle.cc (modified) (1 diff)
• ModulesLinkDef.h (modified) (4 diffs)
• MomentumSmearing.cc (modified) (4 diffs)
• MomentumSmearing.h (modified) (1 diff)
• ParticleDensity.cc (added)
• ParticleDensity.h (added)
• ParticlePropagator.cc (modified) (9 diffs)
• SimpleCalorimeter.cc (modified) (5 diffs)
• SimpleCalorimeter.h (modified) (1 diff)
• TauTagging.cc (modified) (5 diffs)
• TimeOfFlight.cc (added)
• TimeOfFlight.h (added)
• TimeSmearing.cc (modified) (6 diffs)
• TimeSmearing.h (modified) (3 diffs)
• TrackCovariance.cc (added)
• TrackCovariance.h (added)
• TrackSmearing.cc (modified) (8 diffs)
• TreeWriter.cc (modified) (12 diffs)
• TreeWriter.h (modified) (2 diffs)
• TruthVertexFinder.cc (added)
• TruthVertexFinder.h (added)

modules/AngularSmearing.cc

@@ -98 +98 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, eta, phi, e;
+  Double_t pt, eta, phi, e, m;
 
   fItInputArray->Reset();
   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
   {
-    const TLorentzVector &candidatePosition = candidate->Position;
     const TLorentzVector &candidateMomentum = candidate->Momentum;
-    eta = candidatePosition.Eta();
-    phi = candidatePosition.Phi();
+    eta = candidateMomentum.Eta();
+    phi = candidateMomentum.Phi();
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
+    m = candidateMomentum.M();
 
     // apply smearing formula for eta,phi
-
-    eta = gRandom->Gaus(eta, fFormulaEta->Eval(pt, eta, phi, e));
-    phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta, phi, e));
+    eta = gRandom->Gaus(eta, fFormulaEta->Eval(pt, eta, phi, e, candidate));
+    phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta, phi, e, candidate));
 
     if(pt <= 0.0) continue;
@@ -119 +118 @@
     mother = candidate;
     candidate = static_cast<Candidate *>(candidate->Clone());
-    eta = candidateMomentum.Eta();
-    phi = candidateMomentum.Phi();
-    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
+    candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
     candidate->AddCandidate(mother);
 

modules/Calorimeter.cc

@@ -231 +231 @@
   fItParticleInputArray->Reset();
   number = -1;
+  fTowerRmax=0.;
   while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
     ++number;
+
+    // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
+    if (particlePosition.Perp() > fTowerRmax)
+      fTowerRmax=particlePosition.Perp();
 
     pdgCode = TMath::Abs(particle->PID);
@@ -450 +455 @@
 
     fTower->AddCandidate(particle);
+    fTower->Position = position;
   }
 
@@ -461 +467 @@
 {
   Candidate *track, *tower, *mother;
-  Double_t energy, pt, eta, phi;
+  Double_t energy, pt, eta, phi, r;
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalNeutralEnergy, hcalNeutralEnergy;
@@ -511 +517 @@
   for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
   {
-    weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
+    weight = TMath::Power((fTower->ECalEnergyTimePairs[i].first),2);
     sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
     sumWeight += weight;
@@ -517 +523 @@
   }
 
+  // check whether barrel or endcap tower
+  if (fTower->Position.Perp() < fTowerRmax && TMath::Abs(eta) > 0.)
+    r = fTower->Position.Z()/TMath::SinH(eta);
+  else
+    r = fTower->Position.Pt();
+
   if(sumWeight > 0.0)
   {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime / sumWeight);
+    fTower->Position.SetPtEtaPhiE(r, eta, phi, sumWeightedTime / sumWeight);
   }
   else
   {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
+    fTower->Position.SetPtEtaPhiE(r, eta, phi, 999999.9);
   }
 
@@ -559 +571 @@
   if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
   {
-    // create new photon tower
+    // create new photon tower assuming null mass
     tower = static_cast<Candidate *>(fTower->Clone());
     pt = ecalNeutralEnergy / TMath::CosH(eta);
@@ -646 +658 @@
       track = static_cast<Candidate *>(track->Clone());
       track->AddCandidate(mother);
-
       track->Momentum *= rescaleFactor;
+      track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
 
       fEFlowTrackOutputArray->Add(track);

modules/Calorimeter.h

@@ -60 +60 @@
   Double_t fTimingEnergyMin;
   Bool_t fElectronsFromTrack;
+  Double_t fTowerRmax;
 
   Int_t fTowerTrackHits, fTowerPhotonHits;

modules/Delphes.cc

@@ -61 +61 @@
   fFactory(0)
 {
-  TFolder *folder = new TFolder(name, "");
+  TFolder *folder;
+
   fFactory = new DelphesFactory("ObjectFactory");
+
+  folder = new TFolder(name, "");
 
   SetName(name);
@@ -80 +83 @@
   if(folder)
   {
+    gROOT->GetListOfBrowsables()->Remove(folder);
     folder->Clear();
     delete folder;

modules/DenseTrackFilter.cc

@@ -237 +237 @@
 
   Candidate *candidate, *track;
-  Double_t pt, eta, phi;
+  Double_t pt, eta, phi, m;
   Int_t numberOfCandidates;
 
@@ -251 +251 @@
   eta = candidate->Momentum.Eta();
   phi = candidate->Momentum.Phi();
+  m = candidate->Momentum.M();
+
   eta = gRandom->Gaus(eta, fEtaPhiRes);
   phi = gRandom->Gaus(phi, fEtaPhiRes);
-  candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
+  candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
   candidate->AddCandidate(track);
 

modules/DualReadoutCalorimeter.cc

@@ -24 +24 @@
   // This implementation of dual calorimetry relies on several approximations:
   // - If hadronic energy is found in the tower the energy resolution then the full tower enrgy is smeared according to hadronic resolution (pessimistic for (e,n) or (pi+,gamma))
-  // - While e+ vs pi+ (or gamma vs n) separation is in principle possible for single particles (using C/S, PMT timing, lateral shower profile) it is not obvious it can be done overlapping particles. 
+  // - While e+ vs pi+ (or gamma vs n) separation is in principle possible for single particles (using C/S, PMT timing, lateral shower profile) it is not obvious it can be done overlapping particles.
   //   Now we assume that regarless of the number of particle hits per tower we can always distinguish e+ vs pi+, which is probably not true in the case (e+,n) vs (pi+,gamma) without longitudinal segmentation.
 
@@ -63 +63 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-  
+
   fECalResolutionFormula = new DelphesFormula;
   fHCalResolutionFormula = new DelphesFormula;
@@ -75 +75 @@
   fTowerTrackArray = new TObjArray;
   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
-  
+
 }
 
@@ -82 +82 @@
 DualReadoutCalorimeter::~DualReadoutCalorimeter()
 {
-  
+
   if(fECalResolutionFormula) delete fECalResolutionFormula;
   if(fHCalResolutionFormula) delete fHCalResolutionFormula;
@@ -94 +94 @@
   if(fTowerTrackArray) delete fTowerTrackArray;
   if(fItTowerTrackArray) delete fItTowerTrackArray;
- 
+
 }
 
@@ -172 +172 @@
   fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
   fHCalEnergyMin = GetDouble("HCalEnergyMin", 0.0);
+  fEnergyMin = GetDouble("EnergyMin", 0.0);
 
   fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
   fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
+  fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
 
   // switch on or off the dithering of the center of DualReadoutCalorimeter towers
@@ -245 +247 @@
   fItParticleInputArray->Reset();
   number = -1;
+  fTowerRmax=0.;
+
+  //cout<<"--------- new event ---------- "<<endl;
+
   while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
     ++number;
+
+    // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
+    if (particlePosition.Perp() > fTowerRmax)
+      fTowerRmax=particlePosition.Perp();
 
     pdgCode = TMath::Abs(particle->PID);
@@ -377 +387 @@
       fHCalTrackEnergy = 0.0;
       fTrackEnergy = 0.0;
-      
+
       fECalTrackSigma = 0.0;
       fHCalTrackSigma = 0.0;
       fTrackSigma = 0.0;
-      
+
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
+
+      fTowerTime = 0.0;
+      fTowerTimeWeight = 0.0;
 
       fECalTowerTrackArray->Clear();
       fHCalTowerTrackArray->Clear();
       fTowerTrackArray->Clear();
-    
+
     }
 
@@ -412 +425 @@
       }
 
-      
-      /*
-      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
-      {
-        fECalTrackEnergy += ecalEnergy;
-        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());        
-        if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;        
-        else energyGuess = momentum.E();
-
-        fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
-        fECalTowerTrackArray->Add(track);
-      }
-     
-      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
-      {
-        fHCalTrackEnergy += hcalEnergy;
-        hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-        if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
-        else energyGuess = momentum.E();
-
-        fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
-        fHCalTowerTrackArray->Add(track);
-      }
-      
-      // muons
-      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
-      {
-        fEFlowTrackOutputArray->Add(track);
-      }
-      */
-      
       // in Dual Readout we do not care if tracks are ECAL of HCAL
       if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9)
-      { 
+      {
         fTrackEnergy += energy;
-        // this sigma will be used to determine whether neutral excess is significant. We choose the resolution according to bthe higest deposited fraction (in practice had for charged hadrons and em for electrons)      
+        // this sigma will be used to determine whether neutral excess is significant. We choose the resolution according to bthe higest deposited fraction (in practice had for charged hadrons and em for electrons)
         sigma = 0.0;
         if(fHCalTrackFractions[number] > 0)
@@ -453 +435 @@
         else
           sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-          
-        if(sigma/momentum.E() < track->TrackResolution) 
-          energyGuess = ecalEnergy + hcalEnergy;     
+
+        if(sigma/momentum.E() < track->TrackResolution)
+          energyGuess = ecalEnergy + hcalEnergy;
         else
           energyGuess = momentum.E();
-             
+
         fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
         fTowerTrackArray->Add(track);
-      
       }
       else
@@ -478 +459 @@
     position = particle->Position;
 
+
     // fill current tower
     ecalEnergy = momentum.E() * fECalTowerFractions[number];
@@ -485 +467 @@
     fHCalTowerEnergy += hcalEnergy;
 
-    if(ecalEnergy > fTimingEnergyMin && fTower)
-    {
-      if (abs(particle->PID) != 11 || !fElectronsFromTrack)
-      {
-        fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
-      }
-    }
+    // assume combined timing measurements in ECAL/HCAL sections
+    fTowerTime += (ecalEnergy + hcalEnergy) * position.T(); //sigma_t ~ 1/sqrt(E)
+    fTowerTimeWeight += ecalEnergy + hcalEnergy;
 
     fTower->AddCandidate(particle);
+    fTower->Position = position;
   }
 
@@ -506 +485 @@
 
   Candidate *track, *tower, *mother;
-  Double_t energy, pt, eta, phi;
+  Double_t energy, pt, eta, phi, r, time;
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy;
-  
+
   Double_t ecalSigma, hcalSigma, sigma;
   Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma;
 
   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
-  
+
   TLorentzVector momentum;
   TFractionMap::iterator itFractionMap;
@@ -522 +501 @@
   if(!fTower) return;
 
-
-  //if (fHCalTowerEnergy < 30 && fECalTowerEnergy < 30) return;
-  //cout<<"----------- New tower ---------"<<endl;
-
-
-  // here we change behaviour w.r.t to standard calorimeter. Start with worse case scenario. If fHCalTowerEnergy > 0, assume total energy smeared by HCAL resolution.
-  // For example, if overlapping charged pions and photons take hadronic resolution as combined measurement
-
-  // if no hadronic fraction at all, then use ECAL resolution
-
-  //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<", Eta: "<<fTowerEta<<endl;
-
-  // if no hadronic energy, use ECAL resolution 
+  // if no hadronic energy, use ECAL resolution
   if (fHCalTowerEnergy <= fHCalEnergyMin)
   {
     energy = fECalTowerEnergy;
     sigma  = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
-    //cout<<"em energy"<<energy<<", sigma: "<<sigma<<endl;
-  }
-
-  // if hadronic fraction > 0, use HCAL resolution 
+  }
+
+  // if hadronic fraction > 0, use HCAL resolution
   else
   {
     energy = fECalTowerEnergy + fHCalTowerEnergy;
     sigma  = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
-    //cout<<"had energy: "<<energy<<", sigma: "<<sigma<<endl;
   }
 
   energy = LogNormal(energy, sigma);
-  //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
-  
+
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
 
@@ -562 +526 @@
   hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
 
+  time = (fTowerTimeWeight < 1.0E-09) ? 0.0 : fTowerTime / fTowerTimeWeight;
+
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
@@ -568 +534 @@
   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
 
-  //cout<<"Measured energy: "<<energy<<endl;
-
   if(fSmearTowerCenter)
   {
@@ -583 +547 @@
   pt = energy / TMath::CosH(eta);
 
-  // Time calculation for tower
-  fTower->NTimeHits = 0;
-  sumWeightedTime = 0.0;
-  sumWeight = 0.0;
-
-  for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
-  {
-    weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
-    sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
-    sumWeight += weight;
-    fTower->NTimeHits++;
-  }
-
-  if(sumWeight > 0.0)
-  {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
-  }
-  else
-  {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
-  }
+  // check whether barrel or endcap tower
+
+  // endcap
+  if (TMath::Abs(fTower->Position.Pt() - fTowerRmax) > 1.e-06 && TMath::Abs(eta) > 0.){
+    r = fTower->Position.Z()/TMath::SinH(eta);
+  }
+  // barrel
+  else {
+    r = fTower->Position.Pt();
+  }
+
+  fTower->Position.SetPtEtaPhiE(r, eta, phi, time);
+  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+  fTower->L = fTower->Position.Vect().Mag();
 
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
   fTower->Ehad = hcalEnergy;
-
+  fTower->Etrk = fTrackEnergy;
   fTower->Edges[0] = fTowerEdges[0];
   fTower->Edges[1] = fTowerEdges[1];
@@ -624 +581 @@
 
   // fill energy flow candidates
-  
+
   fTrackSigma = TMath::Sqrt(fTrackSigma);
   neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
   neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma + sigma*sigma);
 
-  //cout<<"trackEnergy: "<<fTrackEnergy<<", trackSigma: "<<fTrackSigma<<", Ntracks: "<<fTowerTrackArray->GetEntries()<<endl;
-
-  //cout<<"neutralEnergy: "<<neutralEnergy<<", neutralSigma: "<<neutralSigma<<endl;
-
-   // For now, if neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack !!! -> Creating only photons !! EFlowNeutralHadron collection will be empy!!! TO BE FIXED
   if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
   {
-    
-    //cout<<"significant neutral excess found:"<<endl;
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-    pt =  neutralEnergy / TMath::CosH(eta);
-    //cout<<"Creating tower with Pt, Eta, Phi, Energy: "<<pt<<","<<eta<<","<<phi<<","<<neutralEnergy<<endl;
+    pt = neutralEnergy / TMath::CosH(eta);
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
-    tower->Eem = neutralEnergy;
-    tower->Ehad = 0.0;
-    tower->PID = 22;
-
-    fEFlowPhotonOutputArray->Add(tower);
-
+
+    // if no hadronic energy, use ECAL resolution
+    if (fHCalTowerEnergy <= fHCalEnergyMin)
+    {
+      tower->Eem = neutralEnergy;
+      tower->Ehad = 0.0;
+      tower->PID = 22;
+      fEFlowPhotonOutputArray->Add(tower);
+    }
+    // if hadronic fraction > 0, use HCAL resolution
+    else
+    {
+      tower->Eem = 0;
+      tower->Ehad = neutralEnergy;
+      tower->PID = 130;
+      fEFlowNeutralHadronOutputArray->Add(tower);
+    }
 
     //clone tracks
@@ -654 +614 @@
     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     {
-      //cout<<"looping over tracks"<<endl;
       mother = track;
       track = static_cast<Candidate*>(track->Clone());
@@ -662 +621 @@
   }
 
-  
+
   // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
   else if(fTrackEnergy > 0.0)
@@ -670 +629 @@
     weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
 
-    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 
+    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
     rescaleFactor = bestEnergyEstimate/fTrackEnergy;
 
@@ -676 +635 @@
     fItTowerTrackArray->Reset();
     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-    {  
+    {
       mother = track;
-      track = static_cast<Candidate*>(track->Clone());
+      track = static_cast<Candidate *>(track->Clone());
       track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
+      track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
       fEFlowTrackOutputArray->Add(track);
     }
   }
-  
-
-  /*
-  // fill energy flow candidates
-  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
-  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
-
-  //compute neutral excesses
-  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
-  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
-  
-  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
-  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
-  
-   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
-  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
-  {
-    // create new photon tower
-    tower = static_cast<Candidate*>(fTower->Clone());
-    pt =  ecalNeutralEnergy / TMath::CosH(eta);
-    
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
-    tower->Eem = ecalNeutralEnergy;
-    tower->Ehad = 0.0;
-    tower->PID = 22;
-    
-    fEFlowPhotonOutputArray->Add(tower);
-   
-    //clone tracks
-    fItECalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
-    {
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  
-  }
- 
-  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
-  else if(fECalTrackEnergy > 0.0)
-  {
-    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
-    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
-  
-    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo); 
-    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
-
-    //rescale tracks
-    fItECalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
-    {  
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-
-
-  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
-  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
-  {
-    // create new photon tower
-    tower = static_cast<Candidate*>(fTower->Clone());
-    pt =  hcalNeutralEnergy / TMath::CosH(eta);
-    
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
-    tower->Ehad = hcalNeutralEnergy;
-    tower->Eem = 0.0;
-    
-    fEFlowNeutralHadronOutputArray->Add(tower);
-   
-    //clone tracks
-    fItHCalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
-    {
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  
-  }
- 
-  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
-  else if(fHCalTrackEnergy > 0.0)
-  {
-    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
-    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
-  
-    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo); 
-    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
-
-    //rescale tracks
-    fItHCalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
-    {  
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-
-  */
+
+
 }
 

modules/DualReadoutCalorimeter.h

@@ -60 +60 @@
   Double_t fECalTrackEnergy, fHCalTrackEnergy;
   Double_t fTrackEnergy;
+  Double_t fTowerRmax;
 
   Double_t fTimingEnergyMin;
@@ -77 +78 @@
   Double_t fHCalTrackSigma;
   Double_t fTrackSigma;
+
+  Double_t fTowerTime;
+  Double_t fTowerTimeWeight;
 
   Bool_t fSmearTowerCenter;

modules/Efficiency.cc

@@ -78 +78 @@
   fItInputArray = fInputArray->MakeIterator();
 
+  // switch to compute efficiency based on momentum vector eta, phi
+  fUseMomentumVector = GetBool("UseMomentumVector", false);
+
   // create output array
 
@@ -95 +98 @@
 {
   Candidate *candidate;
-  Double_t pt, eta, phi, e, d0, dz, ctgTheta;
+  Double_t pt, eta, phi, e;
+
 
   fItInputArray->Reset();
@@ -104 +108 @@
     eta = candidatePosition.Eta();
     phi = candidatePosition.Phi();
+
+    if (fUseMomentumVector){
+      eta = candidateMomentum.Eta();
+      phi = candidateMomentum.Phi();
+    }
+
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
-    d0 = candidate->D0;
-    dz = candidate->DZ;
-    ctgTheta = candidate->CtgTheta;
 
     // apply an efficency formula
-    if(gRandom->Uniform() > fFormula->Eval(pt, eta, phi, e, d0, dz, ctgTheta)) continue;
+    if(gRandom->Uniform() > fFormula->Eval(pt, eta, phi, e, candidate)) continue;
 
     fOutputArray->Add(candidate);

modules/Efficiency.h

@@ -53 +53 @@
   TObjArray *fOutputArray; //!
 
+  Double_t fUseMomentumVector; //!
+
   ClassDef(Efficiency, 1)
 };

modules/EnergySmearing.cc

@@ -95 +95 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, energy, eta, phi;
+  Double_t pt, energy, eta, phi, m;
 
   fItInputArray->Reset();
@@ -107 +107 @@
     phi = candidatePosition.Phi();
     energy = candidateMomentum.E();
+    m = candidateMomentum.M();
 
     // apply smearing formula
@@ -117 +118 @@
     eta = candidateMomentum.Eta();
     phi = candidateMomentum.Phi();
-    candidate->Momentum.SetPtEtaPhiE(energy / TMath::CosH(eta), eta, phi, energy);
+    pt = (energy > m) ? TMath::Sqrt(energy*energy - m*m)/TMath::CosH(eta) : 0;
+    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     candidate->TrackResolution = fFormula->Eval(pt, eta, phi, energy) / candidateMomentum.E();
     candidate->AddCandidate(mother);

modules/FastJetFinder.cc

@@ -104 +104 @@
   fJetAlgorithm = GetInt("JetAlgorithm", 6);
   fParameterR = GetDouble("ParameterR", 0.5);
+  fParameterP = GetDouble("ParameterP", -1.0);
 
   fConeRadius = GetDouble("ConeRadius", 0.5);
@@ -247 +248 @@
     fDefinition = new JetDefinition(fValenciaPlugin);
     break;
+  case 10:
+    fDefinition = new JetDefinition(ee_genkt_algorithm,fParameterR,fParameterP);
+    break;
+
   }
 
@@ -316 +321 @@
   Double_t deta, dphi, detaMax, dphiMax;
   Double_t time, timeWeight;
+  Double_t neutralEnergyFraction, chargedEnergyFraction;
+
   Int_t number, ncharged, nneutrals;
   Int_t charge;
@@ -418 +425 @@
     nneutrals = 0;
 
+    neutralEnergyFraction =0.;
+    chargedEnergyFraction =0.;
+
     inputList.clear();
     inputList = sequence->constituents(*itOutputList);
@@ -432 +442 @@
 
       if(constituent->Charge == 0)
+      {
         nneutrals++;
+        neutralEnergyFraction += constituent->Momentum.E();
+      }
       else
+      {
         ncharged++;
-
+        chargedEnergyFraction += constituent->Momentum.E();
+      }
+      
       time += TMath::Sqrt(constituent->Momentum.E()) * (constituent->Position.T());
       timeWeight += TMath::Sqrt(constituent->Momentum.E());
@@ -455 +471 @@
     candidate->NCharged = ncharged;
 
+    candidate->NeutralEnergyFraction = (momentum.E() > 0 ) ? neutralEnergyFraction/momentum.E() : 0.0;
+    candidate->ChargedEnergyFraction = (momentum.E() > 0 ) ? chargedEnergyFraction/momentum.E() : 0.0;
+
     //for exclusive clustering, access y_n,n+1 as exclusive_ymerge (fNJets);
     candidate->ExclYmerge23 = excl_ymerge23;
@@ -470 +489 @@
       fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, fRTrim), fastjet::SelectorPtFractionMin(fPtFracTrim));
       fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
-
-      trimmed_jet = join(trimmed_jet.constituents());
-
+      
       candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
 

modules/FastJetFinder.h

@@ -72 +72 @@
   Int_t fJetAlgorithm;
   Double_t fParameterR;
+  Double_t fParameterP;
 
   Double_t fJetPTMin;

modules/IdentificationMap.cc

@@ -170 +170 @@
       {
         // change PID of particle
+        candidate = static_cast<Candidate *>(candidate->Clone());
         if(pdgCodeOut != 0) candidate->PID = charge * pdgCodeOut;
         fOutputArray->Add(candidate);

modules/JetFakeParticle.cc

@@ -146 +146 @@
   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
   {
-    const TLorentzVector &candidatePosition = candidate->Position;
     const TLorentzVector &candidateMomentum = candidate->Momentum;
-    eta = candidatePosition.Eta();
-    phi = candidatePosition.Phi();
+    eta = candidateMomentum.Eta();
+    phi = candidateMomentum.Phi();
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();

modules/ModulesLinkDef.h

@@ -36 +36 @@
 #include "modules/MomentumSmearing.h"
 #include "modules/TrackSmearing.h"
+#include "modules/TrackCovariance.h"
+#include "modules/ClusterCounting.h"
 #include "modules/ImpactParameterSmearing.h"
 #include "modules/TimeSmearing.h"
+#include "modules/TimeOfFlight.h"
 #include "modules/SimpleCalorimeter.h"
 #include "modules/DenseTrackFilter.h"
@@ -72 +75 @@
 #include "modules/VertexFinder.h"
 #include "modules/VertexFinderDA4D.h"
+#include "modules/DecayFilter.h"
+#include "modules/ParticleDensity.h"
+#include "modules/TruthVertexFinder.h"
 #include "modules/ExampleModule.h"
 #include "modules/LLPFilter.h"
@@ -93 +99 @@
 #pragma link C++ class MomentumSmearing+;
 #pragma link C++ class TrackSmearing+;
+#pragma link C++ class TrackCovariance+;
+#pragma link C++ class ClusterCounting+;
 #pragma link C++ class ImpactParameterSmearing+;
 #pragma link C++ class TimeSmearing+;
+#pragma link C++ class TimeOfFlight+;
 #pragma link C++ class SimpleCalorimeter+;
 #pragma link C++ class DenseTrackFilter+;
@@ -129 +138 @@
 #pragma link C++ class VertexFinder+;
 #pragma link C++ class VertexFinderDA4D+;
+#pragma link C++ class DecayFilter+;
+#pragma link C++ class ParticleDensity+;
+#pragma link C++ class TruthVertexFinder+;
 #pragma link C++ class ExampleModule+;
 #pragma link C++ class LLPFilter+;

modules/MomentumSmearing.cc

@@ -78 +78 @@
   fItInputArray = fInputArray->MakeIterator();
 
+  // switch to compute momentum smearing based on momentum vector eta, phi
+  fUseMomentumVector = GetBool("UseMomentumVector", false);
+
   // create output array
 
@@ -95 +98 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, eta, phi, e, res;
+  Double_t pt, eta, phi, e, m, res;
 
   fItInputArray->Reset();
@@ -104 +107 @@
     eta = candidatePosition.Eta();
     phi = candidatePosition.Phi();
+
+    if (fUseMomentumVector){
+      eta = candidateMomentum.Eta();
+      phi = candidateMomentum.Phi();
+    }
+
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
-    res = fFormula->Eval(pt, eta, phi, e);
+    m = candidateMomentum.M();
+    res = fFormula->Eval(pt, eta, phi, e, candidate);
 
     // apply smearing formula
@@ -121 +131 @@
     eta = candidateMomentum.Eta();
     phi = candidateMomentum.Phi();
-    candidate->Momentum.SetPtEtaPhiE(pt, eta, phi, pt * TMath::CosH(eta));
+    candidate->Momentum.SetPtEtaPhiM(pt, eta, phi, m);
     //candidate->TrackResolution = fFormula->Eval(pt, eta, phi, e);
     candidate->TrackResolution = res;

modules/MomentumSmearing.h

@@ -55 +55 @@
   TObjArray *fOutputArray; //!
 
+  Double_t fUseMomentumVector; //!
+
   ClassDef(MomentumSmearing, 1)
 };

modules/ParticlePropagator.cc

@@ -125 +125 @@
   TLorentzVector particlePosition, particleMomentum, beamSpotPosition;
   Double_t px, py, pz, pt, pt2, e, q;
-  Double_t x, y, z, t, r, phi;
-  Double_t x_c, y_c, r_c, phi_c, phi_0;
-  Double_t x_t, y_t, z_t, r_t;
-  Double_t t1, t2, t3, t4, t5, t6;
-  Double_t t_z, t_r, t_ra, t_rb;
-  Double_t tmp, discr, discr2;
-  Double_t delta, gammam, omega, asinrho;
-  Double_t rcu, rc2, xd, yd, zd;
-  Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
+  Double_t x, y, z, t, r;
+  Double_t x_c, y_c, r_c, phi_0;
+  Double_t x_t, y_t, z_t, r_t, phi_t;
+  Double_t t_r, t_z;
+  Double_t tmp;
+  Double_t gammam, omega;
+  Double_t xd, yd, zd;
+  Double_t l, d0, dz, ctgTheta, alpha;
   Double_t bsx, bsy, bsz;
+  Double_t td, pio, phid, vz;
 
   const Double_t c_light = 2.99792458E8;
 
   if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
+  {
     beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  }
   else
   {
@@ -160 +162 @@
     particlePosition = particle->Position;
     particleMomentum = particle->Momentum;
+
     x = particlePosition.X() * 1.0E-3;
     y = particlePosition.Y() * 1.0E-3;
@@ -206 +209 @@
       // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
       tmp = px * y - py * x;
-      discr2 = pt2 * fRadius2 - tmp * tmp;
-
-      if(discr2 < 0.0)
-      {
-        // no solutions
-        continue;
-      }
-
-      tmp = px * x + py * y;
-      discr = TMath::Sqrt(discr2);
-      t1 = (-tmp + discr) / pt2;
-      t2 = (-tmp - discr) / pt2;
-      t = (t1 < 0.0) ? t2 : t1;
-
-      z_t = z + pz * t;
-      if(TMath::Abs(z_t) > fHalfLength)
-      {
-        t3 = (+fHalfLength - z) / pz;
-        t4 = (-fHalfLength - z) / pz;
-        t = (t3 < 0.0) ? t4 : t3;
-      }
+      t_r = (TMath::Sqrt(pt2 * fRadius2 - tmp * tmp) - px * x - py * y) / pt2;
+
+      t_z = (TMath::Sign(fHalfLength, pz) - z) / pz;
+
+      t = TMath::Min(t_r, t_z);
 
       x_t = x + px * t;
@@ -245 +232 @@
 
       fOutputArray->Add(candidate);
+
       if(TMath::Abs(q) > 1.0E-9)
       {
@@ -267 +255 @@
     {
 
-      // 1.  initial transverse momentum p_{T0}: Part->pt
-      //     initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)
-      //     relativistic gamma: gamma = E/mc^2; gammam = gamma * m
-      //     gyration frequency omega = q/(gamma m) fBz
-      //     helix radius r = p_{T0} / (omega gamma m)
+      // 1. initial transverse momentum p_{T0}: Part->pt
+      //    initial transverse momentum direction phi_0 = -atan(p_{X0} / p_{Y0})
+      //    relativistic gamma: gamma = E / mc^2; gammam = gamma * m
+      //    gyration frequency omega = q * Bz / (gammam)
+      //    helix radius r = p_{T0} / (omega * gammam)
 
       gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
-      omega = q * fBz / (gammam); // omega is here in [89875518/s]
+      omega = q * fBz / gammam; // omega is here in [89875518/s]
       r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
 
@@ -283 +271 @@
       y_c = y - r * TMath::Cos(phi_0);
       r_c = TMath::Hypot(x_c, y_c);
-      phi_c = TMath::ATan2(y_c, x_c);
-      phi = phi_c;
-      if(x_c < 0.0) phi += TMath::Pi();
-
-      rcu = TMath::Abs(r);
-      rc2 = r_c * r_c;
-
-      // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
-      xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c;
-      xd = (rc2 > 0.0) ? xd / rc2 : -999;
-      yd = y_c * (-rcu * r_c + rc2);
-      yd = (rc2 > 0.0) ? yd / rc2 : -999;
-      zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt;
-
-      // use perigee momentum rather than original particle
-      // momentum, since the orignal particle momentum isn't known
-
-      px = TMath::Sign(1.0, r) * pt * (-y_c / r_c);
-      py = TMath::Sign(1.0, r) * pt * (x_c / r_c);
-      etap = particleMomentum.Eta();
-      phip = TMath::ATan2(py, px);
-
-      particleMomentum.SetPtEtaPhiE(pt, etap, phip, particleMomentum.E());
+
+      // time of closest approach
+      td = (phi_0 + TMath::ATan2(x_c, y_c)) / omega;
+
+      // remove all the modulo pi that might have come from the atan
+      pio = TMath::Abs(TMath::Pi() / omega);
+      while(TMath::Abs(td) > 0.5 * pio)
+      {
+        td -= TMath::Sign(1.0, td) * pio;
+      }
+
+      vz = pz * c_light / e;
+
+      // calculate coordinates of closest approach to z axis
+      phid = phi_0 - omega * td;
+      xd = x_c - r * TMath::Sin(phid);
+      yd = y_c + r * TMath::Cos(phid);
+      zd = z + vz * td;
+
+      // momentum at closest approach
+      px = pt * TMath::Cos(phid);
+      py = pt * TMath::Sin(phid);
+
+      particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E());
 
       // calculate additional track parameters (correct for beamspot position)
-
-      d0 = ((x - bsx) * py - (y - bsy) * px) / pt;
-      dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
-      p = particleMomentum.P();
+      d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
+      dz = zd - bsz;
       ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
 
@@ -317 +304 @@
       //    t_r : time to exit from the sides
       //    t_z : time to exit from the front or the back
-      t_r = 0.0; // in [ns]
-      int sign_pz = (pz > 0.0) ? 1 : -1;
-      if(pz == 0.0)
-        t_z = 1.0E99;
-      else
-        t_z = gammam / (pz * 1.0E9 / c_light) * (-z + fHalfLength * sign_pz);
+      t_z = (vz == 0.0) ? 1.0E99 : (TMath::Sign(fHalfLength, pz) - z) / vz;
 
       if(r_c + TMath::Abs(r) < fRadius)
@@ -331 +313 @@
       else
       {
-        asinrho = TMath::ASin((fRadius * fRadius - r_c * r_c - r * r) / (2 * TMath::Abs(r) * r_c));
-        delta = phi_0 - phi;
-        if(delta < -TMath::Pi()) delta += 2 * TMath::Pi();
-        if(delta > TMath::Pi()) delta -= 2 * TMath::Pi();
-        t1 = (delta + asinrho) / omega;
-        t2 = (delta + TMath::Pi() - asinrho) / omega;
-        t3 = (delta + TMath::Pi() + asinrho) / omega;
-        t4 = (delta - asinrho) / omega;
-        t5 = (delta - TMath::Pi() - asinrho) / omega;
-        t6 = (delta - TMath::Pi() + asinrho) / omega;
-
-        if(t1 < 0.0) t1 = 1.0E99;
-        if(t2 < 0.0) t2 = 1.0E99;
-        if(t3 < 0.0) t3 = 1.0E99;
-        if(t4 < 0.0) t4 = 1.0E99;
-        if(t5 < 0.0) t5 = 1.0E99;
-        if(t6 < 0.0) t6 = 1.0E99;
-
-        t_ra = TMath::Min(t1, TMath::Min(t2, t3));
-        t_rb = TMath::Min(t4, TMath::Min(t5, t6));
-        t_r = TMath::Min(t_ra, t_rb);
+        alpha = TMath::ACos((r * r + r_c * r_c - fRadius * fRadius) / (2 * TMath::Abs(r) * r_c));
+        t_r = td + TMath::Abs(alpha / omega);
+
         t = TMath::Min(t_r, t_z);
       }
 
       // 4. position in terms of x(t), y(t), z(t)
-      x_t = x_c + r * TMath::Sin(omega * t - phi_0);
-      y_t = y_c + r * TMath::Cos(omega * t - phi_0);
-      z_t = z + pz * 1.0E9 / c_light / gammam * t;
+      phi_t = phi_0 - omega * t;
+      x_t = x_c - r * TMath::Sin(phi_t);
+      y_t = y_c + r * TMath::Cos(phi_t);
+      z_t = z + vz * t;
       r_t = TMath::Hypot(x_t, y_t);
 
-      // compute path length for an helix
-
-      alpha = pz * 1.0E9 / c_light / gammam;
-      l = t * TMath::Sqrt(alpha * alpha + r * r * omega * omega);
+      // lenght of the path from production to tracker
+      l = t * TMath::Hypot(vz, r * omega);
 
       if(r_t > 0.0)
       {
-
         // store these variables before cloning
         if(particle == candidate)
@@ -374 +336 @@
           particle->D0 = d0 * 1.0E3;
           particle->DZ = dz * 1.0E3;
-          particle->P = p;
+          particle->P = particleMomentum.P();
           particle->PT = pt;
           particle->CtgTheta = ctgTheta;
-          particle->Phi = phip;
+          particle->Phi = particleMomentum.Phi();
         }
 

modules/SimpleCalorimeter.cc

@@ -208 +208 @@
   fItParticleInputArray->Reset();
   number = -1;
+  fTowerRmax=0.;
   while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
     ++number;
+
+    // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
+    if (particlePosition.Perp() > fTowerRmax)
+      fTowerRmax=particlePosition.Perp();
 
     pdgCode = TMath::Abs(particle->PID);
@@ -394 +399 @@
     fTowerEnergy += energy;
 
-    fTowerTime += energy * position.T();
-    fTowerTimeWeight += energy;
+    fTowerTime += energy * energy * position.T(); //sigma_t ~ 1/E
+    fTowerTimeWeight += energy * energy;
 
     fTower->AddCandidate(particle);
+    fTower->Position = position;
+
   }
 
@@ -409 +416 @@
 {
   Candidate *tower, *track, *mother;
-  Double_t energy, neutralEnergy, pt, eta, phi;
+  Double_t energy, neutralEnergy, pt, eta, phi, r;
   Double_t sigma, neutralSigma;
   Double_t time;
@@ -443 +450 @@
   pt = energy / TMath::CosH(eta);
 
-  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
+  // endcap
+  if (TMath::Abs(fTower->Position.Pt() - fTowerRmax) > 1.e-06 && TMath::Abs(eta) > 0.){
+    r = fTower->Position.Z()/TMath::SinH(eta);
+  }
+  // barrel
+  else {
+    r = fTower->Position.Pt();
+  }
+
+  fTower->Position.SetPtEtaPhiE(r, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+  fTower->L = fTower->Position.Vect().Mag();
 
   fTower->Eem = (!fIsEcal) ? 0 : energy;
   fTower->Ehad = (fIsEcal) ? 0 : energy;
+  fTower->Etrk = fTrackEnergy;
 
   fTower->Edges[0] = fTowerEdges[0];
@@ -507 +525 @@
       track = static_cast<Candidate *>(track->Clone());
       track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
+      track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
       fEFlowTrackOutputArray->Add(track);
     }

modules/SimpleCalorimeter.h

@@ -62 +62 @@
   Double_t fTrackTime;
 
+  Double_t fTowerRmax;
+
   Double_t fTowerTimeWeight;
   Double_t fTrackTimeWeight;

modules/TauTagging.cc

@@ -83 +83 @@
     daughter1 = static_cast<Candidate *>(fParticleInputArray->At(i));
     pdgCode = TMath::Abs(daughter1->PID);
-    if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15)
-      return -1;
-    else if(pdgCode == 24)
+    //if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15)
+    //  return -1;
+    if(pdgCode == 24)
     {
       if(daughter1->D1 < 0) return -1;
@@ -96 +96 @@
     }
   }
-
   return 0;
 }
@@ -190 +189 @@
 void TauTagging::Process()
 {
-  Candidate *jet, *tau, *daughter;
+  Candidate *jet, *tau, *daughter, *part;
   TLorentzVector tauMomentum;
   Double_t pt, eta, phi, e, eff;
@@ -204 +203 @@
   // loop over all input jets
   fItJetInputArray->Reset();
+
   while((jet = static_cast<Candidate *>(fItJetInputArray->Next())))
   {
+
     const TLorentzVector &jetMomentum = jet->Momentum;
     pdgCode = 0;
@@ -243 +244 @@
       }
     }
+    
+    // fake electrons and muons
+    
+    if (pdgCode == 0)
+    {
+     
+      Double_t drMin = fDeltaR;   
+      fItPartonInputArray->Reset();
+      while((part = static_cast<Candidate *>(fItPartonInputArray->Next())))
+      {
+        if(TMath::Abs(part->PID) == 11 || TMath::Abs(part->PID) == 13) 
+        {
+            tauMomentum = part->Momentum;
+            if (tauMomentum.Pt() < fClassifier->fPTMin) continue;
+            if (TMath::Abs(tauMomentum.Eta()) > fClassifier->fEtaMax) continue;
+
+            Double_t dr = jetMomentum.DeltaR(tauMomentum);
+            if( dr < drMin)
+            {
+               drMin = dr;
+               pdgCode = TMath::Abs(part->PID);
+               charge = part->Charge;
+            }  
+        }
+      }
+    }
+
     // find an efficency formula
     itEfficiencyMap = fEfficiencyMap.find(pdgCode);

modules/TimeSmearing.cc

@@ -19 +19 @@
 /** \class TimeSmearing
  *
- *  Performs transverse momentum resolution smearing.
+ *  Performs time smearing.
  *
- *  \author P. Demin - UCL, Louvain-la-Neuve
+ *  \author M. Selvaggi - CERN
  *
  */
@@ -49 +49 @@
 
 using namespace std;
-
 //------------------------------------------------------------------------------
 
 TimeSmearing::TimeSmearing() :
-  fItInputArray(0)
+  fItInputArray(0), fResolutionFormula(0)
 {
+        fResolutionFormula = new DelphesFormula;
 }
 
@@ -61 +61 @@
 TimeSmearing::~TimeSmearing()
 {
+        if(fResolutionFormula) delete fResolutionFormula;
 }
 
@@ -69 +70 @@
   // read resolution formula
 
-  fTimeResolution = GetDouble("TimeResolution", 1.0E-10);
-  // import input array
+  // read time resolution formula in seconds
+  fResolutionFormula->Compile(GetString("TimeResolution", "30e-12"));
 
+  // import track input array
   fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
   fItInputArray = fInputArray->MakeIterator();
 
+
   // create output array
-
-  fOutputArray = ExportArray(GetString("OutputArray", "muons"));
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
 }
 
@@ -92 +94 @@
 {
   Candidate *candidate, *mother;
-  Double_t ti, tf_smeared, tf;
+  Double_t tf_smeared, tf;
+  Double_t eta, energy;
+  Double_t timeResolution;
+
   const Double_t c_light = 2.99792458E8;
 
@@ -98 +103 @@
   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
   {
-    const TLorentzVector &candidateInitialPosition = candidate->InitialPosition;
+
     const TLorentzVector &candidateFinalPosition = candidate->Position;
+    const TLorentzVector &candidateMomentum = candidate->Momentum;
 
-    ti = candidateInitialPosition.T() * 1.0E-3 / c_light;
     tf = candidateFinalPosition.T() * 1.0E-3 / c_light;
 
+    eta = candidateMomentum.Eta();
+    energy = candidateMomentum.E();
+
     // apply smearing formula
-    tf_smeared = gRandom->Gaus(tf, fTimeResolution);
-    ti = ti + tf_smeared - tf;
-    tf = tf_smeared;
+    timeResolution = fResolutionFormula->Eval(0.0, eta, 0.0, energy);
+    tf_smeared = gRandom->Gaus(tf, timeResolution);
 
     mother = candidate;
     candidate = static_cast<Candidate *>(candidate->Clone());
-    candidate->InitialPosition.SetT(ti * 1.0E3 * c_light);
-    candidate->Position.SetT(tf * 1.0E3 * c_light);
 
-    candidate->ErrorT = fTimeResolution * 1.0E3 * c_light;
+    candidate->Position.SetT(tf_smeared * 1.0E3 * c_light);
+    candidate->ErrorT = timeResolution * 1.0E3 * c_light;
 
     candidate->AddCandidate(mother);
-
     fOutputArray->Add(candidate);
   }
 }
-
-//------------------------------------------------------------------------------

modules/TimeSmearing.h

@@ -22 +22 @@
 /** \class TimeSmearing
  *
- *  Performs transverse time smearing.
+ *  Performs time smearing.
  *
- *  \author Michele Selvaggi - UCL, Louvain-la-Neuve
+ *  \author Michele Selvaggi - CERN
  *
  */
@@ -32 +32 @@
 class TIterator;
 class TObjArray;
+class DelphesFormula;
 
 class TimeSmearing: public DelphesModule
@@ -44 +45 @@
 
 private:
-  Double_t fTimeResolution;
+
+  DelphesFormula *fResolutionFormula;
+  Int_t fVertexTimeMode;
 
   TIterator *fItInputArray; //!

modules/TrackSmearing.cc

@@ -158 +158 @@
   TLorentzVector beamSpotPosition;
   Candidate *candidate, *mother;
-  Double_t pt, eta, d0, d0Error, trueD0, dz, dzError, trueDZ, p, pError, trueP, ctgTheta, ctgThetaError, trueCtgTheta, phi, phiError, truePhi;
+  Double_t pt, eta, e, m, d0, d0Error, trueD0, dz, dzError, trueDZ, p, pError, trueP, ctgTheta, ctgThetaError, trueCtgTheta, phi, phiError, truePhi;
   Double_t x, y, z, t, px, py, pz, theta;
   Double_t q, r;
@@ -223 +223 @@
     pt = momentum.Pt();
     eta = momentum.Eta();
+    e = momentum.E();
+    m = momentum.M();
 
     d0 = trueD0 = candidate->D0;
@@ -232 +234 @@
 
     if(fUseD0Formula)
-      d0Error = fD0Formula->Eval(pt, eta);
+      d0Error = fD0Formula->Eval(pt, eta, phi, e, candidate);
     else
     {
@@ -247 +249 @@
 
     if(fUseDZFormula)
-      dzError = fDZFormula->Eval(pt, eta);
+      dzError = fDZFormula->Eval(pt, eta, phi, e, candidate);
     else
     {
@@ -262 +264 @@
 
     if(fUsePFormula)
-      pError = fPFormula->Eval(pt, eta) * p;
+      pError = fPFormula->Eval(pt, eta, phi, e, candidate) * p;
     else
     {
@@ -277 +279 @@
 
     if(fUseCtgThetaFormula)
-      ctgThetaError = fCtgThetaFormula->Eval(pt, eta);
+      ctgThetaError = fCtgThetaFormula->Eval(pt, eta, phi, e, candidate);
     else
     {
@@ -292 +294 @@
 
     if(fUsePhiFormula)
-      phiError = fPhiFormula->Eval(pt, eta);
+      phiError = fPhiFormula->Eval(pt, eta, phi, e, candidate);
     else
     {
@@ -331 +333 @@
     candidate->Momentum.SetPy(p * TMath::Sin(phi) * TMath::Sin(theta));
     candidate->Momentum.SetPz(p * TMath::Cos(theta));
-    candidate->Momentum.SetE(candidate->Momentum.Pt() * TMath::CosH(eta));
+    candidate->Momentum.SetE(TMath::Sqrt(p*p + m*m));
     candidate->PT = candidate->Momentum.Pt();
 

modules/TreeWriter.cc

@@ -72 +72 @@
   fClassMap[Track::Class()] = &TreeWriter::ProcessTracks;
   fClassMap[Tower::Class()] = &TreeWriter::ProcessTowers;
+  fClassMap[ParticleFlowCandidate::Class()] = &TreeWriter::ProcessParticleFlowCandidates;
   fClassMap[Photon::Class()] = &TreeWriter::ProcessPhotons;
   fClassMap[Electron::Class()] = &TreeWriter::ProcessElectrons;
@@ -123 +124 @@
     fBranchMap.insert(make_pair(branch, make_pair(itClassMap->second, array)));
   }
+
+  param = GetParam("Info");
+  TString infoName;
+  Double_t infoValue;
+
+  size = param.GetSize();
+  for(i = 0; i < size / 2; ++i)
+  {
+    infoName = param[i * 2].GetString();
+    infoValue = param[i * 2 + 1].GetDouble();
+
+    AddInfo(infoName, infoValue);
+  }
 }
 
@@ -138 +152 @@
   it1.Reset();
   array->Clear();
+
   while((candidate = static_cast<Candidate *>(it1.Next())))
   {
@@ -215 +230 @@
     entry->Pz = momentum.Pz();
 
-    entry->D0 = candidate->D0;
-    entry->DZ = candidate->DZ;
-    entry->P = momentum.P();
-    entry->PT = candidate->PT;
-    entry->CtgTheta = candidate->CtgTheta;
-    entry->Phi = candidate->Phi;
-
     entry->Eta = eta;
     entry->Phi = momentum.Phi();
@@ -322 +330 @@
   Candidate *particle = 0;
   Track *entry = 0;
-  Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
+  Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m;
   const Double_t c_light = 2.99792458E8;
 
@@ -356 +364 @@
 
     entry->D0 = candidate->D0;
-    entry->ErrorD0 = candidate->ErrorD0;
     entry->DZ = candidate->DZ;
-    entry->ErrorDZ = candidate->ErrorDZ;
+    entry->Nclusters = candidate->Nclusters;
+    entry->dNdx = candidate->dNdx;
 
     entry->ErrorP = candidate->ErrorP;
     entry->ErrorPT = candidate->ErrorPT;
+
+    // diagonal covariance matrix terms
+    entry->ErrorD0 = candidate->ErrorD0;
+    entry->ErrorC = candidate->ErrorC;
+    entry->ErrorPhi = candidate->ErrorPhi;
+    entry->ErrorDZ = candidate->ErrorDZ;
     entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
-    entry->ErrorPhi = candidate->ErrorPhi;
+
+    // add some offdiagonal covariance matrix elements
+    entry->ErrorD0Phi          = candidate->TrackCovariance(0,1);
+    entry->ErrorD0C            = candidate->TrackCovariance(0,2);
+    entry->ErrorD0DZ           = candidate->TrackCovariance(0,3);
+    entry->ErrorD0CtgTheta     = candidate->TrackCovariance(0,4);
+    entry->ErrorPhiC           = candidate->TrackCovariance(1,2);
+    entry->ErrorPhiDZ          = candidate->TrackCovariance(1,3);
+    entry->ErrorPhiCtgTheta    = candidate->TrackCovariance(1,4);
+    entry->ErrorCDZ            = candidate->TrackCovariance(2,3);
+    entry->ErrorCCtgTheta      = candidate->TrackCovariance(2,4);
+    entry->ErrorDZCtgTheta     = candidate->TrackCovariance(3,4);
 
     entry->Xd = candidate->Xd;
@@ -374 +399 @@
     p = momentum.P();
     phi = momentum.Phi();
+    m = momentum.M();
     ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
 
@@ -386 +412 @@
     entry->Phi = phi;
     entry->CtgTheta = ctgTheta;
+    entry->C = candidate->C;
+    entry->Mass = m;
 
     particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
-    const TLorentzVector &initialPosition = particle->Position;
+    //const TLorentzVector &initialPosition = particle->Position;
+    const TLorentzVector &initialPosition = candidate->InitialPosition;
 
     entry->X = initialPosition.X();
@@ -394 +423 @@
     entry->Z = initialPosition.Z();
     entry->T = initialPosition.T() * 1.0E-3 / c_light;
+    entry->ErrorT =candidate-> ErrorT * 1.0E-3 / c_light;
 
     entry->Particle = particle;
@@ -435 +465 @@
     entry->Eem = candidate->Eem;
     entry->Ehad = candidate->Ehad;
+    entry->Etrk = candidate->Etrk;
     entry->Edges[0] = candidate->Edges[0];
     entry->Edges[1] = candidate->Edges[1];
@@ -444 +475 @@
 
     FillParticles(candidate, &entry->Particles);
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void TreeWriter::ProcessParticleFlowCandidates(ExRootTreeBranch *branch, TObjArray *array)
+{
+
+  TIter iterator(array);
+  Candidate *candidate = 0;
+  Candidate *particle = 0;
+  ParticleFlowCandidate *entry = 0;
+  Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m;
+  const Double_t c_light = 2.99792458E8;
+
+  // loop over all tracks
+  iterator.Reset();
+  while((candidate = static_cast<Candidate *>(iterator.Next())))
+  {
+    const TLorentzVector &position = candidate->Position;
+
+    cosTheta = TMath::Abs(position.CosTheta());
+    signz = (position.Pz() >= 0.0) ? 1.0 : -1.0;
+    eta = (cosTheta == 1.0 ? signz * 999.9 : position.Eta());
+    rapidity = (cosTheta == 1.0 ? signz * 999.9 : position.Rapidity());
+
+    entry = static_cast<ParticleFlowCandidate *>(branch->NewEntry());
+
+    entry->SetBit(kIsReferenced);
+    entry->SetUniqueID(candidate->GetUniqueID());
+
+    entry->PID = candidate->PID;
+
+    entry->Charge = candidate->Charge;
+
+    entry->EtaOuter = eta;
+    entry->PhiOuter = position.Phi();
+
+    entry->XOuter = position.X();
+    entry->YOuter = position.Y();
+    entry->ZOuter = position.Z();
+    entry->TOuter = position.T() * 1.0E-3 / c_light;
+
+    entry->L = candidate->L;
+
+    entry->D0 = candidate->D0;
+    entry->DZ = candidate->DZ;
+    entry->Nclusters = candidate->Nclusters;
+    entry->dNdx = candidate->dNdx;
+
+    entry->ErrorP = candidate->ErrorP;
+    entry->ErrorPT = candidate->ErrorPT;
+    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
+
+
+    // diagonal covariance matrix terms
+
+    entry->ErrorD0 = candidate->ErrorD0;
+    entry->ErrorC = candidate->ErrorC;
+    entry->ErrorPhi = candidate->ErrorPhi;
+    entry->ErrorDZ = candidate->ErrorDZ;
+    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
+
+    // add some offdiagonal covariance matrix elements
+    entry->ErrorD0Phi          = candidate->TrackCovariance(0,1);
+    entry->ErrorD0C            = candidate->TrackCovariance(0,2);
+    entry->ErrorD0DZ           = candidate->TrackCovariance(0,3);
+    entry->ErrorD0CtgTheta     = candidate->TrackCovariance(0,4);
+    entry->ErrorPhiC           = candidate->TrackCovariance(1,2);
+    entry->ErrorPhiDZ          = candidate->TrackCovariance(1,3);
+    entry->ErrorPhiCtgTheta    = candidate->TrackCovariance(1,4);
+    entry->ErrorCDZ            = candidate->TrackCovariance(2,3);
+    entry->ErrorCCtgTheta      = candidate->TrackCovariance(2,4);
+    entry->ErrorDZCtgTheta     = candidate->TrackCovariance(3,4);
+
+    entry->Xd = candidate->Xd;
+    entry->Yd = candidate->Yd;
+    entry->Zd = candidate->Zd;
+
+    const TLorentzVector &momentum = candidate->Momentum;
+
+    e = momentum.E();
+    pt = momentum.Pt();
+    p = momentum.P();
+    phi = momentum.Phi();
+    m = momentum.M();
+    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
+
+    entry->E = e;
+    entry->P = p;
+    entry->PT = pt;
+    entry->Eta = eta;
+    entry->Phi = phi;
+    entry->CtgTheta = ctgTheta;
+    entry->C = candidate->C;
+    entry->Mass = m;
+
+    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
+    //const TLorentzVector &initialPosition = particle->Position;
+    const TLorentzVector &initialPosition = candidate->InitialPosition;
+
+    entry->X = initialPosition.X();
+    entry->Y = initialPosition.Y();
+    entry->Z = initialPosition.Z();
+    entry->T = initialPosition.T() * 1.0E-3 / c_light;
+    entry->ErrorT = candidate-> ErrorT * 1.0E-3 / c_light;
+
+    entry->VertexIndex = candidate->ClusterIndex;
+
+    entry->Eem = candidate->Eem;
+    entry->Ehad = candidate->Ehad;
+    entry->Etrk = candidate->Etrk;
+    entry->Edges[0] = candidate->Edges[0];
+    entry->Edges[1] = candidate->Edges[1];
+    entry->Edges[2] = candidate->Edges[2];
+    entry->Edges[3] = candidate->Edges[3];
+
+    //entry->T = position.T() * 1.0E-3 / c_light;
+    entry->NTimeHits = candidate->NTimeHits;
+
+    FillParticles(candidate, &entry->Particles);
+
   }
 }
@@ -687 +840 @@
     entry->NCharged = candidate->NCharged;
     entry->NNeutrals = candidate->NNeutrals;
+
+    entry->NeutralEnergyFraction = candidate->NeutralEnergyFraction;
+    entry->ChargedEnergyFraction = candidate->ChargedEnergyFraction;
     entry->Beta = candidate->Beta;
     entry->BetaStar = candidate->BetaStar;

modules/TreeWriter.h

@@ -56 +56 @@
   void ProcessTracks(ExRootTreeBranch *branch, TObjArray *array);
   void ProcessTowers(ExRootTreeBranch *branch, TObjArray *array);
+  void ProcessParticleFlowCandidates(ExRootTreeBranch *branch, TObjArray *array);
   void ProcessPhotons(ExRootTreeBranch *branch, TObjArray *array);
   void ProcessElectrons(ExRootTreeBranch *branch, TObjArray *array);
@@ -79 +80 @@
 #endif
 
-  ClassDef(TreeWriter, 1)
+  ClassDef(TreeWriter, 2)
 };