Changeset 91ef0b8 in git
- Timestamp:
- Mar 7, 2021, 9:48:26 AM (4 years ago)
- Branches:
- master
- Children:
- bd33fce
- Parents:
- 0d65492 (diff), 363e269 (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. - Files:
-
- 1 added
- 1 deleted
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
Makefile
r0d65492 r91ef0b8 660 660 tmp/external/TrackCovariance/TrkUtil.$(ObjSuf): \ 661 661 external/TrackCovariance/TrkUtil.$(SrcSuf) 662 tmp/external/TrackCovariance/VertexFit.$(ObjSuf): \ 663 external/TrackCovariance/VertexFit.$(SrcSuf) 662 664 tmp/modules/AngularSmearing.$(ObjSuf): \ 663 665 modules/AngularSmearing.$(SrcSuf) \ … … 1168 1170 tmp/external/TrackCovariance/SolTrack.$(ObjSuf) \ 1169 1171 tmp/external/TrackCovariance/TrkUtil.$(ObjSuf) \ 1172 tmp/external/TrackCovariance/VertexFit.$(ObjSuf) \ 1170 1173 tmp/modules/AngularSmearing.$(ObjSuf) \ 1171 1174 tmp/modules/BTagging.$(ObjSuf) \ -
README.md
r0d65492 r91ef0b8 1 [![C ircleCI](https://circleci.com/gh/delphes/delphes.svg?style=shield)](https://circleci.com/gh/delphes/delphes) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3735069.svg)](https://doi.org/10.5281/zenodo.3735069)1 [![CI](https://github.com/delphes/delphes/actions/workflows/ci.yml/badge.svg)](https://github.com/delphes/delphes/actions/workflows/ci.yml) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3735069.svg)](https://doi.org/10.5281/zenodo.3735069) 2 2 3 3 Delphes -
cards/delphes_card_IDEAtrkCov.tcl
r0d65492 r91ef0b8 8 8 # michele.selvaggi@cern.ch 9 9 ##################################################################### 10 # 10 11 set B 2.0 12 11 13 ####################################### 12 14 # Order of execution of various modules … … 59 61 TreeWriter 60 62 } 61 62 63 63 64 ################################# … … 80 81 81 82 # magnetic field, in T 82 set Bz 2.083 set Bz $B 83 84 } 84 85 … … 310 311 } 311 312 312 set Bz 2.0 313 } 314 313 set Bz $B 314 } 315 315 316 316 ############## … … 767 767 add Branch MissingET/momentum MissingET MissingET 768 768 add Branch ScalarHT/energy ScalarHT ScalarHT 769 } 769 770 # add Info InfoName InfoValue 771 add Info Bz $B 772 } -
classes/DelphesClasses.cc
r0d65492 r91ef0b8 112 112 TMatrixDSym Cv; 113 113 Cv.ResizeTo(5, 5); 114 114 115 115 // convert diagonal term to original units 116 Cv(0, 0)=TMath::Power(ErrorD0 *1.0E-3, 2.);116 Cv(0, 0)=TMath::Power(ErrorD0, 2.); 117 117 Cv(1, 1)=TMath::Power(ErrorPhi, 2.); 118 Cv(2, 2)=TMath::Power(ErrorC *1.0E3, 2.);119 Cv(3, 3)=TMath::Power(ErrorDZ *1.0E-3, 2.);118 Cv(2, 2)=TMath::Power(ErrorC, 2.); 119 Cv(3, 3)=TMath::Power(ErrorDZ, 2.); 120 120 Cv(4, 4)=TMath::Power(ErrorCtgTheta, 2.); 121 121 … … 171 171 TMatrixDSym Cv; 172 172 Cv.ResizeTo(5, 5); 173 173 174 174 // convert diagonal term to original units 175 Cv(0, 0)=TMath::Power(ErrorD0 *1.0E-3, 2.);175 Cv(0, 0)=TMath::Power(ErrorD0, 2.); 176 176 Cv(1, 1)=TMath::Power(ErrorPhi, 2.); 177 Cv(2, 2)=TMath::Power(ErrorC *1.0E3, 2.);178 Cv(3, 3)=TMath::Power(ErrorDZ *1.0E-3, 2.);177 Cv(2, 2)=TMath::Power(ErrorC, 2.); 178 Cv(3, 3)=TMath::Power(ErrorDZ, 2.); 179 179 Cv(4, 4)=TMath::Power(ErrorCtgTheta, 2.); 180 180 … … 234 234 NNeutrals(0), 235 235 NeutralEnergyFraction(0), // charged energy fraction 236 ChargedEnergyFraction(0), // neutral energy fraction 236 ChargedEnergyFraction(0), // neutral energy fraction 237 237 Beta(0), 238 238 BetaStar(0), -
classes/DelphesModule.cc
r0d65492 r91ef0b8 130 130 //------------------------------------------------------------------------------ 131 131 132 void DelphesModule::AddInfo(const char *name, Double_t value) 133 { 134 stringstream message; 135 if(!fTreeWriter) 136 { 137 fTreeWriter = static_cast<ExRootTreeWriter *>(GetObject("TreeWriter", ExRootTreeWriter::Class())); 138 if(!fTreeWriter) 139 { 140 message << "can't access access tree writer"; 141 throw runtime_error(message.str()); 142 } 143 } 144 fTreeWriter->AddInfo(name, value); 145 } 146 147 //------------------------------------------------------------------------------ 148 132 149 ExRootResult *DelphesModule::GetPlots() 133 150 { -
classes/DelphesModule.h
r0d65492 r91ef0b8 55 55 56 56 ExRootTreeBranch *NewBranch(const char *name, TClass *cl); 57 void AddInfo(const char *name, Double_t value); 57 58 58 59 ExRootResult *GetPlots(); -
examples/Example6.C
r0d65492 r91ef0b8 64 64 // CEntral track over min Pt 65 65 TH1* histPtgenC = new TH1F("h_PtgenC", "Generated Pt - Central", 500, 0., 50.); // pt for central tracks; 66 TH1* histPtobsC = new TH1F("h_PtobsC", "Reconstructed Pt - Central", 500, 0., 50.); // pt for central 66 TH1* histPtobsC = new TH1F("h_PtobsC", "Reconstructed Pt - Central", 500, 0., 50.); // pt for central 67 67 TH1* histCtgenH = new TH1F("h_CtgenH", "Generateded Cotangent", 100, -10., 10.); // cot(theta) for high pt tracks; 68 68 TH1* histCtobsH = new TH1F("h_CtobsH", "Reconstructed Cotangent", 100, -10., 10.); // cot(theta) for high pt tracks … … 77 77 // 78 78 // Get magnetifc field 79 Double_t Bz = 2.0;79 Double_t Bz = treeReader->GetInfo("Bz"); 80 80 81 81 // Loop over all events … … 107 107 histCtobs->Fill(obsCtg); 108 108 // 109 // Get associated generated particle 109 // Get associated generated particle 110 110 GenParticle* gp = (GenParticle*)trk->Particle.GetObject(); 111 111 // … … 162 162 GenParticle* gpart = (GenParticle*)branchGenPart->At(it); 163 163 // 164 // Plot charged particle parameters 164 // Plot charged particle parameters 165 165 // Only final state particles (Status = 1) 166 166 if (gpart->Status == 1 && TMath::Abs(gpart->Charge) > 0) { … … 290 290 histAccCtgH->Draw(); 291 291 Cnv5->cd(3); gPad->SetLogy(1); 292 histPtgenC->Draw(); 292 histPtgenC->Draw(); 293 293 histPtobsC->SetLineColor(kRed); 294 294 histPtobsC->Draw("SAME"); -
examples/ExamplePVtxFit.C
r0d65492 r91ef0b8 1 1 /* 2 2 Example of using vertex fitter class to fit primary vertex 3 assumed to be generated in (0,0,0). 4 To run compiled version: 5 root 6 root> .L examples/LoadPvtxFit.C+ 7 root> LoadPVtxFit() 8 root> ExamplePVtxFit("infile.root", Nevents) 3 assumed to be generated in (0,0,0) 4 */ 9 5 10 March 4, 2021 11 F. Bedeschi, INFN-Pisa, Italy 12 */ 6 #ifdef __CLING__ 7 R__LOAD_LIBRARY(libDelphes) 13 8 14 9 #include "classes/DelphesClasses.h" … … 18 13 #include "external/TrackCovariance/VertexFit.h" 19 14 20 #include <TClonesArray.h> 21 #include <TChain.h> 22 #include <TVectorD.h> 23 #include <TMatrixDSym.h> 24 #include <TH1.h> 25 #include <TCanvas.h> 26 #include <TStyle.h> 15 #endif 16 27 17 28 18 //------------------------------------------------------------------------------ … … 50 40 // 51 41 // Loop over all events 52 Int_t Nev = TMath::Min(Nevent, (Int_t) 42 Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries); 53 43 for (Int_t entry = 0; entry < Nev; ++entry) 54 44 { … … 56 46 treeReader->ReadEntry(entry); 57 47 Int_t Ntr = 0; // # of tracks from primary vertex 58 Int_t NtrG = branchTrack->GetEntries(); // Total # of tracks48 Int_t NtrG = branchTrack->GetEntries(); 59 49 TVectorD** pr = new TVectorD * [NtrG]; 60 50 TMatrixDSym** cv = new TMatrixDSym * [NtrG]; … … 68 58 Track* trk = (Track*)branchTrack->At(it); 69 59 // 70 // Get associated generated particle 60 // Get associated generated particle 71 61 GenParticle* gp = (GenParticle*)trk->Particle.GetObject(); 72 62 // 73 // Position of origin (mm)63 // Position of origin in mm 74 64 Double_t x = gp->X; 75 65 Double_t y = gp->Y; 76 66 Double_t z = gp->Z; 77 67 // 78 // Select group oftracks originating from the primary vertex68 // group tracks originating from the primary vertex 79 69 if (x == 0.0 && y == 0.0) 80 70 { … … 90 80 // 91 81 pr[Ntr] = new TVectorD(obsPar); 92 cv[Ntr] = new TMatrixDSym(TrkUtil::CovToMm(trk->CovarianceMatrix())); 82 //std::cout<<"Cov Matrix:"<<std::endl; 83 //trk->CovarianceMatrix().Print(); 84 cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix()); 93 85 Ntr++; 94 86 } … … 115 107 } 116 108 109 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl; 117 110 // 118 111 // Cleanup … … 125 118 // Show resulting histograms 126 119 // 127 TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex fitpulls", 50, 50, 900, 500);120 TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500); 128 121 Cnv->Divide(2, 2); 129 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111); 122 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111); 130 123 hXpull->Fit("gaus"); hXpull->Draw(); 131 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111); 124 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111); 132 125 hYpull->Fit("gaus"); hYpull->Draw(); 133 126 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111); -
external/ExRootAnalysis/ExRootTreeReader.cc
r0d65492 r91ef0b8 12 12 #include "TBranchElement.h" 13 13 #include "TCanvas.h" 14 #include "TParameter.h" 14 15 #include "TClonesArray.h" 15 16 #include "TH2.h" … … 119 120 //------------------------------------------------------------------------------ 120 121 122 Double_t ExRootTreeReader::GetInfo(const char *name) 123 { 124 TTree *tree = NULL; 125 TParameter<Double_t> *param = NULL; 126 Double_t result = -999.9; 127 if(fChain) tree = fChain->GetTree(); 128 if(tree) param = static_cast<TParameter<Double_t> *>(tree->GetUserInfo()->FindObject(name)); 129 if(param) result = param->GetVal(); 130 return result; 131 } 132 133 //------------------------------------------------------------------------------ 134 121 135 Bool_t ExRootTreeReader::Notify() 122 136 { -
external/ExRootAnalysis/ExRootTreeReader.h
r0d65492 r91ef0b8 29 29 30 30 TClonesArray *UseBranch(const char *branchName); 31 Double_t GetInfo(const char *name); 31 32 32 33 private: -
external/ExRootAnalysis/ExRootTreeWriter.cc
r0d65492 r91ef0b8 11 11 #include "ExRootAnalysis/ExRootTreeBranch.h" 12 12 13 #include "TParameter.h" 13 14 #include "TClonesArray.h" 14 15 #include "TFile.h" … … 48 49 fBranches.insert(branch); 49 50 return branch; 51 } 52 53 //------------------------------------------------------------------------------ 54 55 void ExRootTreeWriter::AddInfo(const char *name, Double_t value) 56 { 57 if(!fTree) fTree = NewTree(); 58 fTree->GetUserInfo()->Add(new TParameter<Double_t>(name, value)); 50 59 } 51 60 -
external/ExRootAnalysis/ExRootTreeWriter.h
r0d65492 r91ef0b8 29 29 30 30 ExRootTreeBranch *NewBranch(const char *name, TClass *cl); 31 void AddInfo(const char *name, Double_t value); 31 32 32 33 void Clear(); -
external/TrackCovariance/TrkUtil.cc
r0d65492 r91ef0b8 1 1 #include "TrkUtil.h" 2 #include <TMath.h> 2 3 #include <iostream> 3 4 … … 31 32 Double_t r2 = x.Perp2(); 32 33 Double_t cross = x(0) * p(1) - x(1) * p(0); 33 Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2);34 Double_t phi0 = atan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi034 Double_t T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2); 35 Double_t phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi0 35 36 Double_t D; // Impact parameter D 36 37 if (pt < 10.0) D = (T - pt) / a; … … 41 42 Par(2) = C; // Store C 42 43 //Longitudinal parameters 43 Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));44 Double_t st = asin(B) / C;44 Double_t B = C * TMath::Sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D)); 45 Double_t st = TMath::ASin(B) / C; 45 46 Double_t ct = p(2) / pt; 46 47 Double_t z0 = x(2) - ct * st; … … 69 70 // 70 71 TVector3 Xval; 71 Xval(0) = -D * sin(phi0);72 Xval(1) = D * cos(phi0);72 Xval(0) = -D * TMath::Sin(phi0); 73 Xval(1) = D * TMath::Cos(phi0); 73 74 Xval(2) = z0; 74 75 // … … 92 93 TVector3 Pval; 93 94 Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C); 94 Pval(0) = pt * cos(phi0);95 Pval(1) = pt * sin(phi0);95 Pval(0) = pt * TMath::Cos(phi0); 96 Pval(1) = pt * TMath::Sin(phi0); 96 97 Pval(2) = pt * ct; 97 98 // … … 114 115 pACTS(1) = 1000 * Par(3); // z0 from m to mm 115 116 pACTS(2) = Par(1); // Phi0 is unchanged 116 pACTS(3) = atan2(1.0, Par(4)); // Theta in [0, pi] range117 pACTS(4) = Par(2) / (b * sqrt(1 + Par(4) * Par(4))); // q/p in GeV117 pACTS(3) = TMath::ATan2(1.0, Par(4)); // Theta in [0, pi] range 118 pACTS(4) = Par(2) / (b * TMath::Sqrt(1 + Par(4) * Par(4))); // q/p in GeV 118 119 pACTS(5) = 0.0; // Time: currently undefined 119 120 // … … 132 133 A(0, 0) = 1000.; // D-D conversion to mm 133 134 A(1, 2) = 1.0; // phi0-phi0 134 A(2, 4) = 1.0 / ( sqrt(1.0 + ct * ct) * b); // q/p-C135 A(2, 4) = 1.0 / (TMath::Sqrt(1.0 + ct * ct) * b); // q/p-C 135 136 A(3, 1) = 1000.; // z0-z0 conversion to mm 136 137 A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta) 137 A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta)138 A(4, 4) = -C * ct / (b * TMath::Power(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta) 138 139 // 139 140 TMatrixDSym Cv = Cov; -
external/TrackCovariance/VertexFit.cc
r0d65492 r91ef0b8 15 15 fVtxCst = kFALSE; 16 16 fxCst.ResizeTo(3); 17 fCovCst.ResizeTo(3, 3);17 fCovCst.ResizeTo(3, 3); 18 18 fXv.ResizeTo(3); 19 fcovXv.ResizeTo(3, 3);19 fcovXv.ResizeTo(3, 3); 20 20 } 21 21 // Parameters and covariances … … 26 26 fVtxCst = kFALSE; 27 27 fxCst.ResizeTo(3); 28 fCovCst.ResizeTo(3, 3);28 fCovCst.ResizeTo(3, 3); 29 29 fXv.ResizeTo(3); 30 fcovXv.ResizeTo(3, 3);30 fcovXv.ResizeTo(3, 3); 31 31 // 32 32 fPar = trkPar; … … 35 35 // 36 36 ffi = new Double_t[Ntr]; // Fit phases 37 fx0i = new TVectorD * [Ntr]; // Track expansion points37 fx0i = new TVectorD * [Ntr]; // Track expansion points 38 38 for (Int_t i = 0; i < Ntr; i++) fx0i[i] = new TVectorD(3); 39 39 fai = new TVectorD * [Ntr]; // dx/dphi 40 40 for (Int_t i = 0; i < Ntr; i++) fai[i] = new TVectorD(3); 41 41 fa2i = new Double_t[Ntr]; // a'Wa 42 fDi = new TMatrixDSym *[Ntr]; // W-WBW42 fDi = new TMatrixDSym * [Ntr]; // W-WBW 43 43 for (Int_t i = 0; i < Ntr; i++) fDi[i] = new TMatrixDSym(3); 44 44 fWi = new TMatrixDSym * [Ntr]; // (ACA')^-1 … … 46 46 fWinvi = new TMatrixDSym * [Ntr]; // ACA' 47 47 for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3); 48 49 48 } 50 49 // ObsTrk list … … 55 54 fVtxCst = kFALSE; 56 55 fxCst.ResizeTo(3); 57 fCovCst.ResizeTo(3, 3);56 fCovCst.ResizeTo(3, 3); 58 57 fXv.ResizeTo(3); 59 fcovXv.ResizeTo(3, 3);60 // 61 fPar = new TVectorD *[Ntr];62 fCov = new TMatrixDSym *[Ntr];58 fcovXv.ResizeTo(3, 3); 59 // 60 fPar = new TVectorD * [Ntr]; 61 fCov = new TMatrixDSym * [Ntr]; 63 62 fChi2List.ResizeTo(Ntr); 64 63 for (Int_t i = 0; i < Ntr; i++) 65 64 { 66 65 fPar[i] = new TVectorD(track[i]->GetObsPar()); 67 //cout << "fPar = " << endl; fPar[i]->Print();68 66 fCov[i] = new TMatrixDSym(track[i]->GetCov()); 69 //cout << "fCov = " << endl; fCov[i]->Print();70 67 } 71 68 // … … 82 79 fWinvi = new TMatrixDSym * [Ntr]; // ACA' 83 80 for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3); 84 85 //cout << "VertexFit constructor executed" << endl;86 81 } 87 82 // … … 89 84 VertexFit::~VertexFit() 90 85 { 91 fxCst.Clear(); 86 fxCst.Clear(); 92 87 fCovCst.Clear(); 93 fXv.Clear(); 94 fcovXv.Clear(); 88 fXv.Clear(); 89 fcovXv.Clear(); 95 90 fChi2List.Clear(); 96 91 // 97 for (Int_t i = 0; i < fNtr; i++)98 { 99 fPar[i]->Clear(); 92 for (Int_t i = 0; i < fNtr; i++) 93 { 94 fPar[i]->Clear(); 100 95 fCov[i]->Clear(); 101 96 // … … 148 143 H(1, 1) = nn; 149 144 TVectorD c(2); 150 c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += 145 c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += n(i) * (y0(i) - x0(i)); 151 146 c(1) = 0; for (Int_t i = 0; i < 3; i++)c(1) += -k(i) * (y0(i) - x0(i)); 152 147 TVectorD smin = (H * c); … … 158 153 Double_t R2 = TMath::Sqrt(Y(0) * Y(0) + Y(1) * Y(1)); 159 154 // 160 return 0.5 *(R1+R2);155 return 0.5 * (R1 + R2); 161 156 } 162 157 Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2) … … 170 165 TMatrixDSym H(2); 171 166 H(0, 0) = -TMath::Cos(p2(1)); 172 H(0, 1) = 167 H(0, 1) = TMath::Cos(p1(1)); 173 168 H(1, 0) = -TMath::Sin(p2(1)); 174 H(1, 1) = 169 H(1, 1) = TMath::Sin(p1(1)); 175 170 Double_t Det = TMath::Sin(p2(1) - p1(1)); 176 171 H *= 1.0 / Det; … … 195 190 if (Ntry > NtryMax) 196 191 { 197 cout << "FastRv: maximum number of iteration reached" <<endl;192 std::cout << "FastRv: maximum number of iteration reached" << std::endl; 198 193 break; 199 194 } … … 208 203 } 209 204 210 TMatrixDSym VertexFit::RegInv3(TMatrixDSym &Smat0)205 TMatrixDSym VertexFit::RegInv3(TMatrixDSym& Smat0) 211 206 { 212 207 // … … 217 212 if (N != 3) 218 213 { 219 cout << "RegInv3 called with matrix size != 3. Abort & return standard inversion." <<endl;214 std::cout << "RegInv3 called with matrix size != 3. Abort & return standard inversion." << std::endl; 220 215 return Smat.Invert(); 221 216 } … … 232 227 for (Int_t j = 0; j < 2; j++)Q(i, j) = RegMat(i, j); 233 228 } 234 Double_t Det = 1 - Q(0, 1) *Q(1, 0);229 Double_t Det = 1 - Q(0, 1) * Q(1, 0); 235 230 TMatrixDSym H(2); 236 231 H = Q; … … 241 236 p(1) = RegMat(1, 2); 242 237 Double_t pHp = H.Similarity(p); 243 Double_t h = pHp -Det;238 Double_t h = pHp - Det; 244 239 // 245 240 TMatrixDSym pp(2); pp.Rank1Update(p); 246 TMatrixDSym F = (h *H) - pp.Similarity(H);241 TMatrixDSym F = (h * H) - pp.Similarity(H); 247 242 F *= 1.0 / Det; 248 TVectorD b = H *p;243 TVectorD b = H * p; 249 244 TMatrixDSym InvReg(3); 250 245 for (Int_t i = 0; i < 2; i++) … … 263 258 else 264 259 { 265 cout << "RegInv3: found negative elements in diagonal. Return standard inversion." << endl; 266 return Smat.Invert(); 260 D.Zero(); 261 for (Int_t i = 0; i < N; i++) D(i, i) = 1.0 / TMath::Sqrt(TMath::Abs(Smat(i, i))); 262 TMatrixDSym RegMat = Smat.Similarity(D); 263 RegMat.Invert(); 264 return RegMat.Similarity(D); 267 265 } 268 266 } … … 294 292 A(2, 0) = 0.0; 295 293 // phi0 296 A(0, 1) = -D *TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);297 A(1, 1) = -D *TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);294 A(0, 1) = -D * TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C); 295 A(1, 1) = -D * TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C); 298 296 A(2, 1) = 0.0; 299 297 // C 300 A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C *C);301 A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C *C);302 A(2, 2) = -ct *phi / (2 * C*C);298 A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C * C); 299 A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C * C); 300 A(2, 2) = -ct * phi / (2 * C * C); 303 301 // z0 304 302 A(0, 3) = 0.0; … … 353 351 Double_t ct = par(4); 354 352 // 355 x0(0) = -D * TMath::Sin(p0);356 x0(1) = D *TMath::Cos(p0);353 x0(0) = -D * TMath::Sin(p0); 354 x0(1) = D * TMath::Cos(p0); 357 355 x0(2) = z0; 358 356 // … … 378 376 x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C); 379 377 x(1) = x0(1) - (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C); 380 x(2) = x0(2) + ct *phi / (2 * C);378 x(2) = x0(2) + ct * phi / (2 * C); 381 379 // 382 380 return x; … … 395 393 // 396 394 // Stored quantities 397 Double_t *fi = new Double_t[fNtr]; // Phases398 TVectorD **x0i = new TVectorD*[fNtr]; // Track expansion point399 TVectorD **ai = new TVectorD*[fNtr]; // dx/dphi400 Double_t *a2i = new Double_t[fNtr]; // a'Wa401 TMatrixDSym **Di = new TMatrixDSym*[fNtr]; // W-WBW402 TMatrixDSym **Wi = new TMatrixDSym*[fNtr]; // (ACA')^-1403 TMatrixDSym **Winvi = new TMatrixDSym*[fNtr]; // ACA'395 Double_t* fi = new Double_t[fNtr]; // Phases 396 TVectorD** x0i = new TVectorD * [fNtr]; // Track expansion point 397 TVectorD** ai = new TVectorD * [fNtr]; // dx/dphi 398 Double_t* a2i = new Double_t[fNtr]; // a'Wa 399 TMatrixDSym** Di = new TMatrixDSym * [fNtr]; // W-WBW 400 TMatrixDSym** Wi = new TMatrixDSym * [fNtr]; // (ACA')^-1 401 TMatrixDSym** Winvi = new TMatrixDSym * [fNtr]; // ACA' 404 402 // 405 403 // vertex radius approximation … … 415 413 // 416 414 // Find track pair with largest phi difference 417 Int_t isel=0; Int_t jsel=0; // selected track indices 418 Double_t dphiMax = 0.0; // Max phi difference 419 for (Int_t i = 0; i < fNtr-1; i++) 415 Int_t isel = 0; Int_t jsel = 0; // selected track indices 416 Double_t dphiMax = 0.0; // Max phi difference 417 418 for (Int_t i = 0; i < fNtr - 1; i++) 420 419 { 421 420 //ObsTrk* ti = tracks[i]; … … 423 422 Double_t phi1 = pari(1); 424 423 425 for (Int_t j = i +1; j < fNtr; j++)424 for (Int_t j = i + 1; j < fNtr; j++) 426 425 { 427 426 //ObsTrk* tj = tracks[j]; … … 438 437 } 439 438 // 440 //441 //ObsTrk* t1 = tracks[isel];442 439 TVectorD p1 = *fPar[isel]; 443 //ObsTrk* t2 = tracks[jsel];444 440 TVectorD p2 = *fPar[jsel]; 445 441 Double_t R = FastRv1(p1, p2); 446 if (R > 1.0) R = Rd; 447 R = 0.9*R + 0.1*Rd; 448 // 449 //cout << "VertexFinder: fast R = " << R << endl; 442 if (R > 1000.0) R = Rd; 443 R = 0.9 * R + 0.1 * Rd; 450 444 // 451 445 // Iteration properties … … 464 458 H.Zero(); // Reset H matrix 465 459 DW1D.Zero(); 466 // 467 // cout << "VertexFinder: start loop on tracks" <<endl;460 // 461 //std::cout << "VertexFinder: start loop on tracks" << std::endl; 468 462 for (Int_t i = 0; i < fNtr; i++) 469 463 { 470 // Get track helix parameters and their covariance matrix 471 //ObsTrk *t = tracks[i]; 464 // Get track helix parameters and their covariance matrix 472 465 TVectorD par = *fPar[i]; 473 TMatrixDSym Cov = *fCov[i]; 466 TMatrixDSym Cov = *fCov[i]; 467 474 468 Double_t fs; 475 469 if (Ntry <= 0) // Initialize all phases on first pass … … 477 471 Double_t D = par(0); 478 472 Double_t C = par(2); 479 Double_t arg = TMath::Max(1.0e-6, (R *R - D*D) / (1 + 2 * C*D));480 fs = 2 * TMath::ASin(C *TMath::Sqrt(arg));473 Double_t arg = TMath::Max(1.0e-6, (R * R - D * D) / (1 + 2 * C * D)); 474 fs = 2 * TMath::ASin(C * TMath::Sqrt(arg)); 481 475 fi[i] = fs; 482 476 } … … 484 478 // Starting values 485 479 // 486 fs = fi[i]; 487 // cout << "VertexFinder: phase fs set" <<endl;480 fs = fi[i]; // Get phase 481 //std::cout << "VertexFinder: phase fs set" << std::endl; 488 482 TVectorD xs = Fill_x(par, fs); 489 // cout << "VertexFinder: position xs set" <<endl;483 //std::cout << "VertexFinder: position xs set" << std::endl; 490 484 x0i[i] = new TVectorD(xs); // Start helix position 491 // cout << "VertexFinder: position x0i stored" <<endl;485 //std::cout << "VertexFinder: position x0i stored" << std::endl; 492 486 // W matrix = (A*C*A')^-1; W^-1 = A*C*A' 493 487 TMatrixD A = Fill_A(par, fs); // A = dx/da = derivatives wrt track parameters 494 // cout << "VertexFinder: derivatives A set" <<endl;488 //std::cout << "VertexFinder: derivatives A set" << std::endl; 495 489 TMatrixDSym Winv = Cov.Similarity(A); // W^-1 = A*C*A' 496 490 Winvi[i] = new TMatrixDSym(Winv); // Store W^-1 matrix 497 // cout << "VertexFinder: Winvi stored" <<endl;491 //std::cout << "VertexFinder: Winvi stored" << std::endl; 498 492 TMatrixDSym W = RegInv3(Winv); // W = (A*C*A')^-1 499 493 Wi[i] = new TMatrixDSym(W); // Store W matrix 500 // cout << "VertexFinder: Wi stored" <<endl;494 //std::cout << "VertexFinder: Wi stored" << std::endl; 501 495 TVectorD a = Fill_a(par, fs); // a = dx/ds = derivatives wrt phase 502 // cout << "VertexFinder: derivatives a set" <<endl;496 //std::cout << "VertexFinder: derivatives a set" << std::endl; 503 497 ai[i] = new TVectorD(a); // Store a 504 // cout << "VertexFinder: derivatives a stored" <<endl;498 //std::cout << "VertexFinder: derivatives a stored" << std::endl; 505 499 Double_t a2 = W.Similarity(a); 506 500 a2i[i] = a2; // Store a2 507 501 // Build D matrix 508 TMatrixDSym B(3); 502 TMatrixDSym B(3); 503 509 504 B.Rank1Update(a, 1.0); 510 505 B *= -1. / a2; 511 506 B.Similarity(W); 512 TMatrixDSym Ds = W +B; // D matrix507 TMatrixDSym Ds = W + B; // D matrix 513 508 Di[i] = new TMatrixDSym(Ds); // Store D matrix 514 // cout << "VertexFinder: matrix Di stored" <<endl;509 //std::cout << "VertexFinder: matrix Di stored" << std::endl; 515 510 TMatrixDSym DsW1Ds = Winv.Similarity(Ds); // Service matrix to calculate covX 516 DW1D += DsW1Ds; 511 DW1D += DsW1Ds; 517 512 // Update hessian 518 513 H += Ds; … … 523 518 // update vertex position 524 519 TMatrixDSym H1 = RegInv3(H); 525 x = H1 *cterm;526 // cout << "VertexFinder: x vertex set" <<endl;520 x = H1 * cterm; 521 //std::cout << "VertexFinder: x vertex set" << std::endl; 527 522 // Update vertex covariance 528 523 covX = DW1D.Similarity(H1); 529 // cout << "VertexFinder: cov vertex set" <<endl;524 //std::cout << "VertexFinder: cov vertex set" << std::endl; 530 525 // Update phases and chi^2 531 526 Chi2 = 0.0; 532 527 for (Int_t i = 0; i < fNtr; i++) 533 528 { 534 TVectorD lambda = (*Di[i]) *(*x0i[i] - x);529 TVectorD lambda = (*Di[i]) * (*x0i[i] - x); 535 530 TMatrixDSym Wm1 = *Winvi[i]; 536 531 fChi2List(i) = Wm1.Similarity(lambda); 537 532 Chi2 += fChi2List(i); 538 533 TVectorD a = *ai[i]; 539 TVectorD b = (*Wi[i]) *(x - *x0i[i]);540 for (Int_t j = 0; j < 3; j++)fi[i] += a(j) *b(j) / a2i[i];534 TVectorD b = (*Wi[i]) * (x - *x0i[i]); 535 for (Int_t j = 0; j < 3; j++)fi[i] += a(j) * b(j) / a2i[i]; 541 536 } 542 537 543 //cout << "VertexFinder: end chi2 calculation" << endl;544 538 // 545 539 TVectorD dx = x - x0; … … 559 553 // 560 554 561 // cout << "VertexFinder: before store intermediate data" <<endl;555 //std::cout << "VertexFinder: before store intermediate data" << std::endl; 562 556 for (Int_t i = 0; i < fNtr; i++) 563 557 { 564 // cout << "VertexFinder: inside store intermediate data" <<endl;565 // cout << "i = " << i << ", fi[i] = " << fi[i] <<endl;566 // cout << "i = " << i << ", ffi[i] = " << ffi[i] <<endl;558 //std::cout << "VertexFinder: inside store intermediate data" << std::endl; 559 //std::cout << "i = " << i << ", fi[i] = " << fi[i] << std::endl; 560 //std::cout << "i = " << i << ", ffi[i] = " << ffi[i] << std::endl; 567 561 ffi[i] = fi[i]; // Fit phases 568 // cout << "VertexFinder: fi stored" <<endl;562 //std::cout << "VertexFinder: fi stored" << std::endl; 569 563 fx0i[i] = x0i[i]; // Track expansion points 570 // cout << "VertexFinder: x0i stored" <<endl;564 //std::cout << "VertexFinder: x0i stored" << std::endl; 571 565 fai[i] = ai[i]; // dx/dphi 572 // cout << "VertexFinder: ai stored" <<endl;566 //std::cout << "VertexFinder: ai stored" << std::endl; 573 567 fa2i[i] = a2i[i]; // a'Wa 574 // cout << "VertexFinder: a2i stored" <<endl;568 //std::cout << "VertexFinder: a2i stored" << std::endl; 575 569 fDi[i] = Di[i]; // W-WBW 576 // cout << "VertexFinder: Di stored" <<endl;570 //std::cout << "VertexFinder: Di stored" << std::endl; 577 571 fWi[i] = Wi[i]; // (ACA')^-1 578 // cout << "VertexFinder: Wi stored" <<endl;572 //std::cout << "VertexFinder: Wi stored" << std::endl; 579 573 fWinvi[i] = Winvi[i]; // ACA' 580 // cout << "VertexFinder: Winvi stored" <<endl;574 //std::cout << "VertexFinder: Winvi stored" << std::endl; 581 575 } 582 // cout << "Iteration " << Ntry << " completed - Before cleanup" <<endl;576 //std::cout << "Iteration " << Ntry << " completed - Before cleanup" << std::endl; 583 577 // 584 578 // Cleanup … … 599 593 } 600 594 601 // cout << "Iteration " << Ntry << " completed - After cleanup" <<endl;595 //std::cout << "Iteration " << Ntry << " completed - After cleanup" << std::endl; 602 596 } 603 597 // 604 598 fVtxDone = kTRUE; // Set fitting completion flag 605 599 // 606 delete[] fi; // Phases 600 delete[] fi; // Phases 607 601 delete[] x0i; // Track expansion point 608 602 delete[] ai; // dx/dphi … … 615 609 TVectorD VertexFit::GetVtx() 616 610 { 617 // cout << "GetVtx: flag set to " << fVtxDone <<endl;611 //std::cout << "GetVtx: flag set to " << fVtxDone << std::endl; 618 612 if (!fVtxDone)VertexFinder(); 619 613 return fXv; … … 641 635 void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov) // Add gaussian vertex constraint 642 636 { 643 cout << "VertexFit::AddVtxConstraint: Not implemented yet" <<endl;637 std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl; 644 638 } 645 639 // 646 640 void VertexFit::AddTrk(TVectorD par, TMatrixDSym Cov) // Add track to input list 647 641 { 648 cout << "VertexFit::AddTrk: Not implemented yet" <<endl;642 std::cout << "VertexFit::AddTrk: Not implemented yet" << std::endl; 649 643 } 650 644 void VertexFit::RemoveTrk(Int_t iTrk) // Remove iTrk track 651 645 { 652 cout << "VertexFit::RemoveTrk: Not implemented yet" <<endl;653 } 646 std::cout << "VertexFit::RemoveTrk: Not implemented yet" << std::endl; 647 } -
modules/TreeWriter.cc
r0d65492 r91ef0b8 123 123 fBranchMap.insert(make_pair(branch, make_pair(itClassMap->second, array))); 124 124 } 125 126 param = GetParam("Info"); 127 TString infoName; 128 Double_t infoValue; 129 130 size = param.GetSize(); 131 for(i = 0; i < size / 2; ++i) 132 { 133 infoName = param[i * 2].GetString(); 134 infoValue = param[i * 2 + 1].GetDouble(); 135 136 AddInfo(infoName, infoValue); 137 } 125 138 } 126 139 … … 138 151 it1.Reset(); 139 152 array->Clear(); 140 153 141 154 while((candidate = static_cast<Candidate *>(it1.Next()))) 142 155 { … … 537 550 entry->CtgTheta = ctgTheta; 538 551 entry->C = candidate->C; 539 552 540 553 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); 541 554 const TLorentzVector &initialPosition = particle->Position;
Note:
See TracChangeset
for help on using the changeset viewer.