Changeset 79a7b3e in git

Timestamp:

Jan 24, 2020, 3:55:17 PM (5 years ago)

Author:

GitHub <noreply@…>

Branches:

Timing

Children:

62764fb

Parents:

364dbe1 (diff), 4ac0049 (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.

git-author:

Michele Selvaggi <michele.selvaggi@…> (01/24/20 15:55:17)

git-committer:

GitHub <noreply@…> (01/24/20 15:55:17)

Message:

Merge pull request #5 from kaanyuxel/master

New Time Smearing Module for Neutral Particles and FCC-hh Card

Files:

• cards/FCC/FCChh_PileUpVtx.tcl (modified) (8 diffs)
• cards/FCC/FCChh_PileUpVtx_new.tcl (added)
• cards/FCC/trackResolutionFCChh.tcl (added)
• modules/ModulesLinkDef.h (modified) (2 diffs)
• modules/PileUpSubtractor4D.cc (added)
• modules/PileUpSubtractor4D.h (added)
• modules/TimeSmearing.cc (modified) (6 diffs)
• modules/TimeSmearing.h (modified) (1 diff)
• modules/TimeSmearingNeutral.cc (added)
• modules/TimeSmearingNeutral.h (added)
• modules/TrackPileUpSubtractor.cc (modified) (4 diffs)
• modules/TrackTimingPileUpSubtractor.cc (modified) (4 diffs)
• modules/TrackTimingPileUpSubtractor.h (modified) (1 diff)
• modules/VertexFinderDA4D.cc (modified) (12 diffs)
• modules/VertexFinderDA4D.h (modified) (1 diff)

cards/FCC/FCChh_PileUpVtx.tcl

@@ -17 +17 @@
 set ExecutionPath {
 
-  BeamSpotFilter
   PileUpMerger
   ParticlePropagator
@@ -33 +32 @@
   TrackMerger
 
-
   TrackSmearing
-  TimeSmearing  
-
-  VertexFinderDA4D  
-
-  TrackTimingPileUpSubtractor  
 
   ECal
@@ -47 +40 @@
   EFlowMerger
   EFlowFilter
+
+  TimeSmearingMIP
+  TimeSmearingPhotons
+  TimeSmearingNH    
+
+  VertexFinderDA4D
+  PileUpSubtractor4D
+
+  HighMassVertexRecover    
 
   PhotonEfficiency
@@ -81 +83 @@
 
   TreeWriter
-}
-
-#######################
-# GenBeamSpotFilter
-# Saves a particle intended to represent the beamspot
-#######################
-
-module BeamSpotFilter BeamSpotFilter {
-    set InputArray Delphes/stableParticles
-    set OutputArray beamSpotParticle
-
 }
 
@@ -308 +299 @@
 
   # from http://mersi.web.cern.ch/mersi/layouts/.private/Baseline_tilted_200_Pixel_1_1_1/index.html
-  source trackResolutionCMS.tcl 
+  source trackResolutionFCChh.tcl 
   # FIXME !!!! we need to add track resolution of FCC-hh baseline detector !!!!!
-}
-
-########################################
-#   Time Smearing
-########################################
-
-module TimeSmearing TimeSmearing {
-  set InputArray TrackSmearing/tracks
-  set OutputArray tracks
-
-  # assume 20 ps resolution for now
-  set TimeResolution 20E-12
-}
-
-##################################
-# Primary vertex reconstruction
-##################################
-
-
-module VertexFinderDA4D VertexFinderDA4D {
-  set InputArray TimeSmearing/tracks
-
-  set OutputArray tracks
-  set VertexOutputArray vertices
-
-  set Verbose 0
-  set MinPT 1.0
-
-  # in mm
-  set VertexSpaceSize 0.5
-
-  # in s
-  set VertexTimeSize 10E-12
-
-  set UseTc 1
-  set BetaMax 0.1
-  set BetaStop 1.0
-  set CoolingFactor 0.8
-  set MaxIterations 100
-
-  # in mm
-  set DzCutOff 40
-  set D0CutOff 30
-
-}
-
-##########################
-# Track pile-up subtractor
-##########################
-
-module TrackTimingPileUpSubtractor TrackTimingPileUpSubtractor {
-# add InputArray InputArray OutputArray
-
-  add InputArray ChargedHadronMomentumSmearing/chargedHadrons
-  add InputArray ElectronMomentumSmearing/electrons
-  add InputArray MuonMomentumSmearing/muons
-  
-  set VertexInputArray VertexFinderDA4D/vertices
-  # assume perfect pile-up subtraction for tracks with |z| > fZVertexResolution
-  # Z vertex resolution in m
-  set ZVertexResolution {0.0001}
 }
 
