Fork me on GitHub

Changeset 7bca620 in git


Ignore:
Timestamp:
Feb 10, 2022, 11:42:45 AM (3 years ago)
Author:
Franco BEDESCHI <bed@…>
Branches:
master
Children:
5c03893
Parents:
b750b0a
Message:

Add feature to force starting radius in vertex fit

Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • examples/ExamplePVtxFind.C

    rb750b0a r7bca620  
    4545        TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD  primary tracks", 41, -0.5, 40.5);
    4646        TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5);
     47        // Chi2 for cuts
     48        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       
    4754        //
    4855        // Loop over all events
    … …  
    5360                treeReader->ReadEntry(entry);
    5461                Int_t NtrG = branchTrack->GetEntries();
    55                 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];
     62                TVectorD** pr = new TVectorD * [NtrG];          // Track Parameters
     63                TMatrixDSym** cv = new TMatrixDSym * [NtrG];    // Track covariances
     64                Bool_t *fprvx = new Bool_t[NtrG];               // Primary vertex flag
    5865                //
    5966                // test Particle branch
    … …  
    8895                                pr[it] = new TVectorD(obsPar);
    8996                                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;
    10697                        //
    10798                        // Find true primary vertex
    … …  
    114105                                Double_t z = gp->Z;
    115106                                Bool_t prim = kTRUE;    // Is primary?
     107                                fprvx[it] = kFALSE;
    116108                                Int_t mp = gp->M1;      // Mother
    117109                                while (mp > 0) {
    … …  
    129121                                if (prim) {             // It's a primary track
    130122                                        Nprim++;
     123                                        fprvx[it] = kTRUE;
    131124                                        xpv = x;        // Store true primary
    132125                                        ypv = y;
    … …  
    136129                        }               // End loop on tracks
    137130                }
    138                 std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;
     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;
    139136                //
    140137                // Find primary vertex
    … …  
    161158                TVectorD** PrSk = new TVectorD * [1];
    162159                TMatrixDSym** CvSk = new TMatrixDSym * [1];
    163                 Double_t MaxChi2 = 2.0+3*2.;
     160                Double_t MaxChi2 = 9.;
    164161                for (Int_t n = 0; n < NtrG; n++) {
    165162                        PrSk[0] = new TVectorD(*pr[n]);
    … …  
    169166                        Double_t Chi2One = Vskim->GetVtxChi2();
    170167                        //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
     168                        if(fprvx[n])hChi2SngP->Fill(Chi2One);
     169                        else hChi2SngN->Fill(Chi2One);
     170                        //
    171171                        if (Chi2One < MaxChi2) {
    172172                                nSkimmed[nSkim] = n;
    … …  
    184184                // Load tracks for primary fit
    185185                Int_t MinTrk = 1;       // Minumum # tracks for vertex fit
     186                std::vector<Int_t> trnum;
    186187                if (nSkim >= MinTrk) {
    187188                        TVectorD** PrFit = new TVectorD * [nSkim];
    … …  
    190191                                PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
    191192                                CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
     193                                trnum.push_back(nSkimmed[n]);
    192194                        }
    193195                        delete[] nSkimmed;
    … …  
    200202                        //
    201203                        // Remove tracks with large chi2
    202                         Double_t MaxChi2Fit = 4.5;
     204                        Double_t MaxChi2Fit = 8.0;
    203205                        Bool_t Done = kFALSE;
    204206                        while (!Done) {
    … …  
    213215                                Double_t Chi2Mx = Chi2L[iMax];
    214216                                //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
     217                                if(fprvx[trnum[iMax]])hChi2MaxP->Fill(Chi2Mx);
     218                                else hChi2MaxN->Fill(Chi2Mx);
    215219                                if (Chi2Mx > MaxChi2Fit && Nfound > 1) {
    216220                                        //std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
    217221                                        Vtx->RemoveTrk(iMax);
     222                                        trnum.erase(trnum.begin() + iMax);
    218223                                        //std::cout << "After remove." << std::endl;
    219224                                        Nfound--;
    … …  
    230235                        if (Nfound >= Nmin) {
    231236                                TVectorD xvtx = Vtx->GetVtx();
    232                                 std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
     237                                //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
    233238                                TMatrixDSym covX = Vtx->GetVtxCov();
    234239                                Double_t Chi2 = Vtx->GetVtxChi2();
    … …  
    256261                        delete[] PrFit;
    257262                        delete[] CvFit;
     263                        delete[] fprvx;
     264                        trnum.clear();
    258265                }
    259266
    … …  
    271278        TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
    272279        Cnv->Divide(2, 2);
    273         Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     280        Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
    274281        hXpull->Fit("gaus"); hXpull->Draw();
    275         Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     282        Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
    276283        hYpull->Fit("gaus"); hYpull->Draw();
    277         Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     284        Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
    278285        hZpull->Fit("gaus"); hZpull->Draw();
    279286        Cnv->cd(4); hChi2->Draw();
    … …  
    288295        CnvN->cd(3);
    289296        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");
    290308}

  • external/TrackCovariance/ObsTrk.cc

    rb750b0a r7bca620  
    7474        fCovILC.ResizeTo(5, 5);
    7575        fGenPar = XPtoPar(fGenX, fGenP, Q);
     76        fGenParMm = ParToMm(fGenPar);
    7677        fGenParACTS = ParToACTS(fGenPar);
    7778        fGenParILC = ParToILC(fGenPar);
    7879        //
    7980        fObsPar = GenToObsPar(fGenPar);
     81        fObsParMm = ParToMm(fObsPar);
    8082        fObsParACTS = ParToACTS(fObsPar);
    8183        fObsParILC = ParToILC(fObsPar);
    … …  
    8385        fObsP = ParToP(fObsPar);
    8486        fObsQ = ParToQ(fObsPar);
     87        fCovMm = CovToMm(fCov);
    8588        fCovACTS = CovToACTS(fObsPar, fCov);
    8689        fCovILC = CovToILC(fCov);

  • external/TrackCovariance/VertexFit.cc

    rb750b0a r7bca620  
    1717{
    1818        fNtr = 0;
    19         fRold = -1.0;
     19        fRstart = -1.0;
    2020        fVtxDone = kFALSE;
    2121        fVtxCst = kFALSE;
    … …  
    3131{
    3232        fNtr = Ntr;
    33         fRold = -1.0;
     33        fRstart = -1.0;
    3434        fVtxDone = kFALSE;
    3535        fVtxCst = kFALSE;
    … …  
    5555{
    5656        fNtr = Ntr;
    57         fRold = -1.0;
     57        fRstart = -1.0;
    5858        fVtxDone = kFALSE;
    5959        fVtxCst = kFALSE;
    … …  
    220220        std::vector<TVectorD*> x0i;                             // Tracks at ma
    221221        std::vector<TVectorD*> ni;                              // Track derivative wrt phase
    222         std::vector<TMatrixDSym*> Ci;                   // Position error matrix at fixed phase
     222        std::vector<TMatrixDSym*> Ci;                           // Position error matrix at fixed phase
    223223        std::vector<TVectorD*> wi;                              // Ci*ni
     224        std::vector<Double_t> s_in;                             // Starting phase
     225        //
     226        //
    224227        //
    225228        // Track loop
    226229        for (Int_t i = 0; i < fNtr; i++)
    227230        {
    228                 Double_t s = 0.;
    229231                TVectorD par = *fPar[i];
    230232                TMatrixDSym Cov = *fCov[i];
    231                 x0i.push_back(new TVectorD(Fill_x0(par)));
     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)));
    232240                ni.push_back(new TVectorD(derXds(par, s)));
    233241                TMatrixD A = derXdPar(par, s);
    … …  
    235243                TMatrixDSym Cinv = RegInv(*Ci[i]);
    236244                wi.push_back(new TVectorD(Cinv * (*ni[i])));
     245                s_in.push_back(s);
    237246        }
    238247        //std::cout << "Vtx init completed. fNtr = "<<fNtr << std::endl;
    … …  
    262271        for (Int_t i = 0; i < fNtr; i++){
    263272                Double_t si = Dot(*wi[i], fXv - (*x0i[i])) / Ci[i]->Similarity(*wi[i]);
    264                 ffi.push_back(si);
     273                ffi.push_back(si+s_in[i]);
    265274                //TVectorD xvi = Fill_x(*fPar[i],si);
    266275                //std::cout << "Fast vertex "<<i<<": xvi = "<<xvi(0)<<", "<<xvi(1)<<", "<<xvi(2)
    … …  
    280289        Ci.clear();
    281290        wi.clear();
     291        s_in.clear();
    282292}
    283293//
    … …  
    402412        //
    403413        fVtxDone = kTRUE;               // Set fit completion flag
    404         fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1));     // Store fit
     414        //fRstart = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1)); // Store fit
    405415        //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2)
    406416        //      << ", after "<<Ntry<<" iterations"<<std::endl;

  • external/TrackCovariance/VertexFit.h

    rb750b0a r7bca620  
    3636        // Results
    3737        Bool_t fVtxDone;                        // Flag vertex fit completed
    38         Double_t fRold;                         // Current value of vertex radius
     38        Double_t fRstart;                       // Starting value of vertex radius (0 = none)
    3939        TVectorD fXv;                           // Found vertex
    4040        TMatrixDSym fcovXv;                     // Vertex covariance
    … …  
    8383        // Handle tracks/constraints
    8484        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
     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
    8788        //
    8889};
Note: See TracChangeset for help on using the changeset viewer.