22 RanGen =
new CLHEP::HepJamesRandom;
64 TH1D* ene =
new TH1D(
"ene",
"generated energy",210,0.,1050.);
65 TH1D* the =
new TH1D(
"the",
"generated theta",90,0.,90.);
66 TH1D*
phi =
new TH1D(
"phi",
"generated phi",120,0.,360.);
67 TH3F* ver =
new TH3F(
"ver",
"Z-X-Y coordinates",100,-25.,25.,20,-10.,10.,40,-10.,10.);
104 double E = 0.;
double Theta = 0.;
double Phi = 0.;
double RxzV = 0.;
double PhiV = 0.;
105 if (
int(
Nsel)%100 == 0)
std::cout <<
" generated " << int(
Nsel) <<
" events" << std::endl;
107 bool notSelected =
true;
109 bool badMomentumGenerated =
true;
110 while (badMomentumGenerated){
120 badMomentumGenerated =
false;
128 bool badVertexGenerated =
true;
129 while (badVertexGenerated){
134 double rotPhi = PhiV +
Pi;
if (rotPhi >
TwoPi) rotPhi -=
TwoPi;
135 double disPhi = std::fabs(rotPhi - Phi);
if (disPhi > Pi) disPhi =
TwoPi - disPhi;
136 if (disPhi < dPhi ||
AcptAllMu) badVertexGenerated =
false;
145 double verMom = absMom*
cos(Theta);
146 double horMom = absMom*
sin(Theta);
147 double Px = horMom*
sin(Phi);
149 double Pz = horMom*
cos(Phi);
150 double Vx = RxzV*
sin(PhiV);
158 double Vz = RxzV*
cos(PhiV);
245 if (
Debug)
std::cout <<
"\nEntered CosmicMuonGenerator::nextMultiEvent()" << std::endl;
246 bool EvtRejected =
true;
247 bool MuInMaxDist =
false;
250 while (EvtRejected) {
258 if (ientry < 0)
return false;
263 std::cout <<
"CosmicMuonGenerator.cc::nextMultiEvent: No more events in file!" << std::endl;
272 " muons in event!" << std::endl;
273 std::cout <<
"trying next event from file" << std::endl;
288 for (
int imu=0; imu<nmuons; ++imu) {
302 for (
int jmu=0; jmu<imu; ++jmu) {
311 if (MuMuDist < MinDist) MinDist = MuMuDist;
320 std::cout <<
"CosmicMuonGenerator.cc: Warning! No muon pair closer than "
321 << 2.*
Target3dRadius/1000. <<
"m MinDist=" << MinDist/1000. <<
"m at surface" << std::endl;
322 std::cout <<
"Fraction of too wide opening angle multi muon events: "
325 std::cout <<
"trying next event from file" << std::endl;
338 std::cout <<
"start trial do loop: MuMuDist=" << MinDist/1000. <<
"[m] Nmuons="
355 if (phi_at < -
Pi) phi_at +=
TwoPi;
356 else if (phi_at >
Pi) phi_at -=
TwoPi;
368 double phi_rel_min = 0.;
369 double phi_rel_max = 0.;
374 for (
int imu=0; imu<nmuons; ++imu) {
380 double phi_rel = phi_mu - phi_at;
381 if (phi_rel < -
Pi) phi_rel +=
TwoPi;
382 else if (phi_rel >
Pi) phi_rel -=
TwoPi;
383 if (phi_rel < phi_rel_min) phi_rel_min = phi_rel;
384 else if (phi_rel > phi_rel_max) phi_rel_max =phi_rel;
394 double JdRxzV_dR_trans = 1.;
395 double JdPhiV_dPhi_trans = 1.;
396 double JdR_trans_sqrt = 1.;
402 double R_min =
std::max(0., R_mu_min);
405 std::cout <<
"CosmicMuonGenerator.cc: Warning! R_at=" << R_at
413 double psR1min = R_min + 0.25*(R_max-R_min);
415 double psR1 = psR1max - psR1min;
417 double psR2min = R_min;
418 double psR2max = R_max;
419 double psR2 = psR2max - psR2min;
423 double psR3 = psR3max - psR3min;
426 double psRall = psR3;
428 double fR1=psR1/psRall, fR2=psR2/psRall, fR3=psR3/psRall;
429 double pR1=0.25, pR2=0.7, pR3=0.05;
433 double psPh1 = 0.5*(phi_rel_max - phi_rel_min);
434 double psPh2 = phi_rel_max - phi_rel_min;
435 double psPh3 =
TwoPi;
436 double psPhall = psPh3;
438 double fPh1=psPh1/psPhall, fPh2=psPh2/psPhall, fPh3=psPh3/psPhall;
440 double pPh1=0.25, pPh2=0.7, pPh3=0.05;
448 int NmuHitTarget = 0;
457 double which_R_class =
RanGen->flat();
458 if (which_R_class < pR1) {
459 RxzV = psR1min + psR1 *
RanGen->flat();
462 else if (which_R_class < pR1+pR2) {
463 RxzV = psR2min + psR2 *
RanGen->flat();
467 RxzV = psR3min + psR3 *
RanGen->flat();
475 double which_phi_class =
RanGen->flat();
476 if (which_phi_class < pPh1) {
477 PhiV = phi_at + phi_rel_min + psPh1 *
RanGen->flat();
478 JdPhiV_dPhi_trans = fPh1/pPh1 *
TwoPi/psPh1;
480 else if (which_phi_class < pPh1+pPh2) {
481 PhiV = phi_at + phi_rel_min + psPh2 *
RanGen->flat();
482 JdPhiV_dPhi_trans = fPh2/pPh2 *
TwoPi/psPh2;
485 PhiV = phi_at + phi_rel_min + psPh3 *
RanGen->flat();
486 JdPhiV_dPhi_trans = fPh3/pPh3 *
TwoPi/psPh3;
491 else if (PhiV >
Pi) PhiV-=
TwoPi;
495 trials += JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans;
508 int NmuHitTargetSphere = 0;
509 for (
int imu=0; imu<nmuons; ++imu) {
520 double theta_sec = atan2(std::fabs(
Py_mu[imu]),pt_sec);
523 double h_prod = r_sec *
tan(theta_sec);
524 if (h_prod + h_sf >
Vy_at)
Vy_at = h_prod + h_sf;
542 double rotPhi = PhiVmu +
Pi;
if (rotPhi > Pi) rotPhi -=
TwoPi;
543 double disPhi = std::fabs(rotPhi - PhiPmu);
if (disPhi > Pi) disPhi =
TwoPi - disPhi;
545 NmuHitTargetSphere++;
553 if (temp > 0) btemp =
true;
583 double Px_sf_this =0., Py_sf_this=0., Pz_sf_this=0.;
586 double Vx_sf_this=0., Vy_sf_this=0., Vz_sf_this=0.;
587 double T0_sf_this=0.;
591 for (
int imu=0; imu<nmuons; ++imu) {
599 if (Id_sf_this == 5) Id_sf_this = -13;
600 else if (Id_sf_this == 6) Id_sf_this = 13;
602 else if (Id_sf_this == 1) Id_sf_this = 22;
603 else if (Id_sf_this == 2) Id_sf_this = -11;
604 else if (Id_sf_this == 3) Id_sf_this = 11;
605 else if (Id_sf_this == 7) Id_sf_this = 111;
606 else if (Id_sf_this == 8) Id_sf_this = 211;
607 else if (Id_sf_this == 9) Id_sf_this = -211;
608 else if (Id_sf_this == 10) Id_sf_this = 130;
609 else if (Id_sf_this == 11) Id_sf_this = 321;
610 else if (Id_sf_this == 12) Id_sf_this = -321;
611 else if (Id_sf_this == 13) Id_sf_this = 2112;
612 else if (Id_sf_this == 14) Id_sf_this = 2212;
613 else if (Id_sf_this == 15) Id_sf_this = -2212;
614 else if (Id_sf_this == 16) Id_sf_this = 310;
615 else if (Id_sf_this == 17) Id_sf_this = 221;
616 else if (Id_sf_this == 18) Id_sf_this = 3122;
619 std::cout <<
"CosmicMuonGenerator.cc: Warning! Muon Id=" << Id_sf_this
620 <<
" from file read in" <<std::endl;
626 Px_sf_this =
Px_mu[imu];
627 Py_sf_this =
Py_mu[imu];
628 Pz_sf_this =
Pz_mu[imu];
630 Vx_sf_this =
Vx_mu[imu];
631 Vy_sf_this =
Vy_mu[imu];
632 Vz_sf_this =
Vz_mu[imu];
635 OneMuoEvt.
create(Id_sf_this, Px_sf_this, Py_sf_this, Pz_sf_this, E_sf_this, MuonMass, Vx_sf_this, Vy_sf_this, Vz_sf_this, T0_sf_this);
645 Id_sf.push_back(Id_sf_this);
646 Px_sf.push_back(Px_sf_this);
647 Py_sf.push_back(Py_sf_this);
648 Pz_sf.push_back(Pz_sf_this);
649 E_sf.push_back(E_sf_this);
651 Vx_sf.push_back(Vx_sf_this);
652 Vy_sf.push_back(Vy_sf_this);
653 Vz_sf.push_back(Vz_sf_this);
654 T0_sf.push_back(T0_sf_this);
678 <<
" without accepting event!" << std::endl;
695 " muons hit target: N(mu=)" << NmuHitTarget << std::endl;
704 double T0_ug_min, T0_ug_max;
705 T0_ug_min = T0_ug_max =
T0_ug[0];
707 double minDeltaT0 = 2*Tbox;
708 for (
unsigned int imu=0; imu<
Id_ug.size(); ++imu) {
709 double T0_this =
T0_ug[imu];
710 if (T0_this < T0_ug_min) T0_ug_min = T0_this;
711 if (T0_this > T0_ug_max) T0_ug_max = T0_this;
715 for (
unsigned int jmu=0; jmu<imu; ++jmu) {
716 if (std::fabs(
T0_ug[imu]-
T0_ug[jmu]) < minDeltaT0) minDeltaT0 = std::fabs(
T0_ug[imu]-
T0_ug[jmu]);
729 double T0_offset, T0diff;
730 for (
int tboxtrial=0; tboxtrial<1000; ++tboxtrial) {
731 T0_offset =
RanGen->flat()*(T0_max -T0_min);
733 T0diff = T0_offset - T0_max;
735 for (
unsigned int imu=0; imu<
Id_ug.size(); ++imu) {
745 if (
Debug)
std::cout <<
"initial T0_at=" <<
T0_at <<
" T0_min=" << T0_min <<
" T0_max=" << T0_max
746 <<
" T0_offset=" << T0_offset;
749 for (
unsigned int imu=0; imu<
Id_ug.size(); ++imu) {
751 T0_sf[imu] += T0diff;
752 T0_ug[imu] += T0diff;
754 std::cout <<
" after: T0_sf[" << imu <<
"]=" <<
T0_sf[imu] <<
" T0_ug=" << T0_ug[imu] << std::endl;
761 EventWeight = JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans
762 / (trials * TboxTrials);
766 <<
" T0_at=" <<
T0_at
768 <<
" EventWeight=" <<
EventWeight <<
" Nmuons=" <<
Id_sf.size() << std::endl;
785 std::cout <<
"*********************************************************" << std::endl;
786 std::cout <<
"*********************************************************" << std::endl;
788 std::cout <<
"*** C O S M I C M U O N S T A T I S T I C S ***" << std::endl;
790 std::cout <<
"*********************************************************" << std::endl;
791 std::cout <<
"*********************************************************" << std::endl;
793 std::cout <<
" number of initial cosmic events: " << int(
Ngen) << std::endl;
794 std::cout <<
" number of actually diced events: " << int(
Ndiced) << std::endl;
795 std::cout <<
" number of generated and accepted events: " << int(
Nsel) << std::endl;
797 std::cout <<
" event selection efficiency: " << selEff*100. <<
"%" << std::endl;
799 std::cout <<
" events with ~100 GeV and 1 - cos(theta) < 1/2pi: " << n100cos << std::endl;
801 std::cout <<
" momentum range: " <<
MinP <<
" ... " <<
MaxP <<
" GeV" << std::endl;
809 std::cout <<
" area of initial cosmics at CMS detector bottom surface: " << area <<
" m^2" << std::endl;
811 std::cout <<
" area of initial cosmics on Surface + PlugWidth: " << area <<
" m^2" << std::endl;
853 std::cout <<
" !!! MinP > 100 GeV. Cannot apply normalisation!" << std::endl;
855 std::cout <<
" !!! Note: generated cosmics exceed parameterisation. No flux calculated!" << std::endl;
858 std::cout <<
" !!! Not enough statistics to apply normalisation (rate=1 +- 1) !!!" << std::endl;
861 std::cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
863 rateErr_syst <<
" (syst) " <<
" muons per second" << std::endl;
865 std::cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
867 std::cout <<
"*********************************************************" << std::endl;
868 std::cout <<
"*********************************************************" << std::endl;
874 std::cout <<
" CMG-ERR: min.energy is out of range (0 GeV ... inf]" << std::endl << std::endl; }
876 std::cout <<
" CMG-ERR: max.energy is out of range (0 GeV ... inf]" << std::endl << std::endl; }
878 std::cout <<
" CMG-ERR: max.energy is not greater than min.energy" << std::endl << std::endl; }
880 std::cout <<
" CMG-ERR: min.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl; }
882 std::cout <<
" CMG-ERR: max.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl; }
884 std::cout <<
" CMG-ERR: max.theta is not greater than min.theta" << std::endl << std::endl; }
886 std::cout <<
" CMG-ERR: min.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl; }
888 std::cout <<
" CMG-ERR: max.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl; }
890 std::cout <<
" CMG-ERR: max.phi is not greater than min.phi" << std::endl << std::endl; }
892 std::cout <<
" CMG-ERR: max.t0 is not greater than min.t0" << std::endl << std::endl; }
894 std::cout <<
" CMG-ERR: E-loss scale factor is out of range [0 ... inf)" << std::endl << std::endl; }
896 std::cout <<
" CMG-ERR: min.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl; }
898 std::cout <<
" CMG-ERR: max.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl; }
900 std::cout <<
" CMG-ERR: max.Enu is not greater than min.Enu" << std::endl << std::endl; }
906 bool goodAngles =
false;
907 bool phiaccepted =
false;
908 bool thetaaccepted =
false;
919 double disPhi = std::fabs(PhiV - Phi);
if (disPhi >
Pi) disPhi =
TwoPi - disPhi;
921 if (disPhi < dPhi) phiaccepted =
true;
923 double ThetaV = asin(RxzV/rVY);
930 std::cout <<
"Rejected phi=" << Phi <<
" PhiV=" << PhiV
931 <<
" dPhi=" << dPhi <<
" disPhi=" << disPhi
933 <<
" Theta=" << Theta << std::endl;
935 if (std::fabs(Theta-ThetaV) < dTheta) thetaaccepted =
true;
936 if (phiaccepted && thetaaccepted) goodAngles =
true;
948 TH2F* disXY =
new TH2F(
"disXY",
"X-Y view",160,-rCMS,rCMS,160,-rCMS,rCMS);
949 TH2F* disZY =
new TH2F(
"disZY",
"Z-Y view",150,-zCMS,zCMS,160,-rCMS,rCMS);
950 gStyle->SetPalette(1,0);
951 gStyle->SetMarkerColor(1);
952 gStyle->SetMarkerSize(1.5);
953 TCanvas *disC =
new TCanvas(
"disC",
"Cosmic Muon Event Display",0,0,800,410);
957 disXY->SetMinimum(log10(
MinP));
958 disXY->SetMaximum(log10(
MaxP));
959 disXY->GetXaxis()->SetLabelSize(0.05);
960 disXY->GetXaxis()->SetTitleSize(0.05);
961 disXY->GetXaxis()->SetTitleOffset(1.0);
962 disXY->GetXaxis()->SetTitle(
"X [m]");
963 disXY->GetYaxis()->SetLabelSize(0.05);
964 disXY->GetYaxis()->SetTitleSize(0.05);
965 disXY->GetYaxis()->SetTitleOffset(0.8);
966 disXY->GetYaxis()->SetTitle(
"Y [m]");
970 disZY->SetMinimum(log10(
MinP));
971 disZY->SetMaximum(log10(
MaxP));
972 disZY->GetXaxis()->SetLabelSize(0.05);
973 disZY->GetXaxis()->SetTitleSize(0.05);
974 disZY->GetXaxis()->SetTitleOffset(1.0);
975 disZY->GetXaxis()->SetTitle(
"Z [m]");
976 disZY->GetYaxis()->SetLabelSize(0.05);
977 disZY->GetYaxis()->SetTitleSize(0.05);
978 disZY->GetYaxis()->SetTitleOffset(0.8);
979 disZY->GetYaxis()->SetTitle(
"Y [m]");
993 TMarker* InteractionPoint =
new TMarker(0.,0.,2);
994 TArc* r8m =
new TArc(0.,0.,(RadiusDet/1000.));
995 TLatex* logEaxis =
new TLatex(); logEaxis->SetTextSize(0.05);
1003 float yStep = disXY->GetYaxis()->GetBinWidth(1);
1004 int NbinY = disXY->GetYaxis()->GetNbins();
1005 for (
int iy=0; iy<NbinY; ++iy){
1009 float rXY =
sqrt(verX*verX + verY*verY)*1000.;
1010 float absZ = std::fabs(verZ)*1000.;
1011 if (rXY < RadiusDet && absZ < Z_DistDet){
1012 disXY->Fill(verX,verY,log10(energy));
1013 disZY->Fill(verZ,verY,log10(energy));
1014 disC->cd(1); disXY->Draw(
"COLZ"); InteractionPoint->Draw(
"SAME"); r8m->Draw(
"SAME");
1015 logEaxis->DrawLatex((0.65*RadiusDet/1000.),(1.08*RadiusDet/1000.),
"log_{10}E(#mu^{#pm})");
1016 disC->cd(2); disZY->Draw(
"COL"); InteractionPoint->Draw(
"SAME");
const double Z[kNumberCalorimeter]
Int_t shower_nParticlesWritten
void setZDistOfTarget(double Z)
void propagate(double ElossScaleFac, double RadiusTarget, double Z_DistTarget, double Z_CentrTarget, bool TrackerOnly, bool MTCCHalf)
void initialize(CLHEP::HepRandomEngine *rng=0)
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 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
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)
std::vector< double > Vz_sf
std::vector< double > Vz_mu
double dPhi(double phi1, double phi2)
const T & max(const T &a, const T &b)
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::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
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)