Fork me on GitHub

Changes in / [56fb0be:dd263e4] in git


Ignore:
Files:
1 deleted
7 edited

Legend:

Unmodified
Added
Removed

  • examples/ExamplePVtxFind.C

    r56fb0be rdd263e4  
    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        
    5447        //
    5548        // Loop over all events
    … …  
    6053                treeReader->ReadEntry(entry);
    6154                Int_t NtrG = branchTrack->GetEntries();
    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
     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];
    6558                //
    6659                // test Particle branch
    … …  
    9588                                pr[it] = new TVectorD(obsPar);
    9689                                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;
    97106                        //
    98107                        // Find true primary vertex
    … …  
    105114                                Double_t z = gp->Z;
    106115                                Bool_t prim = kTRUE;    // Is primary?
    107                                 fprvx[it] = kFALSE;
    108116                                Int_t mp = gp->M1;      // Mother
    109117                                while (mp > 0) {
    … …  
    121129                                if (prim) {             // It's a primary track
    122130                                        Nprim++;
    123                                         fprvx[it] = kTRUE;
    124131                                        xpv = x;        // Store true primary
    125132                                        ypv = y;
    … …  
    129136                        }               // End loop on tracks
    130137                }
    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;
    136139                //
    137140                // Find primary vertex
    … …  
    158161                TVectorD** PrSk = new TVectorD * [1];
    159162                TMatrixDSym** CvSk = new TMatrixDSym * [1];
    160                 Double_t MaxChi2 = 9.;
     163                Double_t MaxChi2 = 2.0+3*2.;
    161164                for (Int_t n = 0; n < NtrG; n++) {
    162165                        PrSk[0] = new TVectorD(*pr[n]);
    … …  
    166169                        Double_t Chi2One = Vskim->GetVtxChi2();
    167170                        //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;
    187186                if (nSkim >= MinTrk) {
    188187                        TVectorD** PrFit = new TVectorD * [nSkim];
    … …  
    191190                                PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
    192191                                CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
    193                                 trnum.push_back(nSkimmed[n]);
    194192                        }
    195193                        delete[] nSkimmed;
    … …  
    202200                        //
    203201                        // Remove tracks with large chi2
    204                         Double_t MaxChi2Fit = 8.0;
     202                        Double_t MaxChi2Fit = 4.5;
    205203                        Bool_t Done = kFALSE;
    206204                        while (!Done) {
    … …  
    215213                                Double_t Chi2Mx = Chi2L[iMax];
    216214                                //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
    217                                 if(fprvx[trnum[iMax]])hChi2MaxP->Fill(Chi2Mx);
    218                                 else hChi2MaxN->Fill(Chi2Mx);
    219215                                if (Chi2Mx > MaxChi2Fit && Nfound > 1) {
    220216                                        //std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
    221217                                        Vtx->RemoveTrk(iMax);
    222                                         trnum.erase(trnum.begin() + iMax);
    223218                                        //std::cout << "After remove." << std::endl;
    224219                                        Nfound--;
    … …  
    235230                        if (Nfound >= Nmin) {
    236231                                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;
    238233                                TMatrixDSym covX = Vtx->GetVtxCov();
    239234                                Double_t Chi2 = Vtx->GetVtxChi2();
    … …  
    261256                        delete[] PrFit;
    262257                        delete[] CvFit;
    263                         delete[] fprvx;
    264                         trnum.clear();
    265258                }
    266259
    … …  
    278271        TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
    279272        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);
    281274        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);
    283276        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);
    285278        hZpull->Fit("gaus"); hZpull->Draw();
    286279        Cnv->cd(4); hChi2->Draw();
    … …  
    295288        CnvN->cd(3);
    296289        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");
    308290}

  • external/TrackCovariance/ObsTrk.cc

    r56fb0be rdd263e4  
    7474        fCovILC.ResizeTo(5, 5);
    7575        fGenPar = XPtoPar(fGenX, fGenP, Q);
    76         fGenParMm = ParToMm(fGenPar);
    7776        fGenParACTS = ParToACTS(fGenPar);
    7877        fGenParILC = ParToILC(fGenPar);
    7978        //
    8079        fObsPar = GenToObsPar(fGenPar);
    81         fObsParMm = ParToMm(fObsPar);
    8280        fObsParACTS = ParToACTS(fObsPar);
    8381        fObsParILC = ParToILC(fObsPar);
    … …  
    8583        fObsP = ParToP(fObsPar);
    8684        fObsQ = ParToQ(fObsPar);
    87         fCovMm = CovToMm(fCov);
    8885        fCovACTS = CovToACTS(fObsPar, fCov);
    8986        fCovILC = CovToILC(fCov);
    … …  
    134131        Double_t ZinNeg = fG->GetZminNeg();
    135132        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         TVector3 fXfirst(Xfirst, Yfirst, Zfirst);
    140133        if (inside)
    141134        {
    … …  
    151144                //std::cout<<"ObsTrk:: outside: x= "<<fGenX(0)<<", y= "<<fGenX(1)
    152145                //                         <<", z= "<<fGenX(2)<<std::endl;
     146                SolTrack* trk = new SolTrack(fGenX, fGenP, fG);
    153147                Bool_t Res = kTRUE; Bool_t MS = kTRUE;
    154148                trk->CovCalc(Res, MS);                                  // Calculate covariance matrix
    155                 Cov = trk->Cov();
    156         }                                       // Track covariance
    157         delete trk;
     149                Cov = trk->Cov();                                       // Track covariance
     150                delete trk;
     151        }
    158152        //
    159153        fCov = Cov;

  • external/TrackCovariance/ObsTrk.h

    r56fb0be rdd263e4  
    4646        TMatrixDSym fCovILC;                    // Covariance of track parameters in ILC format
    4747                                                                        // (d0, phi0, w, z0, tan(lambda))
    48         TVector3 fXfirst;                       // x,y,z of first track hit
    4948        //
    5049        // Service routines
    … …  
    9089        TMatrixDSym GetCovMm()  { return fCov; }        // in mm
    9190        TMatrixDSym GetCovACTS(){ return fCovACTS; }
    92         TMatrixDSym GetCovILC() { return fCovILC; }
    93         // First hit
    94         TVector3 GetFirstHit()  { return fXfirst; }
     91        TMatrixDSym GetCovILC(){ return fCovILC; }
     92        //
    9593};
    9694

  • external/TrackCovariance/SolTrack.cc

    r56fb0be rdd263e4  
    232232        //
    233233        return kmh;
    234 }
    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 hits
    240         Xfirst = 0.;
    241         Yfirst = 0.;
    242         Zfirst = 0.;
    243         Int_t Nmh = nmHit();    // # measurement hits
    244         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 radius           
    255                         dh[i] = TMath::ASin(C() * TMath::Sqrt((rr * rr - D() * D()) / (1. + 2 * C() * D()))) / C();     // Arc length traveled
    256                 }
    257                 //
    258                 Int_t *hord = new Int_t[Nmh];                   // hit order by increasing arc length
    259                 TMath::Sort(Nmh, dh, hord, kFALSE);             // Order by increasing arc length
    260                 iFirst = hord[0];                               // First hit pointer
    261                 Xfirst = Xh[iFirst];
    262                 Yfirst = Yh[iFirst];
    263                 Zfirst = Zh[iFirst];
    264                 iFirstLay = ih[iFirst];
    265                 //
    266                 // Clean
    267                 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 number
    276234}
    277235//

  • external/TrackCovariance/SolTrack.h

    r56fb0be rdd263e4  
    7373        Int_t HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh);
    7474        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 position
    7675        //
    7776        // Track graph

  • external/TrackCovariance/VertexFit.cc

    r56fb0be rdd263e4  
    1717{
    1818        fNtr = 0;
    19         fRstart = -1.0;
     19        fRold = -1.0;
    2020        fVtxDone = kFALSE;
    2121        fVtxCst = kFALSE;
    … …  
    3131{
    3232        fNtr = Ntr;
    33         fRstart = -1.0;
     33        fRold = -1.0;
    3434        fVtxDone = kFALSE;
    3535        fVtxCst = kFALSE;
    … …  
    5555{
    5656        fNtr = Ntr;
    57         fRstart = -1.0;
     57        fRold = -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         //
    227224        //
    228225        // Track loop
    229226        for (Int_t i = 0; i < fNtr; i++)
    230227        {
     228                Double_t s = 0.;
    231229                TVectorD par = *fPar[i];
    232230                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)));
    240232                ni.push_back(new TVectorD(derXds(par, s)));
    241233                TMatrixD A = derXdPar(par, s);
    … …  
    243235                TMatrixDSym Cinv = RegInv(*Ci[i]);
    244236                wi.push_back(new TVectorD(Cinv * (*ni[i])));
    245                 s_in.push_back(s);
    246237        }
    247238        //std::cout << "Vtx init completed. fNtr = "<<fNtr << std::endl;
    … …  
    271262        for (Int_t i = 0; i < fNtr; i++){
    272263                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);
    274265                //TVectorD xvi = Fill_x(*fPar[i],si);
    275266                //std::cout << "Fast vertex "<<i<<": xvi = "<<xvi(0)<<", "<<xvi(1)<<", "<<xvi(2)
    … …  
    289280        Ci.clear();
    290281        wi.clear();
    291         s_in.clear();
    292282}
    293283//
    … …  
    412402        //
    413403        fVtxDone = kTRUE;               // Set fit completion flag
    414         //fRstart = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1)); // Store fit
     404        fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1));     // Store fit
    415405        //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2)
    416406        //      << ", after "<<Ntry<<" iterations"<<std::endl;

  • external/TrackCovariance/VertexFit.h

    r56fb0be rdd263e4  
    3636        // Results
    3737        Bool_t fVtxDone;                        // Flag vertex fit completed
    38         Double_t fRstart;                       // Starting value of vertex radius (0 = none)
     38        Double_t fRold;                         // Current value of vertex radius
    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
    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
    8887        //
    8988};
Note: See TracChangeset for help on using the changeset viewer.