30 RanGen =
new CLHEP::HepJamesRandom;
73 TH1D* ene =
new TH1D(
"ene",
"generated energy", 210, 0., 1050.);
74 TH1D* the =
new TH1D(
"the",
"generated theta", 90, 0., 90.);
75 TH1D*
phi =
new TH1D(
"phi",
"generated phi", 120, 0., 360.);
76 TH3F* ver =
new TH3F(
"ver",
"Z-X-Y coordinates", 100, -25., 25., 20, -10., 10., 40, -10., 10.);
120 if (
int(
Nsel) % 100 == 0)
123 bool notSelected =
true;
124 while (notSelected) {
125 bool badMomentumGenerated =
true;
126 while (badMomentumGenerated) {
135 badMomentumGenerated =
false;
143 bool badVertexGenerated =
true;
144 while (badVertexGenerated) {
151 double rotPhi = PhiV +
Pi;
154 double disPhi = std::fabs(rotPhi -
Phi);
156 disPhi =
TwoPi - disPhi;
158 badVertexGenerated =
false;
168 double verMom = absMom *
cos(Theta);
169 double horMom = absMom *
sin(Theta);
170 double Px = horMom *
sin(
Phi);
172 double Pz = horMom *
cos(
Phi);
173 double Vx = RxzV *
sin(PhiV);
181 double Vz = RxzV *
cos(PhiV);
277 std::cout <<
"\nEntered CosmicMuonGenerator::nextMultiEvent()" << std::endl;
278 bool EvtRejected =
true;
279 bool MuInMaxDist =
false;
282 while (EvtRejected) {
286 std::cout <<
"CosmicMuonGenerator::nextMultiEvent(): SimTree_jentry=" 295 std::cout <<
"CosmicMuonGenerator.cc::nextMultiEvent: No more events in file!" << std::endl;
301 std::cout <<
"CosmicMuonGenerator.cc: Warning! Less than " <<
MultiMuonNmin <<
" muons in event!" << std::endl;
302 std::cout <<
"trying next event from file" << std::endl;
317 for (
int imu = 0; imu <
nmuons; ++imu) {
332 for (
int jmu = 0; jmu < imu; ++jmu) {
342 if (MuMuDist < MinDist)
353 <<
"m MinDist=" << MinDist / 1000. <<
"m at surface" << std::endl;
354 std::cout <<
"Fraction of too wide opening angle multi muon events: " 357 std::cout <<
"trying next event from file" << std::endl;
370 std::cout <<
"start trial do loop: MuMuDist=" << MinDist / 1000. <<
"[m] Nmuons=" <<
nmuons 388 else if (phi_at >
Pi)
401 double phi_rel_min = 0.;
402 double phi_rel_max = 0.;
407 for (
int imu = 0; imu <
nmuons; ++imu) {
414 double phi_mu = atan2(
Px_mu[imu],
Pz_mu[imu]);
415 double phi_rel = phi_mu - phi_at;
418 else if (phi_rel >
Pi)
420 if (phi_rel < phi_rel_min)
421 phi_rel_min = phi_rel;
422 else if (phi_rel > phi_rel_max)
423 phi_rel_max = phi_rel;
430 double JdRxzV_dR_trans = 1.;
431 double JdPhiV_dPhi_trans = 1.;
432 double JdR_trans_sqrt = 1.;
438 double R_min =
std::max(0., R_mu_min);
449 double psR1min = R_min + 0.25 * (R_max - R_min);
451 double psR1 = psR1max - psR1min;
453 double psR2min = R_min;
454 double psR2max = R_max;
455 double psR2 = psR2max - psR2min;
459 double psR3 = psR3max - psR3min;
462 double psRall = psR3;
464 double fR1 = psR1 / psRall, fR2 = psR2 / psRall, fR3 = psR3 / psRall;
465 double pR1 = 0.25, pR2 = 0.7, pR3 = 0.05;
468 double psPh1 = 0.5 * (phi_rel_max - phi_rel_min);
469 double psPh2 = phi_rel_max - phi_rel_min;
470 double psPh3 =
TwoPi;
471 double psPhall = psPh3;
473 double fPh1 = psPh1 / psPhall, fPh2 = psPh2 / psPhall,
474 fPh3 = psPh3 / psPhall;
476 double pPh1 = 0.25, pPh2 = 0.7, pPh3 = 0.05;
483 int NmuHitTarget = 0;
490 double which_R_class =
RanGen->flat();
491 if (which_R_class < pR1) {
492 RxzV = psR1min + psR1 *
RanGen->flat();
494 }
else if (which_R_class < pR1 + pR2) {
495 RxzV = psR2min + psR2 *
RanGen->flat();
498 RxzV = psR3min + psR3 *
RanGen->flat();
506 double which_phi_class =
RanGen->flat();
507 if (which_phi_class < pPh1) {
508 PhiV = phi_at + phi_rel_min + psPh1 *
RanGen->flat();
509 JdPhiV_dPhi_trans = fPh1 / pPh1 *
TwoPi / psPh1;
510 }
else if (which_phi_class < pPh1 + pPh2) {
511 PhiV = phi_at + phi_rel_min + psPh2 *
RanGen->flat();
512 JdPhiV_dPhi_trans = fPh2 / pPh2 *
TwoPi / psPh2;
514 PhiV = phi_at + phi_rel_min + psPh3 *
RanGen->flat();
515 JdPhiV_dPhi_trans = fPh3 / pPh3 *
TwoPi / psPh3;
525 trials += JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans;
538 int NmuHitTargetSphere = 0;
539 for (
int imu = 0; imu <
nmuons; ++imu) {
550 double theta_sec = atan2(std::fabs(
Py_mu[imu]), pt_sec);
552 double h_prod = r_sec *
tan(theta_sec);
553 if (h_prod + h_sf >
Vy_at)
554 Vy_at = h_prod + h_sf;
572 double PhiPmu = atan2(
Px_mu[imu],
Pz_mu[imu]);
573 double PhiVmu = atan2(
Vx_mu[imu],
Vz_mu[imu]);
574 double rotPhi = PhiVmu +
Pi;
577 double disPhi = std::fabs(rotPhi - PhiPmu);
579 disPhi =
TwoPi - disPhi;
581 NmuHitTargetSphere++;
615 double Px_sf_this = 0., Py_sf_this = 0., Pz_sf_this = 0.;
616 double E_sf_this = 0.;
618 double Vx_sf_this = 0., Vy_sf_this = 0., Vz_sf_this = 0.;
619 double T0_sf_this = 0.;
623 for (
int imu = 0; imu <
nmuons; ++imu) {
633 else if (Id_sf_this == 6)
636 else if (Id_sf_this == 1)
638 else if (Id_sf_this == 2)
640 else if (Id_sf_this == 3)
642 else if (Id_sf_this == 7)
644 else if (Id_sf_this == 8)
646 else if (Id_sf_this == 9)
648 else if (Id_sf_this == 10)
650 else if (Id_sf_this == 11)
652 else if (Id_sf_this == 12)
654 else if (Id_sf_this == 13)
656 else if (Id_sf_this == 14)
658 else if (Id_sf_this == 15)
660 else if (Id_sf_this == 16)
662 else if (Id_sf_this == 17)
664 else if (Id_sf_this == 18)
668 std::cout <<
"CosmicMuonGenerator.cc: Warning! Muon Id=" << Id_sf_this <<
" from file read in" << std::endl;
674 Px_sf_this =
Px_mu[imu];
675 Py_sf_this =
Py_mu[imu];
676 Pz_sf_this =
Pz_mu[imu];
678 Vx_sf_this =
Vx_mu[imu];
679 Vy_sf_this =
Vy_mu[imu];
680 Vz_sf_this =
Vz_mu[imu];
699 Id_sf.push_back(Id_sf_this);
700 Px_sf.push_back(Px_sf_this);
701 Py_sf.push_back(Py_sf_this);
702 Pz_sf.push_back(Pz_sf_this);
703 E_sf.push_back(E_sf_this);
705 Vx_sf.push_back(Vx_sf_this);
706 Vy_sf.push_back(Vy_sf_this);
707 Vz_sf.push_back(Vz_sf_this);
708 T0_sf.push_back(T0_sf_this);
730 <<
" without accepting event!" << std::endl;
745 std::cout <<
"CosmicMuonGenerator.cc: Warning! less than " <<
MultiMuonNmin <<
" muons hit target: N(mu=)" 746 << NmuHitTarget << std::endl;
754 double T0_ug_min, T0_ug_max;
755 T0_ug_min = T0_ug_max =
T0_ug[0];
757 double minDeltaT0 = 2 * Tbox;
758 for (
unsigned int imu = 0; imu <
Id_ug.size(); ++imu) {
759 double T0_this =
T0_ug[imu];
760 if (T0_this < T0_ug_min)
762 if (T0_this > T0_ug_max)
765 std::cout <<
"imu=" << imu <<
" T0_this=" << T0_this
767 for (
unsigned int jmu = 0; jmu < imu; ++jmu) {
768 if (std::fabs(
T0_ug[imu] -
T0_ug[jmu]) < minDeltaT0)
769 minDeltaT0 = std::fabs(
T0_ug[imu] -
T0_ug[jmu]);
782 double T0_offset, T0diff;
783 for (
int tboxtrial = 0; tboxtrial < 1000; ++tboxtrial) {
784 T0_offset =
RanGen->flat() * (T0_max - T0_min);
786 T0diff = T0_offset - T0_max;
788 for (
unsigned int imu = 0; imu <
Id_ug.size(); ++imu) {
799 std::cout <<
"initial T0_at=" <<
T0_at <<
" T0_min=" << T0_min <<
" T0_max=" << T0_max
800 <<
" T0_offset=" << T0_offset;
803 std::cout <<
" T0diff=" << T0diff << std::endl;
804 for (
unsigned int imu = 0; imu <
Id_ug.size(); ++imu) {
807 T0_sf[imu] += T0diff;
808 T0_ug[imu] += T0diff;
810 std::cout <<
" after: T0_sf[" << imu <<
"]=" <<
T0_sf[imu] <<
" T0_ug=" <<
T0_ug[imu] << std::endl;
813 std::cout <<
"T0diff=" << T0diff <<
" T0_at=" <<
T0_at << std::endl;
816 EventWeight = JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans / (trials * TboxTrials);
819 std::cout <<
"CosmicMuonGenerator.cc: Theta_at=" <<
Theta_at <<
" phi_at=" << phi_at <<
" Px_at=" <<
Px_at 822 <<
" Nmuons=" <<
Id_sf.size() << std::endl;
834 std::cout <<
"*********************************************************" << std::endl;
835 std::cout <<
"*********************************************************" << std::endl;
837 std::cout <<
"*** C O S M I C M U O N S T A T I S T I C S ***" << std::endl;
839 std::cout <<
"*********************************************************" << std::endl;
840 std::cout <<
"*********************************************************" << std::endl;
842 std::cout <<
" number of initial cosmic events: " <<
int(
Ngen) << std::endl;
844 std::cout <<
" number of generated and accepted events: " <<
int(
Nsel) << std::endl;
846 std::cout <<
" event selection efficiency: " << selEff * 100. <<
"%" << std::endl;
849 std::cout <<
" events with ~100 GeV and 1 - cos(theta) < 1/2pi: " << n100cos << std::endl;
851 std::cout <<
" momentum range: " <<
MinP <<
" ... " <<
MaxP <<
" GeV" << std::endl;
860 std::cout <<
" area of initial cosmics at CMS detector bottom surface: " <<
area <<
" m^2" << std::endl;
862 std::cout <<
" area of initial cosmics on Surface + PlugWidth: " <<
area <<
" m^2" << std::endl;
903 std::cout <<
" !!! MinP > 100 GeV. Cannot apply normalisation!" << std::endl;
905 std::cout <<
" !!! Note: generated cosmics exceed parameterisation. No flux calculated!" << std::endl;
908 std::cout <<
" !!! Not enough statistics to apply normalisation (rate=1 +- 1) !!!" << std::endl;
911 std::cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
914 <<
" muons per second" << std::endl;
918 std::cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
920 std::cout <<
"*********************************************************" << std::endl;
921 std::cout <<
"*********************************************************" << std::endl;
928 std::cout <<
" CMG-ERR: min.energy is out of range (0 GeV ... inf]" << std::endl << std::endl;
932 std::cout <<
" CMG-ERR: max.energy is out of range (0 GeV ... inf]" << std::endl << std::endl;
936 std::cout <<
" CMG-ERR: max.energy is not greater than min.energy" << std::endl << std::endl;
940 std::cout <<
" CMG-ERR: min.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl;
944 std::cout <<
" CMG-ERR: max.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl;
948 std::cout <<
" CMG-ERR: max.theta is not greater than min.theta" << std::endl << std::endl;
952 std::cout <<
" CMG-ERR: min.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl;
956 std::cout <<
" CMG-ERR: max.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl;
960 std::cout <<
" CMG-ERR: max.phi is not greater than min.phi" << std::endl << std::endl;
964 std::cout <<
" CMG-ERR: max.t0 is not greater than min.t0" << std::endl << std::endl;
968 std::cout <<
" CMG-ERR: E-loss scale factor is out of range [0 ... inf)" << std::endl << std::endl;
972 std::cout <<
" CMG-ERR: min.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl;
976 std::cout <<
" CMG-ERR: max.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl;
980 std::cout <<
" CMG-ERR: max.Enu is not greater than min.Enu" << std::endl << std::endl;
986 bool goodAngles =
false;
987 bool phiaccepted =
false;
988 bool thetaaccepted =
false;
1001 double disPhi = std::fabs(PhiV -
Phi);
1003 disPhi =
TwoPi - disPhi;
1010 double ThetaV = asin(RxzV / rVY);
1019 std::cout <<
"Rejected phi=" <<
Phi <<
" PhiV=" << PhiV <<
" dPhi=" <<
dPhi <<
" disPhi=" << disPhi
1020 <<
" RxzV=" << RxzV <<
" Target3dRadius=" <<
Target3dRadius <<
" Theta=" << Theta << std::endl;
1022 if (std::fabs(Theta - ThetaV) < dTheta)
1023 thetaaccepted =
true;
1024 if (phiaccepted && thetaaccepted)
1030 #if ROOT_INTERACTIVE 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);
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]");
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]");
1073 #if ROOT_INTERACTIVE 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);
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));
1105 disXY->Draw(
"COLZ");
1106 InteractionPoint->Draw(
"SAME");
1108 logEaxis->DrawLatex((0.65 * RadiusDet / 1000.), (1.08 * RadiusDet / 1000.),
"log_{10}E(#mu^{#pm})");
1111 InteractionPoint->Draw(
"SAME");
Int_t shower_nParticlesWritten
void setZDistOfTarget(double Z)
void propagate(double ElossScaleFac, double RadiusTarget, double Z_DistTarget, double Z_CentrTarget, bool TrackerOnly, bool MTCCHalf)
std::vector< double > Theta_mu
int events_n100cos(double energy, double theta)
void setMinEnu(double MinEn)
void setTIFOnly_constant(bool TIF)
std::vector< double > Px_ug
void setNuProdAlt(double NuPrdAlt)
void setZCentrOfTarget(double Z)
Double_t particle__Py[kMaxparticle_]
std::vector< double > Py_mu
void setRhoAir(double VarRhoAir)
const double Z_DistTracker
Sin< T >::type sin(const T &t)
int initialize(double, double, double, double, CLHEP::HepRandomEngine *, bool, bool)
void setRandomEngine(CLHEP::HepRandomEngine *v)
void setRadiusOfTarget(double R)
Double_t particle__Time[kMaxparticle_]
void setNumberOfEvents(unsigned int N)
std::vector< double > Py_sf
virtual void Init(TTree *tree)
std::vector< double > E_sf
Double_t particle__Px[kMaxparticle_]
std::vector< double > Pz_mu
const double SpeedOfLight
void initialize(CLHEP::HepRandomEngine *rng=nullptr)
Double_t particle__Pz[kMaxparticle_]
std::vector< double > E_ug
std::string MultiMuonFileName
void setMultiMuonFileFirstEvent(int MultiMuFile1stEvt)
const double NorthCMSzDeltaPhi
virtual Int_t GetEntry(Long64_t entry)
void setRhoPlug(double VarRhoPlug)
std::vector< double > Vz_ug
void setMinPhi(double Phi)
const double SurfaceOfEarth
void setMaxPhi(double Phi)
void setRandomEngine(CLHEP::HepRandomEngine *v)
std::vector< double > Vz_sf
std::vector< double > Vz_mu
std::vector< double > Vx_sf
void setMultiMuonNmin(int MultiMuNmin)
Cos< T >::type cos(const T &t)
int NcloseMultiMuonEvents
std::vector< double > Vy_sf
void setMinTheta(double Theta)
CLHEP::HepRandomEngine * RanGen
Tan< T >::type tan(const T &t)
std::vector< double > Vx_mu
void setMaxEnu(double MaxEn)
void setMultiMuon(bool MultiMu)
void setClayWidth(double ClayLaeyrWidth)
SingleParticleEvent OneMuoEvt
std::vector< double > Vy_mu
std::vector< double > T0_sf
std::vector< double > Px_mu
int MultiMuonFileFirstEvent
std::vector< double > Px_sf
std::vector< double > Vy_ug
void setAcptAllMu(bool AllMu)
Int_t particle__ParticleID[kMaxparticle_]
double momentum_times_charge()
void setPlugVz(double PlugVtz)
void setTIFOnly_linear(bool TIF)
std::pair< OmniClusterRef, TrackingParticleRef > P
std::vector< double > Pz_ug
Double_t particle__x[kMaxparticle_]
void create(int id, double px, double py, double pz, double e, double m, double vx, double vy, double vz, double t0)
void setMultiMuonFileName(std::string MultiMuonFileName)
void setMTCCHalf(bool MTCC)
const double RadiusTracker
int NskippedMultiMuonEvents
void setMinP_CMS(double P)
void setElossScaleFactor(double ElossScaleFact)
void setMaxTheta(double Theta)
unsigned int NumberOfEvents
ALPAKA_FN_ACC ALPAKA_FN_INLINE uint32_t iy(uint32_t id)
void setRhoWall(double VarRhoSWall)
void setPlugVx(double PlugVtx)
void setRhoRock(double VarRhoRock)
int initializeNuMu(double, double, double, double, double, double, double, double, double, CLHEP::HepRandomEngine *)
std::vector< double > Vx_ug
std::vector< double > P_mu
std::vector< double > Py_ug
std::vector< double > T0_ug
Power< A, B >::type pow(const A &a, const B &b)
void setRhoClay(double VarRhoClay)
std::vector< double > Pz_sf
Double_t particle__y[kMaxparticle_]
void setTrackerOnly(bool Tracker)