Fork me on GitHub

Changeset 187fc41 in git


Ignore:
Timestamp:
May 2, 2016, 5:33:35 PM (9 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
1dad056
Parents:
2dc1092
Message:

added BeamSpot array in propagator, and calculation of track parameters

Location:
modules
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • modules/ParticlePropagator.cc

    r2dc1092 r187fc41  
    6666{
    6767}
    68 
    69 
    70 //------------------------------------------------------------------------------
    71 
    72 TLorentzVector ParticlePropagator::BeamSpotPosition(const TObjArray *array)
    73 {
    74   Candidate *candidate;
    75   Bool_t passed = false;
    76   TLorentzVector bs;
    77  
    78   fItInputArray->Reset();
    79   while((candidate = static_cast<Candidate*>(fItInputArray->Next())) && !passed)
    80   {
    81     if(candidate->IsPU == 0) passed = true;
    82     bs = candidate->Position;
    83   }
    84  
    85   return bs;
    86  
    87 }
    88 
    8968
    9069
     
    11392  fItInputArray = fInputArray->MakeIterator();
    11493
     94  // import beamspot
     95 
     96  fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
     97
    11598  // create output arrays
    11699
     
    133116{
    134117  Candidate *candidate, *mother;
    135   TLorentzVector candidatePosition, candidateMomentum;
     118  TLorentzVector candidatePosition, candidateMomentum, beamSpotPosition;
    136119  Double_t px, py, pz, pt, pt2, e, q;
    137120  Double_t x, y, z, t, r, phi;
     
    142125  Double_t tmp, discr, discr2;
    143126  Double_t delta, gammam, omega, asinrho;
    144   Double_t rcu, rc2, dxy, xd, yd, zd;
    145 
     127  Double_t rcu, rc2, xd, yd, zd;
     128  Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
     129  Double_t bsx, bsy, bsz;
     130         
    146131  const Double_t c_light = 2.99792458E8;
    147 
    148   const TLorentzVector &bs = BeamSpotPosition(fInputArray);
    149 
     132 
     133  if (!fBeamSpotInputArray->GetSize () || !fBeamSpotInputArray->At(0))
     134    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
     135  else
     136  {
     137    Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At(0));
     138    beamSpotPosition = beamSpotCandidate.Position;
     139  }
     140 
    150141  fItInputArray->Reset();
    151142  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
     
    156147    y = candidatePosition.Y()*1.0E-3;
    157148    z = candidatePosition.Z()*1.0E-3;
     149   
     150    bsx = beamSpotPosition.X()*1.0E-3;
     151    bsy = beamSpotPosition.Y()*1.0E-3;
     152    bsz = beamSpotPosition.Z()*1.0E-3;
     153   
    158154    q = candidate->Charge;
    159155
     
    205201      y_t = y + py*t;
    206202      z_t = z + pz*t;
     203     
     204      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
    207205
    208206      mother = candidate;
     
    210208
    211209      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
    212 
     210      candidate->L = l*1.0E3;
     211   
    213212      candidate->Momentum = candidateMomentum;
    214213      candidate->AddCandidate(mother);
     
    262261      yd = (rc2 > 0.0) ? yd / rc2 : -999;
    263262      zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
    264 
    265       // calculate impact paramater
    266       dxy = (xd*py - yd*px)/pt;
    267 
     263     
     264      // use perigee momentum rather than original particle
     265      // momentum, since the orignal particle momentum isn't known
     266     
     267      px = TMath::Sign(1.0,r) * pt * (-y_c / r_c);
     268      py = TMath::Sign(1.0,r) * pt * (x_c / r_c);
     269      etap = candidateMomentum.Eta();
     270      phip = TMath::ATan2(py, px);
     271       
     272      candidateMomentum.SetPtEtaPhiE(pt, etap, phip, candidateMomentum.E());
     273     
     274      // calculate additional track parameters (correct for beamspot position)
     275       
     276      d0        = (  (x - bsx) * py - (y - bsy) * px) / pt;
     277      dz        = z - ((x - bsx) * px + (py - bsy) * py) / pt * (pz / pt);
     278      p         = candidateMomentum.P();
     279      ctgTheta  = 1.0 / TMath::Tan (candidateMomentum.Theta ());
     280         
    268281      // 3. time evaluation t = TMath::Min(t_r, t_z)
    269282      //    t_r : time to exit from the sides
     
    311324      r_t = TMath::Hypot(x_t, y_t);
    312325
     326     
     327      // compute path length for an helix
     328     
     329      alpha = pz*1.0E9 / c_light / gammam;
     330      l = t * TMath::Sqrt(alpha*alpha + r*r*omega*omega);
     331         
    313332      if(r_t > 0.0)
    314333      {
     
    319338
    320339        candidate->Momentum = candidateMomentum;
    321         candidate->Dxy = dxy*1.0E3;
    322         candidate->Xd = xd*1.0E3;
     340       
     341        candidate->L =  l*1.0E3;
     342        candidate->D0 = d0*1.0E3;
     343        candidate->DZ = dz*1.0E3;
     344        candidate->P = p;
     345        candidate->PT = pt;
     346        candidate->CtgTheta = ctgTheta;
     347        candidate->Phi = phi;
     348
     349        candidate->Xd = xd*1.0E3;
    323350        candidate->Yd = yd*1.0E3;
    324351        candidate->Zd = zd*1.0E3;
  • modules/ParticlePropagator.h

    r2dc1092 r187fc41  
    5656
    5757  const TObjArray *fInputArray; //!
     58  const TObjArray *fBeamSpotInputArray; //!
    5859
    5960  TObjArray *fOutputArray; //!
     
    6263  TObjArray *fMuonOutputArray; //!
    6364
    64   TLorentzVector BeamSpotPosition(const TObjArray *array);
    65 
    6665  ClassDef(ParticlePropagator, 1)
    6766};
Note: See TracChangeset for help on using the changeset viewer.