Fork me on GitHub

Changes in / [91ef0b8:0d65492] in git


Ignore:
Files:
1 added
1 deleted
15 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r91ef0b8 r0d65492  
    660660tmp/external/TrackCovariance/TrkUtil.$(ObjSuf): \
    661661        external/TrackCovariance/TrkUtil.$(SrcSuf)
    662 tmp/external/TrackCovariance/VertexFit.$(ObjSuf): \
    663         external/TrackCovariance/VertexFit.$(SrcSuf)
    664662tmp/modules/AngularSmearing.$(ObjSuf): \
    665663        modules/AngularSmearing.$(SrcSuf) \
    … …  
    11701168        tmp/external/TrackCovariance/SolTrack.$(ObjSuf) \
    11711169        tmp/external/TrackCovariance/TrkUtil.$(ObjSuf) \
    1172         tmp/external/TrackCovariance/VertexFit.$(ObjSuf) \
    11731170        tmp/modules/AngularSmearing.$(ObjSuf) \
    11741171        tmp/modules/BTagging.$(ObjSuf) \

  • README.md

    r91ef0b8 r0d65492  
    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)
     1[![CircleCI](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)
    22
    33Delphes

  • cards/delphes_card_IDEAtrkCov.tcl

    r91ef0b8 r0d65492  
    88#        michele.selvaggi@cern.ch
    99#####################################################################
    10 
    11 set B 2.0
    12 
     10#
    1311#######################################
    1412# Order of execution of various modules
    … …  
    6159  TreeWriter
    6260}
     61
    6362
    6463#################################
    … …  
    8180
    8281  # magnetic field, in T
    83   set Bz $B
     82  set Bz 2.0
    8483}
    8584
    … …  
    311310    }
    312311
    313     set Bz $B
    314 }
     312    set Bz 2.0
     313}
     314
    315315
    316316##############
    … …  
    767767    add Branch MissingET/momentum MissingET MissingET
    768768    add Branch ScalarHT/energy ScalarHT ScalarHT
    769 
    770     # add Info InfoName InfoValue
    771     add Info Bz $B
    772 }
     769}

  • classes/DelphesClasses.cc

    r91ef0b8 r0d65492  
    112112  TMatrixDSym Cv;
    113113  Cv.ResizeTo(5, 5);
    114 
     114 
    115115  // convert diagonal term to original units
    116   Cv(0, 0)=TMath::Power(ErrorD0, 2.);
     116  Cv(0, 0)=TMath::Power(ErrorD0*1.0E-3, 2.);
    117117  Cv(1, 1)=TMath::Power(ErrorPhi, 2.);
    118   Cv(2, 2)=TMath::Power(ErrorC, 2.);
    119   Cv(3, 3)=TMath::Power(ErrorDZ, 2.);
     118  Cv(2, 2)=TMath::Power(ErrorC*1.0E3, 2.);
     119  Cv(3, 3)=TMath::Power(ErrorDZ*1.0E-3, 2.);
    120120  Cv(4, 4)=TMath::Power(ErrorCtgTheta, 2.);
    121121
    … …  
    171171  TMatrixDSym Cv;
    172172  Cv.ResizeTo(5, 5);
    173 
     173 
    174174  // convert diagonal term to original units
    175   Cv(0, 0)=TMath::Power(ErrorD0, 2.);
     175  Cv(0, 0)=TMath::Power(ErrorD0*1.0E-3, 2.);
    176176  Cv(1, 1)=TMath::Power(ErrorPhi, 2.);
    177   Cv(2, 2)=TMath::Power(ErrorC, 2.);
    178   Cv(3, 3)=TMath::Power(ErrorDZ, 2.);
     177  Cv(2, 2)=TMath::Power(ErrorC*1.0E3, 2.);
     178  Cv(3, 3)=TMath::Power(ErrorDZ*1.0E-3, 2.);
    179179  Cv(4, 4)=TMath::Power(ErrorCtgTheta, 2.);
    180180
    … …  
    234234  NNeutrals(0),
    235235  NeutralEnergyFraction(0),  // charged energy fraction
    236   ChargedEnergyFraction(0),  // neutral energy fraction
     236  ChargedEnergyFraction(0),  // neutral energy fraction 
    237237  Beta(0),
    238238  BetaStar(0),

  • classes/DelphesModule.cc

    r91ef0b8 r0d65492  
    130130//------------------------------------------------------------------------------
    131131
    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 
    149132ExRootResult *DelphesModule::GetPlots()
    150133{

  • classes/DelphesModule.h

    r91ef0b8 r0d65492  
    5555
    5656  ExRootTreeBranch *NewBranch(const char *name, TClass *cl);
    57   void AddInfo(const char *name, Double_t value);
    5857
    5958  ExRootResult *GetPlots();

  • examples/Example6.C

    r91ef0b8 r0d65492  
    6464        // CEntral track over min Pt
    6565        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 
    6767        TH1* histCtgenH = new TH1F("h_CtgenH", "Generateded Cotangent", 100, -10., 10.);                // cot(theta) for high pt tracks;
    6868        TH1* histCtobsH = new TH1F("h_CtobsH", "Reconstructed Cotangent", 100, -10., 10.);      // cot(theta) for high pt tracks
    … …  
    7777        //
    7878        // Get magnetifc field
    79         Double_t Bz = treeReader->GetInfo("Bz");
     79        Double_t Bz = 2.0;
    8080
    8181        // Loop over all events
    … …  
    107107                                histCtobs->Fill(obsCtg);
    108108                                //
    109                                 // Get associated generated particle
     109                                // Get associated generated particle 
    110110                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
    111111                                //
    … …  
    162162                                        GenParticle* gpart = (GenParticle*)branchGenPart->At(it);
    163163                                        //
    164                                         // Plot charged particle parameters
     164                                        // Plot charged particle parameters 
    165165                                        // Only final state particles (Status = 1)
    166166                                        if (gpart->Status == 1 && TMath::Abs(gpart->Charge) > 0) {
    … …  
    290290        histAccCtgH->Draw();
    291291        Cnv5->cd(3); gPad->SetLogy(1);
    292         histPtgenC->Draw();
     292        histPtgenC->Draw(); 
    293293        histPtobsC->SetLineColor(kRed);
    294294        histPtobsC->Draw("SAME");

  • examples/ExamplePVtxFit.C

    r91ef0b8 r0d65492  
    11/*
    22Example of using vertex fitter class to fit primary vertex
    3 assumed to be generated in (0,0,0)
     3assumed to be generated in (0,0,0).
     4To run compiled version:
     5root
     6root> .L examples/LoadPvtxFit.C+
     7root> LoadPVtxFit()
     8root> ExamplePVtxFit("infile.root", Nevents)
     9
     10March 4, 2021
     11F. Bedeschi, INFN-Pisa, Italy
    412*/
    5 
    6 #ifdef __CLING__
    7 R__LOAD_LIBRARY(libDelphes)
    813
    914#include "classes/DelphesClasses.h"
    … …  
    1318#include "external/TrackCovariance/VertexFit.h"
    1419
    15 #endif
    16 
     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>
    1727
    1828//------------------------------------------------------------------------------
    … …  
    4050        //
    4151        // Loop over all events
    42         Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
     52        Int_t Nev = TMath::Min(Nevent, (Int_t) numberOfEntries);
    4353        for (Int_t entry = 0; entry < Nev; ++entry)
    4454        {
    … …  
    4656                treeReader->ReadEntry(entry);
    4757                Int_t Ntr = 0;  // # of tracks from primary vertex
    48                 Int_t NtrG = branchTrack->GetEntries();
     58                Int_t NtrG = branchTrack->GetEntries(); // Total # of tracks
    4959                TVectorD** pr = new TVectorD * [NtrG];
    5060                TMatrixDSym** cv = new TMatrixDSym * [NtrG];
    … …  
    5868                                Track* trk = (Track*)branchTrack->At(it);
    5969                                //
    60                                 // Get associated generated particle
     70                                // Get associated generated particle 
    6171                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
    6272                                //
    63                                 // Position of origin in mm
     73                                // Position of origin (mm)
    6474                                Double_t x = gp->X;
    6575                                Double_t y = gp->Y;
    6676                                Double_t z = gp->Z;
    6777                                //
    68                                 // group tracks originating from the primary vertex
     78                                // Select group of tracks originating from the primary vertex
    6979                                if (x == 0.0 && y == 0.0)
    7080                                {
    … …  
    8090                                        //
    8191                                        pr[Ntr] = new TVectorD(obsPar);
    82                                         //std::cout<<"Cov Matrix:"<<std::endl;
    83                                         //trk->CovarianceMatrix().Print();
    84                                         cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
     92                                        cv[Ntr] = new TMatrixDSym(TrkUtil::CovToMm(trk->CovarianceMatrix()));
    8593                                        Ntr++;
    8694                                }
    … …  
    107115                }
    108116
    109                 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
    110117                //
    111118                // Cleanup
    … …  
    118125        // Show resulting histograms
    119126        //
    120         TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
     127        TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex fit pulls", 50, 50, 900, 500);
    121128        Cnv->Divide(2, 2);
    122         Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     129        Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111); 
    123130        hXpull->Fit("gaus"); hXpull->Draw();
    124         Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     131        Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111); 
    125132        hYpull->Fit("gaus"); hYpull->Draw();
    126133        Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);

  • external/ExRootAnalysis/ExRootTreeReader.cc

    r91ef0b8 r0d65492  
    1212#include "TBranchElement.h"
    1313#include "TCanvas.h"
    14 #include "TParameter.h"
    1514#include "TClonesArray.h"
    1615#include "TH2.h"
    … …  
    120119//------------------------------------------------------------------------------
    121120
    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 
    135121Bool_t ExRootTreeReader::Notify()
    136122{

  • external/ExRootAnalysis/ExRootTreeReader.h

    r91ef0b8 r0d65492  
    2929
    3030  TClonesArray *UseBranch(const char *branchName);
    31   Double_t GetInfo(const char *name);
    3231
    3332private:

  • external/ExRootAnalysis/ExRootTreeWriter.cc

    r91ef0b8 r0d65492  
    1111#include "ExRootAnalysis/ExRootTreeBranch.h"
    1212
    13 #include "TParameter.h"
    1413#include "TClonesArray.h"
    1514#include "TFile.h"
    … …  
    4948  fBranches.insert(branch);
    5049  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));
    5950}
    6051

  • external/ExRootAnalysis/ExRootTreeWriter.h

    r91ef0b8 r0d65492  
    2929
    3030  ExRootTreeBranch *NewBranch(const char *name, TClass *cl);
    31   void AddInfo(const char *name, Double_t value);
    3231
    3332  void Clear();

  • external/TrackCovariance/TrkUtil.cc

    r91ef0b8 r0d65492  
    11#include "TrkUtil.h"
    2 #include <TMath.h>
    32#include <iostream>
    43
    … …  
    3231        Double_t r2 = x.Perp2();
    3332        Double_t cross = x(0) * p(1) - x(1) * p(0);
    34         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
     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);    // Phi0
    3635        Double_t D;                                                     // Impact parameter D
    3736        if (pt < 10.0) D = (T - pt) / a;
    … …  
    4241        Par(2) = C;             // Store C
    4342        //Longitudinal parameters
    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;
     43        Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));
     44        Double_t st = asin(B) / C;
    4645        Double_t ct = p(2) / pt;
    4746        Double_t z0 = x(2) - ct * st;
    … …  
    7069        //
    7170        TVector3 Xval;
    72         Xval(0) = -D * TMath::Sin(phi0);
    73         Xval(1) = D * TMath::Cos(phi0);
     71        Xval(0) = -D * sin(phi0);
     72        Xval(1) = D * cos(phi0);
    7473        Xval(2) = z0;
    7574        //
    … …  
    9392        TVector3 Pval;
    9493        Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C);
    95         Pval(0) = pt * TMath::Cos(phi0);
    96         Pval(1) = pt * TMath::Sin(phi0);
     94        Pval(0) = pt * cos(phi0);
     95        Pval(1) = pt * sin(phi0);
    9796        Pval(2) = pt * ct;
    9897        //
    … …  
    115114        pACTS(1) = 1000 * Par(3);               // z0 from m to mm
    116115        pACTS(2) = Par(1);                      // Phi0 is unchanged
    117         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
     116        pACTS(3) = atan2(1.0, Par(4));          // Theta in [0, pi] range
     117        pACTS(4) = Par(2) / (b * sqrt(1 + Par(4) * Par(4)));            // q/p in GeV
    119118        pACTS(5) = 0.0;                         // Time: currently undefined
    120119        //
    … …  
    133132        A(0, 0) = 1000.;                // D-D  conversion to mm
    134133        A(1, 2) = 1.0;          // phi0-phi0
    135         A(2, 4) = 1.0 / (TMath::Sqrt(1.0 + ct * ct) * b);       // q/p-C
     134        A(2, 4) = 1.0 / (sqrt(1.0 + ct * ct) * b);      // q/p-C
    136135        A(3, 1) = 1000.;                // z0-z0 conversion to mm
    137136        A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta)
    138         A(4, 4) = -C * ct / (b * TMath::Power(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta)
     137        A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta)
    139138        //
    140139        TMatrixDSym Cv = Cov;

  • external/TrackCovariance/VertexFit.cc

    r91ef0b8 r0d65492  
    1515        fVtxCst = kFALSE;
    1616        fxCst.ResizeTo(3);
    17         fCovCst.ResizeTo(3, 3);
     17        fCovCst.ResizeTo(3,3);
    1818        fXv.ResizeTo(3);
    19         fcovXv.ResizeTo(3, 3);
     19        fcovXv.ResizeTo(3,3);
    2020}
    2121// Parameters and covariances
    … …  
    2626        fVtxCst = kFALSE;
    2727        fxCst.ResizeTo(3);
    28         fCovCst.ResizeTo(3, 3);
     28        fCovCst.ResizeTo(3,3);
    2929        fXv.ResizeTo(3);
    30         fcovXv.ResizeTo(3, 3);
     30        fcovXv.ResizeTo(3,3);
    3131        //
    3232        fPar = trkPar;
    … …  
    3535        //
    3636        ffi = new Double_t[Ntr];                                // Fit phases
    37         fx0i = new TVectorD * [Ntr];                    // Track expansion points
     37        fx0i = new TVectorD* [Ntr];                     // Track expansion points
    3838        for (Int_t i = 0; i < Ntr; i++) fx0i[i] = new TVectorD(3);
    3939        fai = new TVectorD * [Ntr];                     // dx/dphi
    4040        for (Int_t i = 0; i < Ntr; i++) fai[i] = new TVectorD(3);
    4141        fa2i = new Double_t[Ntr];                               // a'Wa
    42         fDi = new TMatrixDSym * [Ntr];          // W-WBW
     42        fDi = new TMatrixDSym*[Ntr];            // W-WBW
    4343        for (Int_t i = 0; i < Ntr; i++) fDi[i] = new TMatrixDSym(3);
    4444        fWi = new TMatrixDSym * [Ntr];  // (ACA')^-1
    … …  
    4646        fWinvi = new TMatrixDSym * [Ntr];       // ACA'
    4747        for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3);
     48       
    4849}
    4950// ObsTrk list
    … …  
    5455        fVtxCst = kFALSE;
    5556        fxCst.ResizeTo(3);
    56         fCovCst.ResizeTo(3, 3);
     57        fCovCst.ResizeTo(3,3);
    5758        fXv.ResizeTo(3);
    58         fcovXv.ResizeTo(3, 3);
    59         //
    60         fPar = new TVectorD * [Ntr];
    61         fCov = new TMatrixDSym * [Ntr];
     59        fcovXv.ResizeTo(3,3);
     60        //
     61        fPar = new TVectorD*[Ntr];
     62        fCov = new TMatrixDSym*[Ntr];
    6263        fChi2List.ResizeTo(Ntr);
    6364        for (Int_t i = 0; i < Ntr; i++)
    6465        {
    6566                fPar[i] = new TVectorD(track[i]->GetObsPar());
     67                //cout << "fPar = " << endl; fPar[i]->Print();
    6668                fCov[i] = new TMatrixDSym(track[i]->GetCov());
     69                //cout << "fCov = " << endl; fCov[i]->Print();
    6770        }
    6871        //
    … …  
    7982        fWinvi = new TMatrixDSym * [Ntr];       // ACA'
    8083        for (Int_t i = 0; i < Ntr; i++) fWinvi[i] = new TMatrixDSym(3);
     84
     85        //cout << "VertexFit constructor executed" << endl;
    8186}
    8287//
    … …  
    8489VertexFit::~VertexFit()
    8590{
    86         fxCst.Clear();
     91        fxCst.Clear(); 
    8792        fCovCst.Clear();
    88         fXv.Clear();
    89         fcovXv.Clear();
     93        fXv.Clear();   
     94        fcovXv.Clear(); 
    9095        fChi2List.Clear();
    9196        //
    92         for (Int_t i = 0; i < fNtr; i++)
    93         {
    94                 fPar[i]->Clear();
     97        for (Int_t i= 0; i < fNtr; i++)
     98        {
     99                fPar[i]->Clear(); 
    95100                fCov[i]->Clear();
    96101                //
    … …  
    143148        H(1, 1) = nn;
    144149        TVectorD c(2);
    145         c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += n(i) * (y0(i) - x0(i));
     150        c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) +=  n(i) * (y0(i) - x0(i));
    146151        c(1) = 0; for (Int_t i = 0; i < 3; i++)c(1) += -k(i) * (y0(i) - x0(i));
    147152        TVectorD smin = (H * c);
    … …  
    153158        Double_t R2 = TMath::Sqrt(Y(0) * Y(0) + Y(1) * Y(1));
    154159        //
    155         return 0.5 * (R1 + R2);
     160        return 0.5*(R1+R2);
    156161}
    157162Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2)
    … …  
    165170        TMatrixDSym H(2);
    166171        H(0, 0) = -TMath::Cos(p2(1));
    167         H(0, 1) = TMath::Cos(p1(1));
     172        H(0, 1) =  TMath::Cos(p1(1));
    168173        H(1, 0) = -TMath::Sin(p2(1));
    169         H(1, 1) = TMath::Sin(p1(1));
     174        H(1, 1) =  TMath::Sin(p1(1));
    170175        Double_t Det = TMath::Sin(p2(1) - p1(1));
    171176        H *= 1.0 / Det;
    … …  
    190195                if (Ntry > NtryMax)
    191196                {
    192                         std::cout << "FastRv: maximum number of iteration reached" << std::endl;
     197                        cout << "FastRv: maximum number of iteration reached" << endl;
    193198                        break;
    194199                }
    … …  
    203208}
    204209
    205 TMatrixDSym VertexFit::RegInv3(TMatrixDSym& Smat0)
     210TMatrixDSym VertexFit::RegInv3(TMatrixDSym &Smat0)
    206211{
    207212        //
    … …  
    212217        if (N != 3)
    213218        {
    214                 std::cout << "RegInv3 called with  matrix size != 3. Abort & return standard inversion." << std::endl;
     219                cout << "RegInv3 called with  matrix size != 3. Abort & return standard inversion." << endl;
    215220                return Smat.Invert();
    216221        }
    … …  
    227232                        for (Int_t j = 0; j < 2; j++)Q(i, j) = RegMat(i, j);
    228233                }
    229                 Double_t Det = 1 - Q(0, 1) * Q(1, 0);
     234                Double_t Det = 1 - Q(0, 1)*Q(1, 0);
    230235                TMatrixDSym H(2);
    231236                H = Q;
    … …  
    236241                p(1) = RegMat(1, 2);
    237242                Double_t pHp = H.Similarity(p);
    238                 Double_t h = pHp - Det;
     243                Double_t h = pHp-Det;
    239244                //
    240245                TMatrixDSym pp(2); pp.Rank1Update(p);
    241                 TMatrixDSym F = (h * H) - pp.Similarity(H);
     246                TMatrixDSym F = (h*H) - pp.Similarity(H);
    242247                F *= 1.0 / Det;
    243                 TVectorD b = H * p;
     248                TVectorD b = H*p;
    244249                TMatrixDSym InvReg(3);
    245250                for (Int_t i = 0; i < 2; i++)
    … …  
    258263        else
    259264        {
    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);
     265                cout << "RegInv3: found negative elements in diagonal. Return standard inversion." << endl;
     266                return Smat.Invert();
    265267        }
    266268}
    … …  
    292294        A(2, 0) = 0.0;
    293295        // phi0
    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);
     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);
    296298        A(2, 1) = 0.0;
    297299        // 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);
     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);
    301303        // z0
    302304        A(0, 3) = 0.0;
    … …  
    351353        Double_t ct = par(4);
    352354        //
    353         x0(0) = -D * TMath::Sin(p0);
    354         x0(1) = D * TMath::Cos(p0);
     355        x0(0) = -D *TMath::Sin(p0);
     356        x0(1) = D*TMath::Cos(p0);
    355357        x0(2) = z0;
    356358        //
    … …  
    376378        x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
    377379        x(1) = x0(1) - (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
    378         x(2) = x0(2) + ct * phi / (2 * C);
     380        x(2) = x0(2) + ct*phi / (2 * C);
    379381        //
    380382        return x;
    … …  
    393395        //
    394396        // Stored quantities
    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'
     397        Double_t *fi = new Double_t[fNtr];                              // Phases
     398        TVectorD **x0i = new TVectorD*[fNtr];                   // Track expansion point
     399        TVectorD **ai = new TVectorD*[fNtr];                            // dx/dphi
     400        Double_t *a2i = new Double_t[fNtr];                             // a'Wa
     401        TMatrixDSym **Di = new TMatrixDSym*[fNtr];              // W-WBW
     402        TMatrixDSym **Wi = new TMatrixDSym*[fNtr];              // (ACA')^-1
     403        TMatrixDSym **Winvi = new TMatrixDSym*[fNtr];   // ACA'
    402404        //
    403405        // vertex radius approximation
    … …  
    413415        //
    414416        // Find track pair with largest phi difference
    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++)
     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++)
    419420        {
    420421                //ObsTrk* ti = tracks[i];
    … …  
    422423                Double_t phi1 = pari(1);
    423424
    424                 for (Int_t j = i + 1; j < fNtr; j++)
     425                for (Int_t j = i+1; j < fNtr; j++)
    425426                {
    426427                        //ObsTrk* tj = tracks[j];
    … …  
    437438        }
    438439        //
     440        //
     441        //ObsTrk* t1 = tracks[isel];
    439442        TVectorD p1 = *fPar[isel];
     443        //ObsTrk* t2 = tracks[jsel];
    440444        TVectorD p2 = *fPar[jsel];
    441445        Double_t R = FastRv1(p1, p2);
    442         if (R > 1000.0) R = Rd;
    443         R = 0.9 * R + 0.1 * Rd;
     446        if (R > 1.0) R = Rd;
     447        R = 0.9*R + 0.1*Rd;
     448        //
     449        //cout << "VertexFinder: fast R = " << R << endl;
    444450        //
    445451        // Iteration properties
    … …  
    458464                H.Zero();               // Reset H matrix
    459465                DW1D.Zero();
    460                 //
    461                 //std::cout << "VertexFinder: start loop on tracks" << std::endl;
     466                // 
     467                //cout << "VertexFinder: start loop on tracks" << endl;
    462468                for (Int_t i = 0; i < fNtr; i++)
    463469                {
    464                         // Get track helix parameters and their covariance matrix
     470                        // Get track helix parameters and their covariance matrix
     471                        //ObsTrk *t = tracks[i];
    465472                        TVectorD par = *fPar[i];
    466                         TMatrixDSym Cov = *fCov[i];
    467 
     473                        TMatrixDSym Cov = *fCov[i];
    468474                        Double_t fs;
    469475                        if (Ntry <= 0)  // Initialize all phases on first pass
    … …  
    471477                                Double_t D = par(0);
    472478                                Double_t C = par(2);
    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));
     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));
    475481                                fi[i] = fs;
    476482                        }
    … …  
    478484                        // Starting values
    479485                        //
    480                         fs = fi[i];             // Get phase
    481                         //std::cout << "VertexFinder: phase fs set" << std::endl;
     486                        fs = fi[i];                                                             // Get phase
     487                        //cout << "VertexFinder: phase fs set" << endl;
    482488                        TVectorD xs = Fill_x(par, fs);
    483                         //std::cout << "VertexFinder: position xs set" << std::endl;
     489                        //cout << "VertexFinder: position xs set" << endl;
    484490                        x0i[i] = new TVectorD(xs);                              // Start helix position
    485                         //std::cout << "VertexFinder: position x0i stored" << std::endl;
     491                        //cout << "VertexFinder: position x0i stored" << endl;
    486492                        // W matrix = (A*C*A')^-1; W^-1 = A*C*A'
    487493                        TMatrixD A = Fill_A(par, fs);                   // A = dx/da = derivatives wrt track parameters
    488                         //std::cout << "VertexFinder: derivatives A set" << std::endl;
     494                        //cout << "VertexFinder: derivatives A set" << endl;
    489495                        TMatrixDSym Winv = Cov.Similarity(A);   // W^-1 = A*C*A'
    490496                        Winvi[i] = new TMatrixDSym(Winv);               // Store W^-1 matrix
    491                         //std::cout << "VertexFinder: Winvi stored" << std::endl;
     497                        //cout << "VertexFinder: Winvi stored" << endl;
    492498                        TMatrixDSym W = RegInv3(Winv);                  // W = (A*C*A')^-1
    493499                        Wi[i] = new TMatrixDSym(W);                             // Store W matrix
    494                         //std::cout << "VertexFinder: Wi stored" << std::endl;
     500                        //cout << "VertexFinder: Wi stored" << endl;
    495501                        TVectorD a = Fill_a(par, fs);                   // a = dx/ds = derivatives wrt phase
    496                         //std::cout << "VertexFinder: derivatives a set" << std::endl;
     502                        //cout << "VertexFinder: derivatives a set" << endl;
    497503                        ai[i] = new TVectorD(a);                                // Store a
    498                         //std::cout << "VertexFinder: derivatives a stored" << std::endl;
     504                        //cout << "VertexFinder: derivatives a stored" << endl;
    499505                        Double_t a2 = W.Similarity(a);
    500506                        a2i[i] = a2;                                                    // Store a2
    501507                        // Build D matrix
    502                         TMatrixDSym B(3);
    503 
     508                        TMatrixDSym B(3);
    504509                        B.Rank1Update(a, 1.0);
    505510                        B *= -1. / a2;
    506511                        B.Similarity(W);
    507                         TMatrixDSym Ds = W + B;                                 // D matrix
     512                        TMatrixDSym Ds = W+B;                                   // D matrix
    508513                        Di[i] = new TMatrixDSym(Ds);                    // Store D matrix
    509                         //std::cout << "VertexFinder: matrix Di stored" << std::endl;
     514                        //cout << "VertexFinder: matrix Di stored" << endl;
    510515                        TMatrixDSym DsW1Ds = Winv.Similarity(Ds);       // Service matrix to calculate covX
    511                         DW1D += DsW1Ds;
     516                        DW1D += DsW1Ds;                                                         
    512517                        // Update hessian
    513518                        H += Ds;
    … …  
    518523                // update vertex position
    519524                TMatrixDSym H1 = RegInv3(H);
    520                 x = H1 * cterm;
    521                 //std::cout << "VertexFinder: x vertex set" << std::endl;
     525                x = H1*cterm;
     526                //cout << "VertexFinder: x vertex set" << endl;
    522527                // Update vertex covariance
    523528                covX = DW1D.Similarity(H1);
    524                 //std::cout << "VertexFinder: cov vertex set" << std::endl;
     529                //cout << "VertexFinder: cov vertex set" << endl;
    525530                // Update phases and chi^2
    526531                Chi2 = 0.0;
    527532                for (Int_t i = 0; i < fNtr; i++)
    528533                {
    529                         TVectorD lambda = (*Di[i]) * (*x0i[i] - x);
     534                        TVectorD lambda = (*Di[i])*(*x0i[i] - x);
    530535                        TMatrixDSym Wm1 = *Winvi[i];
    531536                        fChi2List(i) = Wm1.Similarity(lambda);
    532537                        Chi2 += fChi2List(i);
    533538                        TVectorD a = *ai[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];
     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];
    536541                }
    537542
     543                //cout << "VertexFinder: end chi2 calculation" << endl;
    538544                //
    539545                TVectorD dx = x - x0;
    … …  
    553559                //
    554560
    555                 //std::cout << "VertexFinder: before store intermediate data" << std::endl;
     561                //cout << "VertexFinder: before store intermediate data" << endl;
    556562                for (Int_t i = 0; i < fNtr; i++)
    557563                {
    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;
     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;
    561567                        ffi[i] = fi[i];                         // Fit phases
    562                         //std::cout << "VertexFinder: fi stored" << std::endl;
     568                        //cout << "VertexFinder: fi stored" << endl;
    563569                        fx0i[i] = x0i[i];                       // Track expansion points
    564                         //std::cout << "VertexFinder: x0i stored" << std::endl;
     570                        //cout << "VertexFinder: x0i stored" << endl;
    565571                        fai[i] = ai[i];                         // dx/dphi
    566                         //std::cout << "VertexFinder: ai stored" << std::endl;
     572                        //cout << "VertexFinder: ai stored" << endl;
    567573                        fa2i[i] = a2i[i];                       // a'Wa
    568                         //std::cout << "VertexFinder: a2i stored" << std::endl;
     574                        //cout << "VertexFinder: a2i stored" << endl;
    569575                        fDi[i] = Di[i];                         // W-WBW
    570                         //std::cout << "VertexFinder: Di stored" << std::endl;
     576                        //cout << "VertexFinder: Di stored" << endl;
    571577                        fWi[i] = Wi[i];                         // (ACA')^-1
    572                         //std::cout << "VertexFinder: Wi stored" << std::endl;
     578                        //cout << "VertexFinder: Wi stored" << endl;
    573579                        fWinvi[i] = Winvi[i];           // ACA'
    574                         //std::cout << "VertexFinder: Winvi stored" << std::endl;
     580                        //cout << "VertexFinder: Winvi stored" << endl;
    575581                }
    576                 //std::cout << "Iteration " << Ntry << " completed - Before cleanup" << std::endl;
     582                //cout << "Iteration " << Ntry << " completed - Before cleanup" << endl;
    577583                //
    578584                // Cleanup
    … …  
    593599                }
    594600
    595                 //std::cout << "Iteration " << Ntry << " completed - After cleanup" << std::endl;
     601                //cout << "Iteration " << Ntry << " completed - After cleanup" << endl;
    596602        }
    597603        //
    598604        fVtxDone = kTRUE;       // Set fitting completion flag
    599605        //
    600         delete[] fi;            // Phases
     606        delete[] fi;            // Phases 
    601607        delete[] x0i;           // Track expansion point
    602608        delete[] ai;            // dx/dphi
    … …  
    609615TVectorD VertexFit::GetVtx()
    610616{
    611         //std::cout << "GetVtx: flag set to " << fVtxDone << std::endl;
     617        //cout << "GetVtx: flag set to " << fVtxDone << endl;
    612618        if (!fVtxDone)VertexFinder();
    613619        return fXv;
    … …  
    635641void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov)  // Add gaussian vertex constraint
    636642{
    637         std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
     643        cout << "VertexFit::AddVtxConstraint: Not implemented yet" << endl;
    638644}
    639645//
    640646void VertexFit::AddTrk(TVectorD par, TMatrixDSym Cov)                   // Add track to input list
    641647{
    642         std::cout << "VertexFit::AddTrk: Not implemented yet" << std::endl;
     648        cout << "VertexFit::AddTrk: Not implemented yet" << endl;
    643649}
    644650void VertexFit::RemoveTrk(Int_t iTrk)                                                   // Remove iTrk track
    645651{
    646         std::cout << "VertexFit::RemoveTrk: Not implemented yet" << std::endl;
    647 }
     652        cout << "VertexFit::RemoveTrk: Not implemented yet" << endl;
     653}

  • modules/TreeWriter.cc

    r91ef0b8 r0d65492  
    123123    fBranchMap.insert(make_pair(branch, make_pair(itClassMap->second, array)));
    124124  }
    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   }
    138125}
    139126
    … …  
    151138  it1.Reset();
    152139  array->Clear();
    153 
     140 
    154141  while((candidate = static_cast<Candidate *>(it1.Next())))
    155142  {
    … …  
    550537    entry->CtgTheta = ctgTheta;
    551538    entry->C = candidate->C;
    552 
     539   
    553540    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    554541    const TLorentzVector &initialPosition = particle->Position;
Note: See TracChangeset for help on using the changeset viewer.