@@ -539 +469 @@
 }
 
-
 #################
 # Electron filter
@@ -606 +535 @@
 }
 
+########################################
+#   Time Smearing Neutral MIP
+########################################
+
+module TimeSmearing TimeSmearingMIP {
+  set InputArray HCal/eflowTracks
+  set OutputArray tracks
+
+  # assume 30 ps resolution for now
+  set TimeResolution {30E-12}
+}
+
+########################################
+#   Time Smearing Neutral Photons
+########################################
+
+module TimeSmearing TimeSmearingPhotons {
+  set InputArray ECal/eflowPhotons
+  set OutputArray photons
+  set TimeResolution {sqrt(20^2 + 150^2)/energy^2}
+}
+
+########################################
+#   Time Smearing Neutral NeutralHadrons
+########################################
+#
+module TimeSmearing TimeSmearingNH {
+  set InputArray HCal/eflowNeutralHadrons
+  set OutputArray neutralhadrons
+
+  # assume 30 ps resolution for now
+  set TimeResolution {sqrt(20^2 + 150^2)/energy^2}
+}
+
+
+##################################
+# Primary vertex reconstruction
+##################################
+
+
+module VertexFinderDA4D VertexFinderDA4D {
+  set InputArray TimeSmearingMIP/tracks
+
+  set OutputArray tracks
+  set VertexOutputArray vertices
+
+  set Verbose 0
+  set MinPT 1.0
+
+  # in mm
+  set VertexSpaceSize 0.5
+
+  # in s
+  set VertexTimeSize 10E-12
+
+  set UseTc 1
+  set BetaMax 0.1
+  set BetaStop 1.0
+  set CoolingFactor 0.8
+  set MaxIterations 100
+
+  # in mm
+  set DzCutOff 40
+  set D0CutOff 30
+
+}
+
+##########################
+# Track pile-up subtractor
+##########################
+
+module PileUpSubtractor4D PileUpSubtractor4D {
+# add InputArray InputArray OutputArray
+
+  add InputArray TimeSmearingMIP/tracks
+  add InputArray TimeSmearingPhotons/photons
+  add InputArray TimeSmearingNH/neutralhadrons
+
+  set VertexInputArray VertexFinderDA4D/vertices
+
+  set fChargedMinSignificance {3}
+  set fNeutralMinSignificance {3}
+}
+
+######################################
+# Heavy(slow) particles vertex recover
+######################################
+
+module HighMassVertexRecover HighMassVertexRecover {
+
+  set TrackInputArray VertexFinderDA4D/tracks
+  set VertexInputArray VertexFinderDA4D/vertices
+
+  set TrackOutputArray tracks
+  set VertexOutputArray vertices
+
+  set Verbose 0
+
+}
 
 ###################
@@ -1057 +1085 @@
   add Branch ScalarHT/energy ScalarHT ScalarHT
   add Branch VertexFinderDA4D/vertices Vertex4D Vertex
-}
-
+
+  add Branch HighMassVertexRecover/tracks Track Track
+}
+

modules/ModulesLinkDef.h

@@ -56 +56 @@
 #include "modules/PileUpMerger.h"
 #include "modules/JetPileUpSubtractor.h"
-#include "modules/TrackPileUpSubtractor.h"
 #include "modules/TrackTimingPileUpSubtractor.h"
 #include "modules/TaggingParticlesSkimmer.h"
@@ -114 +113 @@
 #pragma link C++ class PileUpMerger+;
 #pragma link C++ class JetPileUpSubtractor+;
-#pragma link C++ class TrackPileUpSubtractor+;
 #pragma link C++ class TrackTimingPileUpSubtractor+;
 #pragma link C++ class TaggingParticlesSkimmer+;

