Fork me on GitHub

Changes in / [90132e0:56af73f] in git


Ignore:
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • DelphesEnv.sh

    r90132e0 r56af73f  
    1 DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )"
     1export PYTHONPATH=`pwd`/python:$PYTHONPATH
     2export LD_LIBRARY_PATH=`pwd`:$LD_LIBRARY_PATH
    23
    3 export DELPHES_HOME="$DIR"
    4 export PYTHONPATH="$DIR/python:${PYTHONPATH}"
    5 export LD_LIBRARY_PATH="$DIR:${LD_LIBRARY_PATH}"
    6 export LIBRARY_PATH="$DIR:${LIBRARY_PATH}"

  • Makefile

    r90132e0 r56af73f  
    232232        classes/DelphesClasses.h \
    233233        classes/DelphesFactory.h \
     234        classes/DelphesHepMC3Reader.h \
    234235        modules/Delphes.h \
    235236        external/ExRootAnalysis/ExRootProgressBar.h \

  • README

    r90132e0 r56af73f  
    1616Finally, we can run Delphes:
    1717
    18    ./DelphesHepMC3
     18   ./DelphesHepMC
    1919
    2020Command line parameters:
    2121
    22    ./DelphesHepMC3 config_file output_file [input_file(s)]
     22   ./DelphesHepMC config_file output_file [input_file(s)]
    2323     config_file - configuration file in Tcl format
    2424     output_file - output file in ROOT format,

  • README.md

    r90132e0 r56af73f  
    3030
    3131```
    32    ./DelphesHepMC3
     32   ./DelphesHepMC
    3333```
    3434
    … …  
    3636
    3737```
    38    ./DelphesHepMC3 config_file output_file [input_file(s)]
     38   ./DelphesHepMC config_file output_file [input_file(s)]
    3939     config_file - configuration file in Tcl format
    4040     output_file - output file in ROOT format,

  • classes/DelphesHepMC3Reader.cc

    r90132e0 r56af73f  
    5656DelphesHepMC3Reader::DelphesHepMC3Reader() :
    5757  fInputFile(0), fBuffer(0), fPDG(0),
    58   fVertexCounter(-2), fParticleCounter(-1)
     58  fVertexCounter(-1), fParticleCounter(-1)
    5959{
    6060  fBuffer = new char[kBufferSize];
    … …  
    8181void DelphesHepMC3Reader::Clear()
    8282{
    83   fWeights.clear();
     83  fWeight.clear();
    8484  fMomentumCoefficient = 1.0;
    8585  fPositionCoefficient = 1.0;
    86   fVertexCounter = -2;
     86  fVertexCounter = -1;
    8787  fParticleCounter = -1;
    88   fVertices.clear();
    89   fParticles.clear();
    90   fInVertexMap.clear();
    91   fOutVertexMap.clear();
    92   fMotherMap.clear();
     88  fVertexMap.clear();
    9389  fDaughterMap.clear();
    9490}
    … …  
    9894bool DelphesHepMC3Reader::EventReady()
    9995{
    100   return (fVertexCounter == -1) && (fParticleCounter == 0);
     96  return (fParticleCounter == 0);
    10197}
    10298
    … …  
    123119    Clear();
    124120
     121    fX = 0.0;
     122    fY = 0.0;
     123    fZ = 0.0;
     124    fT = 0.0;
     125
    125126    rc = bufferStream.ReadInt(fEventNumber)
    126127      && bufferStream.ReadInt(fVertexCounter)
    … …  
    167168    while(bufferStream.ReadDbl(weight))
    168169    {
    169       fWeights.push_back(weight);
     170      fWeight.push_back(weight);
    170171    }
    171172  }
    … …  
    256257  else if(key == 'V')
    257258  {
    258     fParticles.clear();
    259 
    260259    fX = 0.0;
    261260    fY = 0.0;
    … …  
    263262    fT = 0.0;
    264263
     264    fM1 = 0;
     265    fM2 = 0;
     266
    265267    rc = bufferStream.ReadInt(fVertexCode)
    266268      && bufferStream.ReadInt(fVertexStatus);
    … …  
    284286    while(bufferStream.ReadInt(code))
    285287    {
    286       fParticles.push_back(code);
    287       bufferStream.FindChr(',');
    288     }
     288      if(code < fM1 || fM1 == 0) fM1 = code;
     289      if(code > fM2) fM2 = code;
     290      fVertexMap[code] = fVertexCode;
     291    }
     292
     293    if(fM1 == fM2) fM2 = 0;
    289294
    290295    if(bufferStream.FindChr('@'))
    … …  
    302307      }
    303308    }
    304 
    305     AnalyzeVertex(factory, fVertexCode);
    306309  }
    307310  else if(key == 'P' && fParticleCounter > 0)
    … …  
    326329    }
    327330
    328     AnalyzeParticle(factory);
     331    itDaughterMap = fDaughterMap.find(fOutVertexCode);
     332    if(itDaughterMap == fDaughterMap.end())
     333    {
     334      fDaughterMap[fOutVertexCode] = make_pair(fParticleCode, fParticleCode);
     335    }
     336    else
     337    {
     338      itDaughterMap->second.second = fParticleCode;
     339    }
     340
     341    AnalyzeParticle(factory, allParticleOutputArray,
     342      stableParticleOutputArray, partonOutputArray);
    329343  }
    330344
    331345  if(EventReady())
    332346  {
    333     FinalizeParticles(allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
     347    FinalizeParticles(allParticleOutputArray);
    334348  }
    335349
    … …  
    349363  element->ProcessID = fProcessID;
    350364  element->MPI = fMPI;
    351   element->Weight = fWeights.size() > 0 ? fWeights[0] : 1.0;
     365  element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
    352366  element->CrossSection = fCrossSection;
    353367  element->CrossSectionError = fCrossSectionError;
    … …  
    373387{
    374388  Weight *element;
    375   vector<double>::const_iterator itWeights;
    376 
    377   for(itWeights = fWeights.begin(); itWeights != fWeights.end(); ++itWeights)
     389  vector<double>::const_iterator itWeight;
     390
     391  for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
    378392  {
    379393    element = static_cast<Weight *>(branch->NewEntry());
    380394
    381     element->Weight = *itWeights;
    382   }
    383 }
    384 
    385 //---------------------------------------------------------------------------
    386 
    387 void DelphesHepMC3Reader::AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate)
    388 {
    389   int index;
    390   TLorentzVector *position;
    391   TObjArray *array;
    392   vector<int>::iterator itParticle;
    393   map<int, int>::iterator itVertexMap;
    394 
    395   itVertexMap = fOutVertexMap.find(code);
    396   if(itVertexMap == fOutVertexMap.end())
    397   {
    398     --fVertexCounter;
    399 
    400     index = fVertices.size();
    401     fOutVertexMap[code] = index;
    402     if(candidate && code > 0) fInVertexMap[code] = index;
    403 
    404     position = factory->New<TLorentzVector>();
    405     array = factory->NewArray();
    406     position->SetXYZT(0.0, 0.0, 0.0, 0.0);
    407     fVertices.push_back(make_pair(position, array));
    408   }
    409   else
    410   {
    411     index = itVertexMap->second;
    412     position = fVertices[index].first;
    413     array = fVertices[index].second;
    414   }
    415 
    416   if(candidate)
    417   {
    418     array->Add(candidate);
    419   }
    420   else
    421   {
    422     position->SetXYZT(fX, fY, fZ, fT);
    423     for(itParticle = fParticles.begin(); itParticle != fParticles.end(); ++itParticle)
    424     {
    425       fInVertexMap[*itParticle] = index;
    426     }
    427   }
    428 }
    429 
    430 //---------------------------------------------------------------------------
    431 
    432 void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory)
    433 {
    434   Candidate *candidate;
    435 
    436   candidate = factory->NewCandidate();
    437 
    438   candidate->PID = fPID;
    439 
    440   candidate->Status = fParticleStatus;
    441 
    442   candidate->Mass = fMass;
    443 
    444   candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
    445 
    446   candidate->D1 = fParticleCode;
    447 
    448   AnalyzeVertex(factory, fOutVertexCode, candidate);
    449 }
    450 
    451 //---------------------------------------------------------------------------
    452 
    453 void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray,
     395    element->Weight = *itWeight;
     396  }
     397}
     398
     399//---------------------------------------------------------------------------
     400
     401void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory,
     402  TObjArray *allParticleOutputArray,
    454403  TObjArray *stableParticleOutputArray,
    455404  TObjArray *partonOutputArray)
    456405{
    457   TLorentzVector *position;
    458   TObjArray *array;
    459406  Candidate *candidate;
    460407  TParticlePDG *pdgParticle;
    461408  int pdgCode;
     409
     410  candidate = factory->NewCandidate();
     411
     412  candidate->PID = fPID;
     413  pdgCode = TMath::Abs(candidate->PID);
     414
     415  candidate->Status = fParticleStatus;
     416
     417  pdgParticle = fPDG->GetParticle(fPID);
     418  candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
     419  candidate->Mass = fMass;
     420
     421  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
     422  if(fMomentumCoefficient != 1.0)
     423  {
     424    candidate->Momentum *= fMomentumCoefficient;
     425  }
     426
     427  candidate->Position.SetXYZT(fX, fY, fZ, fT);
     428  if(fPositionCoefficient != 1.0)
     429  {
     430    candidate->Position *= fPositionCoefficient;
     431  }
     432
     433  candidate->D1 = -1;
     434  candidate->D2 = -1;
     435
     436  if(fOutVertexCode < 0)
     437  {
     438    candidate->M1 = fM1 - 1;
     439    candidate->M2 = fM2 - 1;
     440  }
     441  else
     442  {
     443    candidate->M1 = fOutVertexCode - 1;
     444    candidate->M2 = -1;
     445  }
     446
     447  allParticleOutputArray->Add(candidate);
     448
     449  if(!pdgParticle) return;
     450
     451  if(fParticleStatus == 1)
     452  {
     453    stableParticleOutputArray->Add(candidate);
     454  }
     455  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
     456  {
     457    partonOutputArray->Add(candidate);
     458  }
     459}
     460
     461//---------------------------------------------------------------------------
     462
     463void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray)
     464{
     465  Candidate *candidate;
    462466  map<int, int >::iterator itVertexMap;
    463   map<int, pair<int, int> >::iterator itMotherMap;
    464467  map<int, pair<int, int> >::iterator itDaughterMap;
    465   int i, j, code, counter;
    466 
    467   counter = 0;
    468   for(i = 0; i < fVertices.size(); ++i)
    469   {
    470     position = fVertices[i].first;
    471     array = fVertices[i].second;
    472 
    473     for(j = 0; j < array->GetEntriesFast(); ++j)
    474     {
    475       candidate = static_cast<Candidate *>(array->At(j));
    476 
    477       candidate->Position = *position;
    478       if(fPositionCoefficient != 1.0)
    479       {
    480         candidate->Position *= fPositionCoefficient;
    481       }
    482 
    483       if(fMomentumCoefficient != 1.0)
    484       {
    485         candidate->Momentum *= fMomentumCoefficient;
    486       }
    487 
    488       candidate->M1 = i;
    489 
    490       itDaughterMap = fDaughterMap.find(i);
    491       if(itDaughterMap == fDaughterMap.end())
    492       {
    493         fDaughterMap[i] = make_pair(counter, counter);
    494       }
    495       else
    496       {
    497         itDaughterMap->second.second = counter;
    498       }
    499 
    500       code = candidate->D1;
    501 
    502       itVertexMap = fInVertexMap.find(code);
    503       if(itVertexMap == fInVertexMap.end())
    504       {
    505         candidate->D1 = -1;
    506       }
    507       else
    508       {
    509         code = itVertexMap->second;
    510 
    511         candidate->D1 = code;
    512 
    513         itMotherMap = fMotherMap.find(code);
    514         if(itMotherMap == fMotherMap.end())
    515         {
    516           fMotherMap[code] = make_pair(counter, -1);
    517         }
    518         else
    519         {
    520           itMotherMap->second.second = counter;
    521         }
    522       }
    523 
    524       allParticleOutputArray->Add(candidate);
    525 
    526       ++counter;
    527 
    528       pdgParticle = fPDG->GetParticle(candidate->PID);
    529 
    530       candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
    531 
    532       if(!pdgParticle) continue;
    533 
    534       pdgCode = TMath::Abs(candidate->PID);
    535 
    536       if(candidate->Status == 1)
    537       {
    538         stableParticleOutputArray->Add(candidate);
    539       }
    540       else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
    541       {
    542         partonOutputArray->Add(candidate);
    543       }
    544     }
    545   }
     468  int i, index;
    546469
    547470  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
    … …  
    549472    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
    550473
    551     itMotherMap = fMotherMap.find(candidate->M1);
    552     if(itMotherMap == fMotherMap.end())
    553     {
    554       candidate->M1 = -1;
    555       candidate->M2 = -1;
    556     }
    557     else
    558     {
    559       candidate->M1 = itMotherMap->second.first;
    560       candidate->M2 = itMotherMap->second.second;
    561     }
    562 
    563     if(candidate->D1 < 0)
     474    index = i + 1;
     475
     476    itVertexMap = fVertexMap.find(index);
     477    if(itVertexMap != fVertexMap.end())
     478    {
     479      index = itVertexMap->second;
     480    }
     481
     482    itDaughterMap = fDaughterMap.find(index);
     483    if(itDaughterMap == fDaughterMap.end())
    564484    {
    565485      candidate->D1 = -1;
    … …  
    568488    else
    569489    {
    570       itDaughterMap = fDaughterMap.find(candidate->D1);
    571       if(itDaughterMap == fDaughterMap.end())
    572       {
    573         candidate->D1 = -1;
    574         candidate->D2 = -1;
    575       }
    576       else
    577       {
    578         candidate->D1 = itDaughterMap->second.first;
    579         candidate->D2 = itDaughterMap->second.second;
    580       }
    581     }
    582   }
    583 }
    584 
    585 //---------------------------------------------------------------------------
     490      candidate->D1 = itDaughterMap->second.first - 1;
     491      candidate->D2 = itDaughterMap->second.second - 1;
     492    }
     493  }
     494}
     495
     496//---------------------------------------------------------------------------

  • classes/DelphesHepMC3Reader.h

    r90132e0 r56af73f  
    3636class TStopwatch;
    3737class TDatabasePDG;
    38 class TLorentzVector;
    3938class ExRootTreeBranch;
    4039class DelphesFactory;
    41 class Candidate;
    4240
    4341class DelphesHepMC3Reader
    … …  
    6361
    6462private:
    65   void AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate = 0);
    66 
    67   void AnalyzeParticle(DelphesFactory *factory);
    68 
    69   void FinalizeParticles(TObjArray *allParticleOutputArray,
     63  void AnalyzeParticle(DelphesFactory *factory,
     64    TObjArray *allParticleOutputArray,
    7065    TObjArray *stableParticleOutputArray,
    7166    TObjArray *partonOutputArray);
     67
     68  void FinalizeParticles(TObjArray *allParticleOutputArray);
    7269
    7370  FILE *fInputFile;
    … …  
    8279  double fMomentumCoefficient, fPositionCoefficient;
    8380
    84   std::vector<double> fWeights;
     81  std::vector<double> fWeight;
    8582
    8683  double fCrossSection, fCrossSectionError;
    … …  
    9592  double fPx, fPy, fPz, fE, fMass;
    9693
    97   std::vector<std::pair<TLorentzVector *, TObjArray *> > fVertices;
    98   std::vector<int> fParticles;
     94  int fM1, fM2;
    9995
    100   std::map<int, int> fInVertexMap;
    101   std::map<int, int> fOutVertexMap;
    102 
    103   std::map<int, std::pair<int, int> > fMotherMap;
     96  std::map<int, int> fVertexMap;
    10497  std::map<int, std::pair<int, int> > fDaughterMap;
    10598};

  • modules/Delphes.cc

    r90132e0 r56af73f  
    6161  fFactory(0)
    6262{
    63   TFolder *folder;
    64 
     63  TFolder *folder = new TFolder(name, "");
    6564  fFactory = new DelphesFactory("ObjectFactory");
    66 
    67   folder = new TFolder(name, "");
    6865
    6966  SetName(name);
    … …  
    8380  if(folder)
    8481  {
    85     gROOT->GetListOfBrowsables()->Remove(folder);
    8682    folder->Clear();
    8783    delete folder;

  • modules/ParticlePropagator.cc

    r90132e0 r56af73f  
    126126  Double_t px, py, pz, pt, pt2, e, q;
    127127  Double_t x, y, z, t, r;
    128   Double_t x_c, y_c, r_c, phi_0;
    129   Double_t x_t, y_t, z_t, r_t, phi_t;
    130   Double_t t_r, t_z;
    131   Double_t tmp;
     128  Double_t x_c, y_c, r_c, phi_c, phi_0;
     129  Double_t x_t, y_t, z_t, r_t;
     130  Double_t t_z, t_r;
     131  Double_t discr;
    132132  Double_t gammam, omega;
    133133  Double_t xd, yd, zd;
    134134  Double_t l, d0, dz, ctgTheta, alpha;
    135135  Double_t bsx, bsy, bsz;
    136   Double_t td, pio, phid, vz;
     136  Double_t rxp, rdp, t_R;
     137  Double_t td, pio, phid, sign_pz, vz;
    137138
    138139  const Double_t c_light = 2.99792458E8;
    139140
    140141  if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
    141   {
    142142    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
    143   }
    144143  else
    145144  {
    … …  
    162161    particlePosition = particle->Position;
    163162    particleMomentum = particle->Momentum;
     163
     164    // Constants
    164165
    165166    x = particlePosition.X() * 1.0E-3;
    … …  
    207208    else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
    208209    {
    209       // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
    210       tmp = px * y - py * x;
    211       t_r = (TMath::Sqrt(pt2 * fRadius2 - tmp * tmp) - px * x - py * y) / pt2;
    212 
    213       t_z = (TMath::Sign(fHalfLength, pz) - z) / pz;
    214 
    215       t = TMath::Min(t_r, t_z);
    216 
    217       x_t = x + px * t;
    218       y_t = y + py * t;
    219       z_t = z + pz * t;
    220 
    221       l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z));
     210
     211      rxp = x*py - y*px;
     212      rdp = x*px + y*py;
     213
     214      discr = fRadius*fRadius*pt*pt - rxp*rxp;
     215
     216      t_R = e * (sqrt(discr) - rdp) / (c_light * pt * pt);
     217      t_z = e * (TMath::Sign(fHalfLengthMax, pz) - z) / ( c_light * pz);
     218
     219      t = TMath::Min(t_R, t_z);
     220
     221      x_t = x + px*t*c_light/e;
     222      y_t = y + py*t*c_light/e;
     223      z_t = z + pz*t*c_light/e;
     224      r_t = TMath::Hypot(x_t, y_t);
     225
     226      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
    222227
    223228      mother = candidate;
    224       candidate = static_cast<Candidate *>(candidate->Clone());
     229      candidate = static_cast<Candidate*>(candidate->Clone());
    225230
    226231      candidate->InitialPosition = particlePosition;
    227       candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3);
    228       candidate->L = l * 1.0E3;
     232      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*c_light*1.0E3);
     233      candidate->L = l*1.0E3;
    229234
    230235      candidate->Momentum = particleMomentum;
    … …  
    255260    {
    256261
    257       // 1. initial transverse momentum p_{T0}: Part->pt
    258       //    initial transverse momentum direction phi_0 = -atan(p_{X0} / p_{Y0})
    259       //    relativistic gamma: gamma = E / mc^2; gammam = gamma * m
    260       //    gyration frequency omega = q * Bz / (gammam)
    261       //    helix radius r = p_{T0} / (omega * gammam)
    262 
    263       gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
    264       omega = q * fBz / gammam; // omega is here in [89875518/s]
    265       r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
     262      // 1.  initial transverse momentum p_{T0}: Part->pt
     263      //     initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)
     264      //     relativistic gamma: gamma = E/mc^2; gammam = gamma * m
     265      //     gyration frequency omega = q/(gamma m) fBz
     266      //     helix radius r = p_{T0} / (omega gamma m)
     267
     268      gammam = e*1.0E9 / (c_light*c_light);      // gammam in [eV/c^2]
     269      omega = q * fBz / (gammam);                // omega is here in [89875518/s]
     270      r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
    266271
    267272      phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
    268273
    269274      // 2. helix axis coordinates
    270       x_c = x + r * TMath::Sin(phi_0);
    271       y_c = y - r * TMath::Cos(phi_0);
     275      x_c = x + r*TMath::Sin(phi_0);
     276      y_c = y - r*TMath::Cos(phi_0);
    272277      r_c = TMath::Hypot(x_c, y_c);
    273 
    274       // time of closest approach
    275       td = (phi_0 + TMath::ATan2(x_c, y_c)) / omega;
    276 
    277       // remove all the modulo pi that might have come from the atan
    278       pio = TMath::Abs(TMath::Pi() / omega);
    279       while(TMath::Abs(td) > 0.5 * pio)
    280       {
    281         td -= TMath::Sign(1.0, td) * pio;
    282       }
    283 
    284       vz = pz * c_light / e;
    285 
    286       // calculate coordinates of closest approach to z axis
    287       phid = phi_0 - omega * td;
    288       xd = x_c - r * TMath::Sin(phid);
    289       yd = y_c + r * TMath::Cos(phid);
    290       zd = z + vz * td;
    291 
    292       // momentum at closest approach
    293       px = pt * TMath::Cos(phid);
    294       py = pt * TMath::Sin(phid);
     278      phi_c = TMath::ATan(y_c/x_c);
     279      if(x_c < 0.0) phi_c -= TMath::Sign(1., phi_c)*TMath::Pi();
     280
     281      //Find the time of closest approach
     282      td = (phi_0 - TMath::ATan(-x_c/y_c))/omega;
     283
     284      //Remove all the modulo pi that might have come from the atan
     285      pio = fabs(TMath::Pi()/omega);
     286      while(fabs(td) > 0.5*pio)
     287      {
     288        td -= TMath::Sign(1., td)*pio;
     289      }
     290
     291      //Compute the coordinate of closed approach to z axis
     292      //if wants wtr beamline need to be changedto re-center with a traslation of the z axis
     293      phid = phi_0 - omega*td;
     294      xd = x_c - r*TMath::Sin(phid);
     295      yd = y_c + r*TMath::Cos(phid);
     296      zd = z + c_light*(pz/e)*td;
     297
     298      //Compute momentum at closest approach (perigee??)
     299      px = pt*TMath::Cos(phid);
     300      py = pt*TMath::Sin(phid);
    295301
    296302      particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E());
    … …  
    299305      d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
    300306      dz = zd - bsz;
    301       ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
     307      ctgTheta  = 1.0 / TMath::Tan (particleMomentum.Theta());
    302308
    303309      // 3. time evaluation t = TMath::Min(t_r, t_z)
    304310      //    t_r : time to exit from the sides
    305311      //    t_z : time to exit from the front or the back
    306       t_z = (vz == 0.0) ? 1.0E99 : (TMath::Sign(fHalfLength, pz) - z) / vz;
    307 
    308       if(r_c + TMath::Abs(r) < fRadius)
    309       {
    310         // helix does not cross the cylinder sides
     312      t = 0;
     313      t_z = 0;
     314      sign_pz = (pz > 0.0) ? 1 : -1;
     315      if(pz == 0.0) t_z = 1.0E99;
     316      else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
     317
     318      if(r_c + TMath::Abs(r)  < fRadius)   // helix does not cross the cylinder sides
     319      {
    311320        t = t_z;
    312321      }
    313322      else
    314323      {
    315         alpha = TMath::ACos((r * r + r_c * r_c - fRadius * fRadius) / (2 * TMath::Abs(r) * r_c));
    316         t_r = td + TMath::Abs(alpha / omega);
     324        alpha = -(fRadius*fRadius - r*r - r_c*r_c)/(2*fabs(r)*r_c);
     325        alpha = fabs(TMath::ACos(alpha));
     326        t_r = td + alpha/fabs(omega);
    317327
    318328        t = TMath::Min(t_r, t_z);
    319329      }
    320330
    321       // 4. position in terms of x(t), y(t), z(t)
    322       phi_t = phi_0 - omega * t;
    323       x_t = x_c - r * TMath::Sin(phi_t);
    324       y_t = y_c + r * TMath::Cos(phi_t);
    325       z_t = z + vz * t;
    326       r_t = TMath::Hypot(x_t, y_t);
    327 
    328       // lenght of the path from production to tracker
    329       l = t * TMath::Hypot(vz, r * omega);
     331      x_t = x_c - r*TMath::Sin(phi_0 - omega*t);
     332      y_t = y_c + r*TMath::Cos(phi_0 - omega*t);
     333      z_t = z + c_light*t*pz/e;
     334      r_t =  TMath::Hypot(x_t, y_t);
     335
     336      // compute path length for an helix
     337      vz = pz*1.0E9 / c_light / gammam;
     338      //lenght of the path from production to tracker
     339      l = t * TMath::Sqrt(vz*vz + r*r*omega*omega);
    330340
    331341      if(r_t > 0.0)
Note: See TracChangeset for help on using the changeset viewer.