Fork me on GitHub

Changeset c1ce3fe in git for modules/PhotonConversions.cc


Ignore:
Timestamp:
Jun 29, 2015, 10:28:05 PM (9 years ago)
Author:
Pavel Demin <pavel.demin@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
ff37d75
Parents:
839deb7
Message:

adapt DelphesCylindricalFormula and PhotonConversions to ROOT 6.04

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/PhotonConversions.cc

    r839deb7 rc1ce3fe  
    5757
    5858PhotonConversions::PhotonConversions() :
    59   fDecayXsec(0), fConversionMap(0), fItInputArray(0)
     59  fItInputArray(0), fConversionMap(0), fDecayXsec(0)
    6060{
    6161  fDecayXsec = new TF1;
     
    7474{
    7575  fRadius = GetDouble("Radius", 1.0);
    76   fRadius2 =  fRadius*fRadius;
    77  
     76  fRadius2 = fRadius*fRadius;
     77
    7878  fHalfLength = GetDouble("HalfLength", 3.0);
    79  
     79
    8080  fEtaMax = GetDouble("EtaMax", 5.0);
    8181  fEtaMin = GetDouble("EtaMin", 2.0);
    82  
    83   fStep = GetDouble("Step", 0.1); // in meters
    84  
    85   fConversionMap->Compile(GetString("ConversionMap", "0."));
    86  
    87   fDecayXsec->Compile("1 - 4/3*x*(1-x)");
    88   fDecayXsec->SetRange(0.0,1.0);
    89  
     82
     83  fStep = GetDouble("Step", 0.1); // in meters
     84
     85  fConversionMap->Compile(GetString("ConversionMap", "0.0"));
     86
     87#if  ROOT_VERSION_CODE < ROOT_VERSION(6,04,00)
     88  fDecayXsec->Compile("1.0 - 4.0/3.0 * x * (1.0 - x)");
     89#else
     90  fDecayXsec->GetFormula()->Compile("1.0 - 4.0/3.0 * x * (1.0 - x)");
     91#endif
     92  fDecayXsec->SetRange(0.0, 1.0);
     93
    9094  // import array with output from filter/classifier module
    9195
     
    111115void PhotonConversions::Process()
    112116{
    113   Candidate *candidate, *mother, *ep, *em;
     117  Candidate *candidate, *ep, *em;
    114118  TLorentzVector candidatePosition, candidateMomentum;
    115119  TVector3 pos_i;
    116   Double_t px, py, pz, pt, pt2, e, eta, phi, q;
     120  Double_t px, py, pz, pt, pt2, e, eta, phi;
    117121  Double_t x, y, z, t;
    118122  Double_t x_t, y_t, z_t, r_t;
     
    123127  Double_t rate, p_conv, x1, x2;
    124128  Bool_t converted;
    125  
    126   const Double_t c_light = 2.99792458E8;
    127  
     129
    128130  fItInputArray->Reset();
    129131  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    130132  {
    131    
    132     if(candidate->PID != 22) fOutputArray->Add(candidate);
     133
     134    if(candidate->PID != 22)
     135    {
     136      fOutputArray->Add(candidate);
     137    }
    133138    else
    134139    {
     
    138143      y = candidatePosition.Y()*1.0E-3;
    139144      z = candidatePosition.Z()*1.0E-3;
    140       q = candidate->Charge;
    141145
    142146      // check that particle position is inside the cylinder
    143       if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
    144       {
    145         continue;
    146       }
     147      if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength) continue;
    147148
    148149      px = candidateMomentum.Px();
     
    154155      phi = candidateMomentum.Phi();
    155156      e = candidateMomentum.E();
    156    
     157
    157158      if(eta < fEtaMin || eta > fEtaMax) continue;
    158    
    159        // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
     159
     160      // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
    160161      tmp = px*y - py*x;
    161162      discr2 = pt2*fRadius2 - tmp*tmp;
     
    182183
    183184      // final position
    184 
    185185      x_t = x + px*t;
    186186      y_t = y + py*t;
    187187      z_t = z + pz*t;
    188      
     188
    189189      r_t = TMath::Sqrt(x_t*x_t + y_t*y_t + z_t*z_t);
    190        
    191        
    192       // here starts conversion code   
     190
     191
     192      // here starts conversion code
    193193      nsteps = Int_t(r_t/fStep);
    194        
     194
    195195      x_i = x;
    196196      y_i = y;
    197197      z_i = z;
    198      
     198
    199199      dt = t/nsteps;
    200      
     200
    201201      converted = false;
    202      
    203       for(i=0; i<nsteps; ++i)
     202
     203      for(i = 0; i < nsteps; ++i)
    204204      {
    205        
    206         x_i += px*dt;
    207         y_i += py*dt;
    208         z_i += pz*dt;
     205        x_i += px*dt;
     206        y_i += py*dt;
     207        z_i += pz*dt;
    209208        pos_i.SetXYZ(x_i,y_i,z_i);
    210        
    211         //convert photon position into cylindrical coordinates, cylindrical r,phi,z !!
    212        
    213         r_i = TMath::Sqrt(x_i*x_i + y_i*y_i);
     209
     210        // convert photon position into cylindrical coordinates, cylindrical r,phi,z !!
     211
     212        r_i = TMath::Sqrt(x_i*x_i + y_i*y_i);
    214213        phi_i = pos_i.Phi();
    215        
     214
    216215        // read conversion rate/meter from card
    217         rate = fConversionMap->Eval(r_i, phi_i, z_i);
    218        
    219         // convert into conversion probability
    220         p_conv = 1 - TMath::Exp(-7.0/9.0*fStep*rate);
    221        
    222            
    223         //case conversion occurs
    224         if(gRandom->Uniform() < p_conv)
    225         {
    226           converted = true;
    227          
    228           // generate x1 and x2, the fraction of the photon energy taken resp. by e+ and e-
    229           x1 = fDecayXsec->GetRandom();
    230           x2 = 1 - x1;
    231        
    232           mother = candidate;
    233           ep = static_cast<Candidate*>(candidate->Clone());
    234           em = static_cast<Candidate*>(candidate->Clone());
    235        
    236           ep->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3);
    237           em->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3);
    238          
    239           ep->Momentum.SetPtEtaPhiE(x1*pt, eta, phi, x1*e);
    240           em->Momentum.SetPtEtaPhiE(x2*pt, eta, phi, x2*e);
    241        
    242           ep->PID = -11;
    243           em->PID = 11;
    244          
    245           ep->Charge = 1.0;
    246           em->Charge = -1.0;
    247        
    248           ep->IsFromConversion = 1;
    249           em->IsFromConversion = 1;
    250          
    251           fOutputArray->Add(em);
    252           fOutputArray->Add(ep);
    253        
    254           break;
    255        
    256         }
    257      
     216        rate = fConversionMap->Eval(r_i, phi_i, z_i);
     217
     218        // convert into conversion probability
     219        p_conv = 1 - TMath::Exp(-7.0/9.0*fStep*rate);
     220
     221        // case conversion occurs
     222        if(gRandom->Uniform() < p_conv)
     223        {
     224          converted = true;
     225
     226          // generate x1 and x2, the fraction of the photon energy taken resp. by e+ and e-
     227          x1 = fDecayXsec->GetRandom();
     228          x2 = 1 - x1;
     229
     230          ep = static_cast<Candidate*>(candidate->Clone());
     231          em = static_cast<Candidate*>(candidate->Clone());
     232
     233          ep->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3);
     234          em->Position.SetXYZT(x_i*1.0E3, y_i*1.0E3, z_i*1.0E3, candidatePosition.T() + nsteps*dt*e*1.0E3);
     235
     236          ep->Momentum.SetPtEtaPhiE(x1*pt, eta, phi, x1*e);
     237          em->Momentum.SetPtEtaPhiE(x2*pt, eta, phi, x2*e);
     238
     239          ep->PID = -11;
     240          em->PID = 11;
     241
     242          ep->Charge = 1.0;
     243          em->Charge = -1.0;
     244
     245          ep->IsFromConversion = 1;
     246          em->IsFromConversion = 1;
     247
     248          fOutputArray->Add(em);
     249          fOutputArray->Add(ep);
     250
     251          break;
     252        }
    258253      }
    259      
    260      if(!converted) fOutputArray->Add(candidate);
    261  
     254      if(!converted) fOutputArray->Add(candidate);
    262255    }
    263256  }
Note: See TracChangeset for help on using the changeset viewer.