File indexing completed on 2024-04-06 12:13:28
0001
0002
0003
0004
0005
0006 #define sim_cxx
0007
0008 #include "GeneratorInterface/CosmicMuonGenerator/interface/CosmicMuonGenerator.h"
0009
0010 void CosmicMuonGenerator::runCMG() {
0011 initialize();
0012 for (unsigned int iGen = 0; iGen < NumberOfEvents; ++iGen) {
0013 nextEvent();
0014 }
0015 terminate();
0016 }
0017
0018 void CosmicMuonGenerator::setRandomEngine(CLHEP::HepRandomEngine* v) {
0019 if (delRanGen)
0020 delete RanGen;
0021 RanGen = v;
0022 delRanGen = false;
0023 Cosmics->setRandomEngine(v);
0024 }
0025
0026 void CosmicMuonGenerator::initialize(CLHEP::HepRandomEngine* rng) {
0027 if (delRanGen)
0028 delete RanGen;
0029 if (!rng) {
0030 RanGen = new CLHEP::HepJamesRandom;
0031 RanGen->setSeed(RanSeed, 0);
0032 delRanGen = true;
0033 } else {
0034 RanGen = rng;
0035 delRanGen = false;
0036 }
0037 checkIn();
0038 if (NumberOfEvents > 0) {
0039
0040 double RadiusTargetEff = RadiusOfTarget;
0041 double Z_DistTargetEff = ZDistOfTarget;
0042
0043 if (TrackerOnly == true) {
0044 RadiusTargetEff = RadiusTracker;
0045 Z_DistTargetEff = Z_DistTracker;
0046 }
0047 Target3dRadius = sqrt(RadiusTargetEff * RadiusTargetEff + Z_DistTargetEff * Z_DistTargetEff) + MinStepSize;
0048 if (Debug)
0049 std::cout << " radius of sphere around target = " << Target3dRadius << " mm" << std::endl;
0050
0051 if (MinTheta > 90. * Deg2Rad)
0052 SurfaceRadius = (RadiusCMS) * (-tan(MinTheta)) + MinStepSize;
0053 else
0054 SurfaceRadius = (SurfaceOfEarth + PlugWidth + RadiusTargetEff) * tan(MaxTheta) + Target3dRadius;
0055 if (Debug)
0056 std::cout << " starting point radius at Surface + PlugWidth = " << SurfaceRadius << " mm" << std::endl;
0057
0058 OneMuoEvt.PlugVx = PlugVx;
0059 OneMuoEvt.PlugVz = PlugVz;
0060 OneMuoEvt.RhoAir = RhoAir;
0061 OneMuoEvt.RhoWall = RhoWall;
0062 OneMuoEvt.RhoRock = RhoRock;
0063 OneMuoEvt.RhoClay = RhoClay;
0064 OneMuoEvt.RhoPlug = RhoPlug;
0065
0066 if (MinTheta >= 90. * Deg2Rad)
0067 Cosmics->initializeNuMu(MinP, MaxP, MinTheta, MaxTheta, MinEnu, MaxEnu, MinPhi, MaxPhi, NuProdAlt, RanGen);
0068 else
0069 Cosmics->initialize(MinP, MaxP, MinTheta, MaxTheta, RanGen, TIFOnly_constant, TIFOnly_linear);
0070
0071 #if ROOT_INTERACTIVE
0072
0073 TH1D* ene = new TH1D("ene", "generated energy", 210, 0., 1050.);
0074 TH1D* the = new TH1D("the", "generated theta", 90, 0., 90.);
0075 TH1D* phi = new TH1D("phi", "generated phi", 120, 0., 360.);
0076 TH3F* ver = new TH3F("ver", "Z-X-Y coordinates", 100, -25., 25., 20, -10., 10., 40, -10., 10.);
0077 #endif
0078 if (EventDisplay)
0079 initEvDis();
0080 std::cout << std::endl;
0081
0082 if (MultiMuon) {
0083 MultiIn = nullptr;
0084
0085 std::cout << "MultiMuonFileName.c_str()=" << MultiMuonFileName.c_str() << std::endl;
0086 MultiIn = new TFile(MultiMuonFileName.c_str());
0087
0088 if (!MultiIn)
0089 std::cout << "MultiMuon=True: MultiMuonFileName='" << MultiMuonFileName.c_str() << "' does not exist"
0090 << std::endl;
0091 else
0092 std::cout << "MultiMuonFile: " << MultiMuonFileName.c_str() << " opened!" << std::endl;
0093
0094 MultiTree = (TTree*)MultiIn->Get("sim");
0095 SimTree = new sim(MultiTree);
0096 SimTree->Init(MultiTree);
0097 SimTreeEntries = SimTree->fChain->GetEntriesFast();
0098 std::cout << "SimTreeEntries=" << SimTreeEntries << std::endl;
0099
0100 if (MultiMuonFileFirstEvent <= 0)
0101 SimTree_jentry = 0;
0102 else
0103 SimTree_jentry = MultiMuonFileFirstEvent - 1;
0104
0105 NcloseMultiMuonEvents = 0;
0106 NskippedMultiMuonEvents = 0;
0107 }
0108
0109 if (!MultiMuon || (MultiMuon && MultiIn))
0110 NotInitialized = false;
0111 }
0112 }
0113
0114 void CosmicMuonGenerator::nextEvent() {
0115 double E = 0.;
0116 double Theta = 0.;
0117 double Phi = 0.;
0118 double RxzV = 0.;
0119 double PhiV = 0.;
0120 if (int(Nsel) % 100 == 0)
0121 std::cout << " generated " << int(Nsel) << " events" << std::endl;
0122
0123 bool notSelected = true;
0124 while (notSelected) {
0125 bool badMomentumGenerated = true;
0126 while (badMomentumGenerated) {
0127 if (MinTheta > 90. * Deg2Rad)
0128 Cosmics->generateNuMu();
0129 else
0130 Cosmics->generate();
0131
0132 E = sqrt(Cosmics->momentum_times_charge() * Cosmics->momentum_times_charge() + MuonMass * MuonMass);
0133 Theta = TMath::ACos(Cosmics->cos_theta());
0134 Ngen += 1.;
0135 badMomentumGenerated = false;
0136 Phi = RanGen->flat() * (MaxPhi - MinPhi) + MinPhi;
0137 }
0138 Norm->events_n100cos(E, Theta);
0139 Ndiced += 1;
0140
0141
0142 double Nver = 0.;
0143 bool badVertexGenerated = true;
0144 while (badVertexGenerated) {
0145 RxzV = sqrt(RanGen->flat()) * SurfaceRadius;
0146 PhiV = RanGen->flat() * TwoPi;
0147
0148 double dPhi = Pi;
0149 if (RxzV > Target3dRadius)
0150 dPhi = asin(Target3dRadius / RxzV);
0151 double rotPhi = PhiV + Pi;
0152 if (rotPhi > TwoPi)
0153 rotPhi -= TwoPi;
0154 double disPhi = std::fabs(rotPhi - Phi);
0155 if (disPhi > Pi)
0156 disPhi = TwoPi - disPhi;
0157 if (disPhi < dPhi || AcptAllMu)
0158 badVertexGenerated = false;
0159 Nver += 1.;
0160 }
0161 Ngen += (Nver - 1.);
0162
0163
0164 int id = 13;
0165 if (Cosmics->momentum_times_charge() > 0.)
0166 id = -13;
0167 double absMom = sqrt(E * E - MuonMass * MuonMass);
0168 double verMom = absMom * cos(Theta);
0169 double horMom = absMom * sin(Theta);
0170 double Px = horMom * sin(Phi);
0171 double Py = -verMom;
0172 double Pz = horMom * cos(Phi);
0173 double Vx = RxzV * sin(PhiV);
0174
0175 double Vy;
0176 if (MinTheta > 90. * Deg2Rad)
0177 Vy = -RadiusCMS;
0178 else
0179 Vy = SurfaceOfEarth + PlugWidth;
0180
0181 double Vz = RxzV * cos(PhiV);
0182 double T0 = (RanGen->flat() * (MaxT0 - MinT0) + MinT0) * SpeedOfLight;
0183
0184 Id_at = id;
0185 Px_at = Px;
0186 Py_at = Py;
0187 Pz_at = Pz;
0188 E_at = E;
0189 Vx_at = Vx;
0190 Vy_at = Vy;
0191 Vz_at = Vz;
0192 T0_at = T0;
0193
0194 OneMuoEvt.create(id, Px, Py, Pz, E, MuonMass, Vx, Vy, Vz, T0);
0195
0196 if (goodOrientation()) {
0197 if (MinTheta > 90. * Deg2Rad)
0198 OneMuoEvt.propagate(0., RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
0199 else
0200 OneMuoEvt.propagate(ElossScaleFactor, RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
0201 }
0202
0203 if ((OneMuoEvt.hitTarget() && sqrt(OneMuoEvt.e() * OneMuoEvt.e() - MuonMass * MuonMass) > MinP_CMS) ||
0204 AcptAllMu == true) {
0205 Nsel += 1.;
0206 notSelected = false;
0207 }
0208 }
0209
0210 EventWeight = 1.;
0211
0212
0213 Id_sf.resize(1);
0214 Px_sf.resize(1);
0215 Py_sf.resize(1);
0216 Pz_sf.resize(1);
0217 E_sf.resize(1);
0218
0219 Vx_sf.resize(1);
0220 Vy_sf.resize(1);
0221 Vz_sf.resize(1);
0222 T0_sf.resize(1);
0223
0224 Id_sf[0] = Id_at;
0225 Px_sf[0] = Px_at;
0226 Py_sf[0] = Py_at;
0227 Pz_sf[0] = Pz_at;
0228 E_sf[0] = E_at;
0229 Vx_sf[0] = Vx_at;
0230 Vy_sf[0] = Vy_at;
0231 Vz_sf[0] = Vz_at;
0232 T0_sf[0] = T0_at;
0233
0234
0235 Id_ug.resize(1);
0236 Px_ug.resize(1);
0237 Py_ug.resize(1);
0238 Pz_ug.resize(1);
0239 E_ug.resize(1);
0240
0241 Vx_ug.resize(1);
0242 Vy_ug.resize(1);
0243 Vz_ug.resize(1);
0244 T0_ug.resize(1);
0245
0246 Id_ug[0] = OneMuoEvt.id();
0247 Px_ug[0] = OneMuoEvt.px();
0248 Py_ug[0] = OneMuoEvt.py();
0249 Pz_ug[0] = OneMuoEvt.pz();
0250 E_ug[0] = OneMuoEvt.e();
0251
0252 Vx_ug[0] = OneMuoEvt.vx();
0253 Vy_ug[0] = OneMuoEvt.vy();
0254 Vz_ug[0] = OneMuoEvt.vz();
0255 T0_ug[0] = OneMuoEvt.t0();
0256
0257
0258 #if ROOT_INTERACTIVE
0259 ene->Fill(OneMuoEvt.e());
0260 the->Fill((OneMuoEvt.theta() * Rad2Deg));
0261 phi->Fill((OneMuoEvt.phi() * Rad2Deg));
0262 ver->Fill((OneMuoEvt.vz() / 1000.), (OneMuoEvt.vx() / 1000.), (OneMuoEvt.vy() / 1000.));
0263 #endif
0264 if (Debug) {
0265 std::cout << "new event" << std::endl;
0266 std::cout << " Px,Py,Pz,E,m = " << OneMuoEvt.px() << ", " << OneMuoEvt.py() << ", " << OneMuoEvt.pz() << ", "
0267 << OneMuoEvt.e() << ", " << OneMuoEvt.m() << " GeV" << std::endl;
0268 std::cout << " Vx,Vy,Vz,t0 = " << OneMuoEvt.vx() << ", " << OneMuoEvt.vy() << ", " << OneMuoEvt.vz() << ", "
0269 << OneMuoEvt.t0() << " mm" << std::endl;
0270 }
0271 if (EventDisplay)
0272 displayEv();
0273 }
0274
0275 bool CosmicMuonGenerator::nextMultiEvent() {
0276 if (Debug)
0277 std::cout << "\nEntered CosmicMuonGenerator::nextMultiEvent()" << std::endl;
0278 bool EvtRejected = true;
0279 bool MuInMaxDist = false;
0280 double MinDist;
0281
0282 while (EvtRejected) {
0283
0284
0285 Long64_t ientry = SimTree->GetEntry(SimTree_jentry);
0286 std::cout << "CosmicMuonGenerator::nextMultiEvent(): SimTree_jentry="
0287 << SimTree_jentry
0288
0289 << " SimTreeEntries=" << SimTreeEntries << std::endl;
0290 if (ientry < 0)
0291 return false;
0292 if (SimTree_jentry < SimTreeEntries) {
0293 SimTree_jentry++;
0294 } else {
0295 std::cout << "CosmicMuonGenerator.cc::nextMultiEvent: No more events in file!" << std::endl;
0296 return false;
0297 }
0298
0299 int nmuons = SimTree->shower_nParticlesWritten;
0300 if (nmuons < MultiMuonNmin) {
0301 std::cout << "CosmicMuonGenerator.cc: Warning! Less than " << MultiMuonNmin << " muons in event!" << std::endl;
0302 std::cout << "trying next event from file" << std::endl;
0303 NskippedMultiMuonEvents++;
0304 continue;
0305 }
0306
0307 Px_mu.resize(nmuons);
0308 Py_mu.resize(nmuons);
0309 Pz_mu.resize(nmuons);
0310 P_mu.resize(nmuons);
0311
0312 MinDist = 99999.e9;
0313 double MuMuDist;
0314 MuInMaxDist = false;
0315
0316 int NmuPmin = 0;
0317 for (int imu = 0; imu < nmuons; ++imu) {
0318 Px_mu[imu] =
0319 -SimTree->particle__Px[imu] * sin(NorthCMSzDeltaPhi) + SimTree->particle__Py[imu] * cos(NorthCMSzDeltaPhi);
0320 Pz_mu[imu] =
0321 SimTree->particle__Px[imu] * cos(NorthCMSzDeltaPhi) + SimTree->particle__Py[imu] * sin(NorthCMSzDeltaPhi);
0322 Py_mu[imu] = -SimTree->particle__Pz[imu];
0323 P_mu[imu] = sqrt(Px_mu[imu] * Px_mu[imu] + Py_mu[imu] * Py_mu[imu] + Pz_mu[imu] * Pz_mu[imu]);
0324
0325 if (P_mu[imu] < MinP_CMS && AcptAllMu == false)
0326 continue;
0327 else if (SimTree->particle__ParticleID[imu] != 5 && SimTree->particle__ParticleID[imu] != 6)
0328 continue;
0329 else
0330 NmuPmin++;
0331
0332 for (int jmu = 0; jmu < imu; ++jmu) {
0333 if (P_mu[jmu] < MinP_CMS && AcptAllMu == false)
0334 continue;
0335 if (SimTree->particle__ParticleID[imu] != 5 && SimTree->particle__ParticleID[imu] != 6)
0336 continue;
0337 MuMuDist = sqrt((SimTree->particle__x[imu] - SimTree->particle__x[jmu]) *
0338 (SimTree->particle__x[imu] - SimTree->particle__x[jmu]) +
0339 (SimTree->particle__y[imu] - SimTree->particle__y[jmu]) *
0340 (SimTree->particle__y[imu] - SimTree->particle__y[jmu])) *
0341 10.;
0342 if (MuMuDist < MinDist)
0343 MinDist = MuMuDist;
0344 if (MuMuDist < 2. * Target3dRadius)
0345 MuInMaxDist = true;
0346 }
0347 }
0348 if (MultiMuonNmin >= 2) {
0349 if (MuInMaxDist) {
0350 NcloseMultiMuonEvents++;
0351 } else {
0352 std::cout << "CosmicMuonGenerator.cc: Warning! No muon pair closer than " << 2. * Target3dRadius / 1000.
0353 << "m MinDist=" << MinDist / 1000. << "m at surface" << std::endl;
0354 std::cout << "Fraction of too wide opening angle multi muon events: "
0355 << 1 - double(NcloseMultiMuonEvents) / SimTree_jentry << std::endl;
0356 std::cout << "NcloseMultiMuonEvents=" << NcloseMultiMuonEvents << std::endl;
0357 std::cout << "trying next event from file" << std::endl;
0358 NskippedMultiMuonEvents++;
0359 continue;
0360 }
0361 }
0362
0363 if (NmuPmin < MultiMuonNmin && AcptAllMu == false) {
0364 NskippedMultiMuonEvents++;
0365 continue;
0366 }
0367
0368 if (Debug)
0369 if (MultiMuonNmin >= 2)
0370 std::cout << "start trial do loop: MuMuDist=" << MinDist / 1000. << "[m] Nmuons=" << nmuons
0371 << " NcloseMultiMuonEvents=" << NcloseMultiMuonEvents
0372 << " NskippedMultiMuonEvents=" << NskippedMultiMuonEvents << std::endl;
0373
0374
0375 Id_at = SimTree->shower_EventID;
0376
0377 double M_at = 0.;
0378
0379 Id_at = 2212;
0380 M_at = 938.272e-3;
0381
0382
0383 E_at = SimTree->shower_Energy;
0384 Theta_at = SimTree->shower_Theta;
0385 double phi_at = SimTree->shower_Phi - NorthCMSzDeltaPhi;
0386 if (phi_at < -Pi)
0387 phi_at += TwoPi;
0388 else if (phi_at > Pi)
0389 phi_at -= TwoPi;
0390 double P_at = sqrt(E_at * E_at - M_at * M_at);
0391
0392
0393 Px_at = P_at * sin(Theta_at) * sin(phi_at);
0394 Py_at = -P_at * cos(Theta_at);
0395 Pz_at = P_at * sin(Theta_at) * cos(phi_at);
0396
0397
0398 double theta_mu_max = Theta_at;
0399 double theta_mu_min = Theta_at;
0400
0401 double phi_rel_min = 0.;
0402 double phi_rel_max = 0.;
0403
0404
0405
0406 Theta_mu.resize(nmuons);
0407 for (int imu = 0; imu < nmuons; ++imu) {
0408 Theta_mu[imu] = acos(-Py_mu[imu] / P_mu[imu]);
0409 if (Theta_mu[imu] > theta_mu_max)
0410 theta_mu_max = Theta_mu[imu];
0411 if (Theta_mu[imu] < theta_mu_min)
0412 theta_mu_min = Theta_mu[imu];
0413
0414 double phi_mu = atan2(Px_mu[imu], Pz_mu[imu]);
0415 double phi_rel = phi_mu - phi_at;
0416 if (phi_rel < -Pi)
0417 phi_rel += TwoPi;
0418 else if (phi_rel > Pi)
0419 phi_rel -= TwoPi;
0420 if (phi_rel < phi_rel_min)
0421 phi_rel_min = phi_rel;
0422 else if (phi_rel > phi_rel_max)
0423 phi_rel_max = phi_rel;
0424 }
0425
0426 double h_sf = SurfaceOfEarth + PlugWidth;
0427
0428 double R_at = h_sf * tan(Theta_at);
0429
0430 double JdRxzV_dR_trans = 1.;
0431 double JdPhiV_dPhi_trans = 1.;
0432 double JdR_trans_sqrt = 1.;
0433
0434
0435 double R_mu_max = (h_sf + Target3dRadius) * tan(theta_mu_max);
0436 double R_max = std::min(SurfaceRadius, R_mu_max);
0437 double R_mu_min = (h_sf - Target3dRadius) * tan(theta_mu_min);
0438 double R_min = std::max(0., R_mu_min);
0439
0440 if (R_at > SurfaceRadius) {
0441 std::cout << "CosmicMuonGenerator.cc: Warning! R_at=" << R_at << " > SurfaceRadius=" << SurfaceRadius
0442 << std::endl;
0443 }
0444
0445
0446
0447
0448
0449 double psR1min = R_min + 0.25 * (R_max - R_min);
0450 double psR1max = std::min(SurfaceRadius, R_max - 0.25 * (R_max - R_min));
0451 double psR1 = psR1max - psR1min;
0452
0453 double psR2min = R_min;
0454 double psR2max = R_max;
0455 double psR2 = psR2max - psR2min;
0456
0457 double psR3min = 0.;
0458 double psR3max = SurfaceRadius;
0459 double psR3 = psR3max - psR3min;
0460
0461
0462 double psRall = psR3;
0463
0464 double fR1 = psR1 / psRall, fR2 = psR2 / psRall, fR3 = psR3 / psRall;
0465 double pR1 = 0.25, pR2 = 0.7, pR3 = 0.05;
0466
0467
0468 double psPh1 = 0.5 * (phi_rel_max - phi_rel_min);
0469 double psPh2 = phi_rel_max - phi_rel_min;
0470 double psPh3 = TwoPi;
0471 double psPhall = psPh3;
0472
0473 double fPh1 = psPh1 / psPhall, fPh2 = psPh2 / psPhall,
0474 fPh3 = psPh3 / psPhall;
0475
0476 double pPh1 = 0.25, pPh2 = 0.7, pPh3 = 0.05;
0477
0478 Trials = 0;
0479 double trials = 0.;
0480 Vx_mu.resize(nmuons);
0481 Vy_mu.resize(nmuons);
0482 Vz_mu.resize(nmuons);
0483 int NmuHitTarget = 0;
0484 while (NmuHitTarget < MultiMuonNmin) {
0485 NmuHitTarget = 0;
0486 double Nver = 0.;
0487
0488
0489 double RxzV;
0490 double which_R_class = RanGen->flat();
0491 if (which_R_class < pR1) {
0492 RxzV = psR1min + psR1 * RanGen->flat();
0493 JdRxzV_dR_trans = fR1 / pR1 * SurfaceRadius / psR1;
0494 } else if (which_R_class < pR1 + pR2) {
0495 RxzV = psR2min + psR2 * RanGen->flat();
0496 JdRxzV_dR_trans = fR2 / pR2 * SurfaceRadius / psR2;
0497 } else {
0498 RxzV = psR3min + psR3 * RanGen->flat();
0499 JdRxzV_dR_trans = fR3 / pR3 * SurfaceRadius / psR3;
0500 }
0501
0502 JdR_trans_sqrt = 2. * RxzV / SurfaceRadius;
0503
0504
0505 double PhiV;
0506 double which_phi_class = RanGen->flat();
0507 if (which_phi_class < pPh1) {
0508 PhiV = phi_at + phi_rel_min + psPh1 * RanGen->flat();
0509 JdPhiV_dPhi_trans = fPh1 / pPh1 * TwoPi / psPh1;
0510 } else if (which_phi_class < pPh1 + pPh2) {
0511 PhiV = phi_at + phi_rel_min + psPh2 * RanGen->flat();
0512 JdPhiV_dPhi_trans = fPh2 / pPh2 * TwoPi / psPh2;
0513 } else {
0514 PhiV = phi_at + phi_rel_min + psPh3 * RanGen->flat();
0515 JdPhiV_dPhi_trans = fPh3 / pPh3 * TwoPi / psPh3;
0516 }
0517
0518
0519 if (PhiV < -Pi)
0520 PhiV += TwoPi;
0521 else if (PhiV > Pi)
0522 PhiV -= TwoPi;
0523
0524 Nver++;
0525 trials += JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans;
0526 Trials++;
0527 if (trials > max_Trials)
0528 break;
0529 Ngen += (Nver - 1.);
0530
0531 Vx_at = RxzV * sin(PhiV);
0532
0533 Vy_at = h_sf;
0534
0535
0536 Vz_at = RxzV * cos(PhiV);
0537
0538 int NmuHitTargetSphere = 0;
0539 for (int imu = 0; imu < nmuons; ++imu) {
0540 Vx_mu[imu] = Vx_at + (-SimTree->particle__x[imu] * sin(NorthCMSzDeltaPhi) +
0541 SimTree->particle__y[imu] * cos(NorthCMSzDeltaPhi)) *
0542 10;
0543 Vy_mu[imu] = h_sf;
0544 Vz_mu[imu] = Vz_at + (SimTree->particle__x[imu] * cos(NorthCMSzDeltaPhi) +
0545 SimTree->particle__y[imu] * sin(NorthCMSzDeltaPhi)) *
0546 10;
0547
0548
0549 double pt_sec = sqrt(Px_mu[imu] * Px_mu[imu] + Pz_mu[imu] * Pz_mu[imu]);
0550 double theta_sec = atan2(std::fabs(Py_mu[imu]), pt_sec);
0551 double r_sec = sqrt((Vx_mu[imu] - Vx_at) * (Vx_mu[imu] - Vx_at) + (Vz_mu[imu] - Vz_at) * (Vz_mu[imu] - Vz_at));
0552 double h_prod = r_sec * tan(theta_sec);
0553 if (h_prod + h_sf > Vy_at)
0554 Vy_at = h_prod + h_sf;
0555
0556
0557 if (SimTree->particle__ParticleID[imu] != 5 && SimTree->particle__ParticleID[imu] != 6)
0558 continue;
0559
0560 if (P_mu[imu] < MinP_CMS && AcptAllMu == false)
0561 continue;
0562
0563
0564 double Vxz_mu = sqrt(Vx_mu[imu] * Vx_mu[imu] + Vz_mu[imu] * Vz_mu[imu]);
0565 theta_mu_max = atan((Vxz_mu + Target3dRadius) / (h_sf - Target3dRadius));
0566 theta_mu_min = atan((Vxz_mu - Target3dRadius) / (h_sf + Target3dRadius));
0567 if (Theta_mu[imu] > theta_mu_min && Theta_mu[imu] < theta_mu_max) {
0568
0569 double dPhi = Pi;
0570 if (Vxz_mu > Target3dRadius)
0571 dPhi = asin(Target3dRadius / Vxz_mu);
0572 double PhiPmu = atan2(Px_mu[imu], Pz_mu[imu]);
0573 double PhiVmu = atan2(Vx_mu[imu], Vz_mu[imu]);
0574 double rotPhi = PhiVmu + Pi;
0575 if (rotPhi > Pi)
0576 rotPhi -= TwoPi;
0577 double disPhi = std::fabs(rotPhi - PhiPmu);
0578 if (disPhi > Pi)
0579 disPhi = TwoPi - disPhi;
0580 if (disPhi < dPhi) {
0581 NmuHitTargetSphere++;
0582 }
0583 }
0584
0585 }
0586
0587 if (NmuHitTargetSphere < MultiMuonNmin)
0588 continue;
0589
0590
0591 Id_sf.clear();
0592 Px_sf.clear();
0593 Py_sf.clear();
0594 Pz_sf.clear();
0595 E_sf.clear();
0596
0597 Vx_sf.clear();
0598 Vy_sf.clear();
0599 Vz_sf.clear();
0600 T0_sf.clear();
0601
0602
0603 Id_ug.clear();
0604 Px_ug.clear();
0605 Py_ug.clear();
0606 Pz_ug.clear();
0607 E_ug.clear();
0608
0609 Vx_ug.clear();
0610 Vy_ug.clear();
0611 Vz_ug.clear();
0612 T0_ug.clear();
0613
0614 int Id_sf_this = 0;
0615 double Px_sf_this = 0., Py_sf_this = 0., Pz_sf_this = 0.;
0616 double E_sf_this = 0.;
0617
0618 double Vx_sf_this = 0., Vy_sf_this = 0., Vz_sf_this = 0.;
0619 double T0_sf_this = 0.;
0620
0621 T0_at = SimTree->shower_GH_t0 * SpeedOfLight;
0622
0623 for (int imu = 0; imu < nmuons; ++imu) {
0624 if (P_mu[imu] < MinP_CMS && AcptAllMu == false)
0625 continue;
0626
0627 if (SimTree->particle__ParticleID[imu] != 5 && SimTree->particle__ParticleID[imu] != 6)
0628 continue;
0629
0630 Id_sf_this = SimTree->particle__ParticleID[imu];
0631 if (Id_sf_this == 5)
0632 Id_sf_this = -13;
0633 else if (Id_sf_this == 6)
0634 Id_sf_this = 13;
0635
0636 else if (Id_sf_this == 1)
0637 Id_sf_this = 22;
0638 else if (Id_sf_this == 2)
0639 Id_sf_this = -11;
0640 else if (Id_sf_this == 3)
0641 Id_sf_this = 11;
0642 else if (Id_sf_this == 7)
0643 Id_sf_this = 111;
0644 else if (Id_sf_this == 8)
0645 Id_sf_this = 211;
0646 else if (Id_sf_this == 9)
0647 Id_sf_this = -211;
0648 else if (Id_sf_this == 10)
0649 Id_sf_this = 130;
0650 else if (Id_sf_this == 11)
0651 Id_sf_this = 321;
0652 else if (Id_sf_this == 12)
0653 Id_sf_this = -321;
0654 else if (Id_sf_this == 13)
0655 Id_sf_this = 2112;
0656 else if (Id_sf_this == 14)
0657 Id_sf_this = 2212;
0658 else if (Id_sf_this == 15)
0659 Id_sf_this = -2212;
0660 else if (Id_sf_this == 16)
0661 Id_sf_this = 310;
0662 else if (Id_sf_this == 17)
0663 Id_sf_this = 221;
0664 else if (Id_sf_this == 18)
0665 Id_sf_this = 3122;
0666
0667 else {
0668 std::cout << "CosmicMuonGenerator.cc: Warning! Muon Id=" << Id_sf_this << " from file read in" << std::endl;
0669 Id_sf_this = 99999;
0670 }
0671
0672 T0_sf_this = SimTree->particle__Time[imu] * SpeedOfLight;
0673
0674 Px_sf_this = Px_mu[imu];
0675 Py_sf_this = Py_mu[imu];
0676 Pz_sf_this = Pz_mu[imu];
0677 E_sf_this = sqrt(P_mu[imu] * P_mu[imu] + MuonMass * MuonMass);
0678 Vx_sf_this = Vx_mu[imu];
0679 Vy_sf_this = Vy_mu[imu];
0680 Vz_sf_this = Vz_mu[imu];
0681
0682 OneMuoEvt.create(Id_sf_this,
0683 Px_sf_this,
0684 Py_sf_this,
0685 Pz_sf_this,
0686 E_sf_this,
0687 MuonMass,
0688 Vx_sf_this,
0689 Vy_sf_this,
0690 Vz_sf_this,
0691 T0_sf_this);
0692
0693 if (goodOrientation()) {
0694 OneMuoEvt.propagate(ElossScaleFactor, RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
0695 }
0696
0697 if ((OneMuoEvt.hitTarget() && sqrt(OneMuoEvt.e() * OneMuoEvt.e() - MuonMass * MuonMass) > MinP_CMS) ||
0698 AcptAllMu == true) {
0699 Id_sf.push_back(Id_sf_this);
0700 Px_sf.push_back(Px_sf_this);
0701 Py_sf.push_back(Py_sf_this);
0702 Pz_sf.push_back(Pz_sf_this);
0703 E_sf.push_back(E_sf_this);
0704
0705 Vx_sf.push_back(Vx_sf_this);
0706 Vy_sf.push_back(Vy_sf_this);
0707 Vz_sf.push_back(Vz_sf_this);
0708 T0_sf.push_back(T0_sf_this);
0709
0710
0711 Id_ug.push_back(OneMuoEvt.id());
0712 Px_ug.push_back(OneMuoEvt.px());
0713 Py_ug.push_back(OneMuoEvt.py());
0714 Pz_ug.push_back(OneMuoEvt.pz());
0715 E_ug.push_back(OneMuoEvt.e());
0716
0717 Vx_ug.push_back(OneMuoEvt.vx());
0718 Vy_ug.push_back(OneMuoEvt.vy());
0719 Vz_ug.push_back(OneMuoEvt.vz());
0720 T0_ug.push_back(OneMuoEvt.t0());
0721
0722 NmuHitTarget++;
0723 }
0724 }
0725
0726 }
0727
0728 if (trials > max_Trials) {
0729 std::cout << "CosmicMuonGenerator.cc: Warning! trials reach max_trials=" << max_Trials
0730 << " without accepting event!" << std::endl;
0731 if (Debug) {
0732 std::cout << " N(mu)=" << Id_ug.size();
0733 if (!Id_ug.empty())
0734 std::cout << " E[0]=" << E_ug[0] << " theta="
0735 << acos(-Py_ug[0] / sqrt(Px_ug[0] * Px_ug[0] + Py_ug[0] * Py_ug[0] + Pz_ug[0] * Pz_ug[0]))
0736 << " R_xz=" << sqrt(Vx_sf[0] * Vx_sf[0] + Vy_sf[0] * Vy_sf[0]);
0737 std::cout << " Theta_at=" << Theta_at << std::endl;
0738 }
0739 std::cout << "Unweighted int num of Trials = " << Trials << std::endl;
0740 std::cout << "trying next event (" << SimTree_jentry << ") from file" << std::endl;
0741 NskippedMultiMuonEvents++;
0742 continue;
0743 } else {
0744 if (NmuHitTarget < MultiMuonNmin) {
0745 std::cout << "CosmicMuonGenerator.cc: Warning! less than " << MultiMuonNmin << " muons hit target: N(mu=)"
0746 << NmuHitTarget << std::endl;
0747 std::cout << "trying next event (" << SimTree_jentry << ") from file" << std::endl;
0748 NskippedMultiMuonEvents++;
0749 continue;
0750 } else {
0751
0752
0753
0754 double T0_ug_min, T0_ug_max;
0755 T0_ug_min = T0_ug_max = T0_ug[0];
0756 double Tbox = (MaxT0 - MinT0) * SpeedOfLight;
0757 double minDeltaT0 = 2 * Tbox;
0758 for (unsigned int imu = 0; imu < Id_ug.size(); ++imu) {
0759 double T0_this = T0_ug[imu];
0760 if (T0_this < T0_ug_min)
0761 T0_ug_min = T0_this;
0762 if (T0_this > T0_ug_max)
0763 T0_ug_max = T0_this;
0764 if (Debug)
0765 std::cout << "imu=" << imu << " T0_this=" << T0_this
0766 << " P=" << sqrt(pow(Px_ug[imu], 2) + pow(Py_ug[imu], 2) + pow(Pz_ug[imu], 2)) << std::endl;
0767 for (unsigned int jmu = 0; jmu < imu; ++jmu) {
0768 if (std::fabs(T0_ug[imu] - T0_ug[jmu]) < minDeltaT0)
0769 minDeltaT0 = std::fabs(T0_ug[imu] - T0_ug[jmu]);
0770 }
0771 }
0772
0773 if (int(Id_ug.size()) >= MultiMuonNmin && MultiMuonNmin >= 2 && minDeltaT0 > Tbox)
0774 continue;
0775
0776 double T0_min = T0_ug_min + MinT0 * SpeedOfLight;
0777 double T0_max = T0_ug_max + MaxT0 * SpeedOfLight;
0778
0779
0780 int TboxTrials = 0;
0781 int NmuInTbox;
0782 double T0_offset, T0diff;
0783 for (int tboxtrial = 0; tboxtrial < 1000; ++tboxtrial) {
0784 T0_offset = RanGen->flat() * (T0_max - T0_min);
0785 TboxTrials++;
0786 T0diff = T0_offset - T0_max;
0787 NmuInTbox = 0;
0788 for (unsigned int imu = 0; imu < Id_ug.size(); ++imu) {
0789 if (T0_ug[imu] + T0diff > MinT0 * SpeedOfLight && T0_ug[imu] + T0diff < MaxT0 * SpeedOfLight)
0790 NmuInTbox++;
0791 }
0792 if (NmuInTbox >= MultiMuonNmin)
0793 break;
0794 }
0795 if (NmuInTbox < MultiMuonNmin)
0796 continue;
0797
0798 if (Debug)
0799 std::cout << "initial T0_at=" << T0_at << " T0_min=" << T0_min << " T0_max=" << T0_max
0800 << " T0_offset=" << T0_offset;
0801 T0_at += T0diff;
0802 if (Debug)
0803 std::cout << " T0diff=" << T0diff << std::endl;
0804 for (unsigned int imu = 0; imu < Id_ug.size(); ++imu) {
0805 if (Debug)
0806 std::cout << "before: T0_sf[" << imu << "]=" << T0_sf[imu] << " T0_ug=" << T0_ug[imu];
0807 T0_sf[imu] += T0diff;
0808 T0_ug[imu] += T0diff;
0809 if (Debug)
0810 std::cout << " after: T0_sf[" << imu << "]=" << T0_sf[imu] << " T0_ug=" << T0_ug[imu] << std::endl;
0811 }
0812 if (Debug)
0813 std::cout << "T0diff=" << T0diff << " T0_at=" << T0_at << std::endl;
0814
0815 Nsel += 1;
0816 EventWeight = JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans / (trials * TboxTrials);
0817 EvtRejected = false;
0818 if (Debug)
0819 std::cout << "CosmicMuonGenerator.cc: Theta_at=" << Theta_at << " phi_at=" << phi_at << " Px_at=" << Px_at
0820 << " Py_at=" << Py_at << " Pz_at=" << Pz_at << " T0_at=" << T0_at << " Vx_at=" << Vx_at
0821 << " Vy_at=" << Vy_at << " Vz_at=" << Vz_at << " EventWeight=" << EventWeight
0822 << " Nmuons=" << Id_sf.size() << std::endl;
0823 }
0824 }
0825
0826 }
0827
0828 return true;
0829 }
0830
0831 void CosmicMuonGenerator::terminate() {
0832 if (NumberOfEvents > 0) {
0833 std::cout << std::endl;
0834 std::cout << "*********************************************************" << std::endl;
0835 std::cout << "*********************************************************" << std::endl;
0836 std::cout << "*** ***" << std::endl;
0837 std::cout << "*** C O S M I C M U O N S T A T I S T I C S ***" << std::endl;
0838 std::cout << "*** ***" << std::endl;
0839 std::cout << "*********************************************************" << std::endl;
0840 std::cout << "*********************************************************" << std::endl;
0841 std::cout << std::endl;
0842 std::cout << " number of initial cosmic events: " << int(Ngen) << std::endl;
0843 std::cout << " number of actually diced events: " << int(Ndiced) << std::endl;
0844 std::cout << " number of generated and accepted events: " << int(Nsel) << std::endl;
0845 double selEff = Nsel / Ngen;
0846 std::cout << " event selection efficiency: " << selEff * 100. << "%" << std::endl;
0847 int n100cos =
0848 Norm->events_n100cos(0., 0.);
0849 std::cout << " events with ~100 GeV and 1 - cos(theta) < 1/2pi: " << n100cos << std::endl;
0850 std::cout << std::endl;
0851 std::cout << " momentum range: " << MinP << " ... " << MaxP << " GeV" << std::endl;
0852 std::cout << " theta range: " << MinTheta * Rad2Deg << " ... " << MaxTheta * Rad2Deg << " deg"
0853 << std::endl;
0854 std::cout << " phi range: " << MinPhi * Rad2Deg << " ... " << MaxPhi * Rad2Deg << " deg" << std::endl;
0855 std::cout << " time range: " << MinT0 << " ... " << MaxT0 << " ns" << std::endl;
0856 std::cout << " energy loss: " << ElossScaleFactor * 100. << "%" << std::endl;
0857 std::cout << std::endl;
0858 double area = 1.e-6 * Pi * SurfaceRadius * SurfaceRadius;
0859 if (MinTheta > 90. * Deg2Rad)
0860 std::cout << " area of initial cosmics at CMS detector bottom surface: " << area << " m^2" << std::endl;
0861 else
0862 std::cout << " area of initial cosmics on Surface + PlugWidth: " << area << " m^2" << std::endl;
0863 std::cout << " depth of CMS detector (from Surface): " << SurfaceOfEarth / 1000 << " m" << std::endl;
0864
0865
0866
0867
0868 if ((n100cos > 0 && MaxTheta < 84.26 * Deg2Rad) || MinTheta > 90. * Deg2Rad) {
0869
0870
0871
0872
0873
0874
0875
0876 if (MinTheta > 90. * Deg2Rad) {
0877 double Omega = (cos(MinTheta) - cos(MaxTheta)) * (MaxPhi - MinPhi);
0878
0879 EventRate = (Ndiced * 3.e-13) * Omega * 4. * RadiusOfTarget * ZDistOfTarget * 1.e-2 *
0880 selEff;
0881 rateErr_stat = EventRate / sqrt((double)Ndiced);
0882 rateErr_syst = EventRate / 3.e-13 * 1.0e-13;
0883 } else {
0884 EventRate = (Ndiced * Norm->norm(n100cos)) * area * selEff;
0885 rateErr_stat = EventRate / sqrt((double)n100cos);
0886 rateErr_syst = EventRate / 2.63e-3 * 0.06e-3;
0887 }
0888
0889
0890 if (MaxTheta < 0.572) {
0891 double spacean = 2. * Pi * (1. - cos(MaxTheta));
0892 EventRate = (Ndiced * Norm->norm(n100cos)) * area * selEff * spacean;
0893 rateErr_stat = EventRate / sqrt((double)n100cos);
0894 rateErr_syst = EventRate / 2.63e-3 * 0.06e-3;
0895 }
0896
0897 } else {
0898 EventRate = Nsel;
0899 rateErr_stat = Nsel;
0900 rateErr_syst = Nsel;
0901 std::cout << std::endl;
0902 if (MinP > 100.)
0903 std::cout << " !!! MinP > 100 GeV. Cannot apply normalisation!" << std::endl;
0904 else if (MaxTheta > 84.26 * Deg2Rad)
0905 std::cout << " !!! Note: generated cosmics exceed parameterisation. No flux calculated!" << std::endl;
0906
0907 else
0908 std::cout << " !!! Not enough statistics to apply normalisation (rate=1 +- 1) !!!" << std::endl;
0909 }
0910
0911 std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0912 std::cout << " rate is " << EventRate << " +-" << rateErr_stat << " (stat) "
0913 << "+-" << rateErr_syst << " (syst) "
0914 << " muons per second" << std::endl;
0915 if (EventRate != 0)
0916 std::cout << " number of events corresponds to " << Nsel / EventRate << " s"
0917 << std::endl;
0918 std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0919 std::cout << std::endl;
0920 std::cout << "*********************************************************" << std::endl;
0921 std::cout << "*********************************************************" << std::endl;
0922 }
0923 }
0924
0925 void CosmicMuonGenerator::checkIn() {
0926 if (MinP < 0.) {
0927 NumberOfEvents = 0;
0928 std::cout << " CMG-ERR: min.energy is out of range (0 GeV ... inf]" << std::endl << std::endl;
0929 }
0930 if (MaxP < 0.) {
0931 NumberOfEvents = 0;
0932 std::cout << " CMG-ERR: max.energy is out of range (0 GeV ... inf]" << std::endl << std::endl;
0933 }
0934 if (MaxP <= MinP) {
0935 NumberOfEvents = 0;
0936 std::cout << " CMG-ERR: max.energy is not greater than min.energy" << std::endl << std::endl;
0937 }
0938 if (MinTheta < 0.) {
0939 NumberOfEvents = 0;
0940 std::cout << " CMG-ERR: min.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl;
0941 }
0942 if (MaxTheta < 0.) {
0943 NumberOfEvents = 0;
0944 std::cout << " CMG-ERR: max.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl;
0945 }
0946 if (MaxTheta <= MinTheta) {
0947 NumberOfEvents = 0;
0948 std::cout << " CMG-ERR: max.theta is not greater than min.theta" << std::endl << std::endl;
0949 }
0950 if (MinPhi < 0.) {
0951 NumberOfEvents = 0;
0952 std::cout << " CMG-ERR: min.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl;
0953 }
0954 if (MaxPhi < 0.) {
0955 NumberOfEvents = 0;
0956 std::cout << " CMG-ERR: max.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl;
0957 }
0958 if (MaxPhi <= MinPhi) {
0959 NumberOfEvents = 0;
0960 std::cout << " CMG-ERR: max.phi is not greater than min.phi" << std::endl << std::endl;
0961 }
0962 if (MaxT0 <= MinT0) {
0963 NumberOfEvents = 0;
0964 std::cout << " CMG-ERR: max.t0 is not greater than min.t0" << std::endl << std::endl;
0965 }
0966 if (ElossScaleFactor < 0.) {
0967 NumberOfEvents = 0;
0968 std::cout << " CMG-ERR: E-loss scale factor is out of range [0 ... inf)" << std::endl << std::endl;
0969 }
0970 if (MinEnu < 0.) {
0971 NumberOfEvents = 0;
0972 std::cout << " CMG-ERR: min.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl;
0973 }
0974 if (MaxEnu < 0.) {
0975 NumberOfEvents = 0;
0976 std::cout << " CMG-ERR: max.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl;
0977 }
0978 if (MaxEnu <= MinEnu) {
0979 NumberOfEvents = 0;
0980 std::cout << " CMG-ERR: max.Enu is not greater than min.Enu" << std::endl << std::endl;
0981 }
0982 }
0983
0984 bool CosmicMuonGenerator::goodOrientation() {
0985
0986 bool goodAngles = false;
0987 bool phiaccepted = false;
0988 bool thetaaccepted = false;
0989 double RxzV = sqrt(OneMuoEvt.vx() * OneMuoEvt.vx() + OneMuoEvt.vz() * OneMuoEvt.vz());
0990
0991 double rVY;
0992 if (MinTheta > 90. * Deg2Rad)
0993 rVY = -sqrt(RxzV * RxzV + RadiusCMS * RadiusCMS);
0994 else
0995 rVY = sqrt(RxzV * RxzV + (SurfaceOfEarth + PlugWidth) * (SurfaceOfEarth + PlugWidth));
0996
0997 double Phi = OneMuoEvt.phi();
0998 double PhiV = atan2(OneMuoEvt.vx(), OneMuoEvt.vz()) + Pi;
0999 if (PhiV > TwoPi)
1000 PhiV -= TwoPi;
1001 double disPhi = std::fabs(PhiV - Phi);
1002 if (disPhi > Pi)
1003 disPhi = TwoPi - disPhi;
1004 double dPhi = Pi;
1005 if (RxzV > Target3dRadius)
1006 dPhi = asin(Target3dRadius / RxzV);
1007 if (disPhi < dPhi)
1008 phiaccepted = true;
1009 double Theta = OneMuoEvt.theta();
1010 double ThetaV = asin(RxzV / rVY);
1011 double dTheta = Pi;
1012 if (std::fabs(rVY) > Target3dRadius)
1013 dTheta = asin(Target3dRadius / std::fabs(rVY));
1014
1015
1016
1017 if (!phiaccepted && RxzV < Target3dRadius)
1018
1019 std::cout << "Rejected phi=" << Phi << " PhiV=" << PhiV << " dPhi=" << dPhi << " disPhi=" << disPhi
1020 << " RxzV=" << RxzV << " Target3dRadius=" << Target3dRadius << " Theta=" << Theta << std::endl;
1021
1022 if (std::fabs(Theta - ThetaV) < dTheta)
1023 thetaaccepted = true;
1024 if (phiaccepted && thetaaccepted)
1025 goodAngles = true;
1026 return goodAngles;
1027 }
1028
1029 void CosmicMuonGenerator::initEvDis() {
1030 #if ROOT_INTERACTIVE
1031 float rCMS = RadiusCMS / 1000.;
1032 float zCMS = Z_DistCMS / 1000.;
1033 if (TrackerOnly == true) {
1034 rCMS = RadiusTracker / 1000.;
1035 zCMS = Z_DistTracker / 1000.;
1036 }
1037 TH2F* disXY = new TH2F("disXY", "X-Y view", 160, -rCMS, rCMS, 160, -rCMS, rCMS);
1038 TH2F* disZY = new TH2F("disZY", "Z-Y view", 150, -zCMS, zCMS, 160, -rCMS, rCMS);
1039 gStyle->SetPalette(1, 0);
1040 gStyle->SetMarkerColor(1);
1041 gStyle->SetMarkerSize(1.5);
1042 TCanvas* disC = new TCanvas("disC", "Cosmic Muon Event Display", 0, 0, 800, 410);
1043 disC->Divide(2, 1);
1044 disC->cd(1);
1045 gPad->SetTicks(1, 1);
1046 disXY->SetMinimum(log10(MinP));
1047 disXY->SetMaximum(log10(MaxP));
1048 disXY->GetXaxis()->SetLabelSize(0.05);
1049 disXY->GetXaxis()->SetTitleSize(0.05);
1050 disXY->GetXaxis()->SetTitleOffset(1.0);
1051 disXY->GetXaxis()->SetTitle("X [m]");
1052 disXY->GetYaxis()->SetLabelSize(0.05);
1053 disXY->GetYaxis()->SetTitleSize(0.05);
1054 disXY->GetYaxis()->SetTitleOffset(0.8);
1055 disXY->GetYaxis()->SetTitle("Y [m]");
1056 disC->cd(2);
1057 gPad->SetGrid(1, 1);
1058 gPad->SetTicks(1, 1);
1059 disZY->SetMinimum(log10(MinP));
1060 disZY->SetMaximum(log10(MaxP));
1061 disZY->GetXaxis()->SetLabelSize(0.05);
1062 disZY->GetXaxis()->SetTitleSize(0.05);
1063 disZY->GetXaxis()->SetTitleOffset(1.0);
1064 disZY->GetXaxis()->SetTitle("Z [m]");
1065 disZY->GetYaxis()->SetLabelSize(0.05);
1066 disZY->GetYaxis()->SetTitleSize(0.05);
1067 disZY->GetYaxis()->SetTitleOffset(0.8);
1068 disZY->GetYaxis()->SetTitle("Y [m]");
1069 #endif
1070 }
1071
1072 void CosmicMuonGenerator::displayEv() {
1073 #if ROOT_INTERACTIVE
1074 double RadiusDet = RadiusCMS;
1075 double Z_DistDet = Z_DistCMS;
1076 if (TrackerOnly == true) {
1077 RadiusDet = RadiusTracker;
1078 Z_DistDet = Z_DistTracker;
1079 }
1080 disXY->Reset();
1081 disZY->Reset();
1082 TMarker* InteractionPoint = new TMarker(0., 0., 2);
1083 TArc* r8m = new TArc(0., 0., (RadiusDet / 1000.));
1084 TLatex* logEaxis = new TLatex();
1085 logEaxis->SetTextSize(0.05);
1086 float energy = float(OneMuoEvt.e());
1087 float verX = float(OneMuoEvt.vx() / 1000.);
1088 float verY = float(OneMuoEvt.vy() / 1000.);
1089 float verZ = float(OneMuoEvt.vz() / 1000.);
1090 float dirX = float(OneMuoEvt.px()) / std::fabs(OneMuoEvt.py());
1091 float dirY = float(OneMuoEvt.py()) / std::fabs(OneMuoEvt.py());
1092 float dirZ = float(OneMuoEvt.pz()) / std::fabs(OneMuoEvt.py());
1093 float yStep = disXY->GetYaxis()->GetBinWidth(1);
1094 int NbinY = disXY->GetYaxis()->GetNbins();
1095 for (int iy = 0; iy < NbinY; ++iy) {
1096 verX += dirX * yStep;
1097 verY += dirY * yStep;
1098 verZ += dirZ * yStep;
1099 float rXY = sqrt(verX * verX + verY * verY) * 1000.;
1100 float absZ = std::fabs(verZ) * 1000.;
1101 if (rXY < RadiusDet && absZ < Z_DistDet) {
1102 disXY->Fill(verX, verY, log10(energy));
1103 disZY->Fill(verZ, verY, log10(energy));
1104 disC->cd(1);
1105 disXY->Draw("COLZ");
1106 InteractionPoint->Draw("SAME");
1107 r8m->Draw("SAME");
1108 logEaxis->DrawLatex((0.65 * RadiusDet / 1000.), (1.08 * RadiusDet / 1000.), "log_{10}E(#mu^{#pm})");
1109 disC->cd(2);
1110 disZY->Draw("COL");
1111 InteractionPoint->Draw("SAME");
1112 gPad->Update();
1113 }
1114 }
1115 #endif
1116 }
1117
1118 void CosmicMuonGenerator::setNumberOfEvents(unsigned int N) {
1119 if (NotInitialized)
1120 NumberOfEvents = N;
1121 }
1122
1123 void CosmicMuonGenerator::setRanSeed(int N) {
1124 if (NotInitialized)
1125 RanSeed = N;
1126 }
1127
1128 void CosmicMuonGenerator::setMinP(double P) {
1129 if (NotInitialized)
1130 MinP = P;
1131 }
1132
1133 void CosmicMuonGenerator::setMinP_CMS(double P) {
1134 if (NotInitialized)
1135 MinP_CMS = P;
1136 }
1137
1138 void CosmicMuonGenerator::setMaxP(double P) {
1139 if (NotInitialized)
1140 MaxP = P;
1141 }
1142
1143 void CosmicMuonGenerator::setMinTheta(double Theta) {
1144 if (NotInitialized)
1145 MinTheta = Theta * Deg2Rad;
1146 }
1147
1148 void CosmicMuonGenerator::setMaxTheta(double Theta) {
1149 if (NotInitialized)
1150 MaxTheta = Theta * Deg2Rad;
1151 }
1152
1153 void CosmicMuonGenerator::setMinPhi(double Phi) {
1154 if (NotInitialized)
1155 MinPhi = Phi * Deg2Rad;
1156 }
1157
1158 void CosmicMuonGenerator::setMaxPhi(double Phi) {
1159 if (NotInitialized)
1160 MaxPhi = Phi * Deg2Rad;
1161 }
1162
1163 void CosmicMuonGenerator::setMinT0(double T0) {
1164 if (NotInitialized)
1165 MinT0 = T0;
1166 }
1167
1168 void CosmicMuonGenerator::setMaxT0(double T0) {
1169 if (NotInitialized)
1170 MaxT0 = T0;
1171 }
1172
1173 void CosmicMuonGenerator::setElossScaleFactor(double ElossScaleFact) {
1174 if (NotInitialized)
1175 ElossScaleFactor = ElossScaleFact;
1176 }
1177
1178 void CosmicMuonGenerator::setRadiusOfTarget(double R) {
1179 if (NotInitialized)
1180 RadiusOfTarget = R;
1181 }
1182
1183 void CosmicMuonGenerator::setZDistOfTarget(double Z) {
1184 if (NotInitialized)
1185 ZDistOfTarget = Z;
1186 }
1187
1188 void CosmicMuonGenerator::setZCentrOfTarget(double Z) {
1189 if (NotInitialized)
1190 ZCentrOfTarget = Z;
1191 }
1192
1193 void CosmicMuonGenerator::setTrackerOnly(bool Tracker) {
1194 if (NotInitialized)
1195 TrackerOnly = Tracker;
1196 }
1197
1198 void CosmicMuonGenerator::setMultiMuon(bool MultiMu) {
1199 if (NotInitialized)
1200 MultiMuon = MultiMu;
1201 }
1202 void CosmicMuonGenerator::setMultiMuonFileName(std::string MultiMuFile) {
1203 if (NotInitialized)
1204 MultiMuonFileName = MultiMuFile;
1205 }
1206 void CosmicMuonGenerator::setMultiMuonFileFirstEvent(int MultiMuFile1stEvt) {
1207 if (NotInitialized)
1208 MultiMuonFileFirstEvent = MultiMuFile1stEvt;
1209 }
1210 void CosmicMuonGenerator::setMultiMuonNmin(int MultiMuNmin) {
1211 if (NotInitialized)
1212 MultiMuonNmin = MultiMuNmin;
1213 }
1214
1215 void CosmicMuonGenerator::setTIFOnly_constant(bool TIF) {
1216 if (NotInitialized)
1217 TIFOnly_constant = TIF;
1218 }
1219
1220 void CosmicMuonGenerator::setTIFOnly_linear(bool TIF) {
1221 if (NotInitialized)
1222 TIFOnly_linear = TIF;
1223 }
1224 void CosmicMuonGenerator::setMTCCHalf(bool MTCC) {
1225 if (NotInitialized)
1226 MTCCHalf = MTCC;
1227 }
1228
1229 void CosmicMuonGenerator::setPlugVx(double PlugVtx) {
1230 if (NotInitialized)
1231 PlugVx = PlugVtx;
1232 }
1233 void CosmicMuonGenerator::setPlugVz(double PlugVtz) {
1234 if (NotInitialized)
1235 PlugVz = PlugVtz;
1236 }
1237 void CosmicMuonGenerator::setRhoAir(double VarRhoAir) {
1238 if (NotInitialized)
1239 RhoAir = VarRhoAir;
1240 }
1241 void CosmicMuonGenerator::setRhoWall(double VarRhoWall) {
1242 if (NotInitialized)
1243 RhoWall = VarRhoWall;
1244 }
1245 void CosmicMuonGenerator::setRhoRock(double VarRhoRock) {
1246 if (NotInitialized)
1247 RhoRock = VarRhoRock;
1248 }
1249 void CosmicMuonGenerator::setRhoClay(double VarRhoClay) {
1250 if (NotInitialized)
1251 RhoClay = VarRhoClay;
1252 }
1253 void CosmicMuonGenerator::setRhoPlug(double VarRhoPlug) {
1254 if (NotInitialized)
1255 RhoPlug = VarRhoPlug;
1256 }
1257 void CosmicMuonGenerator::setClayWidth(double ClayLayerWidth) {
1258 if (NotInitialized)
1259 ClayWidth = ClayLayerWidth;
1260 }
1261
1262 void CosmicMuonGenerator::setMinEnu(double MinEn) {
1263 if (NotInitialized)
1264 MinEnu = MinEn;
1265 }
1266 void CosmicMuonGenerator::setMaxEnu(double MaxEn) {
1267 if (NotInitialized)
1268 MaxEnu = MaxEn;
1269 }
1270 void CosmicMuonGenerator::setNuProdAlt(double NuPrdAlt) {
1271 if (NotInitialized)
1272 NuProdAlt = NuPrdAlt;
1273 }
1274
1275 double CosmicMuonGenerator::getRate() { return EventRate; }
1276
1277 void CosmicMuonGenerator::setAcptAllMu(bool AllMu) {
1278 if (NotInitialized)
1279 AcptAllMu = AllMu;
1280 }