Fork me on GitHub

Ignore:
Timestamp:
Apr 12, 2021, 6:33:23 PM (4 years ago)
Author:
Franco BEDESCHI <bed@…>
Branches:
master
Children:
c5696dd
Parents:
3cfe61d
Message:

Vertex fit add/remove tracks - TrkUtil has cluster counting info

File:
1 edited

Legend:

Unmodified
Added
Removed
  • external/TrackCovariance/TrkUtil.cc

    r3cfe61d r82db145  
    11#include "TrkUtil.h"
    2 #include <TMath.h>
    32#include <iostream>
    43
     
    76{
    87        fBz = Bz;
     8        fGasSel = 0;                            // Default is He-Isobuthane (90-10)
     9        fRmin = 0.0;                            // Lower                DCH radius
     10        fRmax = 0.0;                            // Higher       DCH radius
     11        fZmin = 0.0;                            // Lower                DCH z
     12        fZmax = 0.0;                            // Higher       DCH z
    913}
    1014TrkUtil::TrkUtil()
    1115{
    1216        fBz = 0.0;
     17        fGasSel = 0;                            // Default is He-Isobuthane (90-10)
     18        fRmin = 0.0;                            // Lower                DCH radius
     19        fRmax = 0.0;                            // Higher       DCH radius
     20        fZmin = 0.0;                            // Lower                DCH z
     21        fZmax = 0.0;                            // Higher       DCH z
    1322}
    1423//
     
    1726{
    1827        fBz = 0.0;
     28        fGasSel = 0;                            // Default is He-Isobuthane (90-10)
     29        fRmin = 0.0;                            // Lower                DCH radius
     30        fRmax = 0.0;                            // Higher       DCH radius
     31        fZmin = 0.0;                            // Lower                DCH z
     32        fZmax = 0.0;                            // Higher       DCH z
    1933}
    2034//
     
    2943        Double_t pt = p.Pt();
    3044        Double_t C = a / (2 * pt);                      // Half curvature
    31         //cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << endl;
     45        //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl;
    3246        Double_t r2 = x.Perp2();
    3347        Double_t cross = x(0) * p(1) - x(1) * p(0);
     
    7993TVector3 TrkUtil::ParToP(TVectorD Par)
    8094{
    81         if (fBz == 0.0)
    82 std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl;
     95        if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl;
    8396        //
    8497        return ParToP(Par,fBz);
     
    113126        Double_t b = -cSpeed() * fBz / 2.;
    114127        pACTS(0) = 1000 * Par(0);               // D from m to mm
    115         pACTS(1) = 1000 * Par(3);               // z0 from m to mm
     128        pACTS(1) = 1000 * Par(3);       // z0 from m to mm
    116129        pACTS(2) = Par(1);                      // Phi0 is unchanged
    117130        pACTS(3) = TMath::ATan2(1.0, Par(4));           // Theta in [0, pi] range
     
    136149        A(3, 1) = 1000.;                // z0-z0 conversion to mm
    137150        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)
     151        A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta)
    139152        //
    140153        TMatrixDSym Cv = Cov;
     
    218231        return Cmm;
    219232}
     233//
     234// Setup chamber volume
     235void TrkUtil::SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax)
     236{
     237        fRmin = Rmin;                           // Lower                DCH radius
     238        fRmax = Rmax;                           // Higher       DCH radius
     239        fZmin = Zmin;                           // Lower                DCH z
     240        fZmax = Zmax;                           // Higher       DCH z
     241}
     242//
     243// Get Trakck length inside DCH volume
     244Double_t TrkUtil::TrkLen(TVectorD Par)
     245{
     246        Double_t tLength = 0.0;
     247        // Check if geometry is initialized
     248        if (fZmin == 0.0 && fZmax == 0.0)
     249        {
     250                // No geometry set so send a warning and return 0
     251                std::cout << "TrkUtil::TrkLen() called without a DCH volume defined" << std::endl;
     252        }
     253        else
     254        {
     255                //******************************************************************
     256                // Determine the track length inside the chamber   ****
     257                //******************************************************************
     258                //
     259                // Track pararameters
     260                Double_t D = Par(0);            // Transverse impact parameter
     261                Double_t phi0 = Par(1);         // Transverse direction at minimum approach
     262                Double_t C = Par(2);            // Half curvature
     263                Double_t z0 = Par(3);           // Z at minimum approach
     264                Double_t ct = Par(4);           // cot(theta)
     265                //std::cout << "TrkUtil:: parameters: D= " << D << ", phi0= " << phi0
     266                //      << ", C= " << C << ", z0= " << z0 << ", ct= " << ct << std::endl;
     267                //
     268                // Track length per unit phase change
     269                Double_t Scale = TMath::Sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C));
     270                //
     271                // Find intersections with chamber boundaries
     272                //
     273                Double_t phRin = 0.0;                   // phase of inner cylinder
     274                Double_t phRin2= 0.0;                   // phase of inner cylinder intersection (2nd branch)
     275                Double_t phRhi = 0.0;                   // phase of outer cylinder intersection
     276                Double_t phZmn = 0.0;                   // phase of left wall intersection
     277                Double_t phZmx = 0.0;                   // phase of right wall intersection
     278                //  ... with inner cylinder
     279                Double_t Rtop = TMath::Abs((1.0 + C*D) / C);
     280
     281                if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment ***
     282                {
     283                        Double_t ph = 2 * TMath::ASin(C*TMath::Sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D)));
     284                        Double_t z = z0 + ct*ph / (2.0*C);
     285
     286                        //std::cout << "Rin intersection: ph = " << ph<<", z= "<<z << std::endl;
     287
     288                        if (z < fZmax && z > fZmin)     phRin = TMath::Abs(ph); // Intersection inside chamber volume   
     289                        //
     290                        // Include second branch of loopers
     291                        Double_t ph2 = TMath::TwoPi() - TMath::Abs(ph);
     292                        if (ph < 0)ph2 = -ph2;
     293                        z = z0 + ct * ph2 / (2.0 * C);
     294                        if (z < fZmax && z > fZmin)     phRin2 = TMath::Abs(ph2);       // Intersection inside chamber volume
     295                }
     296                //  ... with outer cylinder
     297                if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment ***
     298                {
     299                        Double_t ph = 2 * TMath::ASin(C*TMath::Sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D)));
     300                        Double_t z = z0 + ct*ph / (2.0*C);
     301                        if (z < fZmax && z > fZmin)     phRhi = TMath::Abs(ph); // Intersection inside chamber volume   
     302                }
     303                //  ... with left wall
     304                Double_t Zdir = (fZmin - z0) / ct;
     305                if (Zdir > 0.0)
     306                {
     307                        Double_t ph = 2.0*C*Zdir;
     308                        Double_t Rint = TMath::Sqrt(D*D + (1.0 + 2.0*C*D)*pow(TMath::Sin(ph / 2), 2) / (C*C));
     309                        if (Rint < fRmax && Rint > fRmin)       phZmn = TMath::Abs(ph); // Intersection inside chamber volume   
     310                }
     311                //  ... with right wall
     312                Zdir = (fZmax - z0) / ct;
     313                if (Zdir > 0.0)
     314                {
     315                        Double_t ph = 2.0*C*Zdir;
     316                        Double_t Rint = TMath::Sqrt(D*D + (1.0 + 2.0*C*D)*pow(TMath::Sin(ph / 2), 2) / (C*C));
     317                        if (Rint < fRmax && Rint > fRmin)       phZmx = TMath::Abs(ph); // Intersection inside chamber volume   
     318                }
     319                //
     320                // Order phases and keep the lowest two non-zero ones
     321                //
     322                const Int_t Nint = 5;
     323                Double_t dPhase = 0.0;  // Phase difference between two close intersections
     324                Double_t ph_arr[Nint] = { phRin, phRin2, phRhi, phZmn, phZmx };
     325                Int_t srtind[Nint];
     326                TMath::Sort(Nint, ph_arr, srtind, kFALSE);
     327                Int_t iPos = -1;                // First element > 0
     328                for (Int_t i = 0; i < Nint; i++)
     329                {
     330                        if (ph_arr[srtind[i]] <= 0.0) iPos = i;
     331                }
     332
     333                if (iPos < Nint - 2)
     334                {
     335                        dPhase = ph_arr[srtind[iPos + 2]] - ph_arr[srtind[iPos + 1]];
     336                        tLength = dPhase*Scale;
     337                }
     338        }
     339        return tLength;
     340}
     341//
     342// Return number of ionization clusters
     343Bool_t TrkUtil::IonClusters(Double_t &Ncl, Double_t mass, TVectorD Par)
     344{
     345        //
     346        // Units are meters/Tesla/GeV
     347        //
     348        Ncl = 0.0;
     349        Bool_t Signal = kFALSE;
     350        Double_t tLen = 0;
     351        // Check if geometry is initialized
     352        if (fZmin == 0.0 && fZmax == 0.0)
     353        {
     354                // No geometry set so send a warning and return 0
     355                std::cout << "TrkUtil::IonClusters() called without a volume defined" << std::endl;
     356        }
     357        else tLen = TrkLen(Par);
     358
     359        //******************************************************************
     360        // Now get the number of clusters                       ****
     361        //******************************************************************
     362        //
     363        Double_t muClu = 0.0;   // mean number of clusters
     364        Double_t bg = 0.0;              // beta*gamma
     365        Ncl = 0.0;
     366        if (tLen > 0.0)
     367        {
     368                Signal = kTRUE;
     369                //
     370                // Find beta*gamma
     371                if (fBz == 0.0)
     372                {
     373                        Signal = kFALSE;
     374                        std::cout << "TrkUtil::IonClusters: Please set Bz!!!" << std::endl;
     375                }
     376                else
     377                {
     378                        TVector3 p = ParToP(Par);
     379                        bg = p.Mag() / mass;
     380                        muClu = Nclusters(bg)*tLen;                             // Avg. number of clusters
     381
     382                        Ncl = gRandom->PoissonD(muClu);                 // Actual number of clusters
     383                }
     384
     385        }
     386//
     387        return Signal;
     388}
     389//
     390//
     391Double_t TrkUtil::Nclusters(Double_t begam)
     392{
     393        Int_t Opt = fGasSel;
     394        Double_t Nclu = Nclusters(begam, Opt);
     395        //
     396        return Nclu;
     397}
     398//
     399Double_t TrkUtil::Nclusters(Double_t begam, Int_t Opt) {
     400        //
     401        // Opt = 0: He 90 - Isobutane 10
     402        //     = 1: pure He
     403        //     = 2: Argon 50 - Ethane 50
     404        //     = 3: pure Argon
     405        //
     406        //
     407        std::vector<double> bg{ 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,
     408        12., 15., 20., 50., 100., 200., 500., 1000. };
     409        // He 90 - Isobutane 10
     410        std::vector<double> ncl_He_Iso{ 42.94, 23.6,18.97,12.98,12.2,12.13,
     411        12.24,12.73,13.03,13.29,13.63,14.08,15.56,16.43,16.8,16.95,16.98 };
     412        //
     413        // pure He
     414        std::vector<double> ncl_He{ 11.79,6.5,5.23,3.59,3.38,3.37,3.4,3.54,3.63,
     415        3.7,3.8,3.92,4.33,4.61,4.78,4.87,4.89 };
     416        //
     417        // Argon 50 - Ethane 50
     418        std::vector<double> ncl_Ar_Eth{ 130.04,71.55,57.56,39.44,37.08,36.9,
     419        37.25,38.76,39.68,40.49,41.53,42.91,46.8,48.09,48.59,48.85,48.93 };
     420        //
     421        // pure Argon
     422        std::vector<double> ncl_Ar{ 88.69,48.93,39.41,27.09,25.51,25.43,25.69,
     423        26.78,27.44,28.02,28.77,29.78,32.67,33.75,34.24,34.57,34.68 };
     424        //
     425        Int_t nPoints = (Int_t)bg.size();
     426        bg.push_back(10000.);
     427        std::vector<double> ncl;
     428        switch (Opt)
     429        {
     430        case 0: ncl = ncl_He_Iso;                       // He-Isobutane
     431                break;
     432        case 1: ncl = ncl_He;                           // pure He
     433                break;
     434        case 2: ncl = ncl_Ar_Eth;                       // Argon - Ethane
     435                break;
     436        case 3: ncl = ncl_Ar;                           // pure Argon
     437                break;
     438        }
     439        ncl.push_back(ncl[nPoints - 1]);
     440        Int_t ilow = 0;
     441        while (begam > bg[ilow])ilow++;
     442        ilow--;
     443        //std::cout << "ilow= " << ilow << ", low = " << bg[ilow] << ", val = " << begam
     444        //      << ", high = " << bg[ilow + 1] << std::endl;
     445        //
     446        Int_t ind[3] = { ilow, ilow + 1, ilow + 2 };
     447        TVectorD y(3);
     448        for (Int_t i = 0; i < 3; i++)y(i) = ncl[ind[i]];
     449        TVectorD x(3);
     450        for (Int_t i = 0; i < 3; i++)x(i) = bg[ind[i]];
     451        TMatrixD Xval(3, 3);
     452        for (Int_t i = 0; i < 3; i++)Xval(i, 0) = 1.0;
     453        for (Int_t i = 0; i < 3; i++)Xval(i, 1) = x(i);
     454        for (Int_t i = 0; i < 3; i++)Xval(i, 2) = x(i) * x(i);
     455        //std::cout << "Xval:" << std::endl; Xval.Print();
     456        Xval.Invert();
     457        TVectorD coeff = Xval * y;
     458        Double_t interp = coeff[0] + coeff[1] * begam + coeff[2] * begam * begam;
     459        //std::cout << "val1= (" <<x(0)<<", "<< y(0) << "), val2= ("
     460        //      <<x(1)<<", "<< y(1) << "), val3= ("
     461        //      <<x(2)<<", "<< y(2)
     462        //      << "), result= (" <<begam<<", "<< interp<<")" << std::endl;
     463        //
     464        if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << std::endl;
     465        if (begam < bg[0]) interp = 0.0;
     466        //std::cout << "bg= " << begam << ", Opt= " << Opt <<", interp = "<<interp<< std::endl;
     467        return 100*interp;
     468}
     469//
     470Double_t TrkUtil::funcNcl(Double_t *xp, Double_t *par){
     471        Double_t bg = xp[0];
     472        return Nclusters(bg);
     473}
     474//
     475void TrkUtil::SetGasMix(Int_t Opt)
     476{
     477        if (Opt < 0 || Opt > 3)
     478        {
     479                std::cout << "TrkUtil::SetGasMix Gas option not allowed. No action."
     480                        << std::endl;
     481        }
     482        else fGasSel = Opt;
     483}
Note: See TracChangeset for help on using the changeset viewer.