modules/TimeSmearing.cc

@@ -53 +53 @@
 
 TimeSmearing::TimeSmearing() :
-  fFormula(0), fItInputArray(0)
+  fItInputArray(0), fFormula(0)
 {
-  fFormula = new DelphesFormula;
+        fFormula = new DelphesFormula;
 }
 
@@ -61 +61 @@
 
 TimeSmearing::~TimeSmearing()
-{
-  if(fFormula) delete fFormula;
+{  
+        if(fFormula) delete fFormula;
 }
 
@@ -69 +69 @@
 void TimeSmearing::Init()
 {
-  // read resolution formula
-
-  fFormula->Compile(GetString("TimeResolution", "1.0"));
+  // read resolution formula    
+  fFormula->Compile(GetString("TimeResolution", "0.001"));
 
   // import input array
@@ -77 +76 @@
   fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
   fItInputArray = fInputArray->MakeIterator();
-
   // create output array
 
@@ -95 +93 @@
 {
   Candidate *candidate, *mother;
-  Double_t ti, tf_smeared, tf, timeResolution;
-  Double_t pt, eta, phi, e, d0, dz, ctgTheta;
+  Double_t ti, tf_smeared, tf;
+  Double_t pt, eta, phi, e, p, l;
+  Double_t beta_particle;
+  Double_t timeResolution;
 
 
@@ -114 +114 @@
     pt = candidateMomentum.Pt();
     e = candidateMomentum.E();
-    d0 = candidate->D0;
-    dz = candidate->DZ;
-    ctgTheta = candidate->CtgTheta;
-
-    // apply smearing formula
+    p = candidateMomentum.P();
+    beta_particle = p/e;
+    l = candidate->L;
+    timeResolution = fFormula->Eval(e);
     
-    timeResolution = fFormula->Eval(pt, eta, phi, e, d0, dz, ctgTheta);
-    if(fabs(candidate->Position.Eta())<fEtaMax)
+    if(candidate->Charge != 0)
     { 
       tf_smeared = tf + timeResolution*gRandom->Gaus(0, 1);
+      mother = candidate;
+      candidate = static_cast<Candidate*>(candidate->Clone());  // I am not sure that we need these lines !!!
+      candidate->AddCandidate(mother);
+      candidate->InitialPosition.SetT((100+ti)*1.0E3*c_light);  // Dummy Value, correct value will be computed by VertexFinderDA4D
+      candidate->Position.SetT(tf_smeared*1.0E3*c_light);
+      candidate->ErrorT = timeResolution*1.0E3*c_light;
+      fOutputArray->Add(candidate);
     }
-    else continue;
-
-    // double beta_particle = candidate->Momentum.P()/candidate->Momentum.E();
-    // ti = tf_smeared - candidate->Ld*1.0E-3/(c_light*beta_particle);
-
-    mother = candidate;
-    candidate = static_cast<Candidate*>(candidate->Clone());
-    candidate->AddCandidate(mother);
-    candidate->InitialPosition.SetT((100+ti)*1.0E3*c_light);
-    candidate->Position.SetT(tf_smeared*1.0E3*c_light);
-    candidate->ErrorT = timeResolution*1.0E3*c_light;
-
-    fOutputArray->Add(candidate);
+    else
+    {
+      ti = timeResolution - l*1.0E3/(c_light*beta_particle);   
+      candidate->InitialPosition.SetT(ti);
+      candidate->ErrorT = timeResolution*1.0E3*c_light;         // Do we need to sum with 100 like in upside ?
+      fOutputArray->Add(candidate);
+    }
+   
   }
 }

modules/TimeSmearing.h

@@ -45 +45 @@
 
 private:
-  DelphesFormula *fFormula; //!
+  DelphesFormula *fFormula;
   Double_t fEtaMax;
 

modules/TrackPileUpSubtractor.cc

@@ -117 +117 @@
 void TrackPileUpSubtractor::Process()
 {
-  Candidate *candidate, *particle, *particleTest;
+  Candidate *candidate, *particle;
   map<TIterator *, TObjArray *>::iterator itInputMap;
   TIterator *iterator;
@@ -123 +123 @@
   Double_t z, zvtx = 0;
   Double_t pt, eta, phi, e;
-  Double_t sumPT2 = 0;
-  Double_t sumSquare = 0;
-  int counter = 0;
-  Double_t tempPTSquare = 0;
-  Double_t tempZVertex = 0; 
+
   // find z position of primary vertex
 
-  cout << " ---------- NEW EVENT --------- " << endl; 
-
   fItVertexInputArray->Reset();
-  cout << " NUMBER OF VERTICES : " << fVertexInputArray->GetEntriesFast() << endl;
   while((candidate = static_cast<Candidate *>(fItVertexInputArray->Next())))
   {
-      cout << " ---------- NEW VERTEX --------- " << candidate->IsPU << endl; 
-      tempZVertex = candidate->Position.Z();
-      tempPTSquare = candidate->SumPT2;
-      if(tempPTSquare > sumSquare)
-      {
-        sumSquare = tempPTSquare;
-        zvtx = tempZVertex;
-        cout << " Sum Square : " << sumSquare << " Z of Vertex : " << zvtx << "  Is PileUp : " << candidate->IsPU << endl;
-      } 
-      
-      /*
-      for(int i = 0; i < candidate->GetCandidates()->GetEntriesFast(); i++)
-      {
-        particleTest = static_cast<Candidate *>(candidate->GetCandidates()->At(i));
-        TLorentzVector &candidateMomentumNew = particleTest->Momentum;
-        if(candidateMomentumNew.Pt() > 1.0)
-          sumSquare += (candidateMomentumNew.Pt()*candidateMomentumNew.Pt());
-        cout << candidateMomentumNew.Pt() << "  " << candidate->SumPT2 << " counter : " << candidate->GetCandidates()->GetEntriesFast() << endl;
-      }  
-      */
+    if(!candidate->IsPU)
+    {
+      zvtx = candidate->Position.Z();
+    }
   }
 
@@ -169 +146 @@
     {
       particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
-      const TLorentzVector &candidateMomentum = candidate->Momentum;
+      const TLorentzVector &candidateMomentum = particle->Momentum;
 
       eta = candidateMomentum.Eta();
@@ -176 +153 @@
       e = candidateMomentum.E();
 
-      z = candidate->Position.Z();
-      counter ++; 
-      //sum = pt*pt;
-      if(pt > 1.0)
-        sumPT2 += pt*pt;
+      z = particle->Position.Z();
 
       // apply pile-up subtraction
       // assume perfect pile-up subtraction for tracks outside fZVertexResolution
-      // cout << particle->IsRecoPU << " ParticleSumSquare : " << sumPT2 << " VertexResult : " << sumSquare << " Added SumSquare : " << zvtx << endl;
-      //cout << pt << "  " << candidate->IsPU << "  " <<  TMath::Abs(z - zvtx) << "   " << sumPT2 << "   " << sumPT2 << endl;
 
-      
-      if(particle->Charge != 0 && TMath::Abs(z - zvtx) > fFormula->Eval(pt, eta, phi, e) * 1.0e3)
+      if(candidate->Charge != 0 && candidate->IsPU && TMath::Abs(z - zvtx) > fFormula->Eval(pt, eta, phi, e) * 1.0e3)
       {
         candidate->IsRecoPU = 1;
-        //cout <<  TMath::Abs(z - zvtx) << "   " << fFormula->Eval(pt, eta, phi, e) * 1.0e3 << endl;
       }
       else

modules/TrackTimingPileUpSubtractor.cc

@@ -75 +75 @@
 
   // read resolution formula in m
-  fFormula->Compile(GetString("ZVertexResolution", "0.001"));
+  fChargedMinSignificance = GetDouble("ChargedMinSignificance", 3);
+  fNeutralMinSignificance = GetDouble("NeutralMinSignificance", 3);
 
   fPTMin = GetDouble("PTMin", 0.);
@@ -128 +129 @@
   Double_t tempPTSquare = 0;
   Double_t pt, eta, phi, e;
-  Double_t distance = 0;
+  Double_t distanceCharged, distanceNeutral = 0;
 
   // find z position of primary vertex
@@ -143 +144 @@
       tvtx = candidate->Position.T();
       tvtx_err = candidate->PositionError.T();
-      //cout << " initial : " << candidate->InitialPosition.T() << " final : " << candidate->Position.T() << endl;
     } 
   }
@@ -167 +167 @@
       z = particle->Position.Z();
       z_err = particle->PositionError.Z();
-      t = particle->Position.T();
+      t = particle->InitialPosition.T();
       t_err = particle->PositionError.T();
 
-      // apply pile-up subtraction
-      distance = pow((zvtx - z),2)/pow((zvtx_err - z_err),2) + pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
-      //cout << " t : " << tvtx << "  t(particle)" << t << endl;
-      // here I calculated distance using Z and T of selected vertex (highest sum Pt square) and particles
-      // however z_err of vertices is gives 0 because of using CMS trackResolutionCMS.tcl (in that formula, there is limitation on |eta| < 2.5)
-      // thats why I used TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0
+      distanceCharged = pow((zvtx - z),2)/pow((zvtx_err - z_err),2) + pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
+      distanceNeutral = pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
 
-      if(candidate->Charge != 0 && TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0)
+      if(candidate->Charge != 0 && distanceCharged < fChargedMinSignificance)
       {
         candidate->IsRecoPU = 1;
       }
+      else if(candidate->Charge == 0 && distanceNeutral < fNeutralMinSignificance)
+      {
+        candidate->IsRecoPU = 1;
+      }  
       else
       {

modules/TrackTimingPileUpSubtractor.h

@@ -49 +49 @@
   DelphesFormula *fFormula; //!
 
+  Double_t fChargedMinSignificance;
+  Double_t fNeutralMinSignificance;
+
   Double_t fPTMin;
 

modules/VertexFinderDA4D.cc

@@ -275 +275 @@
   }
 
-  if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Ast");
-
   if( fVerbose > 2){cout << "Cool down untill reaching the temperature to finish increasing the number of vertexes" << endl;}
 
@@ -290 +288 @@
       delta2 = update(beta, tks, vtx, rho0);
 
-      if( fVerbose > 10 ) plot_status(beta, vtx, tks, niter, "Bup");
       if (fVerbose > 3)
       {
@@ -312 +309 @@
       n_it++;
 
-      if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");
     }
 
@@ -320 +316 @@
     {
       split(beta, vtx, tks);
-      if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Asp");
     }
     else
@@ -355 +350 @@
     }
     while (delta2 > 0.3*fD2UpdateLim &&  niter < fMaxIterations);
