Fork me on GitHub

Changeset df408d2 in git


Ignore:
Timestamp:
Mar 2, 2021, 3:14:09 PM (4 years ago)
Author:
Franco BEDESCHI <bed@…>
Branches:
master
Children:
92b8d11
Parents:
d3165fa
Message:

Fixed compatibility with ROOT5/minor changes to TrkUtil

Files:
3 edited

Legend:

Unmodified
Added
Removed

  • examples/Example6.C

    rd3165fa rdf408d2  
    5151        // Track parameter pulls
    5252        TH1* histDpull = new TH1F("h_Dpull", "Pull impact parameter", 100, -10., 10.);
    53         TH1* histP0pull = new TH1F("h_P0pull", "Pull phi_{0}", 100, -10., 10.);
     53        TH1* histP0pull = new TH1F("h_P0pull", "Pull #phi_{0}", 100, -10., 10.);
    5454        TH1* histCpull = new TH1F("h_Cpull", "Pull half curvature", 100, -10., 10.);
    55         TH1* histZ0pull = new TH1F("h_Z0pull", "Pull Z0", 100, -10., 10.);
     55        TH1* histZ0pull = new TH1F("h_Z0pull", "Pull Z_{0}", 100, -10., 10.);
    5656        TH1* histCtpull = new TH1F("h_Ctpull", "Pull cotg(#theta)", 100, -10., 10.);
    5757        //
    … …  
    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
    … …  
    7171        TH1* histAccPt = new TH1F("h_AccPt", "Pt acceptance", 500, 0., 50.);
    7272        // Tracks over min pt
    73         TH1* histAccCtgH = new TH1F("h_AccCtgH", "Cotangent acceptance (high pt)", 500, -10., 10.);
     73        TH1* histAccCtgH = new TH1F("h_AccCtgH", "Cotangent acceptance (high pt)", 100, -10., 10.);
    7474        // tracks above min angle
    7575        TH1* histAccPtC = new TH1F("h_AccPtC", "Pt acceptance (central tracks)", 500, 0., 50.);
    7676
    7777        //
    78         // Get magnetic field
     78        // Get magnetifc field
    7979        Double_t Bz = 2.0;
    8080
    … …  
    106106                                histZ0obs->Fill(obsZ0);
    107107                                histCtobs->Fill(obsCtg);
    108 
     108                                //
     109                                // Get associated generated particle
    109110                                GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
    110 
     111                                //
    111112                                // Position of origin in meters
    112113                                Double_t x = 1.0e-3 * gp->X;
    … …  
    139140                                Double_t pullZ0  = (obsZ0  - genZ0)  / trk->ErrorDZ;
    140141                                Double_t pullCtg = (obsCtg - genCtg) / trk->ErrorCtgTheta;
    141 
    142                                 //
    143                                 histDpull->Fill(pullD0);
     142                                //
     143                                histDpull ->Fill(pullD0);
    144144                                histP0pull->Fill(pullPhi);
    145                                 histCpull->Fill(pullC);
     145                                histCpull ->Fill(pullC);
    146146                                histZ0pull->Fill(pullZ0);
    147147                                histCtpull->Fill(pullCtg);
    … …  
    149149                                // Acceptance plots
    150150                                Double_t obsPt = trk->PT;
     151                                histPtobs->Fill(obsPt);
    151152                                if (TMath::Abs(genCtg) < maxCtAcc)histPtobsC->Fill(obsPt); // pt for central tracks
    152                                 if (genPt > minPtAcc)histCtobsH->Fill(obsCtg);                          // cot(theta) for high pt tracks
     153                                if (obsPt > minPtAcc)histCtobsH->Fill(obsCtg);                          // cot(theta) for high pt tracks
    153154                        }               // End loop on tracks
    154155
    … …  
    161162                                        GenParticle* gpart = (GenParticle*)branchGenPart->At(it);
    162163                                        //
    163                                         // Plot charged particle parameters
     164                                        // Plot charged particle parameters 
    164165                                        // Only final state particles (Status = 1)
    165166                                        if (gpart->Status == 1 && TMath::Abs(gpart->Charge) > 0) {
    166167                                                //
    167                                                 // Position of origin
    168                                                 Double_t x = gpart->X;
    169                                                 Double_t y = gpart->Y;
    170                                                 Double_t z = gpart->Z;
     168                                                // Position of origin in meters
     169                                                Double_t x = 1.0e-3 * gpart->X;
     170                                                Double_t y = 1.0e-3 * gpart->Y;
     171                                                Double_t z = 1.0e-3 * gpart->Z;
    171172                                                TVector3 xv(x, y, z);
    172173                                                //
    … …  
    179180                                                // Original parameters
    180181                                                Double_t Q = gpart->Charge;
    181                                                 TVectorD genPar = TrkUtil::XPtoPar(xv, pv, Q, Bz);      // Get parameters from position and momentum
     182                                                TVectorD gPar = TrkUtil::XPtoPar(xv, pv, Q, Bz);        // Get parameters from position and momentum
     183                                                TVectorD genPar = TrkUtil::ParToMm(gPar);                       // Convert to mm
    182184                                                // Unpack parameters
    183                                                 Double_t genD0 = genPar(0);
     185                                                Double_t genD0  = genPar(0);
    184186                                                Double_t genPhi = genPar(1);
    185                                                 Double_t genC = genPar(2);
    186                                                 Double_t genZ0 = genPar(3);
     187                                                Double_t genC   = genPar(2);
     188                                                Double_t genZ0  = genPar(3);
    187189                                                Double_t genCtg = genPar(4);
     190                                                Double_t genPt  = pv.Pt();
    188191                                                // Fill histograms
    189                                                 histDgen->Fill(genD0);
     192                                                histDgen ->Fill(genD0);
    190193                                                histP0gen->Fill(genPhi);
    191                                                 histCgen->Fill(genC);
     194                                                histCgen ->Fill(genC);
    192195                                                histZ0gen->Fill(genZ0);
    193196                                                histCtgen->Fill(genCtg);
    194197                                                //
    195198                                                // Acceptance plots
    196                                                 if (TMath::Abs(genCtg) < maxCtAcc)histPtgenC->Fill(gpart->PT);
    197                                                 if (gpart->PT > minPtAcc)histCtgenH->Fill(genCtg);
     199                                                histPtgen->Fill(genPt);
     200                                                if (TMath::Abs(genCtg) < maxCtAcc)histPtgenC->Fill(genPt);
     201                                                if (genPt > minPtAcc)histCtgenH->Fill(genCtg);
    198202                                        }
    199203                                }
    … …  
    222226        TCanvas* Cnv2 = new TCanvas("Cnv2", "Delphes observed track pulls", 150, 150, 900, 500);
    223227        Cnv2->Divide(3, 2);
    224         Cnv2->cd(1); histDpull->Draw();
    225         Cnv2->cd(2); histP0pull->Draw();
    226         Cnv2->cd(3); histCpull->Draw();
    227         Cnv2->cd(4); histZ0pull->Draw();
    228         Cnv2->cd(5); histCtpull->Draw();
     228        Cnv2->cd(1); gPad->SetLogy(1);
     229        histDpull->Draw();
     230        Cnv2->cd(2); gPad->SetLogy(1);
     231        histP0pull->Draw();
     232        Cnv2->cd(3); gPad->SetLogy(1);
     233        histCpull->Draw();
     234        Cnv2->cd(4); gPad->SetLogy(1);
     235        histZ0pull->Draw();
     236        Cnv2->cd(5); gPad->SetLogy(1);
     237        histCtpull->Draw();
    229238        //
    230239        // Acceptance plots
    231240        //
    232         TCanvas* Cnv4 = new TCanvas("Cnv4", "Delphes acceptance", 350, 350, 500, 500);
     241        TCanvas* Cnv4 = new TCanvas("Cnv4", "Delphes acceptance", 200, 200, 500, 500);
    233242        Cnv4->Divide(2, 2);
    234243        Cnv4->cd(1);
    235244        histCtgen->Draw();
     245        histCtobs->SetLineColor(kRed);
    236246        histCtobs->Draw("SAME");
    237247        Cnv4->cd(2);
    … …  
    246256        }
    247257        histAccCtg->Draw();
    248         Cnv4->cd(3);
     258        Cnv4->cd(3); gPad->SetLogy(1);
    249259        histPtgen->Draw();
     260        histPtobs->SetLineColor(kRed);
    250261        histPtobs->Draw("SAME");
    251262        Cnv4->cd(4);
    … …  
    261272        histAccPt->Draw();
    262273        //
    263         TCanvas* Cnv5 = new TCanvas("Cnv5", "Delphes acceptance (constrained)", 400, 400, 500, 500);
     274        TCanvas* Cnv5 = new TCanvas("Cnv5", "Delphes acceptance (constrained)", 250, 250, 500, 500);
    264275        Cnv5->Divide(2, 2);
    265276        Cnv5->cd(1);
    266277        histCtgenH->Draw();
     278        histCtobsH->SetLineColor(kRed);
    267279        histCtobsH->Draw("SAME");
    268280        Cnv5->cd(2);
    … …  
    277289        }
    278290        histAccCtgH->Draw();
    279         Cnv5->cd(3);
    280         histPtgenC->Draw();
     291        Cnv5->cd(3); gPad->SetLogy(1);
     292        histPtgenC->Draw();
     293        histPtobsC->SetLineColor(kRed);
    281294        histPtobsC->Draw("SAME");
    282295        Cnv5->cd(4);

  • external/TrackCovariance/TrkUtil.cc

    rd3165fa rdf408d2  
    11#include "TrkUtil.h"
     2#include <iostream>
    23
    34// Constructor
    … …  
    3031        Double_t r2 = x.Perp2();
    3132        Double_t cross = x(0) * p(1) - x(1) * p(0);
    32         Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2);
    33         Double_t phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T);     // Phi0
     33        Double_t T = TMath::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
    3435        Double_t D;                                                     // Impact parameter D
    3536        if (pt < 10.0) D = (T - pt) / a;
    … …  
    4041        Par(2) = C;             // Store C
    4142        //Longitudinal parameters
    42         Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));
    43         Double_t st = TMath::ASin(B) / C;
     43        Double_t B = C * TMath::Sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));
     44        Double_t st = asin(B) / C;
    4445        Double_t ct = p(2) / pt;
    4546        Double_t z0 = x(2) - ct * st;
    … …  
    6869        //
    6970        TVector3 Xval;
    70         Xval(0) = -D * TMath::Sin(phi0);
    71         Xval(1) = D * TMath::Cos(phi0);
     71        Xval(0) = -D * sin(phi0);
     72        Xval(1) = D * cos(phi0);
    7273        Xval(2) = z0;
    7374        //
    … …  
    7677//
    7778TVector3 TrkUtil::ParToP(TVectorD Par)
     79{
     80        if (fBz == 0.0)
     81std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl;
     82        //
     83        return ParToP(Par,fBz);
     84}
     85//
     86TVector3 TrkUtil::ParToP(TVectorD Par, Double_t Bz)
    7887{
    7988        Double_t C = Par(2);
    … …  
    8291        //
    8392        TVector3 Pval;
    84         Double_t pt = fBz * cSpeed() / TMath::Abs(2 * C);
    85         Pval(0) = pt * TMath::Cos(phi0);
    86         Pval(1) = pt * TMath::Sin(phi0);
     93        Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C);
     94        Pval(0) = pt * cos(phi0);
     95        Pval(1) = pt * sin(phi0);
    8796        Pval(2) = pt * ct;
    8897        //
    … …  
    103112        Double_t b = -cSpeed() * fBz / 2.;
    104113        pACTS(0) = 1000 * Par(0);               // D from m to mm
    105         pACTS(1) = 1000 * Par(3);       // z0 from m to mm
     114        pACTS(1) = 1000 * Par(3);               // z0 from m to mm
    106115        pACTS(2) = Par(1);                      // Phi0 is unchanged
    107         pACTS(3) = TMath::ATan2(1.0, Par(4));           // Theta in [0, pi] range
    108         pACTS(4) = Par(2) / (b * 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 * TMath::Sqrt(1 + Par(4) * Par(4)));             // q/p in GeV
    109118        pACTS(5) = 0.0;                         // Time: currently undefined
    110119        //
    … …  
    123132        A(0, 0) = 1000.;                // D-D  conversion to mm
    124133        A(1, 2) = 1.0;          // phi0-phi0
    125         A(2, 4) = 1.0 / (sqrt(1.0 + ct * ct) * b);      // q/p-C
     134        A(2, 4) = 1.0 / (TMath::Sqrt(1.0 + ct * ct) * b);       // q/p-C
    126135        A(3, 1) = 1000.;                // z0-z0 conversion to mm
    127136        A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta)

  • external/TrackCovariance/TrkUtil.h

    rd3165fa rdf408d2  
    2020        void SetBfield(Double_t Bz) { fBz = Bz; }
    2121        TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q);
    22         TVector3 ParToX(TVectorD Par);
    2322        TVector3 ParToP(TVectorD Par);
    24         Double_t ParToQ(TVectorD Par);
    2523        //
    2624        // Conversion to ACTS parametrization
    … …  
    4644        static Double_t cSpeed()
    4745        {
    48                 return TMath::C()*1.0e-9;       // Reduced speed of light
     46                Double_t c = 2.99792458e8;      // speed of light m/sec
     47                //return TMath::C()*1.0e-9;     // Incompatible with root5
     48                return c*1.0e-9;                // Reduced speed of light       
    4949        }
    5050        //
    … …  
    5252        //
    5353        static TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q, Double_t Bz);
     54        static TVector3 ParToX(TVectorD Par);                           // position of minimum distance from z axis
     55        static TVector3 ParToP(TVectorD Par, Double_t Bz);      // Get Momentum from track parameters
     56        static Double_t ParToQ(TVectorD Par);                           // Get track charge
    5457        //
    5558        // Conversion from meters to mm
Note: See TracChangeset for help on using the changeset viewer.