Changes in / [4564cba:54ec182] in git
- Files:
-
- 11 deleted
- 17 edited
-
Makefile (modified) (10 diffs)
-
cards/delphes_card_CMS_CSCCluster.tcl (deleted)
-
classes/ClassesLinkDef.h (modified) (2 diffs)
-
classes/DelphesClasses.cc (modified) (8 diffs)
-
classes/DelphesClasses.h (modified) (6 diffs)
-
classes/DelphesCscClusterFormula.cc (deleted)
-
classes/DelphesCscClusterFormula.h (deleted)
-
classes/DelphesHepMC2Reader.cc (modified) (3 diffs)
-
classes/DelphesHepMC3Reader.cc (modified) (3 diffs)
-
examples/ExampleCscCluster.py (deleted)
-
examples/ExampleKsVtxFit.C (deleted)
-
examples/ExamplePVtxFind.C (modified) (16 diffs)
-
external/TrackCovariance/ObsTrk.cc (modified) (5 diffs)
-
external/TrackCovariance/ObsTrk.h (modified) (2 diffs)
-
external/TrackCovariance/SolTrack.cc (modified) (10 diffs)
-
external/TrackCovariance/SolTrack.h (modified) (1 diff)
-
external/TrackCovariance/VertexFit.cc (modified) (8 diffs)
-
external/TrackCovariance/VertexFit.h (modified) (2 diffs)
-
modules/CscClusterEfficiency.cc (deleted)
-
modules/CscClusterEfficiency.h (deleted)
-
modules/CscClusterId.cc (deleted)
-
modules/CscClusterId.h (deleted)
-
modules/LLPFilter.cc (deleted)
-
modules/LLPFilter.h (deleted)
-
modules/ModulesLinkDef.h (modified) (2 diffs)
-
modules/TrackCovariance.cc (modified) (3 diffs)
-
modules/TreeWriter.cc (modified) (10 diffs)
-
modules/TreeWriter.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Makefile
r4564cba r54ec182 464 464 modules/ParticleDensity.h \ 465 465 modules/TruthVertexFinder.h \ 466 modules/ExampleModule.h \ 467 modules/LLPFilter.h \ 468 modules/CscClusterEfficiency.h \ 469 modules/CscClusterId.h 466 modules/ExampleModule.h 470 467 tmp/modules/ModulesDict$(PcmSuf): \ 471 468 tmp/modules/ModulesDict.$(SrcSuf) … … 521 518 classes/DelphesFactory.h \ 522 519 classes/SortableObject.h 523 tmp/classes/DelphesCscClusterFormula.$(ObjSuf): \524 classes/DelphesCscClusterFormula.$(SrcSuf) \525 classes/DelphesCscClusterFormula.h \526 classes/DelphesClasses.h527 520 tmp/classes/DelphesCylindricalFormula.$(ObjSuf): \ 528 521 classes/DelphesCylindricalFormula.$(SrcSuf) \ … … 749 742 external/ExRootAnalysis/ExRootFilter.h \ 750 743 external/ExRootAnalysis/ExRootResult.h 751 tmp/modules/CscClusterEfficiency.$(ObjSuf): \752 modules/CscClusterEfficiency.$(SrcSuf) \753 modules/CscClusterEfficiency.h \754 classes/DelphesClasses.h \755 classes/DelphesFactory.h \756 classes/DelphesCscClusterFormula.h \757 external/ExRootAnalysis/ExRootClassifier.h \758 external/ExRootAnalysis/ExRootFilter.h \759 external/ExRootAnalysis/ExRootResult.h760 tmp/modules/CscClusterId.$(ObjSuf): \761 modules/CscClusterId.$(SrcSuf) \762 modules/CscClusterId.h \763 classes/DelphesClasses.h \764 classes/DelphesFactory.h \765 classes/DelphesCscClusterFormula.h \766 external/ExRootAnalysis/ExRootClassifier.h \767 external/ExRootAnalysis/ExRootFilter.h \768 external/ExRootAnalysis/ExRootResult.h769 744 tmp/modules/DecayFilter.$(ObjSuf): \ 770 745 modules/DecayFilter.$(SrcSuf) \ … … 907 882 external/ExRootAnalysis/ExRootFilter.h \ 908 883 external/ExRootAnalysis/ExRootResult.h 909 tmp/modules/LLPFilter.$(ObjSuf): \910 modules/LLPFilter.$(SrcSuf) \911 modules/LLPFilter.h \912 classes/DelphesClasses.h \913 classes/DelphesFactory.h \914 classes/DelphesFormula.h \915 external/ExRootAnalysis/ExRootClassifier.h \916 external/ExRootAnalysis/ExRootFilter.h \917 external/ExRootAnalysis/ExRootResult.h918 884 tmp/modules/LeptonDressing.$(ObjSuf): \ 919 885 modules/LeptonDressing.$(SrcSuf) \ … … 1197 1163 DELPHES_OBJ += \ 1198 1164 tmp/classes/DelphesClasses.$(ObjSuf) \ 1199 tmp/classes/DelphesCscClusterFormula.$(ObjSuf) \1200 1165 tmp/classes/DelphesCylindricalFormula.$(ObjSuf) \ 1201 1166 tmp/classes/DelphesFactory.$(ObjSuf) \ … … 1262 1227 tmp/modules/ClusterCounting.$(ObjSuf) \ 1263 1228 tmp/modules/ConstituentFilter.$(ObjSuf) \ 1264 tmp/modules/CscClusterEfficiency.$(ObjSuf) \1265 tmp/modules/CscClusterId.$(ObjSuf) \1266 1229 tmp/modules/DecayFilter.$(ObjSuf) \ 1267 1230 tmp/modules/Delphes.$(ObjSuf) \ … … 1279 1242 tmp/modules/JetFlavorAssociation.$(ObjSuf) \ 1280 1243 tmp/modules/JetPileUpSubtractor.$(ObjSuf) \ 1281 tmp/modules/LLPFilter.$(ObjSuf) \1282 1244 tmp/modules/LeptonDressing.$(ObjSuf) \ 1283 1245 tmp/modules/Merger.$(ObjSuf) \ … … 1966 1928 @touch $@ 1967 1929 1968 modules/LLPFilter.h: \1969 classes/DelphesModule.h1970 @touch $@1971 1972 1930 external/fastjet/internal/MinHeap.hh: \ 1973 1931 external/fastjet/internal/base.hh … … 2299 2257 @touch $@ 2300 2258 2301 modules/CscClusterEfficiency.h: \2302 classes/DelphesModule.h2303 @touch $@2304 2305 2259 external/fastjet/PseudoJetStructureBase.hh: \ 2306 2260 external/fastjet/internal/base.hh … … 2427 2381 external/fastjet/config.h: \ 2428 2382 external/fastjet/config_win.h 2429 @touch $@2430 2431 modules/CscClusterId.h: \2432 classes/DelphesModule.h2433 2383 @touch $@ 2434 2384 -
classes/ClassesLinkDef.h
r4564cba r54ec182 57 57 #pragma link C++ class Electron+; 58 58 #pragma link C++ class Muon+; 59 #pragma link C++ class CscCluster+;60 61 59 #pragma link C++ class Jet+; 62 60 #pragma link C++ class Track+; … … 68 66 69 67 #endif 68 -
classes/DelphesClasses.cc
r4564cba r54ec182 35 35 CompBase *Electron::fgCompare = CompPT<Electron>::Instance(); 36 36 CompBase *Muon::fgCompare = CompPT<Muon>::Instance(); 37 38 37 CompBase *Jet::fgCompare = CompPT<Jet>::Instance(); 39 38 CompBase *Track::fgCompare = CompPT<Track>::Instance(); … … 43 42 CompBase *Vertex::fgCompare = CompSumPT2<Vertex>::Instance(); 44 43 CompBase *Candidate::fgCompare = CompMomentumPt<Candidate>::Instance(); 45 CompBase *CscCluster::fgCompare = CompE<CscCluster>::Instance();46 44 47 45 //------------------------------------------------------------------------------ … … 220 218 Position(0.0, 0.0, 0.0, 0.0), 221 219 InitialPosition(0.0, 0.0, 0.0, 0.0), 222 DecayPosition(0.0, 0.0, 0.0, 0.0),223 220 PositionError(0.0, 0.0, 0.0, 0.0), 224 221 Area(0.0, 0.0, 0.0, 0.0), … … 233 230 Phi(0), ErrorPhi(0), 234 231 Xd(0), Yd(0), Zd(0), 235 XFirstHit(0), YFirstHit(0), ZFirstHit(0),236 232 Nclusters(0.0), 237 233 dNdx(0.0), … … 392 388 object.Position = Position; 393 389 object.InitialPosition = InitialPosition; 394 object.DecayPosition = DecayPosition;395 390 object.PositionError = PositionError; 396 391 object.Area = Area; … … 414 409 object.Yd = Yd; 415 410 object.Zd = Zd; 416 object.XFirstHit = XFirstHit;417 object.YFirstHit = YFirstHit;418 object.ZFirstHit = ZFirstHit;419 411 object.Nclusters = Nclusters; 420 412 object.dNdx = dNdx; … … 531 523 Position.SetXYZT(0.0, 0.0, 0.0, 0.0); 532 524 InitialPosition.SetXYZT(0.0, 0.0, 0.0, 0.0); 533 DecayPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);534 525 Area.SetXYZT(0.0, 0.0, 0.0, 0.0); 535 526 TrackCovariance.Zero(); … … 553 544 Yd = 0.0; 554 545 Zd = 0.0; 555 XFirstHit = 0.0;556 YFirstHit = 0.0;557 ZFirstHit = 0.0;558 546 Nclusters = 0.0; 559 547 dNdx = 0.0; -
classes/DelphesClasses.h
r4564cba r54ec182 158 158 Float_t Z; // particle vertex position (z component) | hepevt.vhep[number][2] 159 159 160 Float_t decayX;161 Float_t decayY;162 Float_t decayZ;163 Float_t decayT;164 165 160 static CompBase *fgCompare; //! 166 161 const CompBase *GetCompare() const { return fgCompare; } … … 462 457 Float_t Yd; // Y coordinate of point of closest approach to vertex 463 458 Float_t Zd; // Z coordinate of point of closest approach to vertex 464 465 Float_t XFirstHit; // X coordinate of point of closest approach to vertex466 Float_t YFirstHit; // Y coordinate of point of closest approach to vertex467 Float_t ZFirstHit; // Z coordinate of point of closest approach to vertex468 459 469 460 Float_t L; // track path length … … 574 565 Float_t Zd; // Z coordinate of point of closest approach to vertex 575 566 576 Float_t XFirstHit; // X coordinate of point of closest approach to vertex577 Float_t YFirstHit; // Y coordinate of point of closest approach to vertex578 Float_t ZFirstHit; // Z coordinate of point of closest approach to vertex579 580 567 Float_t L; // track path length 581 568 Float_t D0; // track transverse impact parameter … … 651 638 ClassDef(HectorHit, 1) 652 639 }; 653 //---------------------------------------------------------------------------654 655 class CscCluster: public SortableObject656 {657 public:658 Float_t Eta; // eta of LLP659 Float_t Phi; // phi of LLP660 Float_t PT; // pt of LLP661 Float_t Px;// px of LLP662 Float_t Py;// py of LLP663 Float_t Pz;// pz of LLP664 Float_t E; // E of LLP665 Float_t Ehad; // had energy of LLP666 Float_t Eem; // em energy of LLP667 Float_t pid; // LLP pid668 Float_t T; // LLP decay time-photon travel time669 Float_t X; // LLP decay x670 Float_t Y; // LLP decay y671 Float_t Z; // LLP decay z672 Float_t R; // LLP decay z673 Float_t beta; // LLP beta674 Float_t ctau; //LLP ctau675 676 677 static CompBase *fgCompare; //!678 const CompBase *GetCompare() const { return fgCompare; }679 680 ClassDef(CscCluster, 5)681 };682 640 683 641 //--------------------------------------------------------------------------- … … 724 682 Float_t DeltaPhi; 725 683 726 TLorentzVector Momentum, Position, InitialPosition, PositionError, DecayPosition,Area;684 TLorentzVector Momentum, Position, InitialPosition, PositionError, Area; 727 685 728 686 Float_t L; // path length … … 750 708 Float_t Zd; 751 709 752 Float_t XFirstHit;753 Float_t YFirstHit;754 Float_t ZFirstHit;755 756 710 // tracking resolution 757 711 -
classes/DelphesHepMC2Reader.cc
r4564cba r54ec182 438 438 { 439 439 Candidate *candidate; 440 Candidate *candidateDaughter;441 440 map<int, pair<int, int> >::iterator itMotherMap; 442 441 map<int, pair<int, int> >::iterator itDaughterMap; … … 447 446 candidate = static_cast<Candidate *>(allParticleOutputArray->At(i)); 448 447 449 450 448 if(candidate->M1 > 0) 451 449 { … … 479 477 candidate->D1 = -1; 480 478 candidate->D2 = -1; 481 const TLorentzVector &decayPosition = candidate->Position; 482 candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position 483 } 479 } 484 480 else 485 481 { 486 482 candidate->D1 = itDaughterMap->second.first; 487 483 candidate->D2 = itDaughterMap->second.second; 488 candidateDaughter = static_cast<Candidate *>(allParticleOutputArray->At(candidate->D1)); 489 const TLorentzVector &decayPosition = candidateDaughter->Position; 490 candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position 491 492 493 } 494 } 495 } 496 } 497 498 //--------------------------------------------------------------------------- 484 } 485 } 486 } 487 } 488 489 //--------------------------------------------------------------------------- -
classes/DelphesHepMC3Reader.cc
r4564cba r54ec182 458 458 TObjArray *array; 459 459 Candidate *candidate; 460 Candidate *candidateDaughter;461 460 TParticlePDG *pdgParticle; 462 461 int pdgCode; … … 574 573 candidate->D1 = -1; 575 574 candidate->D2 = -1; 576 const TLorentzVector &decayPosition = candidate->Position;577 candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position578 575 } 579 576 else … … 581 578 candidate->D1 = itDaughterMap->second.first; 582 579 candidate->D2 = itDaughterMap->second.second; 583 candidateDaughter = static_cast<Candidate *>(allParticleOutputArray->At(candidate->D1)); 584 const TLorentzVector &decayPosition = candidateDaughter->Position; 585 candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position 586 } 587 } 588 } 589 } 590 591 //--------------------------------------------------------------------------- 580 } 581 } 582 } 583 } 584 585 //--------------------------------------------------------------------------- -
examples/ExamplePVtxFind.C
r4564cba r54ec182 45 45 TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD primary tracks", 41, -0.5, 40.5); 46 46 TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5); 47 // Chi2 for cuts48 TH1D* hChi2SngP = new TH1D("hChi2SngP", "#chi^2 for single tracks", 200, 0., 50.);49 TH1D* hChi2MaxP = new TH1D("hChi2MaxP", "#chi^2 max contribution", 200, 0., 50.);50 TH1D* hChi2SngN = new TH1D("hChi2SngN", "#chi^2 for single tracks", 200, 0., 50.);51 TH1D* hChi2MaxN = new TH1D("hChi2MaxN", "#chi^2 max contribution", 200, 0., 50.);52 53 54 47 // 55 48 // Loop over all events … … 60 53 treeReader->ReadEntry(entry); 61 54 Int_t NtrG = branchTrack->GetEntries(); 62 TVectorD** pr = new TVectorD * [NtrG]; // Track Parameters63 T MatrixDSym** cv = new TMatrixDSym * [NtrG]; // Track covariances64 Bool_t *fprvx = new Bool_t[NtrG]; // Primary vertex flag55 if(entry%500 ==0)std::cout << "Event "<<entry<<" opened containing " << NtrG << " tracks" << std::endl; 56 TVectorD** pr = new TVectorD * [NtrG]; 57 TMatrixDSym** cv = new TMatrixDSym * [NtrG]; 65 58 // 66 59 // test Particle branch … … 95 88 pr[it] = new TVectorD(obsPar); 96 89 cv[it] = new TMatrixDSym(trk->CovarianceMatrix()); 90 // Check covariance 91 /* 92 TMatrixDSymEigen Eign(*cv[it]); 93 TVectorD lambda = Eign.GetEigenValues(); 94 Bool_t good = kTRUE; 95 for(Int_t i=0; i<5; i++) if(lambda(i)<0.0) good = kFALSE; 96 if(!good){ 97 std::cout<<"Cov["<<it<<"]"<<" eigenvalues: " 98 <<lambda(0)<<", "<<lambda(1)<<", "<<lambda(2)<<", " 99 <<lambda(3)<<", "<<lambda(4)<<std::endl; 100 TMatrixDSym Ctmp = *cv[it]; 101 std::cout<<"Error D (dir - cov) "<<trk->ErrorD0 102 <<" - "<<TMath::Sqrt(Ctmp(0,0))<<std::endl; 103 } 104 */ 105 //std::cout << "Track loaded it= " << it << std::endl; 97 106 // 98 107 // Find true primary vertex … … 105 114 Double_t z = gp->Z; 106 115 Bool_t prim = kTRUE; // Is primary? 107 fprvx[it] = kFALSE;108 116 Int_t mp = gp->M1; // Mother 109 117 while (mp > 0) { … … 121 129 if (prim) { // It's a primary track 122 130 Nprim++; 123 fprvx[it] = kTRUE;124 131 xpv = x; // Store true primary 125 132 ypv = y; … … 129 136 } // End loop on tracks 130 137 } 131 if(entry%500 ==0){ 132 std::cout << "Event "<<entry<<" opened containing " << NtrG << " / "<< Nprim 133 << " Total / primary tracks"<< std::endl; 134 } 135 //std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl; 138 std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl; 136 139 // 137 140 // Find primary vertex … … 158 161 TVectorD** PrSk = new TVectorD * [1]; 159 162 TMatrixDSym** CvSk = new TMatrixDSym * [1]; 160 Double_t MaxChi2 = 9.;163 Double_t MaxChi2 = 2.0+3*2.; 161 164 for (Int_t n = 0; n < NtrG; n++) { 162 165 PrSk[0] = new TVectorD(*pr[n]); … … 166 169 Double_t Chi2One = Vskim->GetVtxChi2(); 167 170 //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl; 168 if(fprvx[n])hChi2SngP->Fill(Chi2One);169 else hChi2SngN->Fill(Chi2One);170 //171 171 if (Chi2One < MaxChi2) { 172 172 nSkimmed[nSkim] = n; … … 184 184 // Load tracks for primary fit 185 185 Int_t MinTrk = 1; // Minumum # tracks for vertex fit 186 std::vector<Int_t> trnum;187 186 if (nSkim >= MinTrk) { 188 187 TVectorD** PrFit = new TVectorD * [nSkim]; … … 191 190 PrFit[n] = new TVectorD(*pr[nSkimmed[n]]); 192 191 CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]); 193 trnum.push_back(nSkimmed[n]);194 192 } 195 193 delete[] nSkimmed; … … 202 200 // 203 201 // Remove tracks with large chi2 204 Double_t MaxChi2Fit = 8.0;202 Double_t MaxChi2Fit = 4.5; 205 203 Bool_t Done = kFALSE; 206 204 while (!Done) { … … 215 213 Double_t Chi2Mx = Chi2L[iMax]; 216 214 //std::cout << "Chi2Mx "<<Chi2Mx << std::endl; 217 if(fprvx[trnum[iMax]])hChi2MaxP->Fill(Chi2Mx);218 else hChi2MaxN->Fill(Chi2Mx);219 215 if (Chi2Mx > MaxChi2Fit && Nfound > 1) { 220 216 //std::cout << "Before remove. Nfound = "<<Nfound << std::endl; 221 217 Vtx->RemoveTrk(iMax); 222 trnum.erase(trnum.begin() + iMax);223 218 //std::cout << "After remove." << std::endl; 224 219 Nfound--; … … 235 230 if (Nfound >= Nmin) { 236 231 TVectorD xvtx = Vtx->GetVtx(); 237 //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;232 std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl; 238 233 TMatrixDSym covX = Vtx->GetVtxCov(); 239 234 Double_t Chi2 = Vtx->GetVtxChi2(); … … 261 256 delete[] PrFit; 262 257 delete[] CvFit; 263 delete[] fprvx;264 trnum.clear();265 258 } 266 259 … … 278 271 TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500); 279 272 Cnv->Divide(2, 2); 280 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);273 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111); 281 274 hXpull->Fit("gaus"); hXpull->Draw(); 282 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);275 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111); 283 276 hYpull->Fit("gaus"); hYpull->Draw(); 284 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);277 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111); 285 278 hZpull->Fit("gaus"); hZpull->Draw(); 286 279 Cnv->cd(4); hChi2->Draw(); … … 295 288 CnvN->cd(3); 296 289 hTrDiff->Draw(); 297 //298 TCanvas* CnvCh = new TCanvas("CnvCh", "#chi^2", 200, 200, 900, 500);299 CnvCh->Divide(2,1);300 CnvCh->cd(1);301 hChi2SngP->Draw();302 hChi2SngN->SetLineColor(kRed);303 hChi2SngN->Draw("SAME");304 CnvCh->cd(2);305 hChi2MaxP->Draw();306 hChi2MaxN->SetLineColor(kRed);307 hChi2MaxN->Draw("SAME");308 290 } -
external/TrackCovariance/ObsTrk.cc
r4564cba r54ec182 74 74 fCovILC.ResizeTo(5, 5); 75 75 fGenPar = XPtoPar(fGenX, fGenP, Q); 76 fGenParMm = ParToMm(fGenPar);77 76 fGenParACTS = ParToACTS(fGenPar); 78 77 fGenParILC = ParToILC(fGenPar); 79 78 // 80 79 fObsPar = GenToObsPar(fGenPar); 81 fObsParMm = ParToMm(fObsPar);82 80 fObsParACTS = ParToACTS(fObsPar); 83 81 fObsParILC = ParToILC(fObsPar); … … 85 83 fObsP = ParToP(fObsPar); 86 84 fObsQ = ParToQ(fObsPar); 87 fCovMm = CovToMm(fCov);88 85 fCovACTS = CovToACTS(fObsPar, fCov); 89 86 fCovILC = CovToILC(fCov); … … 120 117 //if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl; 121 118 Double_t minAn = fGC->GetMinAng(); 122 //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 119 //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 123 120 // << " is below grid range of " << minAn << std::endl; 124 121 Double_t maxAn = fGC->GetMaxAng(); … … 134 131 Double_t ZinNeg = fG->GetZminNeg(); 135 132 Bool_t inside = TrkUtil::IsInside(fGenX, Rin, ZinNeg, ZinPos); // Check if in inner box 136 SolTrack* trk = new SolTrack(fGenX, fGenP, fG);137 Double_t Xfirst, Yfirst, Zfirst;138 Int_t iLay = trk->FirstHit(Xfirst, Yfirst, Zfirst);139 fXfirst = TVector3(Xfirst, Yfirst, Zfirst);140 //std::cout<<"obs trk: "<<Xfirst<<","<<Yfirst<<","<<Zfirst<<std::endl;141 142 133 if (inside) 143 134 { … … 153 144 //std::cout<<"ObsTrk:: outside: x= "<<fGenX(0)<<", y= "<<fGenX(1) 154 145 // <<", z= "<<fGenX(2)<<std::endl; 146 SolTrack* trk = new SolTrack(fGenX, fGenP, fG); 155 147 Bool_t Res = kTRUE; Bool_t MS = kTRUE; 156 148 trk->CovCalc(Res, MS); // Calculate covariance matrix 157 Cov = trk->Cov(); 158 } // Track covariance159 delete trk;149 Cov = trk->Cov(); // Track covariance 150 delete trk; 151 } 160 152 // 161 153 fCov = Cov; -
external/TrackCovariance/ObsTrk.h
r4564cba r54ec182 46 46 TMatrixDSym fCovILC; // Covariance of track parameters in ILC format 47 47 // (d0, phi0, w, z0, tan(lambda)) 48 TVector3 fXfirst; // x,y,z of first track hit49 48 // 50 49 // Service routines … … 90 89 TMatrixDSym GetCovMm() { return fCov; } // in mm 91 90 TMatrixDSym GetCovACTS(){ return fCovACTS; } 92 TMatrixDSym GetCovILC() { return fCovILC; } 93 // First hit 94 TVector3 GetFirstHit() { return fXfirst; } 91 TMatrixDSym GetCovILC(){ return fCovILC; } 92 // 95 93 }; 96 94 -
external/TrackCovariance/SolTrack.cc
r4564cba r54ec182 176 176 { 177 177 // 178 // Return lists of hits associated to a track including all scattering layers. 178 // Return lists of hits associated to a track including all scattering layers. 179 179 // Return value is the total number of measurement hits 180 180 // kmh = total number of measurement layers hit for given type … … 190 190 { 191 191 Double_t R; Double_t phi; Double_t zz; 192 if (HitLayer(i, R, phi, zz)) 192 if (HitLayer(i, R, phi, zz)) 193 193 { 194 194 zhh[kh] = zz; … … 207 207 { 208 208 // 209 // Return lists of hits associated to a track for all measurement layers. 209 // Return lists of hits associated to a track for all measurement layers. 210 210 // Return value is the total number of measurement hits 211 211 // kmh = total number of measurement layers hit for given type … … 234 234 } 235 235 // 236 Int_t SolTrack::FirstHit(Double_t &Xfirst, Double_t &Yfirst, Double_t &Zfirst)237 {238 Int_t iFirst = -1;239 Int_t iFirstLay = -1; // Default return with no hits240 Xfirst = 0.;241 Yfirst = 0.;242 Zfirst = 0.;243 Int_t Nmh = nmHit(); // # measurement hits244 if(Nmh > 0){245 Int_t *ih = new Int_t [Nmh];246 Double_t *Xh = new Double_t[Nmh];247 Double_t *Yh = new Double_t[Nmh];248 Double_t *Zh = new Double_t[Nmh];249 Double_t *dh = new Double_t[Nmh];250 //251 Int_t n = HitListXYZ(ih, Xh, Yh, Zh);252 //253 for(Int_t i=0; i<Nmh; i++){254 Double_t rr = TMath::Sqrt(Xh[i]*Xh[i]+Yh[i]*Yh[i]); // Hit radius255 dh[i] = TMath::ASin(C() * TMath::Sqrt((rr * rr - D() * D()) / (1. + 2 * C() * D()))) / C(); // Arc length traveled256 }257 //258 Int_t *hord = new Int_t[Nmh]; // hit order by increasing arc length259 TMath::Sort(Nmh, dh, hord, kFALSE); // Order by increasing arc length260 iFirst = hord[0]; // First hit pointer261 Xfirst = Xh[iFirst];262 Yfirst = Yh[iFirst];263 Zfirst = Zh[iFirst];264 iFirstLay = ih[iFirst];265 //266 // Clean267 delete [] ih;268 delete [] Xh;269 delete [] Yh;270 delete [] Zh;271 delete [] dh;272 delete [] hord;273 }274 //275 return iFirstLay; // Return first hit layer number276 }277 //278 236 // Track plot 279 237 // … … 327 285 // 328 286 // 329 // Input flags: 287 // Input flags: 330 288 // Res = .TRUE. turn on resolution effects/Use standard resolutions 331 289 // .FALSE. set all resolutions to 0 … … 410 368 nz = 1.0; 411 369 } 412 Double_t corr = TMath::Abs(pxi*nx + pyi * ny + pzi * nz) / p(); 370 Double_t corr = TMath::Abs(pxi*nx + pyi * ny + pzi * nz) / p(); 413 371 Double_t Rlf = fG->lTh(i) / (corr*fG->lX0(i)); // Rad. length fraction 414 372 thms[ii] = 0.0136*TMath::Sqrt(Rlf)*(1.0 + 0.038*TMath::Log(Rlf)) / p(); // MS angle … … 438 396 Int_t ityp = fG->lTyp(i); // Layer type Barrel or Z 439 397 Int_t nmeai = fG->lND(i); // # measurements in layer 440 398 441 399 if (fG->isMeasure(i)) 442 400 { … … 470 428 Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2); // C derivative 471 429 Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3); // z0 derivative 472 Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4); // cot(theta) derivative 430 Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4); // cot(theta) derivative 473 431 } 474 432 if (ityp == 2) // Z type layer (Measure R-phi at const. Z) … … 531 489 { 532 490 Double_t strk = 0; 533 if (nmk + 1 == 1) strk = fG->lStU(k); // Stereo angle upper 491 if (nmk + 1 == 1) strk = fG->lStU(k); // Stereo angle upper 534 492 if (nmk + 1 == 2) strk = fG->lStL(k); // Stereo angle lower 535 493 //if (im == km && Res) Sm(im, km) += sig*sig; // Detector resolution on diagonal … … 542 500 // 543 501 // Loop on all layers below for MS contributions 544 for (Int_t jj = 0; jj < kk; jj++) 502 for (Int_t jj = 0; jj < kk; jj++) 545 503 { 546 504 Double_t di = dik(ii, jj); -
external/TrackCovariance/SolTrack.h
r4564cba r54ec182 73 73 Int_t HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh); 74 74 Int_t HitListXYZ(Int_t *&ihh, Double_t *&Xh, Double_t *&Yh, Double_t *&Zh); 75 Int_t FirstHit(Double_t &Xfirst, Double_t &Yfirst, Double_t &Zfirst); // First hit position76 75 // 77 76 // Track graph -
external/TrackCovariance/VertexFit.cc
r4564cba r54ec182 17 17 { 18 18 fNtr = 0; 19 fR start= -1.0;19 fRold = -1.0; 20 20 fVtxDone = kFALSE; 21 21 fVtxCst = kFALSE; … … 31 31 { 32 32 fNtr = Ntr; 33 fR start= -1.0;33 fRold = -1.0; 34 34 fVtxDone = kFALSE; 35 35 fVtxCst = kFALSE; … … 55 55 { 56 56 fNtr = Ntr; 57 fR start= -1.0;57 fRold = -1.0; 58 58 fVtxDone = kFALSE; 59 59 fVtxCst = kFALSE; … … 220 220 std::vector<TVectorD*> x0i; // Tracks at ma 221 221 std::vector<TVectorD*> ni; // Track derivative wrt phase 222 std::vector<TMatrixDSym*> Ci; 222 std::vector<TMatrixDSym*> Ci; // Position error matrix at fixed phase 223 223 std::vector<TVectorD*> wi; // Ci*ni 224 std::vector<Double_t> s_in; // Starting phase225 //226 //227 224 // 228 225 // Track loop 229 226 for (Int_t i = 0; i < fNtr; i++) 230 227 { 228 Double_t s = 0.; 231 229 TVectorD par = *fPar[i]; 232 230 TMatrixDSym Cov = *fCov[i]; 233 Double_t s = 0.; 234 // Case when starting radius is provided 235 if(fRstart > TMath::Abs(par(0))){ 236 s = 2.*TMath::ASin(par(2)*TMath::Sqrt((fRstart*fRstart-par(0)*par(0))/(1.+2.*par(2)*par(0)))); 237 } 238 // 239 x0i.push_back(new TVectorD(Fill_x(par, s))); 231 x0i.push_back(new TVectorD(Fill_x0(par))); 240 232 ni.push_back(new TVectorD(derXds(par, s))); 241 233 TMatrixD A = derXdPar(par, s); … … 243 235 TMatrixDSym Cinv = RegInv(*Ci[i]); 244 236 wi.push_back(new TVectorD(Cinv * (*ni[i]))); 245 s_in.push_back(s);246 237 } 247 238 //std::cout << "Vtx init completed. fNtr = "<<fNtr << std::endl; … … 271 262 for (Int_t i = 0; i < fNtr; i++){ 272 263 Double_t si = Dot(*wi[i], fXv - (*x0i[i])) / Ci[i]->Similarity(*wi[i]); 273 ffi.push_back(si +s_in[i]);264 ffi.push_back(si); 274 265 //TVectorD xvi = Fill_x(*fPar[i],si); 275 266 //std::cout << "Fast vertex "<<i<<": xvi = "<<xvi(0)<<", "<<xvi(1)<<", "<<xvi(2) … … 289 280 Ci.clear(); 290 281 wi.clear(); 291 s_in.clear();292 282 } 293 283 // … … 412 402 // 413 403 fVtxDone = kTRUE; // Set fit completion flag 414 //fRstart= TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1)); // Store fit404 fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1)); // Store fit 415 405 //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2) 416 406 // << ", after "<<Ntry<<" iterations"<<std::endl; -
external/TrackCovariance/VertexFit.h
r4564cba r54ec182 36 36 // Results 37 37 Bool_t fVtxDone; // Flag vertex fit completed 38 Double_t fR start; // Starting value of vertex radius (0 = none)38 Double_t fRold; // Current value of vertex radius 39 39 TVectorD fXv; // Found vertex 40 40 TMatrixDSym fcovXv; // Vertex covariance … … 83 83 // Handle tracks/constraints 84 84 void AddVtxConstraint(TVectorD xv, TMatrixDSym cov); // Add gaussian vertex constraint 85 void AddTrk(TVectorD *par, TMatrixDSym *Cov); // Add track to input list 86 void RemoveTrk(Int_t iTrk); // Remove iTrk track 87 void SetStartR(Double_t R) { fRstart = R; }; // Set starting radius 85 void AddTrk(TVectorD *par, TMatrixDSym *Cov); // Add track to input list 86 void RemoveTrk(Int_t iTrk); // Remove iTrk track 88 87 // 89 88 }; -
modules/ModulesLinkDef.h
r4564cba r54ec182 79 79 #include "modules/TruthVertexFinder.h" 80 80 #include "modules/ExampleModule.h" 81 #include "modules/LLPFilter.h"82 #include "modules/CscClusterEfficiency.h"83 #include "modules/CscClusterId.h"84 81 85 82 #ifdef __CINT__ … … 142 139 #pragma link C++ class TruthVertexFinder+; 143 140 #pragma link C++ class ExampleModule+; 144 #pragma link C++ class LLPFilter+;145 #pragma link C++ class CscClusterEfficiency+;146 #pragma link C++ class CscClusterId+;147 141 148 142 #endif -
modules/TrackCovariance.cc
r4564cba r54ec182 121 121 const TLorentzVector &candidateMomentum = particle->Momentum; 122 122 123 124 123 Bool_t inside = TrkUtil::IsInside(candidatePosition.Vect(), Rin, ZinNeg, ZinPos); // Check if in inner box 125 124 Bool_t Accept = kTRUE; … … 131 130 132 131 ObsTrk track(candidatePosition.Vect(), candidateMomentum.Vect(), candidate->Charge, fCovariance, fGeometry); 133 134 132 135 133 mother = candidate; … … 152 150 candidate->Yd = track.GetObsX().Y()*1e03; 153 151 candidate->Zd = track.GetObsX().Z()*1e03; 154 155 candidate->XFirstHit = track.GetFirstHit().X()*1e03;156 candidate->YFirstHit = track.GetFirstHit().Y()*1e03;157 candidate->ZFirstHit = track.GetFirstHit().Z()*1e03;158 152 159 153 candidate->D0 = track.GetObsPar()[0]*1e03; -
modules/TreeWriter.cc
r4564cba r54ec182 45 45 #include "TString.h" 46 46 47 #include <set>48 47 #include <algorithm> 49 48 #include <iostream> … … 77 76 fClassMap[Electron::Class()] = &TreeWriter::ProcessElectrons; 78 77 fClassMap[Muon::Class()] = &TreeWriter::ProcessMuons; 79 fClassMap[CscCluster::Class()] = &TreeWriter::ProcessCscCluster;80 78 fClassMap[Jet::Class()] = &TreeWriter::ProcessJets; 81 79 fClassMap[MissingET::Class()] = &TreeWriter::ProcessMissingET; … … 151 149 { 152 150 TIter it1(candidate->GetCandidates()); 153 set<Candidate *> s;154 set<Candidate *>::iterator it3;155 151 it1.Reset(); 156 s.clear();157 152 array->Clear(); 158 153 … … 164 159 if(candidate->GetCandidates()->GetEntriesFast() == 0) 165 160 { 166 s.insert(candidate);161 array->Add(candidate); 167 162 continue; 168 163 } … … 172 167 if(candidate->GetCandidates()->GetEntriesFast() == 0) 173 168 { 174 s.insert(candidate);169 array->Add(candidate); 175 170 continue; 176 171 } … … 180 175 while((candidate = static_cast<Candidate *>(it2.Next()))) 181 176 { 182 candidate = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); 183 if(candidate->GetCandidates()->GetEntriesFast() == 0) 184 { 185 s.insert(candidate); 186 } 177 array->Add(candidate->GetCandidates()->At(0)); 187 178 } 188 }189 190 for(it3 = s.begin(); it3 != s.end(); ++it3)191 {192 array->Add(*it3);193 179 } 194 180 } … … 252 238 entry->Z = position.Z(); 253 239 entry->T = position.T() * 1.0E-3 / c_light; 254 255 240 } 256 241 } … … 399 384 entry->Zd = candidate->Zd; 400 385 401 entry->XFirstHit = candidate->XFirstHit;402 entry->YFirstHit = candidate->YFirstHit;403 entry->ZFirstHit = candidate->ZFirstHit;404 405 386 const TLorentzVector &momentum = candidate->Momentum; 406 387 … … 563 544 entry->Zd = candidate->Zd; 564 545 565 entry->XFirstHit = candidate->XFirstHit;566 entry->YFirstHit = candidate->YFirstHit;567 entry->ZFirstHit = candidate->ZFirstHit;568 569 546 const TLorentzVector &momentum = candidate->Momentum; 570 547 … … 909 886 } 910 887 } 911 //------------------------------------------------------------------------------912 913 void TreeWriter::ProcessCscCluster(ExRootTreeBranch *branch, TObjArray *array)914 {915 TIter iterator(array);916 Candidate *candidate = 0;917 CscCluster *entry = 0;918 Double_t pt, signPz, cosTheta, eta, rapidity;919 920 const Double_t c_light = 2.99792458E8; // in unit of m/s921 922 array->Sort();923 924 925 // loop over all clusters926 iterator.Reset();927 while((candidate = static_cast<Candidate *>(iterator.Next())))928 {929 const TLorentzVector &momentum = candidate->Momentum;930 const TLorentzVector &position = candidate->DecayPosition;931 932 pt = momentum.Pt();933 cosTheta = TMath::Abs(momentum.CosTheta());934 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;935 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());936 937 entry = static_cast<CscCluster *>(branch->NewEntry());938 939 entry->SetBit(kIsReferenced);940 entry->SetUniqueID(candidate->GetUniqueID());941 942 entry->Eta = eta;943 entry->Phi = momentum.Phi();944 945 entry->PT = momentum.Pt(); // pt of LLP946 entry->Px = momentum.Px();// px of LLP947 entry->Py = momentum.Py();// py of LLP948 entry->Pz = momentum.Pz();// pz of LLP949 entry->E = momentum.E(); // E of LLP950 entry->pid = candidate->PID; // LLP pid951 entry->Eem = candidate->Eem; // LLP Eem952 entry->Ehad = candidate->Ehad; // LLP Ehad953 Double_t beta = momentum.P()/momentum.E();954 Double_t gamma = 1.0/sqrt(1-beta*beta);955 Double_t decayDistance = sqrt(pow(position.X(),2)+pow(position.Y(),2)+pow(position.Z(),2)); // mm956 entry->beta = beta; // LLP pid957 entry->ctau = decayDistance/(beta * gamma); // LLP travel time in its rest frame958 entry->T = decayDistance*(1./beta-1)* 1.0E-3/c_light*1e9; // ns959 entry->X = position.X(); // LLP decay x960 entry->Y = position.Y(); // LLP decay y961 entry->Z = position.Z(); // LLP decay z962 }963 }964 888 965 889 //------------------------------------------------------------------------------ -
modules/TreeWriter.h
r4564cba r54ec182 60 60 void ProcessElectrons(ExRootTreeBranch *branch, TObjArray *array); 61 61 void ProcessMuons(ExRootTreeBranch *branch, TObjArray *array); 62 void ProcessCscCluster(ExRootTreeBranch *branch, TObjArray *array);63 62 void ProcessTauJets(ExRootTreeBranch *branch, TObjArray *array); 64 63 void ProcessJets(ExRootTreeBranch *branch, TObjArray *array); … … 68 67 void ProcessWeight(ExRootTreeBranch *branch, TObjArray *array); 69 68 void ProcessHectorHit(ExRootTreeBranch *branch, TObjArray *array); 70 71 69 72 70 #if !defined(__CINT__) && !defined(__CLING__)
Note:
See TracChangeset
for help on using the changeset viewer.
