Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/PhotonConversions.cc

    r341014c rbdfe7f2  
    1717 */
    1818
     19
    1920/** \class PhotonConversions
    2021 *
     
    2930
    3031#include "classes/DelphesClasses.h"
     32#include "classes/DelphesFactory.h"
    3133#include "classes/DelphesCylindricalFormula.h"
    32 #include "classes/DelphesFactory.h"
    33 
     34
     35#include "ExRootAnalysis/ExRootResult.h"
     36#include "ExRootAnalysis/ExRootFilter.h"
    3437#include "ExRootAnalysis/ExRootClassifier.h"
    35 #include "ExRootAnalysis/ExRootFilter.h"
    36 #include "ExRootAnalysis/ExRootResult.h"
    37 
     38
     39#include "TMath.h"
     40#include "TF1.h"
     41#include "TString.h"
     42#include "TFormula.h"
     43#include "TRandom3.h"
     44#include "TObjArray.h"
    3845#include "TDatabasePDG.h"
    39 #include "TF1.h"
    40 #include "TFormula.h"
    4146#include "TLorentzVector.h"
    42 #include "TMath.h"
    43 #include "TObjArray.h"
    44 #include "TRandom3.h"
    45 #include "TString.h"
    4647#include "TVector3.h"
    4748
    4849#include <algorithm>
     50#include <stdexcept>
    4951#include <iostream>
    5052#include <sstream>
    51 #include <stdexcept>
    5253
    5354using namespace std;
     
    5859  fItInputArray(0), fConversionMap(0), fDecayXsec(0)
    5960{
    60   fDecayXsec = new TF1("decayXsec", "1.0 - 4.0/3.0 * x * (1.0 - x)", 0.0, 1.0);
     61  fDecayXsec = new TF1("decayXsec","1.0 - 4.0/3.0 * x * (1.0 - x)", 0.0, 1.0);
    6162  fConversionMap = new DelphesCylindricalFormula;
    6263}
     
    7374{
    7475  fRadius = GetDouble("Radius", 1.0);
    75   fRadius2 = fRadius * fRadius;
     76  fRadius2 = fRadius*fRadius;
    7677
    7778  fHalfLength = GetDouble("HalfLength", 3.0);
     
    121122
    122123  fItInputArray->Reset();
    123   while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
     124  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    124125  {
    125126
     
    132133      candidatePosition = candidate->Position;
    133134      candidateMomentum = candidate->Momentum;
    134       x = candidatePosition.X() * 1.0E-3;
    135       y = candidatePosition.Y() * 1.0E-3;
    136       z = candidatePosition.Z() * 1.0E-3;
     135      x = candidatePosition.X()*1.0E-3;
     136      y = candidatePosition.Y()*1.0E-3;
     137      z = candidatePosition.Z()*1.0E-3;
    137138
    138139      // check that particle position is inside the cylinder
     
    151152
    152153      // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
    153       tmp = px * y - py * x;
    154       discr2 = pt2 * fRadius2 - tmp * tmp;
     154      tmp = px*y - py*x;
     155      discr2 = pt2*fRadius2 - tmp*tmp;
    155156
    156157      if(discr2 < 0.0)
     
    160161      }
    161162
    162       tmp = px * x + py * y;
     163      tmp = px*x + py*y;
    163164      discr = TMath::Sqrt(discr2);
    164       t1 = (-tmp + discr) / pt2;
    165       t2 = (-tmp - discr) / pt2;
     165      t1 = (-tmp + discr)/pt2;
     166      t2 = (-tmp - discr)/pt2;
    166167      t = (t1 < 0.0) ? t2 : t1;
    167168
    168       z_t = z + pz * t;
     169      z_t = z + pz*t;
    169170      if(TMath::Abs(z_t) > fHalfLength)
    170171      {
     
    175176
    176177      // final position
    177       x_t = x + px * t;
    178       y_t = y + py * t;
    179       z_t = z + pz * t;
    180 
    181       r_t = TMath::Sqrt(x_t * x_t + y_t * y_t + z_t * z_t);
     178      x_t = x + px*t;
     179      y_t = y + py*t;
     180      z_t = z + pz*t;
     181
     182      r_t = TMath::Sqrt(x_t*x_t + y_t*y_t + z_t*z_t);
     183
    182184
    183185      // here starts conversion code
    184       nsteps = Int_t(r_t / fStep);
     186      nsteps = Int_t(r_t/fStep);
    185187
    186188      x_i = x;
     
    188190      z_i = z;
    189191
    190       dt = t / nsteps;
     192      dt = t/nsteps;
    191193
    192194      converted = false;
     
    194196      for(i = 0; i < nsteps; ++i)
    195197      {
    196         x_i += px * dt;
    197         y_i += py * dt;
    198         z_i += pz * dt;
    199         pos_i.SetXYZ(x_i, y_i, z_i);
     198        x_i += px*dt;
     199        y_i += py*dt;
     200        z_i += pz*dt;
     201        pos_i.SetXYZ(x_i,y_i,z_i);
    200202
    201203        // convert photon position into cylindrical coordinates, cylindrical r,phi,z !!
    202204
    203         r_i = TMath::Sqrt(x_i * x_i + y_i * y_i);
     205        r_i = TMath::Sqrt(x_i*x_i + y_i*y_i);
    204206        phi_i = pos_i.Phi();
    205207
     
    208210
    209211        // convert into conversion probability
    210         p_conv = 1 - TMath::Exp(-7.0 / 9.0 * fStep * rate);
     212        p_conv = 1 - TMath::Exp(-7.0/9.0*fStep*rate);
    211213
    212214        // case conversion occurs
     
    219221          x2 = 1 - x1;
    220222
    221           ep = static_cast<Candidate *>(candidate->Clone());
    222           em = static_cast<Candidate *>(candidate->Clone());
    223 
    224           ep->Position.SetXYZT(x_i * 1.0E3, y_i * 1.0E3, z_i * 1.0E3, candidatePosition.T() + nsteps * dt * e * 1.0E3);
    225           em->Position.SetXYZT(x_i * 1.0E3, y_i * 1.0E3, z_i * 1.0E3, candidatePosition.T() + nsteps * dt * e * 1.0E3);
    226 
    227           ep->Momentum.SetPtEtaPhiE(x1 * pt, eta, phi, x1 * e);
    228           em->Momentum.SetPtEtaPhiE(x2 * pt, eta, phi, x2 * e);
     223          ep = static_cast<Candidate*>(candidate->Clone());
     224          em = static_cast<Candidate*>(candidate->Clone());
     225
     226          ep->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3);
     227          em->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3);
     228
     229          ep->Momentum.SetPtEtaPhiE(x1*pt, eta, phi, x1*e);
     230          em->Momentum.SetPtEtaPhiE(x2*pt, eta, phi, x2*e);
    229231
    230232          ep->PID = -11;
     
    249251
    250252//------------------------------------------------------------------------------
     253
Note: See TracChangeset for help on using the changeset viewer.