|
| 1 | +// Copyright CERN and copyright holders of ALICE O2. This software is |
| 2 | +// distributed under the terms of the GNU General Public License v3 (GPL |
| 3 | +// Version 3), copied verbatim in the file "COPYING". |
| 4 | +// |
| 5 | +// See http://alice-o2.web.cern.ch/license for full licensing information. |
| 6 | +// |
| 7 | +// In applying this license CERN does not waive the privileges and immunities |
| 8 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 9 | +// or submit itself to any jurisdiction. |
| 10 | + |
| 11 | +// Generator for muons according to kinematic parametrizations at ALICE |
| 12 | +// (not at the surface). |
| 13 | + |
| 14 | +// Modifications for O2: [email protected] |
| 15 | + |
| 16 | +#include <TParticle.h> |
| 17 | +#include <TRandom.h> |
| 18 | +#include <TVirtualMC.h> |
| 19 | +#include <TGeoGlobalMagField.h> |
| 20 | +#include "GeneratorCosmics.h" |
| 21 | + |
| 22 | +//----------------------------------------------------------------------------- |
| 23 | +GeneratorCosmics::GeneratorCosmics() : TGenerator("GeneratorCosmics", "GeneratorCosmics") |
| 24 | +{ |
| 25 | +} |
| 26 | + |
| 27 | +//----------------------------------------------------------------------------- |
| 28 | +void GeneratorCosmics::GenerateEvent() |
| 29 | +{ |
| 30 | + // |
| 31 | + // Generate muon(s) |
| 32 | + constexpr int MuMinusPDG = 13, MuPlusPDG = -13; |
| 33 | + constexpr float MuMass = 0.1056583; |
| 34 | + // |
| 35 | + if (!mFieldIsSet) { |
| 36 | + auto fld = TGeoGlobalMagField::Instance()->GetField(); |
| 37 | + if (TVirtualMC::GetMC() && TVirtualMC::GetMC()->GetMagField()) { |
| 38 | + fld = TVirtualMC::GetMC()->GetMagField(); |
| 39 | + } |
| 40 | + if (fld) { |
| 41 | + double r[3] = {0, 0, 0}, b[3] = {0, 0, 0}; |
| 42 | + fld->Field(r, b); |
| 43 | + setBkG(b[2]); |
| 44 | + } |
| 45 | + else { |
| 46 | + throw std::runtime_error("Failed to fetch magnetic field"); |
| 47 | + } |
| 48 | + } |
| 49 | + |
| 50 | + fParticles->Clear(); |
| 51 | + int npart = 0; |
| 52 | + std::unique_ptr<TParticle> part; |
| 53 | + // |
| 54 | + while (npart < mNPart) { // until needed numbe of muons generated |
| 55 | + int trials = 0; |
| 56 | + |
| 57 | + do { // until particle passes all selections |
| 58 | + if (++trials > mMaxTrials) { |
| 59 | + throw std::runtime_error("max. trials reached"); |
| 60 | + } |
| 61 | + int pdg = gRandom->Rndm() < MuMinusFraction ? MuMinusPDG : MuPlusPDG; // mu- : mu+ |
| 62 | + float r[3] = {0.f, mROrigin, 0.f}, p[3], ptot = 0, pt = 0; |
| 63 | + |
| 64 | + if (mParam == GenParamType::ParamMI) { |
| 65 | + ptot = mGenFun->GetRandom(); |
| 66 | + p[1] = -ptot; |
| 67 | + if (gRandom->Rndm() > 0.9) { |
| 68 | + p[0] = gRandom->Gaus(0.0, 0.4) * ptot; |
| 69 | + p[2] = gRandom->Gaus(0.0, 0.4) * ptot; |
| 70 | + } else { |
| 71 | + p[0] = gRandom->Gaus(0.0, 0.2) * ptot; |
| 72 | + p[2] = gRandom->Gaus(0.0, 0.2) * ptot; |
| 73 | + } |
| 74 | + ptot = std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]); |
| 75 | + pt = std::sqrt(p[0] * p[0] + p[1] * p[1]); |
| 76 | + } else { |
| 77 | + ptot = mGenFun->GetRandom(); |
| 78 | + float theta = 0, phi = 0; |
| 79 | + do { |
| 80 | + theta = gRandom->Gaus(0.5 * PiConst, 0.42); |
| 81 | + } while (std::abs(theta - 0.5 * PiConst) > mMaxAngleWRTVertical); |
| 82 | + do { |
| 83 | + phi = gRandom->Gaus(-0.5 * PiConst, 0.42); |
| 84 | + } while (std::abs(phi + 0.5 * PiConst) > mMaxAngleWRTVertical); |
| 85 | + |
| 86 | + pt = ptot * std::sin(theta); |
| 87 | + p[0] = pt * std::cos(phi); |
| 88 | + p[1] = pt * std::sin(phi); |
| 89 | + p[2] = ptot * std::cos(theta); |
| 90 | + } |
| 91 | + |
| 92 | + if (ptot < mPMin || ptot > mPMax || std::acos(std::abs(p[1]) / ptot) > mMaxAngleWRTVertical) { |
| 93 | + continue; |
| 94 | + } |
| 95 | + |
| 96 | + // estimate max deflection |
| 97 | + float xpos = 999, zpos = 999; |
| 98 | + if (!getXZatOrigin(xpos, zpos, r, p, -pdg)) { |
| 99 | + continue; |
| 100 | + } |
| 101 | + auto slpX = p[0] / p[1], slpZ = p[2] / p[1]; |
| 102 | + xpos += mROrigin * slpX; // get max bending |
| 103 | + zpos += mROrigin * slpZ; // only, w/o slopes contribution |
| 104 | + float xmin = -mXAcc, xmax = mXAcc, zmin = -mZAcc, zmax = mZAcc; |
| 105 | + if (xpos < 0) { |
| 106 | + xmax -= xpos; |
| 107 | + } else { |
| 108 | + xmin -= xpos; |
| 109 | + } |
| 110 | + if (zpos < 0) { |
| 111 | + zmax -= zpos; |
| 112 | + } else { |
| 113 | + zmin -= zpos; |
| 114 | + } |
| 115 | + |
| 116 | + r[0] = mROrigin * slpX + xmin + gRandom->Rndm() * (xmax - xmin); |
| 117 | + r[2] = mROrigin * slpZ + zmin + gRandom->Rndm() * (zmax - zmin); |
| 118 | + |
| 119 | + // propagate to fixed Y=mROrigin plane to fixed radius in field free region: solve quadratic equation of circle - line intersection |
| 120 | + auto a = slpX * slpX + 1, xred = r[0] - r[1] * slpX, b = xred * slpX, det = b * b - a * (xred * xred - mROrigin * mROrigin); |
| 121 | + if (det < 0.) { |
| 122 | + continue; |
| 123 | + } |
| 124 | + r[1] = (-b + std::sqrt(det)) / a; |
| 125 | + r[0] += (r[1] - mROrigin) * slpX; |
| 126 | + r[2] += (r[1] - mROrigin) * p[2] / p[1]; |
| 127 | + |
| 128 | + // check trigger condition |
| 129 | + if (!getXZatOrigin(xpos, zpos, r, p, -pdg) || std::abs(xpos) > mXAcc || std::abs(zpos) > mZAcc) { |
| 130 | + continue; |
| 131 | + } |
| 132 | + |
| 133 | + auto etot = std::sqrt(MuMass * MuMass + ptot * ptot); |
| 134 | + part = std::make_unique<TParticle>(pdg, 1, -1, -1, -1, -1, p[0], p[1], p[2], etot, r[0], r[1], r[2], 0); |
| 135 | + break; |
| 136 | + } while (1); |
| 137 | + |
| 138 | + fParticles->Add(part.get()); |
| 139 | + npart++; |
| 140 | + } |
| 141 | +} |
| 142 | + |
| 143 | +//----------------------------------------------------------------------------- |
| 144 | +void GeneratorCosmics::Init() |
| 145 | +{ |
| 146 | + // |
| 147 | + // Initialisation, check consistency of selected ranges |
| 148 | + // |
| 149 | + switch (mParam) { |
| 150 | + case GenParamType::ParamMI: |
| 151 | + mGenFun.reset(new TF1("genFun", "exp(-x/30.)", mPMin, mPMax * 2)); |
| 152 | + break; |
| 153 | + case GenParamType::ParamACORDE: |
| 154 | + mGenFun.reset(new TF1("genFun", "x/(1.+(x/12.8)*(x/12.8))^1.96", mPMin, mPMax)); |
| 155 | + break; |
| 156 | + case GenParamType::ParamTPC: |
| 157 | + mGenFun.reset(new TF1("genFun", "x/(1.+(x/3.)*(x/3.))^1.", mPMin, mPMax)); |
| 158 | + break; |
| 159 | + } |
| 160 | + |
| 161 | + printf("Cosmics generator configuration:\n"); |
| 162 | + printf("Parameterization type: %d with %e < p < %e\n", int(mParam), mPMin, mPMax); |
| 163 | + printf("Tracks created at R=%.2f and requested to have |X|<%.2f and |Z|<%.2f at Y=0\n", mROrigin, mXAcc, mZAcc); |
| 164 | + printf("Magnetic field %f\n", mBkG); |
| 165 | + return; |
| 166 | +} |
| 167 | + |
| 168 | +void GeneratorCosmics::setMaxAngleWRTVertical(float max) |
| 169 | +{ |
| 170 | + if (max < 1e-6) { |
| 171 | + max = 1e-6; |
| 172 | + } |
| 173 | + if (max > 90.) { |
| 174 | + throw std::runtime_error("angle must be in ]0 : 90] range"); |
| 175 | + } |
| 176 | + mMaxAngleWRTVertical = max * PiConst / 180.; |
| 177 | +} |
| 178 | + |
| 179 | +bool GeneratorCosmics::getXZatOrigin(float& xpos, float& zpos, const float r[3], const float p[3], int q) const |
| 180 | +{ |
| 181 | + // get position of the track at Y=0 in the lab frame |
| 182 | + constexpr float B2C = -0.299792458e-3, Almost0 = 1e-9, Almost1 = 1 - Almost0; |
| 183 | + // work explicitly in alpha = -pi/2 frame, i.e. track param X = -r[1] |
| 184 | + auto pt = std::sqrt(p[0] * p[0] + p[1] * p[1]), q2pt = q > 0 ? 1.f / pt : -1.f / pt; |
| 185 | + auto phi = std::atan2(p[0], -p[1]), snp = std::sin(phi), tgl = p[2] / pt; |
| 186 | + auto crv = q2pt * mBkG * B2C, x2r = crv * r[1], f1 = snp, f2 = f1 + x2r; |
| 187 | + if (std::abs(f1) >= Almost1 || std::abs(f2) >= Almost1 || std::abs(q2pt) < Almost0) { |
| 188 | + return false; |
| 189 | + } |
| 190 | + auto r1 = std::sqrt((1. - f1) * (1. + f1)), r2 = std::sqrt((1. - f2) * (1. + f2)); |
| 191 | + if (std::abs(r1) < Almost0 || std::abs(r2) < Almost0) { |
| 192 | + return false; |
| 193 | + } |
| 194 | + auto dy2dx = (f1 + f2) / (r1 + r2); |
| 195 | + xpos = r[0] + r[1] * dy2dx; // print dy2dx, it is wrong? |
| 196 | + if (std::abs(x2r) < 0.05) { |
| 197 | + zpos = r[2] + r[1] * (r2 + f2 * dy2dx) * tgl; |
| 198 | + } else { |
| 199 | + auto rot = std::asin(r1 * f2 - r2 * f1); // more economic version from Yura. |
| 200 | + if (f1 * f1 + f2 * f2 > 1 && f1 * f2 < 0) { // special cases of large rotations or large abs angles |
| 201 | + rot = f2 > 0 ? PiConst - rot : -PiConst - rot; |
| 202 | + } |
| 203 | + zpos = r[2] + tgl / crv * rot; |
| 204 | + } |
| 205 | + return true; |
| 206 | +} |
| 207 | + |
| 208 | +//______________________________________________________________________________ |
| 209 | +int GeneratorCosmics::ImportParticles(TClonesArray *particles, Option_t *option) |
| 210 | +{ |
| 211 | + if (particles == 0) return 0; |
| 212 | + TClonesArray &clonesParticles = *particles; |
| 213 | + clonesParticles.Clear(); |
| 214 | + Int_t numpart = fParticles->GetEntries(); |
| 215 | + for (int i = 0; i<numpart; i++) { |
| 216 | + TParticle *particle = (TParticle *)fParticles->At(i); |
| 217 | + new(clonesParticles[i]) TParticle(*particle); |
| 218 | + } |
| 219 | + return numpart; |
| 220 | +} |
| 221 | + |
| 222 | +//______________________________________________________________________________ |
| 223 | +void GeneratorCosmics::setBkG(float b) |
| 224 | +{ |
| 225 | + mBkG = b; |
| 226 | + mFieldIsSet = true; |
| 227 | + printf("Setting field to %f kG\n", mBkG); |
| 228 | + return; |
| 229 | +} |
0 commit comments