Fork me on GitHub

Changeset 341014c in git for modules/ParticlePropagator.cc


Ignore:
Timestamp:
Feb 12, 2019, 9:29:17 PM (6 years ago)
Author:
Pavel Demin <pavel-demin@…>
Branches:
ImprovedOutputFile, Timing, llp, master
Children:
6455202
Parents:
45e58be
Message:

apply .clang-format to all .h, .cc and .cpp files

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/ParticlePropagator.cc

    r45e58be r341014c  
    1717 */
    1818
    19 
    2019/** \class ParticlePropagator
    2120 *
     
    3534#include "classes/DelphesFormula.h"
    3635
     36#include "ExRootAnalysis/ExRootClassifier.h"
     37#include "ExRootAnalysis/ExRootFilter.h"
    3738#include "ExRootAnalysis/ExRootResult.h"
    38 #include "ExRootAnalysis/ExRootFilter.h"
    39 #include "ExRootAnalysis/ExRootClassifier.h"
    40 
     39
     40#include "TDatabasePDG.h"
     41#include "TFormula.h"
     42#include "TLorentzVector.h"
    4143#include "TMath.h"
     44#include "TObjArray.h"
     45#include "TRandom3.h"
    4246#include "TString.h"
    43 #include "TFormula.h"
    44 #include "TRandom3.h"
    45 #include "TObjArray.h"
    46 #include "TDatabasePDG.h"
    47 #include "TLorentzVector.h"
    4847
    4948#include <algorithm>
    50 #include <stdexcept>
    5149#include <iostream>
    5250#include <sstream>
     51#include <stdexcept>
    5352
    5453using namespace std;
     
    6766}
    6867
    69 
    7068//------------------------------------------------------------------------------
    7169
     
    7371{
    7472  fRadius = GetDouble("Radius", 1.0);
    75   fRadius2 = fRadius*fRadius;
     73  fRadius2 = fRadius * fRadius;
    7674  fHalfLength = GetDouble("HalfLength", 3.0);
    7775  fBz = GetDouble("Bz", 0.0);
     
    140138  const Double_t c_light = 2.99792458E8;
    141139
    142   if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0)
     140  if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
    143141    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
    144142  else
    145143  {
    146     Candidate &beamSpotCandidate = *((Candidate *) fBeamSpotInputArray->At(0));
     144    Candidate &beamSpotCandidate = *((Candidate *)fBeamSpotInputArray->At(0));
    147145    beamSpotPosition = beamSpotCandidate.Position;
    148146  }
    149147
    150148  fItInputArray->Reset();
    151   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
     149  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    152150  {
    153151    if(candidate->GetCandidates()->GetEntriesFast() == 0)
     
    157155    else
    158156    {
    159       particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
     157      particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    160158    }
    161159
    162160    particlePosition = particle->Position;
    163161    particleMomentum = particle->Momentum;
    164     x = particlePosition.X()*1.0E-3;
    165     y = particlePosition.Y()*1.0E-3;
    166     z = particlePosition.Z()*1.0E-3;
    167 
    168     bsx = beamSpotPosition.X()*1.0E-3;
    169     bsy = beamSpotPosition.Y()*1.0E-3;
    170     bsz = beamSpotPosition.Z()*1.0E-3;
     162    x = particlePosition.X() * 1.0E-3;
     163    y = particlePosition.Y() * 1.0E-3;
     164    z = particlePosition.Z() * 1.0E-3;
     165
     166    bsx = beamSpotPosition.X() * 1.0E-3;
     167    bsy = beamSpotPosition.Y() * 1.0E-3;
     168    bsz = beamSpotPosition.Z() * 1.0E-3;
    171169
    172170    q = particle->Charge;
     
    193191    {
    194192      mother = candidate;
    195       candidate = static_cast<Candidate*>(candidate->Clone());
     193      candidate = static_cast<Candidate *>(candidate->Clone());
    196194
    197195      candidate->InitialPosition = particlePosition;
     
    207205    {
    208206      // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
    209       tmp = px*y - py*x;
    210       discr2 = pt2*fRadius2 - tmp*tmp;
     207      tmp = px * y - py * x;
     208      discr2 = pt2 * fRadius2 - tmp * tmp;
    211209
    212210      if(discr2 < 0.0)
     
    216214      }
    217215
    218       tmp = px*x + py*y;
     216      tmp = px * x + py * y;
    219217      discr = TMath::Sqrt(discr2);
    220       t1 = (-tmp + discr)/pt2;
    221       t2 = (-tmp - discr)/pt2;
     218      t1 = (-tmp + discr) / pt2;
     219      t2 = (-tmp - discr) / pt2;
    222220      t = (t1 < 0.0) ? t2 : t1;
    223221
    224       z_t = z + pz*t;
     222      z_t = z + pz * t;
    225223      if(TMath::Abs(z_t) > fHalfLength)
    226224      {
     
    230228      }
    231229
    232       x_t = x + px*t;
    233       y_t = y + py*t;
    234       z_t = z + pz*t;
    235 
    236       l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
     230      x_t = x + px * t;
     231      y_t = y + py * t;
     232      z_t = z + pz * t;
     233
     234      l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z));
    237235
    238236      mother = candidate;
    239       candidate = static_cast<Candidate*>(candidate->Clone());
     237      candidate = static_cast<Candidate *>(candidate->Clone());
    240238
    241239      candidate->InitialPosition = particlePosition;
    242       candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*e*1.0E3);
    243       candidate->L = l*1.0E3;
     240      candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3);
     241      candidate->L = l * 1.0E3;
    244242
    245243      candidate->Momentum = particleMomentum;
     
    251249        switch(TMath::Abs(candidate->PID))
    252250        {
    253           case 11:
    254             fElectronOutputArray->Add(candidate);
    255             break;
    256           case 13:
    257             fMuonOutputArray->Add(candidate);
    258             break;
    259           default:
    260             fChargedHadronOutputArray->Add(candidate);
     251        case 11:
     252          fElectronOutputArray->Add(candidate);
     253          break;
     254        case 13:
     255          fMuonOutputArray->Add(candidate);
     256          break;
     257        default:
     258          fChargedHadronOutputArray->Add(candidate);
    261259        }
    262260      }
    263261      else
    264262      {
    265          fNeutralOutputArray->Add(candidate);
     263        fNeutralOutputArray->Add(candidate);
    266264      }
    267265    }
     
    275273      //     helix radius r = p_{T0} / (omega gamma m)
    276274
    277       gammam = e*1.0E9 / (c_light*c_light);      // gammam in [eV/c^2]
    278       omega = q * fBz / (gammam);                // omega is here in [89875518/s]
    279       r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
     275      gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
     276      omega = q * fBz / (gammam); // omega is here in [89875518/s]
     277      r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
    280278
    281279      phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
    282280
    283281      // 2. helix axis coordinates
    284       x_c = x + r*TMath::Sin(phi_0);
    285       y_c = y - r*TMath::Cos(phi_0);
     282      x_c = x + r * TMath::Sin(phi_0);
     283      y_c = y - r * TMath::Cos(phi_0);
    286284      r_c = TMath::Hypot(x_c, y_c);
    287285      phi_c = TMath::ATan2(y_c, x_c);
     
    290288
    291289      rcu = TMath::Abs(r);
    292       rc2 = r_c*r_c;
     290      rc2 = r_c * r_c;
    293291
    294292      // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
    295       xd = x_c*x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;
     293      xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c;
    296294      xd = (rc2 > 0.0) ? xd / rc2 : -999;
    297       yd = y_c*(-rcu*r_c + rc2);
     295      yd = y_c * (-rcu * r_c + rc2);
    298296      yd = (rc2 > 0.0) ? yd / rc2 : -999;
    299       zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
     297      zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt;
    300298
    301299      // use perigee momentum rather than original particle
     
    311309      // calculate additional track parameters (correct for beamspot position)
    312310
    313       d0        = ((x - bsx) * py - (y - bsy) * px) / pt;
    314       dz        = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
    315       p         = particleMomentum.P();
    316       ctgTheta  = 1.0 / TMath::Tan (particleMomentum.Theta());
    317 
     311      d0 = ((x - bsx) * py - (y - bsy) * px) / pt;
     312      dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
     313      p = particleMomentum.P();
     314      ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
    318315
    319316      // 3. time evaluation t = TMath::Min(t_r, t_z)
     
    322319      t_r = 0.0; // in [ns]
    323320      int sign_pz = (pz > 0.0) ? 1 : -1;
    324       if(pz == 0.0) t_z = 1.0E99;
    325       else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
    326 
    327       if(r_c + TMath::Abs(r)  < fRadius)
     321      if(pz == 0.0)
     322        t_z = 1.0E99;
     323      else
     324        t_z = gammam / (pz * 1.0E9 / c_light) * (-z + fHalfLength * sign_pz);
     325
     326      if(r_c + TMath::Abs(r) < fRadius)
    328327      {
    329328        // helix does not cross the cylinder sides
     
    332331      else
    333332      {
    334         asinrho = TMath::ASin((fRadius*fRadius - r_c*r_c - r*r) / (2*TMath::Abs(r)*r_c));
     333        asinrho = TMath::ASin((fRadius * fRadius - r_c * r_c - r * r) / (2 * TMath::Abs(r) * r_c));
    335334        delta = phi_0 - phi;
    336         if(delta <-TMath::Pi()) delta += 2*TMath::Pi();
    337         if(delta > TMath::Pi()) delta -= 2*TMath::Pi();
     335        if(delta < -TMath::Pi()) delta += 2 * TMath::Pi();
     336        if(delta > TMath::Pi()) delta -= 2 * TMath::Pi();
    338337        t1 = (delta + asinrho) / omega;
    339338        t2 = (delta + TMath::Pi() - asinrho) / omega;
     
    359358      x_t = x_c + r * TMath::Sin(omega * t - phi_0);
    360359      y_t = y_c + r * TMath::Cos(omega * t - phi_0);
    361       z_t = z + pz*1.0E9 / c_light / gammam * t;
     360      z_t = z + pz * 1.0E9 / c_light / gammam * t;
    362361      r_t = TMath::Hypot(x_t, y_t);
    363362
    364 
    365363      // compute path length for an helix
    366364
    367       alpha = pz*1.0E9 / c_light / gammam;
    368       l = t * TMath::Sqrt(alpha*alpha + r*r*omega*omega);
     365      alpha = pz * 1.0E9 / c_light / gammam;
     366      l = t * TMath::Sqrt(alpha * alpha + r * r * omega * omega);
    369367
    370368      if(r_t > 0.0)
     
    374372        if(particle == candidate)
    375373        {
    376           particle->D0 = d0*1.0E3;
    377           particle->DZ = dz*1.0E3;
     374          particle->D0 = d0 * 1.0E3;
     375          particle->DZ = dz * 1.0E3;
    378376          particle->P = p;
    379377          particle->PT = pt;
     
    383381
    384382        mother = candidate;
    385         candidate = static_cast<Candidate*>(candidate->Clone());
     383        candidate = static_cast<Candidate *>(candidate->Clone());
    386384
    387385        candidate->InitialPosition = particlePosition;
    388         candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*c_light*1.0E3);
     386        candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * c_light * 1.0E3);
    389387
    390388        candidate->Momentum = particleMomentum;
    391389
    392         candidate->L = l*1.0E3;
    393 
    394         candidate->Xd = xd*1.0E3;
    395         candidate->Yd = yd*1.0E3;
    396         candidate->Zd = zd*1.0E3;
     390        candidate->L = l * 1.0E3;
     391
     392        candidate->Xd = xd * 1.0E3;
     393        candidate->Yd = yd * 1.0E3;
     394        candidate->Zd = zd * 1.0E3;
    397395
    398396        candidate->AddCandidate(mother);
     
    401399        switch(TMath::Abs(candidate->PID))
    402400        {
    403           case 11:
    404             fElectronOutputArray->Add(candidate);
    405             break;
    406           case 13:
    407             fMuonOutputArray->Add(candidate);
    408             break;
    409           default:
    410             fChargedHadronOutputArray->Add(candidate);
     401        case 11:
     402          fElectronOutputArray->Add(candidate);
     403          break;
     404        case 13:
     405          fMuonOutputArray->Add(candidate);
     406          break;
     407        default:
     408          fChargedHadronOutputArray->Add(candidate);
    411409        }
    412410      }
     
    416414
    417415//------------------------------------------------------------------------------
    418 
Note: See TracChangeset for help on using the changeset viewer.