Fork me on GitHub

Changeset 79a7b3e in git


Ignore:
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:
6 added
9 edited

Legend:

Unmodified
Added
Removed
  • cards/FCC/FCChh_PileUpVtx.tcl

    r364dbe1 r79a7b3e  
    1717set ExecutionPath {
    1818
    19   BeamSpotFilter
    2019  PileUpMerger
    2120  ParticlePropagator
     
    3332  TrackMerger
    3433
    35 
    3634  TrackSmearing
    37   TimeSmearing 
    38 
    39   VertexFinderDA4D 
    40 
    41   TrackTimingPileUpSubtractor 
    4235
    4336  ECal
     
    4740  EFlowMerger
    4841  EFlowFilter
     42
     43  TimeSmearingMIP
     44  TimeSmearingPhotons
     45  TimeSmearingNH   
     46
     47  VertexFinderDA4D
     48  PileUpSubtractor4D
     49
     50  HighMassVertexRecover   
    4951
    5052  PhotonEfficiency
     
    8183
    8284  TreeWriter
    83 }
    84 
    85 #######################
    86 # GenBeamSpotFilter
    87 # Saves a particle intended to represent the beamspot
    88 #######################
    89 
    90 module BeamSpotFilter BeamSpotFilter {
    91     set InputArray Delphes/stableParticles
    92     set OutputArray beamSpotParticle
    93 
    9485}
    9586
     
    308299
    309300  # from http://mersi.web.cern.ch/mersi/layouts/.private/Baseline_tilted_200_Pixel_1_1_1/index.html
    310   source trackResolutionCMS.tcl
     301  source trackResolutionFCChh.tcl
    311302  # FIXME !!!! we need to add track resolution of FCC-hh baseline detector !!!!!
    312 }
    313 
    314 ########################################
    315 #   Time Smearing
    316 ########################################
    317 
    318 module TimeSmearing TimeSmearing {
    319   set InputArray TrackSmearing/tracks
    320   set OutputArray tracks
    321 
    322   # assume 20 ps resolution for now
    323   set TimeResolution 20E-12
    324 }
    325 
    326 ##################################
    327 # Primary vertex reconstruction
    328 ##################################
    329 
    330 
    331 module VertexFinderDA4D VertexFinderDA4D {
    332   set InputArray TimeSmearing/tracks
    333 
    334   set OutputArray tracks
    335   set VertexOutputArray vertices
    336 
    337   set Verbose 0
    338   set MinPT 1.0
    339 
    340   # in mm
    341   set VertexSpaceSize 0.5
    342 
    343   # in s
    344   set VertexTimeSize 10E-12
    345 
    346   set UseTc 1
    347   set BetaMax 0.1
    348   set BetaStop 1.0
    349   set CoolingFactor 0.8
    350   set MaxIterations 100
    351 
    352   # in mm
    353   set DzCutOff 40
    354   set D0CutOff 30
    355 
    356 }
    357 
    358 ##########################
    359 # Track pile-up subtractor
    360 ##########################
    361 
    362 module TrackTimingPileUpSubtractor TrackTimingPileUpSubtractor {
    363 # add InputArray InputArray OutputArray
    364 
    365   add InputArray ChargedHadronMomentumSmearing/chargedHadrons
    366   add InputArray ElectronMomentumSmearing/electrons
    367   add InputArray MuonMomentumSmearing/muons
    368  
    369   set VertexInputArray VertexFinderDA4D/vertices
    370   # assume perfect pile-up subtraction for tracks with |z| > fZVertexResolution
    371   # Z vertex resolution in m
    372   set ZVertexResolution {0.0001}
    373303}
    374304
     
    539469}
    540470
    541 
    542471#################
    543472# Electron filter
     
    606535}
    607536
     537########################################
     538#   Time Smearing Neutral MIP
     539########################################
     540
     541module TimeSmearing TimeSmearingMIP {
     542  set InputArray HCal/eflowTracks
     543  set OutputArray tracks
     544
     545  # assume 30 ps resolution for now
     546  set TimeResolution {30E-12}
     547}
     548
     549########################################
     550#   Time Smearing Neutral Photons
     551########################################
     552
     553module TimeSmearing TimeSmearingPhotons {
     554  set InputArray ECal/eflowPhotons
     555  set OutputArray photons
     556  set TimeResolution {sqrt(20^2 + 150^2)/energy^2}
     557}
     558
     559########################################
     560#   Time Smearing Neutral NeutralHadrons
     561########################################
     562#
     563module TimeSmearing TimeSmearingNH {
     564  set InputArray HCal/eflowNeutralHadrons
     565  set OutputArray neutralhadrons
     566
     567  # assume 30 ps resolution for now
     568  set TimeResolution {sqrt(20^2 + 150^2)/energy^2}
     569}
     570
     571
     572##################################
     573# Primary vertex reconstruction
     574##################################
     575
     576
     577module VertexFinderDA4D VertexFinderDA4D {
     578  set InputArray TimeSmearingMIP/tracks
     579
     580  set OutputArray tracks
     581  set VertexOutputArray vertices
     582
     583  set Verbose 0
     584  set MinPT 1.0
     585
     586  # in mm
     587  set VertexSpaceSize 0.5
     588
     589  # in s
     590  set VertexTimeSize 10E-12
     591
     592  set UseTc 1
     593  set BetaMax 0.1
     594  set BetaStop 1.0
     595  set CoolingFactor 0.8
     596  set MaxIterations 100
     597
     598  # in mm
     599  set DzCutOff 40
     600  set D0CutOff 30
     601
     602}
     603
     604##########################
     605# Track pile-up subtractor
     606##########################
     607
     608module PileUpSubtractor4D PileUpSubtractor4D {
     609# add InputArray InputArray OutputArray
     610
     611  add InputArray TimeSmearingMIP/tracks
     612  add InputArray TimeSmearingPhotons/photons
     613  add InputArray TimeSmearingNH/neutralhadrons
     614
     615  set VertexInputArray VertexFinderDA4D/vertices
     616
     617  set fChargedMinSignificance {3}
     618  set fNeutralMinSignificance {3}
     619}
     620
     621######################################
     622# Heavy(slow) particles vertex recover
     623######################################
     624
     625module HighMassVertexRecover HighMassVertexRecover {
     626
     627  set TrackInputArray VertexFinderDA4D/tracks
     628  set VertexInputArray VertexFinderDA4D/vertices
     629
     630  set TrackOutputArray tracks
     631  set VertexOutputArray vertices
     632
     633  set Verbose 0
     634
     635}
    608636
    609637###################
     
    10571085  add Branch ScalarHT/energy ScalarHT ScalarHT
    10581086  add Branch VertexFinderDA4D/vertices Vertex4D Vertex
    1059 }
    1060 
     1087
     1088  add Branch HighMassVertexRecover/tracks Track Track
     1089}
     1090
  • modules/ModulesLinkDef.h

    r364dbe1 r79a7b3e  
    5656#include "modules/PileUpMerger.h"
    5757#include "modules/JetPileUpSubtractor.h"
    58 #include "modules/TrackPileUpSubtractor.h"
    5958#include "modules/TrackTimingPileUpSubtractor.h"
    6059#include "modules/TaggingParticlesSkimmer.h"
     
    114113#pragma link C++ class PileUpMerger+;
    115114#pragma link C++ class JetPileUpSubtractor+;
    116 #pragma link C++ class TrackPileUpSubtractor+;
    117115#pragma link C++ class TrackTimingPileUpSubtractor+;
    118116#pragma link C++ class TaggingParticlesSkimmer+;
  • modules/TimeSmearing.cc

    r364dbe1 r79a7b3e  
    5353
    5454TimeSmearing::TimeSmearing() :
    55   fFormula(0), fItInputArray(0)
     55  fItInputArray(0), fFormula(0)
    5656{
    57   fFormula = new DelphesFormula;
     57        fFormula = new DelphesFormula;
    5858}
    5959
     
    6161
    6262TimeSmearing::~TimeSmearing()
    63 {
    64   if(fFormula) delete fFormula;
     63{ 
     64        if(fFormula) delete fFormula;
    6565}
    6666
     
    6969void TimeSmearing::Init()
    7070{
    71   // read resolution formula
    72 
    73   fFormula->Compile(GetString("TimeResolution", "1.0"));
     71  // read resolution formula   
     72  fFormula->Compile(GetString("TimeResolution", "0.001"));
    7473
    7574  // import input array
     
    7776  fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons"));
    7877  fItInputArray = fInputArray->MakeIterator();
    79 
    8078  // create output array
    8179
     
    9593{
    9694  Candidate *candidate, *mother;
    97   Double_t ti, tf_smeared, tf, timeResolution;
    98   Double_t pt, eta, phi, e, d0, dz, ctgTheta;
     95  Double_t ti, tf_smeared, tf;
     96  Double_t pt, eta, phi, e, p, l;
     97  Double_t beta_particle;
     98  Double_t timeResolution;
    9999
    100100
     
    114114    pt = candidateMomentum.Pt();
    115115    e = candidateMomentum.E();
    116     d0 = candidate->D0;
    117     dz = candidate->DZ;
    118     ctgTheta = candidate->CtgTheta;
    119 
    120     // apply smearing formula
     116    p = candidateMomentum.P();
     117    beta_particle = p/e;
     118    l = candidate->L;
     119    timeResolution = fFormula->Eval(e);
    121120   
    122     timeResolution = fFormula->Eval(pt, eta, phi, e, d0, dz, ctgTheta);
    123     if(fabs(candidate->Position.Eta())<fEtaMax)
     121    if(candidate->Charge != 0)
    124122    {
    125123      tf_smeared = tf + timeResolution*gRandom->Gaus(0, 1);
     124      mother = candidate;
     125      candidate = static_cast<Candidate*>(candidate->Clone());  // I am not sure that we need these lines !!!
     126      candidate->AddCandidate(mother);
     127      candidate->InitialPosition.SetT((100+ti)*1.0E3*c_light);  // Dummy Value, correct value will be computed by VertexFinderDA4D
     128      candidate->Position.SetT(tf_smeared*1.0E3*c_light);
     129      candidate->ErrorT = timeResolution*1.0E3*c_light;
     130      fOutputArray->Add(candidate);
    126131    }
    127     else continue;
    128 
    129     // double beta_particle = candidate->Momentum.P()/candidate->Momentum.E();
    130     // ti = tf_smeared - candidate->Ld*1.0E-3/(c_light*beta_particle);
    131 
    132     mother = candidate;
    133     candidate = static_cast<Candidate*>(candidate->Clone());
    134     candidate->AddCandidate(mother);
    135     candidate->InitialPosition.SetT((100+ti)*1.0E3*c_light);
    136     candidate->Position.SetT(tf_smeared*1.0E3*c_light);
    137     candidate->ErrorT = timeResolution*1.0E3*c_light;
    138 
    139     fOutputArray->Add(candidate);
     132    else
     133    {
     134      ti = timeResolution - l*1.0E3/(c_light*beta_particle);   
     135      candidate->InitialPosition.SetT(ti);
     136      candidate->ErrorT = timeResolution*1.0E3*c_light;         // Do we need to sum with 100 like in upside ?
     137      fOutputArray->Add(candidate);
     138    }
     139   
    140140  }
    141141}
  • modules/TimeSmearing.h

    r364dbe1 r79a7b3e  
    4545
    4646private:
    47   DelphesFormula *fFormula; //!
     47  DelphesFormula *fFormula;
    4848  Double_t fEtaMax;
    4949
  • modules/TrackPileUpSubtractor.cc

    r364dbe1 r79a7b3e  
    117117void TrackPileUpSubtractor::Process()
    118118{
    119   Candidate *candidate, *particle, *particleTest;
     119  Candidate *candidate, *particle;
    120120  map<TIterator *, TObjArray *>::iterator itInputMap;
    121121  TIterator *iterator;
     
    123123  Double_t z, zvtx = 0;
    124124  Double_t pt, eta, phi, e;
    125   Double_t sumPT2 = 0;
    126   Double_t sumSquare = 0;
    127   int counter = 0;
    128   Double_t tempPTSquare = 0;
    129   Double_t tempZVertex = 0;
     125
    130126  // find z position of primary vertex
    131127
    132   cout << " ---------- NEW EVENT --------- " << endl;
    133 
    134128  fItVertexInputArray->Reset();
    135   cout << " NUMBER OF VERTICES : " << fVertexInputArray->GetEntriesFast() << endl;
    136129  while((candidate = static_cast<Candidate *>(fItVertexInputArray->Next())))
    137130  {
    138       cout << " ---------- NEW VERTEX --------- " << candidate->IsPU << endl;
    139       tempZVertex = candidate->Position.Z();
    140       tempPTSquare = candidate->SumPT2;
    141       if(tempPTSquare > sumSquare)
    142       {
    143         sumSquare = tempPTSquare;
    144         zvtx = tempZVertex;
    145         cout << " Sum Square : " << sumSquare << " Z of Vertex : " << zvtx << "  Is PileUp : " << candidate->IsPU << endl;
    146       }
    147      
    148       /*
    149       for(int i = 0; i < candidate->GetCandidates()->GetEntriesFast(); i++)
    150       {
    151         particleTest = static_cast<Candidate *>(candidate->GetCandidates()->At(i));
    152         TLorentzVector &candidateMomentumNew = particleTest->Momentum;
    153         if(candidateMomentumNew.Pt() > 1.0)
    154           sumSquare += (candidateMomentumNew.Pt()*candidateMomentumNew.Pt());
    155         cout << candidateMomentumNew.Pt() << "  " << candidate->SumPT2 << " counter : " << candidate->GetCandidates()->GetEntriesFast() << endl;
    156       } 
    157       */
     131    if(!candidate->IsPU)
     132    {
     133      zvtx = candidate->Position.Z();
     134    }
    158135  }
    159136
     
    169146    {
    170147      particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    171       const TLorentzVector &candidateMomentum = candidate->Momentum;
     148      const TLorentzVector &candidateMomentum = particle->Momentum;
    172149
    173150      eta = candidateMomentum.Eta();
     
    176153      e = candidateMomentum.E();
    177154
    178       z = candidate->Position.Z();
    179       counter ++;
    180       //sum = pt*pt;
    181       if(pt > 1.0)
    182         sumPT2 += pt*pt;
     155      z = particle->Position.Z();
    183156
    184157      // apply pile-up subtraction
    185158      // assume perfect pile-up subtraction for tracks outside fZVertexResolution
    186       // cout << particle->IsRecoPU << " ParticleSumSquare : " << sumPT2 << " VertexResult : " << sumSquare << " Added SumSquare : " << zvtx << endl;
    187       //cout << pt << "  " << candidate->IsPU << "  " <<  TMath::Abs(z - zvtx) << "   " << sumPT2 << "   " << sumPT2 << endl;
    188159
    189      
    190       if(particle->Charge != 0 && TMath::Abs(z - zvtx) > fFormula->Eval(pt, eta, phi, e) * 1.0e3)
     160      if(candidate->Charge != 0 && candidate->IsPU && TMath::Abs(z - zvtx) > fFormula->Eval(pt, eta, phi, e) * 1.0e3)
    191161      {
    192162        candidate->IsRecoPU = 1;
    193         //cout <<  TMath::Abs(z - zvtx) << "   " << fFormula->Eval(pt, eta, phi, e) * 1.0e3 << endl;
    194163      }
    195164      else
  • modules/TrackTimingPileUpSubtractor.cc

    r364dbe1 r79a7b3e  
    7575
    7676  // read resolution formula in m
    77   fFormula->Compile(GetString("ZVertexResolution", "0.001"));
     77  fChargedMinSignificance = GetDouble("ChargedMinSignificance", 3);
     78  fNeutralMinSignificance = GetDouble("NeutralMinSignificance", 3);
    7879
    7980  fPTMin = GetDouble("PTMin", 0.);
     
    128129  Double_t tempPTSquare = 0;
    129130  Double_t pt, eta, phi, e;
    130   Double_t distance = 0;
     131  Double_t distanceCharged, distanceNeutral = 0;
    131132
    132133  // find z position of primary vertex
     
    143144      tvtx = candidate->Position.T();
    144145      tvtx_err = candidate->PositionError.T();
    145       //cout << " initial : " << candidate->InitialPosition.T() << " final : " << candidate->Position.T() << endl;
    146146    }
    147147  }
     
    167167      z = particle->Position.Z();
    168168      z_err = particle->PositionError.Z();
    169       t = particle->Position.T();
     169      t = particle->InitialPosition.T();
    170170      t_err = particle->PositionError.T();
    171171
    172       // apply pile-up subtraction
    173       distance = pow((zvtx - z),2)/pow((zvtx_err - z_err),2) + pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
    174       //cout << " t : " << tvtx << "  t(particle)" << t << endl;
    175       // here I calculated distance using Z and T of selected vertex (highest sum Pt square) and particles
    176       // however z_err of vertices is gives 0 because of using CMS trackResolutionCMS.tcl (in that formula, there is limitation on |eta| < 2.5)
    177       // thats why I used TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0
     172      distanceCharged = pow((zvtx - z),2)/pow((zvtx_err - z_err),2) + pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
     173      distanceNeutral = pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
    178174
    179       if(candidate->Charge != 0 && TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0)
     175      if(candidate->Charge != 0 && distanceCharged < fChargedMinSignificance)
    180176      {
    181177        candidate->IsRecoPU = 1;
    182178      }
     179      else if(candidate->Charge == 0 && distanceNeutral < fNeutralMinSignificance)
     180      {
     181        candidate->IsRecoPU = 1;
     182      } 
    183183      else
    184184      {
  • modules/TrackTimingPileUpSubtractor.h

    r364dbe1 r79a7b3e  
    4949  DelphesFormula *fFormula; //!
    5050
     51  Double_t fChargedMinSignificance;
     52  Double_t fNeutralMinSignificance;
     53
    5154  Double_t fPTMin;
    5255
  • modules/VertexFinderDA4D.cc

    r364dbe1 r79a7b3e  
    275275  }
    276276
    277   if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Ast");
    278 
    279277  if( fVerbose > 2){cout << "Cool down untill reaching the temperature to finish increasing the number of vertexes" << endl;}
    280278
     
    290288      delta2 = update(beta, tks, vtx, rho0);
    291289
    292       if( fVerbose > 10 ) plot_status(beta, vtx, tks, niter, "Bup");
    293290      if (fVerbose > 3)
    294291      {
     
    312309      n_it++;
    313310
    314       if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");
    315311    }
    316312
     
    320316    {
    321317      split(beta, vtx, tks);
    322       if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Asp");
    323318    }
    324319    else
     
    355350    }
    356351    while (delta2 > 0.3*fD2UpdateLim &&  niter < fMaxIterations);
    357     if( fVerbose > 10 ) plot_status(beta, vtx, tks, f, "Dadout");
    358352  }
    359353
     
    373367        }
    374368        while (delta2 > fD2UpdateLim &&  niter < fMaxIterations);
    375         if( fVerbose > 10 ) plot_status(beta, vtx, tks, i_pu, Form("Eprg%d",min_trk));
    376369        i_pu++;
    377370      }
     
    390383      n_it++;
    391384
    392       if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");
    393385    }
    394386  } while( beta < fBetaPurge );
     
    404396      delta2 = update(beta, tks, vtx, rho0);
    405397      niter++;
    406       if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Bup");
    407398    }
    408399    while (delta2 > 0.3*fD2UpdateLim &&  niter < fMaxIterations);
     
    425416    candidate->ClusterIndex = k;
    426417    candidate->Position.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);
     418    candidate->InitialPosition.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);   
    427419    candidate->PositionError.SetXYZT(0.0, 0.0, fVertexZSize , fVertexTSize*1E-9*c_light);
    428420    candidate->SumPT2 = 0;
     
    483475    fTrackOutputArray->Add(tks.tt[i]);
    484476  }
    485 
    486   if(fVerbose > 10) plot_status_end(vtx, tks);
    487477
    488478}
     
    877867        {
    878868          continue;
    879           // plot_split_crush(zn, tn, vtx, tks, k);
    880869          // throw std::invalid_argument( "0 division" );
    881870        }
     
    10421031  }
    10431032}
    1044 
    1045 
    1046 // -----------------------------------------------------------------------------
    1047 // Plot status
    1048 void VertexFinderDA4D::plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it, const char* flag)
    1049 {
    1050   vector<int> vtx_color = {2,4,8,1,5,6,9,14,46,3};
    1051   while(vtx.getSize() > vtx_color.size()) vtx_color.push_back(40);
    1052 
    1053   vector<double> t_PV, dt_PV, z_PV, dz_PV;
    1054   vector<double> t_PU, dt_PU, z_PU, dz_PU;
    1055 
    1056   double ETot = 0;
    1057   vector<double> pk_exp_mBetaE = Compute_pk_exp_mBetaE(beta, vtx, tks, 0);
    1058 
    1059   for(unsigned int i = 0; i < tks.getSize(); i++)
    1060   {
    1061     for(unsigned int k = 0; k < vtx.getSize(); k++)
    1062     {
    1063       unsigned int idx = k*tks.getSize() + i;
    1064       if(pk_exp_mBetaE[idx] == 0) continue;
    1065 
    1066       double p_ygx = pk_exp_mBetaE[idx] / tks.Z[i];
    1067 
    1068       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]);
    1069     }
    1070 
    1071     if(tks.tt[i]->IsPU)
    1072     {
    1073       t_PU.push_back(tks.t[i]);
    1074       dt_PU.push_back(1./tks.dt_o[i]);
    1075       z_PU.push_back(tks.z[i]);
    1076       dz_PU.push_back(1./tks.dz_o[i]);
    1077     }
    1078     else
    1079     {
    1080       t_PV.push_back(tks.t[i]);
    1081       dt_PV.push_back(1./tks.dt_o[i]);
    1082       z_PV.push_back(tks.z[i]);
    1083       dz_PV.push_back(1./tks.dz_o[i]);
    1084     }
    1085   }
    1086 
    1087 
    1088   ETot /= tks.sum_w;
    1089   fEnergy_rec.push_back(ETot);
    1090   fBeta_rec.push_back(beta);
    1091   fNvtx_rec.push_back(vtx.getSize());
    1092 
    1093   double t_min = TMath::Min(  TMath::MinElement(t_PV.size(), &t_PV[0]), TMath::MinElement(t_PU.size(), &t_PU[0])  );
    1094   t_min = TMath::Min(t_min, TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - fVertexTSize;
    1095   double t_max = TMath::Max(  TMath::MaxElement(t_PV.size(), &t_PV[0]), TMath::MaxElement(t_PU.size(), &t_PU[0])  );
    1096   t_max = TMath::Max(t_max, TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + fVertexTSize;
    1097 
    1098   double z_min = TMath::Min(  TMath::MinElement(z_PV.size(), &z_PV[0]), TMath::MinElement(z_PU.size(), &z_PU[0])  );
    1099   z_min = TMath::Min(z_min, TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
    1100   double z_max = TMath::Max(  TMath::MaxElement(z_PV.size(), &z_PV[0]), TMath::MaxElement(z_PU.size(), &z_PU[0])  );
    1101   z_max = TMath::Max(z_max, TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
    1102 
    1103   auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
    1104 
    1105   TGraphErrors* gr_PVtks = new TGraphErrors(t_PV.size(), &t_PV[0], &z_PV[0], &dt_PV[0], &dz_PV[0]);
    1106   gr_PVtks->SetTitle(Form("Clustering space - #beta = %.6f", beta));
    1107   gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
    1108   gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
    1109   gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
    1110   gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
    1111   gr_PVtks->SetMarkerStyle(4);
    1112   gr_PVtks->SetMarkerColor(8);
    1113   gr_PVtks->SetLineColor(8);
    1114   gr_PVtks->Draw("APE1");
    1115 
    1116   TGraphErrors* gr_PUtks = new TGraphErrors(t_PU.size(), &t_PU[0], &z_PU[0], &dt_PU[0], &dz_PU[0]);
    1117   gr_PUtks->SetMarkerStyle(3);
    1118   gr_PUtks->Draw("PE1");
    1119 
    1120   TGraph* gr_vtx = new TGraph(vtx.getSize(), &(vtx.t[0]), &(vtx.z[0]));
    1121   gr_vtx->SetMarkerStyle(28);
    1122   gr_vtx->SetMarkerColor(2);
    1123   gr_vtx->SetMarkerSize(2.);
    1124   gr_vtx->Draw("PE1");
    1125 
    1126   fItInputGenVtx->Reset();
    1127   TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
    1128   Candidate *candidate;
    1129   unsigned int k = 0;
    1130   while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
    1131   {
    1132     gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1133     k++;
    1134   }
    1135   gr_genvtx->SetMarkerStyle(33);
    1136   gr_genvtx->SetMarkerColor(6);
    1137   gr_genvtx->SetMarkerSize(2.);
    1138   gr_genvtx->Draw("PE1");
    1139 
    1140   // auto leg = new TLegend(0.1, 0.1);
    1141   // leg->AddEntry(gr_PVtks, "PV tks", "ep");
    1142   // leg->AddEntry(gr_PUtks, "PU tks", "ep");
    1143   // leg->AddEntry(gr_vtx, "Cluster center", "p");
    1144   // leg->Draw();
    1145 
    1146   c_2Dspace->SetGrid();
    1147   c_2Dspace->SaveAs(fFigFolderPath + Form("/c_2Dspace_beta%010.0f-%s%d.png", 1E7*beta, flag, n_it));
    1148 
    1149   delete c_2Dspace;
    1150 }
    1151 
    1152 // -----------------------------------------------------------------------------
    1153 // Plot status at the end
    1154 void VertexFinderDA4D::plot_status_end(vertex_t &vtx, tracks_t &tks)
    1155 {
    1156   unsigned int nv = vtx.getSize();
    1157 
    1158   // Define colors in a meaningfull way
    1159   vector<int> MyPalette(nv);
    1160 
    1161   const int Number = 3;
    1162   double Red[Number]    = { 1.00, 0.00, 0.00};
    1163   double Green[Number]  = { 0.00, 1.00, 0.00};
    1164   double Blue[Number]   = { 1.00, 0.00, 1.00};
    1165   double Length[Number] = { 0.00, 0.50, 1.00 };
    1166   int FI = TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nv);
    1167   for (unsigned int i=0;i<nv;i++) MyPalette[i] = FI+i;
    1168 
    1169   TCanvas * c_out = new TCanvas("c_out", "c_out", 800, 600);
    1170   double t_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 2*fVertexTSize;
    1171   double t_max = TMath::Max(TMath::MaxElement(tks.getSize(), &tks.t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 2*fVertexTSize;
    1172 
    1173   double z_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 15;
    1174   double z_max = TMath::Max( TMath::MaxElement(tks.getSize(), &tks.z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 15;
    1175 
    1176   // Draw tracks
    1177   for(unsigned int i = 0; i < tks.getSize(); i++)
    1178   {
    1179     double dt[] = {1./tks.dt_o[i]};
    1180     double dz[] = {1./tks.dz_o[i]};
    1181     TGraphErrors* gr = new TGraphErrors(1, &(tks.t[i]), &(tks.z[i]), dt, dz);
    1182 
    1183     gr->SetNameTitle(Form("gr%d",i), Form("gr%d",i));
    1184 
    1185     int marker = tks.tt[i]->IsPU? 1 : 4;
    1186     gr->SetMarkerStyle(marker);
    1187 
    1188     int idx = tks.tt[i]->ClusterIndex;
    1189     int color = idx>=0 ? MyPalette[idx] : 13;
    1190     gr->SetMarkerColor(color);
    1191     gr->SetLineColor(color);
    1192 
    1193     int line_style = idx>=0 ? 1 : 3;
    1194     gr->SetLineStyle(line_style);
    1195 
    1196     if(i==0)
    1197     {
    1198       gr->SetTitle(Form("Clustering space - Tot Vertexes = %d", nv));
    1199       gr->GetXaxis()->SetTitle("t CA [ps]");
    1200       gr->GetXaxis()->SetLimits(t_min, t_max);
    1201       gr->GetYaxis()->SetTitle("z CA [mm]");
    1202       gr->GetYaxis()->SetRangeUser(z_min, z_max);
    1203       gr->Draw("APE1");
    1204     }
    1205     else gr->Draw("PE1");
    1206   }
    1207 
    1208   // Draw vertices
    1209   for(unsigned int k = 0; k < vtx.getSize(); k++)
    1210   {
    1211     TGraph* gr = new TGraph(1, &(vtx.t[k]), &(vtx.z[k]));
    1212 
    1213     gr->SetNameTitle(Form("grv%d",k), Form("grv%d",k));
    1214 
    1215     gr->SetMarkerStyle(41);
    1216     gr->SetMarkerSize(2.);
    1217     gr->SetMarkerColor(MyPalette[k]);
    1218 
    1219     gr->Draw("P");
    1220   }
    1221 
    1222   fItInputGenVtx->Reset();
    1223   TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
    1224   TGraph* gr_genPV = new TGraph(1);
    1225   Candidate *candidate;
    1226   unsigned int k = 0;
    1227   while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
    1228   {
    1229     if(k == 0 ) {
    1230       gr_genPV->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1231     }
    1232     else gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1233 
    1234     k++;
    1235   }
    1236   gr_genvtx->SetMarkerStyle(20);
    1237   gr_genvtx->SetMarkerColorAlpha(kBlack, 0.8);
    1238   gr_genvtx->SetMarkerSize(.8);
    1239   gr_genvtx->Draw("PE1");
    1240   gr_genPV->SetMarkerStyle(33);
    1241   gr_genPV->SetMarkerColorAlpha(kBlack, 1);
    1242   gr_genPV->SetMarkerSize(2.5);
    1243   gr_genPV->Draw("PE1");
    1244 
    1245   // auto note =  new TLatex();
    1246   // note->DrawLatexNDC(0.5, 0.8, Form("#splitline{Vertexes Reco = %d }{Vertexes gen = %d}", vtx.getSize(), k) );
    1247 
    1248   c_out->SetGrid();
    1249   c_out->SaveAs(fFigFolderPath + Form("/c_final.root"));
    1250   delete c_out;
    1251 }
    1252 
    1253 // -----------------------------------------------------------------------------
    1254 // Plot splitting
    1255 void VertexFinderDA4D::plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx)
    1256 {
    1257   vector<double> t, dt, z, dz;
    1258 
    1259   for(unsigned int i = 0; i < tks.getSize(); i++)
    1260   {
    1261       t.push_back(tks.t[i]);
    1262       dt.push_back(1./tks.dt_o[i]);
    1263       z.push_back(tks.z[i]);
    1264       dz.push_back(1./tks.dz_o[i]);
    1265   }
    1266 
    1267 
    1268   double t_min = TMath::Min(TMath::MinElement(t.size(), &t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 50;
    1269   double t_max = TMath::Max(TMath::MaxElement(t.size(), &t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 50;
    1270 
    1271   double z_min = TMath::Min(TMath::MinElement(z.size(), &z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
    1272   double z_max = TMath::Max(TMath::MaxElement(z.size(), &z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
    1273 
    1274   auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
    1275 
    1276   TGraphErrors* gr_PVtks = new TGraphErrors(t.size(), &t[0], &z[0], &dt[0], &dz[0]);
    1277   gr_PVtks->SetTitle(Form("Clustering space"));
    1278   gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
    1279   gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
    1280   gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
    1281   gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
    1282   gr_PVtks->SetMarkerStyle(4);
    1283   gr_PVtks->SetMarkerColor(1);
    1284   gr_PVtks->SetLineColor(1);
    1285   gr_PVtks->Draw("APE1");
    1286 
    1287   TGraph* gr_vtx = new TGraph(1, &(vtx.t[i_vtx]), &(vtx.z[i_vtx]));
    1288   gr_vtx->SetMarkerStyle(28);
    1289   gr_vtx->SetMarkerColor(2);
    1290   gr_vtx->SetMarkerSize(2.);
    1291   gr_vtx->Draw("PE1");
    1292 
    1293   double t_pos[] = {vtx.t[i_vtx], vtx.t[i_vtx]+100};
    1294   double t_neg[] = {vtx.t[i_vtx], vtx.t[i_vtx]-100};
    1295   double z_pos[] = {vtx.z[i_vtx], vtx.z[i_vtx]+(zn/tn)*100};
    1296   double z_neg[] = {vtx.z[i_vtx], vtx.z[i_vtx]-(zn/tn)*100};
    1297 
    1298   TGraph* gr_pos = new TGraph(2, &t_pos[0], &z_pos[0]);
    1299   gr_pos->SetLineColor(8);
    1300   gr_pos->SetMarkerColor(8);
    1301   gr_pos->Draw("PL");
    1302   TGraph* gr_neg = new TGraph(2, &t_neg[0], &z_neg[0]);
    1303   gr_neg->SetLineColor(4);
    1304   gr_neg->SetMarkerColor(4);
    1305   gr_neg->Draw("PL");
    1306 
    1307 
    1308   c_2Dspace->SetGrid();
    1309   c_2Dspace->SaveAs(fFigFolderPath + Form("/crush_splitting.png"));
    1310 
    1311   delete c_2Dspace;
    1312 }
  • modules/VertexFinderDA4D.h

    r364dbe1 r79a7b3e  
    234234    std::vector<double> Compute_pk_exp_mBetaE(double beta, vertex_t &vtx, tracks_t &tks, double Z_init);
    235235
    236     // Plot status of tracks and Vertices
    237     void plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it = 0, const char* flag ="");
    238 
    239     // Plot status at the end of the fitting
    240     void plot_status_end(vertex_t &vtx, tracks_t &tks);
    241 
    242     // Plot Crush
    243     void plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx);
    244 
    245236
    246237  private:
Note: See TracChangeset for help on using the changeset viewer.