Fork me on GitHub

Changeset 7bca620 in git for examples


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

Add feature to force starting radius in vertex fit

Location:
examples
Files:
1 added
1 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}
Note: See TracChangeset for help on using the changeset viewer.