|
| 1 | +R__LOAD_LIBRARY(libStarlib.so) |
| 2 | +R__ADD_INCLUDE_PATH($STARlight_ROOT/include) |
| 3 | + |
| 4 | +#include "randomgenerator.h" |
| 5 | +#include "upcXevent.h" |
| 6 | +#include "starlight.h" |
| 7 | +#include "inputParameters.h" |
| 8 | + |
| 9 | +// usage: o2-sim -g external --configKeyValues 'GeneratorExternal.fileName=GeneratorStarlight.C;GeneratorExternal.funcName=GeneratorStarlight("kCohJpsiToMu")' |
| 10 | + |
| 11 | +namespace o2 |
| 12 | +{ |
| 13 | +namespace eventgen |
| 14 | +{ |
| 15 | +class GeneratorStarlight_class : public Generator |
| 16 | +{ |
| 17 | + public: |
| 18 | + GeneratorStarlight_class(){}; |
| 19 | + ~GeneratorStarlight_class() = default; |
| 20 | + void selectConfiguration(std::string val) { mSelectedConfiguration = val; }; |
| 21 | + void setCollisionSystem(float energyCM, int beam1Z, int beam1A, int beam2Z, int beam2A) {eCM = energyCM; projZ=beam1Z; projA=beam1A; targZ=beam2Z; targA=beam2A;}; |
| 22 | + bool setParameter(std::string line) { |
| 23 | + if (not mInputParameters.setParameter(line)){ |
| 24 | + std::cout << " --- [ERROR] cannot set parameter: " << line << std::endl; |
| 25 | + return false; |
| 26 | + } |
| 27 | + return true; |
| 28 | + } |
| 29 | + int getPdgMother(){return mPdgMother;} |
| 30 | + |
| 31 | + bool Init() override |
| 32 | + { |
| 33 | + Generator::Init(); |
| 34 | + |
| 35 | + float beam1energy = TMath::Sqrt(Double_t(projZ)/projA*targA/targZ)*eCM/2; |
| 36 | + float beam2energy = TMath::Sqrt(Double_t(projA)/projZ*targZ/targA)*eCM/2; |
| 37 | + float gamma1 = beam1energy/0.938272; |
| 38 | + float gamma2 = beam2energy/0.938272; |
| 39 | + float rapMax = 4.1 + 0.5*(TMath::ACosH(gamma2)-TMath::ACosH(gamma1)); |
| 40 | + |
| 41 | + const struct SLConfig { |
| 42 | + const char* name; |
| 43 | + int prod_mode; |
| 44 | + int prod_pid; |
| 45 | + int nwbins; |
| 46 | + float wmin; |
| 47 | + float wmax; |
| 48 | + float dy; |
| 49 | + int pdg_mother; |
| 50 | + bool decay_EvtGen; |
| 51 | + } slConfig[] = { |
| 52 | + {"kTwoGammaToMuLow", 1, 13, 292, 0.4, 15.0, 0.01, -1, 0 }, // from 0.4 to 15 GeV |
| 53 | + {"kTwoGammaToElLow", 1, 11, 292, 0.4, 15.0, 0.01, -1, 0 }, // from 0.4 to 15 GeV |
| 54 | + {"kTwoGammaToMuMedium", 1, 13, 264, 1.8, 15.0, 0.01, -1, 0 }, // from 1.8 to 15 GeV |
| 55 | + {"kTwoGammaToElMedium", 1, 11, 264, 1.8, 15.0, 0.01, -1, 0 }, // from 1.8 to 15 GeV |
| 56 | + {"kTwoGammaToMuHigh", 1, 13, 220, 4.0, 15.0, 0.01, -1, 0 }, // from 4.0 to 15 GeV |
| 57 | + {"kTwoGammaToElHigh", 1, 11, 220, 4.0, 15.0, 0.01, -1, 0 }, // from 4.0 to 15 GeV |
| 58 | + {"kTwoGammaToRhoRho", 1, 33, 20, -1.0, -1.0, 0.01, -1, 0 }, // |
| 59 | + {"kTwoGammaToF2", 1, 225, 20, -1.0, -1.0, 0.01, -1, 0 }, // |
| 60 | + {"kCohRhoToPi", 3, 113, 1200, -1.0, -1.0, 0.02, 113, 0 }, // |
| 61 | + {"kCohRhoToElEl", 3, 113011, 1200, -1.0, -1.0, 0.02, 113, 0 }, // |
| 62 | + {"kCohRhoToMuMu", 3, 113013, 1200, -1.0, -1.0, 0.02, 113, 0 }, // |
| 63 | + {"kCohRhoToPiWithCont", 3, 913, 1200, -1.0, -1.0, 0.02, -1, 0 }, // |
| 64 | + {"kCohRhoToPiFlat", 3, 113, 1, -1.0, 2.5, 0.02, 113, 0 }, // |
| 65 | + {"kCohPhiToKa", 2, 333, 20, -1.0, -1.0, 0.01, 333, 0 }, // |
| 66 | + {"kDirectPhiToKaKa", 3, 933, 20, -1.0, -1.0, 0.01, 333, 0 }, // |
| 67 | + {"kCohOmegaTo2Pi", 2, 223, 20, -1.0, -1.0, 0.01, 223, 0 }, // |
| 68 | + {"kCohOmegaTo3Pi", 2, 223, 20, -1.0, -1.0, 0.01, 223, 0 }, // |
| 69 | + {"kCohOmegaToPiPiPi", 2, 223211111, 20, -1.0, -1.0, 0.01, 333, 0 }, // |
| 70 | + {"kCohJpsiToMu", 2, 443013, 20, -1.0, -1.0, 0.01, 443, 0 }, // |
| 71 | + {"kCohJpsiToEl", 2, 443011, 20, -1.0, -1.0, 0.01, 443, 0 }, // |
| 72 | + {"kCohJpsiToElRad", 2, 443011, 20, -1.0, -1.0, 0.01, 443, 1 }, // |
| 73 | + {"kCohJpsiToProton", 2, 4432212, 20, -1.0, -1.0, 0.01, 443, 0 }, // |
| 74 | + {"kCohPsi2sToMu", 2, 444013, 20, -1.0, -1.0, 0.01, 100443, 0 }, // |
| 75 | + {"kCohPsi2sToEl", 2, 444011, 20, -1.0, -1.0, 0.01, 100443, 0 }, // |
| 76 | + {"kCohPsi2sToMuPi", 2, 444013, 20, -1.0, -1.0, 0.01, 100443, 1 }, // |
| 77 | + {"kCohPsi2sToElPi", 2, 444011, 20, -1.0, -1.0, 0.01, 100443, 1 }, // |
| 78 | + {"kCohUpsilonToMu", 2, 553013, 20, -1.0, -1.0, 0.01, 553, 0 }, // |
| 79 | + {"kCohUpsilonToEl", 2, 553011, 20, -1.0, -1.0, 0.01, 553, 0 }, // |
| 80 | + {"kIncohRhoToPi", 4, 113, 1200, -1.0, -1.0, 0.02, 113, 0 }, // |
| 81 | + {"kIncohRhoToElEl", 4, 113011, 1200, -1.0, -1.0, 0.02, 113, 0 }, // |
| 82 | + {"kIncohRhoToMuMu", 4, 113013, 1200, -1.0, -1.0, 0.02, 113, 0 }, // |
| 83 | + {"kIncohRhoToPiWithCont",4, 913, 1200, -1.0, -1.0, 0.02, -1, 0 }, // |
| 84 | + {"kIncohRhoToPiFlat", 4, 113, 1, -1.0, 2.5, 0.02, 113, 0 }, // |
| 85 | + {"kIncohPhiToKa", 4, 333, 20, -1.0, -1.0, 0.01, 333, 0 }, // |
| 86 | + {"kIncohOmegaTo2Pi", 4, 223, 20, -1.0, -1.0, 0.01, 223, 0 }, // |
| 87 | + {"kIncohOmegaTo3Pi", 4, 223, 20, -1.0, -1.0, 0.01, 223, 0 }, // |
| 88 | + {"kIncohOmegaToPiPiPi", 4, 223211111, 20, -1.0, -1.0, 0.01, 223, 0 }, // |
| 89 | + {"kIncohJpsiToMu", 4, 443013, 20, -1.0, -1.0, 0.01, 443, 0 }, // |
| 90 | + {"kIncohJpsiToEl", 4, 443011, 20, -1.0, -1.0, 0.01, 443, 0 }, // |
| 91 | + {"kIncohJpsiToElRad", 4, 443011, 20, -1.0, -1.0, 0.01, 443, 1 }, // |
| 92 | + {"kIncohJpsiToProton", 4, 4432212, 20, -1.0, -1.0, 0.01, 443, 0 }, // |
| 93 | + {"kIncohJpsiToLLbar", 4, 4433122, 20, -1.0, -1.0, 0.01, 443, 0 }, // |
| 94 | + {"kIncohPsi2sToMu", 4, 444013, 20, -1.0, -1.0, 0.01, 100443, 0 }, // |
| 95 | + {"kIncohPsi2sToEl", 4, 444011, 20, -1.0, -1.0, 0.01, 100443, 0 }, // |
| 96 | + {"kIncohPsi2sToMuPi", 4, 444013, 20, -1.0, -1.0, 0.01, 100443, 1 }, // |
| 97 | + {"kIncohPsi2sToElPi", 4, 444011, 20, -1.0, -1.0, 0.01, 100443, 1 }, // |
| 98 | + {"kIncohUpsilonToMu", 4, 553013, 20, -1.0, -1.0, 0.01, 553, 0 }, // |
| 99 | + {"kIncohUpsilonToEl", 4, 553011, 20, -1.0, -1.0, 0.01, 553, 0 }, // |
| 100 | + }; |
| 101 | + |
| 102 | + const int nProcess = sizeof(slConfig)/sizeof(SLConfig); |
| 103 | + int idx = -1; |
| 104 | + for (int i=0; i<nProcess; ++i) { |
| 105 | + if (mSelectedConfiguration.compare(slConfig[i].name) == 0) { |
| 106 | + idx = i; |
| 107 | + break; |
| 108 | + } |
| 109 | + } |
| 110 | + |
| 111 | + if (idx == -1) { |
| 112 | + std::cout << "STARLIGHT process "<< mSelectedConfiguration <<" is not supported" << std::endl; |
| 113 | + return false; |
| 114 | + } |
| 115 | + |
| 116 | + mPdgMother = slConfig[idx].pdg_mother; |
| 117 | + mDecayEvtGen = slConfig[idx].decay_EvtGen; |
| 118 | + |
| 119 | + uint random_seed; |
| 120 | + unsigned long long int random_value = 0; |
| 121 | + ifstream urandom("/dev/urandom", ios::in|ios::binary); |
| 122 | + urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed)); |
| 123 | + |
| 124 | + setParameter(Form("BEAM_1_Z = %3i #Z of target",targZ)); |
| 125 | + setParameter(Form("BEAM_1_A = %3i #A of target",targA)); |
| 126 | + setParameter(Form("BEAM_2_Z = %3i #Z of projectile",projZ)); |
| 127 | + setParameter(Form("BEAM_2_A = %3i #A of projectile",projA)); |
| 128 | + setParameter(Form("BEAM_1_GAMMA = %6.1f #Gamma of the target",gamma1)); |
| 129 | + setParameter(Form("BEAM_2_GAMMA = %6.1f #Gamma of the projectile",gamma2)); |
| 130 | + setParameter(Form("W_MAX = %.1f #Max value of w",slConfig[idx].wmax)); |
| 131 | + setParameter(Form("W_MIN = %.1f #Min value of w",slConfig[idx].wmin)); |
| 132 | + setParameter(Form("W_N_BINS = %3i #Bins i w",slConfig[idx].nwbins)); |
| 133 | + setParameter(Form("RAP_MAX = %.2f #max y",rapMax)); |
| 134 | + setParameter(Form("RAP_N_BINS = %.0f #Bins i y",rapMax*2./slConfig[idx].dy)); |
| 135 | + setParameter("CUT_PT = 0 #Cut in pT? 0 = (no, 1 = yes)"); |
| 136 | + setParameter("PT_MIN = 0 #Minimum pT in GeV"); |
| 137 | + setParameter("PT_MAX = 10 #Maximum pT in GeV"); |
| 138 | + setParameter("CUT_ETA = 0 #Cut in pseudorapidity? (0 = no, 1 = yes)"); |
| 139 | + setParameter("ETA_MIN = -5 #Minimum pseudorapidity"); |
| 140 | + setParameter("ETA_MAX = 5 #Maximum pseudorapidity"); |
| 141 | + setParameter(Form("PROD_MODE = %i #gg or gP switch (1 = 2-photon, 2 = coherent vector meson (narrow), 3 = coherent vector meson (wide), # 4 = incoherent vector meson, 5 = A+A DPMJet single, 6 = A+A DPMJet double, 7 = p+A DPMJet single, 8 = p+A Pythia single )",slConfig[idx].prod_mode)); |
| 142 | + setParameter(Form("PROD_PID = %6i #Channel of interest (not relevant for photonuclear processes)",slConfig[idx].prod_pid)); |
| 143 | + setParameter(Form("RND_SEED = %i #Random number seed", random_seed)); |
| 144 | + setParameter("BREAKUP_MODE = 5 #Controls the nuclear breakup"); |
| 145 | + setParameter("INTERFERENCE = 0 #Interference (0 = off, 1 = on)"); |
| 146 | + setParameter("IF_STRENGTH = 1. #% of interfernce (0.0 - 0.1)"); |
| 147 | + setParameter("INT_PT_MAX = 0.24 #Maximum pt considered, when interference is turned on"); |
| 148 | + setParameter("INT_PT_N_BINS = 120 #Number of pt bins when interference is turned on"); |
| 149 | + setParameter("XSEC_METHOD = 1 # Set to 0 to use old method for calculating gamma-gamma luminosity"); //CM |
| 150 | + setParameter("BSLOPE_DEFINITION = 0"); // using default slope |
| 151 | + setParameter("BSLOPE_VALUE = 4.0"); // default slope value |
| 152 | + setParameter("PRINT_VM = 0"); // print cross sections and fluxes vs rapidity in stdout for VM photoproduction processes |
| 153 | + |
| 154 | + if (not mInputParameters.init()) { |
| 155 | + std::cout << "InitStarLight parameter initialization has failed" << std::endl; |
| 156 | + return false; |
| 157 | + } |
| 158 | + |
| 159 | + mStarLight = new starlight; |
| 160 | + mStarLight->setInputParameters(&mInputParameters); |
| 161 | + mRandomGenerator.SetSeed(mInputParameters.randomSeed()); |
| 162 | + mStarLight->setRandomGenerator(&mRandomGenerator); |
| 163 | + return mStarLight->init(); |
| 164 | + |
| 165 | + }; |
| 166 | + |
| 167 | + bool generateEvent() override { |
| 168 | + |
| 169 | + if (!mStarLight) { |
| 170 | + std::cout <<"GenerateEvent: StarLight class/object not properly constructed"<<std::endl; |
| 171 | + return false; |
| 172 | + } |
| 173 | + |
| 174 | + mEvent = mStarLight->produceEvent(); |
| 175 | + // boost event to the experiment CM frame |
| 176 | + mEvent.boost(0.5*(TMath::ACosH(mInputParameters.beam1LorentzGamma()) - TMath::ACosH(mInputParameters.beam2LorentzGamma()))); |
| 177 | + |
| 178 | + return true; |
| 179 | + |
| 180 | + }; |
| 181 | + |
| 182 | + // at importParticles we add particles to the output particle vector |
| 183 | + // according to the selected configuration |
| 184 | + bool importParticles() override |
| 185 | + { |
| 186 | + int nVtx(0); |
| 187 | + float vtx(0), vty(0), vtz(0), vtt(0); |
| 188 | + const std::vector<vector3>* slVtx(mEvent.getVertices()); |
| 189 | + if (slVtx == 0) { // not vertex assume 0,0,0,0; |
| 190 | + vtx = vty = vtz = vtt = 0.0; |
| 191 | + } else { // a vertex exits |
| 192 | + slVtx = mEvent.getVertices(); |
| 193 | + nVtx = slVtx->size(); |
| 194 | + } // end if |
| 195 | + |
| 196 | + const std::vector<starlightParticle>* slPartArr(mEvent.getParticles()); |
| 197 | + const int npart(mEvent.getParticles()->size()); |
| 198 | + |
| 199 | + if(mPdgMother != -1){ //Reconstruct mother particle for VM processes |
| 200 | + TLorentzVector lmoth; |
| 201 | + TLorentzVector ldaug; |
| 202 | + for(int ipart=0;ipart<npart;ipart++) { |
| 203 | + const starlightParticle* slPart(&(slPartArr->at(ipart))); |
| 204 | + ldaug.SetPxPyPzE(slPart->GetPx(), slPart->GetPy(), slPart->GetPz(), slPart->GetE()); |
| 205 | + lmoth += ldaug; |
| 206 | + } |
| 207 | + TParticle particle(mPdgMother, |
| 208 | + 11, |
| 209 | + -1, |
| 210 | + -1, |
| 211 | + 1, |
| 212 | + npart, |
| 213 | + lmoth.Px(), |
| 214 | + lmoth.Py(), |
| 215 | + lmoth.Pz(), |
| 216 | + lmoth.E(), |
| 217 | + 0,0,0,0); |
| 218 | + //particle.Print(); |
| 219 | + mParticles.push_back(particle); |
| 220 | + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(mParticles.back(), 11); |
| 221 | + } |
| 222 | + if(!mDecayEvtGen){ // Don't import daughters in case of external decayer |
| 223 | + for(int ipart=0;ipart<npart;ipart++) { |
| 224 | + const starlightParticle* slPart(&(slPartArr->at(ipart))); |
| 225 | + if (nVtx < 1) { // No verticies |
| 226 | + vtx = vty = vtz = vtt = 0.0; |
| 227 | + } else { |
| 228 | + vtx = (slVtx->at((ipart < nVtx ? ipart : 0))).X(); |
| 229 | + vty = (slVtx->at((ipart < nVtx ? ipart : 0))).Y(); |
| 230 | + vtz = (slVtx->at((ipart < nVtx ? ipart : 0))).Z(); |
| 231 | + vtt = 0.0; // no time given. |
| 232 | + } // end if |
| 233 | + TParticle particle(slPart->getPdgCode(), |
| 234 | + 1, |
| 235 | + 0, |
| 236 | + -1, |
| 237 | + slPart->getFirstDaughter(), |
| 238 | + slPart->getLastDaughter(), |
| 239 | + slPart->GetPx(), |
| 240 | + slPart->GetPy(), |
| 241 | + slPart->GetPz(), |
| 242 | + slPart->GetE(), |
| 243 | + vtx,vty,vtz,vtt); |
| 244 | + //particle.Print(); |
| 245 | + mParticles.push_back(particle); |
| 246 | + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(mParticles.back(), 1); |
| 247 | + } |
| 248 | + } |
| 249 | + return true; |
| 250 | + } |
| 251 | + |
| 252 | + private: |
| 253 | + starlight *mStarLight = 0x0; |
| 254 | + inputParameters mInputParameters; // simulation input information. |
| 255 | + randomGenerator mRandomGenerator; // STARLIGHT's own random generator |
| 256 | + upcXEvent mEvent; // object holding STARlight simulated event. |
| 257 | + std::string mSelectedConfiguration = ""; |
| 258 | + int mPdgMother = -1; |
| 259 | + bool mDecayEvtGen = 0; |
| 260 | + float eCM = 5020; //CMS energy |
| 261 | + int projA=208; //Beam |
| 262 | + int targA=208; |
| 263 | + int projZ=82; |
| 264 | + int targZ=82; |
| 265 | + |
| 266 | + }; |
| 267 | + |
| 268 | +} // namespace eventgen |
| 269 | +} // namespace o2 |
| 270 | + |
| 271 | + |
| 272 | +FairGenerator* |
| 273 | + GeneratorStarlight(std::string configuration = "empty",float energyCM = 5020, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208) |
| 274 | +{ |
| 275 | + auto gen = new o2::eventgen::GeneratorStarlight_class(); |
| 276 | + gen->selectConfiguration(configuration); |
| 277 | + gen->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A); |
| 278 | + return gen; |
| 279 | +} |
| 280 | + |
| 281 | + |
| 282 | + |
| 283 | + |
0 commit comments