Changeset 79a7b3e in git
- Timestamp:
- Jan 24, 2020, 3:55:17 PM (5 years ago)
- 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)
- Files:
-
- 6 added
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
cards/FCC/FCChh_PileUpVtx.tcl
r364dbe1 r79a7b3e 17 17 set ExecutionPath { 18 18 19 BeamSpotFilter20 19 PileUpMerger 21 20 ParticlePropagator … … 33 32 TrackMerger 34 33 35 36 34 TrackSmearing 37 TimeSmearing38 39 VertexFinderDA4D40 41 TrackTimingPileUpSubtractor42 35 43 36 ECal … … 47 40 EFlowMerger 48 41 EFlowFilter 42 43 TimeSmearingMIP 44 TimeSmearingPhotons 45 TimeSmearingNH 46 47 VertexFinderDA4D 48 PileUpSubtractor4D 49 50 HighMassVertexRecover 49 51 50 52 PhotonEfficiency … … 81 83 82 84 TreeWriter 83 }84 85 #######################86 # GenBeamSpotFilter87 # Saves a particle intended to represent the beamspot88 #######################89 90 module BeamSpotFilter BeamSpotFilter {91 set InputArray Delphes/stableParticles92 set OutputArray beamSpotParticle93 94 85 } 95 86 … … 308 299 309 300 # from http://mersi.web.cern.ch/mersi/layouts/.private/Baseline_tilted_200_Pixel_1_1_1/index.html 310 source trackResolution CMS.tcl301 source trackResolutionFCChh.tcl 311 302 # FIXME !!!! we need to add track resolution of FCC-hh baseline detector !!!!! 312 }313 314 ########################################315 # Time Smearing316 ########################################317 318 module TimeSmearing TimeSmearing {319 set InputArray TrackSmearing/tracks320 set OutputArray tracks321 322 # assume 20 ps resolution for now323 set TimeResolution 20E-12324 }325 326 ##################################327 # Primary vertex reconstruction328 ##################################329 330 331 module VertexFinderDA4D VertexFinderDA4D {332 set InputArray TimeSmearing/tracks333 334 set OutputArray tracks335 set VertexOutputArray vertices336 337 set Verbose 0338 set MinPT 1.0339 340 # in mm341 set VertexSpaceSize 0.5342 343 # in s344 set VertexTimeSize 10E-12345 346 set UseTc 1347 set BetaMax 0.1348 set BetaStop 1.0349 set CoolingFactor 0.8350 set MaxIterations 100351 352 # in mm353 set DzCutOff 40354 set D0CutOff 30355 356 }357 358 ##########################359 # Track pile-up subtractor360 ##########################361 362 module TrackTimingPileUpSubtractor TrackTimingPileUpSubtractor {363 # add InputArray InputArray OutputArray364 365 add InputArray ChargedHadronMomentumSmearing/chargedHadrons366 add InputArray ElectronMomentumSmearing/electrons367 add InputArray MuonMomentumSmearing/muons368 369 set VertexInputArray VertexFinderDA4D/vertices370 # assume perfect pile-up subtraction for tracks with |z| > fZVertexResolution371 # Z vertex resolution in m372 set ZVertexResolution {0.0001}373 303 } 374 304 … … 539 469 } 540 470 541 542 471 ################# 543 472 # Electron filter … … 606 535 } 607 536 537 ######################################## 538 # Time Smearing Neutral MIP 539 ######################################## 540 541 module 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 553 module 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 # 563 module 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 577 module 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 608 module 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 625 module 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 } 608 636 609 637 ################### … … 1057 1085 add Branch ScalarHT/energy ScalarHT ScalarHT 1058 1086 add Branch VertexFinderDA4D/vertices Vertex4D Vertex 1059 } 1060 1087 1088 add Branch HighMassVertexRecover/tracks Track Track 1089 } 1090 -
modules/ModulesLinkDef.h
r364dbe1 r79a7b3e 56 56 #include "modules/PileUpMerger.h" 57 57 #include "modules/JetPileUpSubtractor.h" 58 #include "modules/TrackPileUpSubtractor.h"59 58 #include "modules/TrackTimingPileUpSubtractor.h" 60 59 #include "modules/TaggingParticlesSkimmer.h" … … 114 113 #pragma link C++ class PileUpMerger+; 115 114 #pragma link C++ class JetPileUpSubtractor+; 116 #pragma link C++ class TrackPileUpSubtractor+;117 115 #pragma link C++ class TrackTimingPileUpSubtractor+; 118 116 #pragma link C++ class TaggingParticlesSkimmer+; -
modules/TimeSmearing.cc
r364dbe1 r79a7b3e 53 53 54 54 TimeSmearing::TimeSmearing() : 55 f Formula(0), fItInputArray(0)55 fItInputArray(0), fFormula(0) 56 56 { 57 57 fFormula = new DelphesFormula; 58 58 } 59 59 … … 61 61 62 62 TimeSmearing::~TimeSmearing() 63 { 64 63 { 64 if(fFormula) delete fFormula; 65 65 } 66 66 … … 69 69 void TimeSmearing::Init() 70 70 { 71 // read resolution formula 72 73 fFormula->Compile(GetString("TimeResolution", "1.0")); 71 // read resolution formula 72 fFormula->Compile(GetString("TimeResolution", "0.001")); 74 73 75 74 // import input array … … 77 76 fInputArray = ImportArray(GetString("InputArray", "MuonMomentumSmearing/muons")); 78 77 fItInputArray = fInputArray->MakeIterator(); 79 80 78 // create output array 81 79 … … 95 93 { 96 94 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; 99 99 100 100 … … 114 114 pt = candidateMomentum.Pt(); 115 115 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); 121 120 122 timeResolution = fFormula->Eval(pt, eta, phi, e, d0, dz, ctgTheta); 123 if(fabs(candidate->Position.Eta())<fEtaMax) 121 if(candidate->Charge != 0) 124 122 { 125 123 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); 126 131 } 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 140 140 } 141 141 } -
modules/TimeSmearing.h
r364dbe1 r79a7b3e 45 45 46 46 private: 47 DelphesFormula *fFormula; //!47 DelphesFormula *fFormula; 48 48 Double_t fEtaMax; 49 49 -
modules/TrackPileUpSubtractor.cc
r364dbe1 r79a7b3e 117 117 void TrackPileUpSubtractor::Process() 118 118 { 119 Candidate *candidate, *particle , *particleTest;119 Candidate *candidate, *particle; 120 120 map<TIterator *, TObjArray *>::iterator itInputMap; 121 121 TIterator *iterator; … … 123 123 Double_t z, zvtx = 0; 124 124 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 130 126 // find z position of primary vertex 131 127 132 cout << " ---------- NEW EVENT --------- " << endl;133 134 128 fItVertexInputArray->Reset(); 135 cout << " NUMBER OF VERTICES : " << fVertexInputArray->GetEntriesFast() << endl;136 129 while((candidate = static_cast<Candidate *>(fItVertexInputArray->Next()))) 137 130 { 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 } 158 135 } 159 136 … … 169 146 { 170 147 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); 171 const TLorentzVector &candidateMomentum = candidate->Momentum;148 const TLorentzVector &candidateMomentum = particle->Momentum; 172 149 173 150 eta = candidateMomentum.Eta(); … … 176 153 e = candidateMomentum.E(); 177 154 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(); 183 156 184 157 // apply pile-up subtraction 185 158 // 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;188 159 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) 191 161 { 192 162 candidate->IsRecoPU = 1; 193 //cout << TMath::Abs(z - zvtx) << " " << fFormula->Eval(pt, eta, phi, e) * 1.0e3 << endl;194 163 } 195 164 else -
modules/TrackTimingPileUpSubtractor.cc
r364dbe1 r79a7b3e 75 75 76 76 // read resolution formula in m 77 fFormula->Compile(GetString("ZVertexResolution", "0.001")); 77 fChargedMinSignificance = GetDouble("ChargedMinSignificance", 3); 78 fNeutralMinSignificance = GetDouble("NeutralMinSignificance", 3); 78 79 79 80 fPTMin = GetDouble("PTMin", 0.); … … 128 129 Double_t tempPTSquare = 0; 129 130 Double_t pt, eta, phi, e; 130 Double_t distance = 0;131 Double_t distanceCharged, distanceNeutral = 0; 131 132 132 133 // find z position of primary vertex … … 143 144 tvtx = candidate->Position.T(); 144 145 tvtx_err = candidate->PositionError.T(); 145 //cout << " initial : " << candidate->InitialPosition.T() << " final : " << candidate->Position.T() << endl;146 146 } 147 147 } … … 167 167 z = particle->Position.Z(); 168 168 z_err = particle->PositionError.Z(); 169 t = particle-> Position.T();169 t = particle->InitialPosition.T(); 170 170 t_err = particle->PositionError.T(); 171 171 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); 178 174 179 if(candidate->Charge != 0 && TMath::Abs(z - zvtx) < 0.005 && TMath::Abs(t - tvtx) < 5.0)175 if(candidate->Charge != 0 && distanceCharged < fChargedMinSignificance) 180 176 { 181 177 candidate->IsRecoPU = 1; 182 178 } 179 else if(candidate->Charge == 0 && distanceNeutral < fNeutralMinSignificance) 180 { 181 candidate->IsRecoPU = 1; 182 } 183 183 else 184 184 { -
modules/TrackTimingPileUpSubtractor.h
r364dbe1 r79a7b3e 49 49 DelphesFormula *fFormula; //! 50 50 51 Double_t fChargedMinSignificance; 52 Double_t fNeutralMinSignificance; 53 51 54 Double_t fPTMin; 52 55 -
modules/VertexFinderDA4D.cc
r364dbe1 r79a7b3e 275 275 } 276 276 277 if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Ast");278 279 277 if( fVerbose > 2){cout << "Cool down untill reaching the temperature to finish increasing the number of vertexes" << endl;} 280 278 … … 290 288 delta2 = update(beta, tks, vtx, rho0); 291 289 292 if( fVerbose > 10 ) plot_status(beta, vtx, tks, niter, "Bup");293 290 if (fVerbose > 3) 294 291 { … … 312 309 n_it++; 313 310 314 if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");315 311 } 316 312 … … 320 316 { 321 317 split(beta, vtx, tks); 322 if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Asp");323 318 } 324 319 else … … 355 350 } 356 351 while (delta2 > 0.3*fD2UpdateLim && niter < fMaxIterations); 357 if( fVerbose > 10 ) plot_status(beta, vtx, tks, f, "Dadout");358 352 } 359 353 … … 373 367 } 374 368 while (delta2 > fD2UpdateLim && niter < fMaxIterations); 375 if( fVerbose > 10 ) plot_status(beta, vtx, tks, i_pu, Form("Eprg%d",min_trk));376 369 i_pu++; 377 370 } … … 390 383 n_it++; 391 384 392 if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");393 385 } 394 386 } while( beta < fBetaPurge ); … … 404 396 delta2 = update(beta, tks, vtx, rho0); 405 397 niter++; 406 if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Bup");407 398 } 408 399 while (delta2 > 0.3*fD2UpdateLim && niter < fMaxIterations); … … 425 416 candidate->ClusterIndex = k; 426 417 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); 427 419 candidate->PositionError.SetXYZT(0.0, 0.0, fVertexZSize , fVertexTSize*1E-9*c_light); 428 420 candidate->SumPT2 = 0; … … 483 475 fTrackOutputArray->Add(tks.tt[i]); 484 476 } 485 486 if(fVerbose > 10) plot_status_end(vtx, tks);487 477 488 478 } … … 877 867 { 878 868 continue; 879 // plot_split_crush(zn, tn, vtx, tks, k);880 869 // throw std::invalid_argument( "0 division" ); 881 870 } … … 1042 1031 } 1043 1032 } 1044 1045 1046 // -----------------------------------------------------------------------------1047 // Plot status1048 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 else1079 {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 end1154 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 way1159 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 tracks1177 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 vertices1209 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 splitting1255 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 234 234 std::vector<double> Compute_pk_exp_mBetaE(double beta, vertex_t &vtx, tracks_t &tks, double Z_init); 235 235 236 // Plot status of tracks and Vertices237 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 fitting240 void plot_status_end(vertex_t &vtx, tracks_t &tks);241 242 // Plot Crush243 void plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx);244 245 236 246 237 private:
Note:
See TracChangeset
for help on using the changeset viewer.