-    if( fVerbose > 10 ) plot_status(beta, vtx, tks, f, "Dadout");
   }
 
@@ -373 +367 @@
         }
         while (delta2 > fD2UpdateLim &&  niter < fMaxIterations);
-        if( fVerbose > 10 ) plot_status(beta, vtx, tks, i_pu, Form("Eprg%d",min_trk));
         i_pu++;
       }
@@ -390 +383 @@
       n_it++;
 
-      if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");
     }
   } while( beta < fBetaPurge );
@@ -404 +396 @@
       delta2 = update(beta, tks, vtx, rho0);
       niter++;
-      if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Bup");
     }
     while (delta2 > 0.3*fD2UpdateLim &&  niter < fMaxIterations);
@@ -425 +416 @@
     candidate->ClusterIndex = k;
     candidate->Position.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);
+    candidate->InitialPosition.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);    
     candidate->PositionError.SetXYZT(0.0, 0.0, fVertexZSize , fVertexTSize*1E-9*c_light);
     candidate->SumPT2 = 0;
@@ -483 +475 @@
     fTrackOutputArray->Add(tks.tt[i]);
   }
-
-  if(fVerbose > 10) plot_status_end(vtx, tks);
 
 }
@@ -877 +867 @@
         {
           continue;
-          // plot_split_crush(zn, tn, vtx, tks, k);
           // throw std::invalid_argument( "0 division" );
         }
@@ -1042 +1031 @@
   }
 }
-
-
-// -----------------------------------------------------------------------------
-// Plot status
-void VertexFinderDA4D::plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it, const char* flag)
-{
-  vector<int> vtx_color = {2,4,8,1,5,6,9,14,46,3};
-  while(vtx.getSize() > vtx_color.size()) vtx_color.push_back(40);
-
-  vector<double> t_PV, dt_PV, z_PV, dz_PV;
-  vector<double> t_PU, dt_PU, z_PU, dz_PU;
-
-  double ETot = 0;
-  vector<double> pk_exp_mBetaE = Compute_pk_exp_mBetaE(beta, vtx, tks, 0);
-
-  for(unsigned int i = 0; i < tks.getSize(); i++)
-  {
-    for(unsigned int k = 0; k < vtx.getSize(); k++)
-    {
-      unsigned int idx = k*tks.getSize() + i;
-      if(pk_exp_mBetaE[idx] == 0) continue;
-
-      double p_ygx = pk_exp_mBetaE[idx] / tks.Z[i];
-
-      ETot += tks.w[i] * p_ygx * Energy(tks.z[i], vtx.z[k], tks.dz2_o[i], tks.t[i], vtx.t[k], tks.dt2_o[i]);
-    }
-
-    if(tks.tt[i]->IsPU)
-    {
-      t_PU.push_back(tks.t[i]);
-      dt_PU.push_back(1./tks.dt_o[i]);
-      z_PU.push_back(tks.z[i]);
-      dz_PU.push_back(1./tks.dz_o[i]);
-    }
-    else
-    {
-      t_PV.push_back(tks.t[i]);
-      dt_PV.push_back(1./tks.dt_o[i]);
-      z_PV.push_back(tks.z[i]);
-      dz_PV.push_back(1./tks.dz_o[i]);
-    }
-  }
-
-
-  ETot /= tks.sum_w;
-  fEnergy_rec.push_back(ETot);
-  fBeta_rec.push_back(beta);
-  fNvtx_rec.push_back(vtx.getSize());
-
-  double t_min = TMath::Min(  TMath::MinElement(t_PV.size(), &t_PV[0]), TMath::MinElement(t_PU.size(), &t_PU[0])  );
-  t_min = TMath::Min(t_min, TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - fVertexTSize;
-  double t_max = TMath::Max(  TMath::MaxElement(t_PV.size(), &t_PV[0]), TMath::MaxElement(t_PU.size(), &t_PU[0])  );
-  t_max = TMath::Max(t_max, TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + fVertexTSize;
-
-  double z_min = TMath::Min(  TMath::MinElement(z_PV.size(), &z_PV[0]), TMath::MinElement(z_PU.size(), &z_PU[0])  );
-  z_min = TMath::Min(z_min, TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
-  double z_max = TMath::Max(  TMath::MaxElement(z_PV.size(), &z_PV[0]), TMath::MaxElement(z_PU.size(), &z_PU[0])  );
-  z_max = TMath::Max(z_max, TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
-
-  auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
-
-  TGraphErrors* gr_PVtks = new TGraphErrors(t_PV.size(), &t_PV[0], &z_PV[0], &dt_PV[0], &dz_PV[0]);
-  gr_PVtks->SetTitle(Form("Clustering space - #beta = %.6f", beta));
-  gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
-  gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
-  gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
-  gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
-  gr_PVtks->SetMarkerStyle(4);
-  gr_PVtks->SetMarkerColor(8);
-  gr_PVtks->SetLineColor(8);
-  gr_PVtks->Draw("APE1");
-
-  TGraphErrors* gr_PUtks = new TGraphErrors(t_PU.size(), &t_PU[0], &z_PU[0], &dt_PU[0], &dz_PU[0]);
-  gr_PUtks->SetMarkerStyle(3);
-  gr_PUtks->Draw("PE1");
-
-  TGraph* gr_vtx = new TGraph(vtx.getSize(), &(vtx.t[0]), &(vtx.z[0]));
-  gr_vtx->SetMarkerStyle(28);
-  gr_vtx->SetMarkerColor(2);
-  gr_vtx->SetMarkerSize(2.);
-  gr_vtx->Draw("PE1");
-
-  fItInputGenVtx->Reset();
-  TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
-  Candidate *candidate;
-  unsigned int k = 0;
-  while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
-  {
-    gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
-    k++;
-  }
-  gr_genvtx->SetMarkerStyle(33);
-  gr_genvtx->SetMarkerColor(6);
-  gr_genvtx->SetMarkerSize(2.);
-  gr_genvtx->Draw("PE1");
-
-  // auto leg = new TLegend(0.1, 0.1);
-  // leg->AddEntry(gr_PVtks, "PV tks", "ep");
-  // leg->AddEntry(gr_PUtks, "PU tks", "ep");
-  // leg->AddEntry(gr_vtx, "Cluster center", "p");
-  // leg->Draw();
-
-  c_2Dspace->SetGrid();
-  c_2Dspace->SaveAs(fFigFolderPath + Form("/c_2Dspace_beta%010.0f-%s%d.png", 1E7*beta, flag, n_it));
-
-  delete c_2Dspace;
-}
-
-// -----------------------------------------------------------------------------
-// Plot status at the end
-void VertexFinderDA4D::plot_status_end(vertex_t &vtx, tracks_t &tks)
-{
-  unsigned int nv = vtx.getSize();
-
-  // Define colors in a meaningfull way
-  vector<int> MyPalette(nv);
-
-  const int Number = 3;
-  double Red[Number]    = { 1.00, 0.00, 0.00};
-  double Green[Number]  = { 0.00, 1.00, 0.00};
-  double Blue[Number]   = { 1.00, 0.00, 1.00};
-  double Length[Number] = { 0.00, 0.50, 1.00 };
-  int FI = TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nv);
-  for (unsigned int i=0;i<nv;i++) MyPalette[i] = FI+i;
-
-  TCanvas * c_out = new TCanvas("c_out", "c_out", 800, 600);
-  double t_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 2*fVertexTSize;
-  double t_max = TMath::Max(TMath::MaxElement(tks.getSize(), &tks.t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 2*fVertexTSize;
-
-  double z_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 15;
-  double z_max = TMath::Max( TMath::MaxElement(tks.getSize(), &tks.z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 15;
-
-  // Draw tracks
-  for(unsigned int i = 0; i < tks.getSize(); i++)
-  {
-    double dt[] = {1./tks.dt_o[i]};
-    double dz[] = {1./tks.dz_o[i]};
-    TGraphErrors* gr = new TGraphErrors(1, &(tks.t[i]), &(tks.z[i]), dt, dz);
-
-    gr->SetNameTitle(Form("gr%d",i), Form("gr%d",i));
-
-    int marker = tks.tt[i]->IsPU? 1 : 4;
-    gr->SetMarkerStyle(marker);
-
-    int idx = tks.tt[i]->ClusterIndex;
-    int color = idx>=0 ? MyPalette[idx] : 13;
-    gr->SetMarkerColor(color);
-    gr->SetLineColor(color);
-
-    int line_style = idx>=0 ? 1 : 3;
-    gr->SetLineStyle(line_style);
-
-    if(i==0)
-    {
-      gr->SetTitle(Form("Clustering space - Tot Vertexes = %d", nv));
-      gr->GetXaxis()->SetTitle("t CA [ps]");
-      gr->GetXaxis()->SetLimits(t_min, t_max);
-      gr->GetYaxis()->SetTitle("z CA [mm]");
-      gr->GetYaxis()->SetRangeUser(z_min, z_max);
-      gr->Draw("APE1");
-    }
-    else gr->Draw("PE1");
-  }
-
-  // Draw vertices
-  for(unsigned int k = 0; k < vtx.getSize(); k++)
-  {
-    TGraph* gr = new TGraph(1, &(vtx.t[k]), &(vtx.z[k]));
-
-    gr->SetNameTitle(Form("grv%d",k), Form("grv%d",k));
-
-    gr->SetMarkerStyle(41);
-    gr->SetMarkerSize(2.);
-    gr->SetMarkerColor(MyPalette[k]);
-
-    gr->Draw("P");
-  }
-
-  fItInputGenVtx->Reset();
-  TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
-  TGraph* gr_genPV = new TGraph(1);
-  Candidate *candidate;
-  unsigned int k = 0;
-  while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
-  {
-    if(k == 0 ) {
-      gr_genPV->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
-    }
-    else gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
-
-    k++;
-  }
-  gr_genvtx->SetMarkerStyle(20);
-  gr_genvtx->SetMarkerColorAlpha(kBlack, 0.8);
-  gr_genvtx->SetMarkerSize(.8);
-  gr_genvtx->Draw("PE1");
-  gr_genPV->SetMarkerStyle(33);
-  gr_genPV->SetMarkerColorAlpha(kBlack, 1);
-  gr_genPV->SetMarkerSize(2.5);
-  gr_genPV->Draw("PE1");
-
-  // auto note =  new TLatex();
-  // note->DrawLatexNDC(0.5, 0.8, Form("#splitline{Vertexes Reco = %d }{Vertexes gen = %d}", vtx.getSize(), k) );
-
-  c_out->SetGrid();
-  c_out->SaveAs(fFigFolderPath + Form("/c_final.root"));
-  delete c_out;
-}
-
-// -----------------------------------------------------------------------------
-// Plot splitting
-void VertexFinderDA4D::plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx)
-{
-  vector<double> t, dt, z, dz;
-
-  for(unsigned int i = 0; i < tks.getSize(); i++)
-  {
-      t.push_back(tks.t[i]);
-      dt.push_back(1./tks.dt_o[i]);
-      z.push_back(tks.z[i]);
-      dz.push_back(1./tks.dz_o[i]);
-  }
-
-
-  double t_min = TMath::Min(TMath::MinElement(t.size(), &t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 50;
-  double t_max = TMath::Max(TMath::MaxElement(t.size(), &t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 50;
-
-  double z_min = TMath::Min(TMath::MinElement(z.size(), &z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
-  double z_max = TMath::Max(TMath::MaxElement(z.size(), &z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
-
-  auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
-
-  TGraphErrors* gr_PVtks = new TGraphErrors(t.size(), &t[0], &z[0], &dt[0], &dz[0]);
-  gr_PVtks->SetTitle(Form("Clustering space"));
-  gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
-  gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
-  gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
-  gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
-  gr_PVtks->SetMarkerStyle(4);
-  gr_PVtks->SetMarkerColor(1);
-  gr_PVtks->SetLineColor(1);
-  gr_PVtks->Draw("APE1");
-
-  TGraph* gr_vtx = new TGraph(1, &(vtx.t[i_vtx]), &(vtx.z[i_vtx]));
-  gr_vtx->SetMarkerStyle(28);
-  gr_vtx->SetMarkerColor(2);
-  gr_vtx->SetMarkerSize(2.);
-  gr_vtx->Draw("PE1");
-
-  double t_pos[] = {vtx.t[i_vtx], vtx.t[i_vtx]+100};
-  double t_neg[] = {vtx.t[i_vtx], vtx.t[i_vtx]-100};
-  double z_pos[] = {vtx.z[i_vtx], vtx.z[i_vtx]+(zn/tn)*100};
-  double z_neg[] = {vtx.z[i_vtx], vtx.z[i_vtx]-(zn/tn)*100};
-
-  TGraph* gr_pos = new TGraph(2, &t_pos[0], &z_pos[0]);
-  gr_pos->SetLineColor(8);
-  gr_pos->SetMarkerColor(8);
-  gr_pos->Draw("PL");
-  TGraph* gr_neg = new TGraph(2, &t_neg[0], &z_neg[0]);
-  gr_neg->SetLineColor(4);
-  gr_neg->SetMarkerColor(4);
-  gr_neg->Draw("PL");
-
-
-  c_2Dspace->SetGrid();
-  c_2Dspace->SaveAs(fFigFolderPath + Form("/crush_splitting.png"));
-
-  delete c_2Dspace;
-}

modules/VertexFinderDA4D.h

@@ -234 +234 @@
     std::vector<double> Compute_pk_exp_mBetaE(double beta, vertex_t &vtx, tracks_t &tks, double Z_init);
 
-    // Plot status of tracks and Vertices
-    void plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it = 0, const char* flag ="");
-
-    // Plot status at the end of the fitting
-    void plot_status_end(vertex_t &vtx, tracks_t &tks);
-
-    // Plot Crush
-    void plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx);
-
 
   private: