CMS 3D CMS Logo

CosmicMuonGenerator.cc
Go to the documentation of this file.
1 // modified by P. Biallass 29.03.2006 to implement new cosmic generator (CMSCGEN.cc) and new normalization of flux (CMSCGENnorm.cc)
4 // 04.12.2008 sonne: replaced Min/MaxE by Min/MaxP to get cos_sf/ug scripts working again
5 // 20.04.2009 sonne: Implemented mechanism to read in multi muon events and propagate each muon
6 #define sim_cxx
7 
9 
11  initialize();
12  for (unsigned int iGen = 0; iGen < NumberOfEvents; ++iGen) {
13  nextEvent();
14  }
15  terminate();
16 }
17 
18 void CosmicMuonGenerator::setRandomEngine(CLHEP::HepRandomEngine* v) {
19  if (delRanGen)
20  delete RanGen;
21  RanGen = v;
22  delRanGen = false;
24 }
25 
26 void CosmicMuonGenerator::initialize(CLHEP::HepRandomEngine* rng) {
27  if (delRanGen)
28  delete RanGen;
29  if (!rng) {
30  RanGen = new CLHEP::HepJamesRandom;
31  RanGen->setSeed(RanSeed, 0); //set seed for Random Generator (seed can be controled by config-file)
32  delRanGen = true;
33  } else {
34  RanGen = rng;
35  delRanGen = false;
36  }
37  checkIn();
38  if (NumberOfEvents > 0) {
39  // set up "surface geometry" dimensions
40  double RadiusTargetEff = RadiusOfTarget; //get this from cfg-file
41  double Z_DistTargetEff = ZDistOfTarget; //get this from cfg-file
42  //double Z_CentrTargetEff = ZCentrOfTarget; //get this from cfg-file
43  if (TrackerOnly == true) {
44  RadiusTargetEff = RadiusTracker;
45  Z_DistTargetEff = Z_DistTracker;
46  }
47  Target3dRadius = sqrt(RadiusTargetEff * RadiusTargetEff + Z_DistTargetEff * Z_DistTargetEff) + MinStepSize;
48  if (Debug)
49  std::cout << " radius of sphere around target = " << Target3dRadius << " mm" << std::endl;
50 
51  if (MinTheta > 90. * Deg2Rad) //upgoing muons from neutrinos
53  else
54  SurfaceRadius = (SurfaceOfEarth + PlugWidth + RadiusTargetEff) * tan(MaxTheta) + Target3dRadius;
55  if (Debug)
56  std::cout << " starting point radius at Surface + PlugWidth = " << SurfaceRadius << " mm" << std::endl;
57 
65  //set energy and angle limits for CMSCGEN, give same seed as above
66  if (MinTheta >= 90. * Deg2Rad) //upgoing muons from neutrinos
68  else
70 
71 #if ROOT_INTERACTIVE
72  // book histos
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.);
77 #endif
78  if (EventDisplay)
79  initEvDis();
80  std::cout << std::endl;
81 
82  if (MultiMuon) {
83  MultiIn = nullptr;
84 
85  std::cout << "MultiMuonFileName.c_str()=" << MultiMuonFileName.c_str() << std::endl;
86  MultiIn = new TFile(MultiMuonFileName.c_str());
87 
88  if (!MultiIn)
89  std::cout << "MultiMuon=True: MultiMuonFileName='" << MultiMuonFileName.c_str() << "' does not exist"
90  << std::endl;
91  else
92  std::cout << "MultiMuonFile: " << MultiMuonFileName.c_str() << " opened!" << std::endl;
93  //MultiTree = (TTree*) gDirectory->Get("sim");
94  MultiTree = (TTree*)MultiIn->Get("sim");
95  SimTree = new sim(MultiTree);
97  SimTreeEntries = SimTree->fChain->GetEntriesFast();
98  std::cout << "SimTreeEntries=" << SimTreeEntries << std::endl;
99 
100  if (MultiMuonFileFirstEvent <= 0)
101  SimTree_jentry = 0;
102  else
103  SimTree_jentry = MultiMuonFileFirstEvent - 1; //1=1st evt (SimTree_jentry=0)
104 
107  }
108 
109  if (!MultiMuon || (MultiMuon && MultiIn))
110  NotInitialized = false;
111  }
112 }
113 
115  double E = 0.;
116  double Theta = 0.;
117  double Phi = 0.;
118  double RxzV = 0.;
119  double PhiV = 0.;
120  if (int(Nsel) % 100 == 0)
121  std::cout << " generated " << int(Nsel) << " events" << std::endl;
122  // generate cosmic (E,theta,phi)
123  bool notSelected = true;
124  while (notSelected) {
125  bool badMomentumGenerated = true;
126  while (badMomentumGenerated) {
127  if (MinTheta > 90. * Deg2Rad) //upgoing muons from neutrinos
129  else
130  Cosmics->generate(); //dice one event now
131 
133  Theta = TMath::ACos(Cosmics->cos_theta()); //angle has to be in RAD here
134  Ngen += 1.; //count number of initial cosmic events (in surface area), vertices will be added later
135  badMomentumGenerated = false;
136  Phi = RanGen->flat() * (MaxPhi - MinPhi) + MinPhi;
137  }
138  Norm->events_n100cos(E, Theta); //test if this muon is in normalization range
139  Ndiced += 1; //one more cosmic is diced
140 
141  // generate vertex
142  double Nver = 0.;
143  bool badVertexGenerated = true;
144  while (badVertexGenerated) {
145  RxzV = sqrt(RanGen->flat()) * SurfaceRadius;
146  PhiV = RanGen->flat() * TwoPi;
147  // check phi range (for a sphere with Target3dRadius around the target)
148  double dPhi = Pi;
149  if (RxzV > Target3dRadius)
150  dPhi = asin(Target3dRadius / RxzV);
151  double rotPhi = PhiV + Pi;
152  if (rotPhi > TwoPi)
153  rotPhi -= TwoPi;
154  double disPhi = std::fabs(rotPhi - Phi);
155  if (disPhi > Pi)
156  disPhi = TwoPi - disPhi;
157  if (disPhi < dPhi || AcptAllMu)
158  badVertexGenerated = false;
159  Nver += 1.;
160  }
161  Ngen += (Nver - 1.); //add number of generated vertices to initial cosmic events
162 
163  // complete event at surface
164  int id = 13; // mu-
165  if (Cosmics->momentum_times_charge() > 0.)
166  id = -13; // mu+
167  double absMom = sqrt(E * E - MuonMass * MuonMass);
168  double verMom = absMom * cos(Theta);
169  double horMom = absMom * sin(Theta);
170  double Px = horMom * sin(Phi); // [GeV/c]
171  double Py = -verMom; // [GeV/c]
172  double Pz = horMom * cos(Phi); // [GeV/c]
173  double Vx = RxzV * sin(PhiV); // [mm]
174 
175  double Vy;
176  if (MinTheta > 90. * Deg2Rad) //upgoing muons from neutrinos
177  Vy = -RadiusCMS;
178  else
179  Vy = SurfaceOfEarth + PlugWidth; // [mm]
180 
181  double Vz = RxzV * cos(PhiV); // [mm]
182  double T0 = (RanGen->flat() * (MaxT0 - MinT0) + MinT0) * SpeedOfLight; // [mm/c];
183 
184  Id_at = id;
185  Px_at = Px;
186  Py_at = Py;
187  Pz_at = Pz;
188  E_at = E; //M_at = MuonMass;
189  Vx_at = Vx;
190  Vy_at = Vy;
191  Vz_at = Vz;
192  T0_at = T0;
193 
194  OneMuoEvt.create(id, Px, Py, Pz, E, MuonMass, Vx, Vy, Vz, T0);
195  // if angles are ok, propagate to target
196  if (goodOrientation()) {
197  if (MinTheta > 90. * Deg2Rad) //upgoing muons from neutrinos
199  else
201  }
202 
204  AcptAllMu == true) {
205  Nsel += 1.; //count number of generated and accepted events
206  notSelected = false;
207  }
208  }
209 
210  EventWeight = 1.;
211 
212  //just one outgoing particle at SurFace
213  Id_sf.resize(1);
214  Px_sf.resize(1);
215  Py_sf.resize(1);
216  Pz_sf.resize(1);
217  E_sf.resize(1);
218  //M_sf.resize(1);
219  Vx_sf.resize(1);
220  Vy_sf.resize(1);
221  Vz_sf.resize(1);
222  T0_sf.resize(1);
223 
224  Id_sf[0] = Id_at;
225  Px_sf[0] = Px_at;
226  Py_sf[0] = Py_at;
227  Pz_sf[0] = Pz_at;
228  E_sf[0] = E_at; //M_fs[0] = MuonMass;
229  Vx_sf[0] = Vx_at;
230  Vy_sf[0] = Vy_at;
231  Vz_sf[0] = Vz_at;
232  T0_sf[0] = T0_at;
233 
234  //just one particle at UnderGround
235  Id_ug.resize(1);
236  Px_ug.resize(1);
237  Py_ug.resize(1);
238  Pz_ug.resize(1);
239  E_ug.resize(1);
240  //M_ug.resize(1);
241  Vx_ug.resize(1);
242  Vy_ug.resize(1);
243  Vz_ug.resize(1);
244  T0_ug.resize(1);
245 
246  Id_ug[0] = OneMuoEvt.id();
247  Px_ug[0] = OneMuoEvt.px();
248  Py_ug[0] = OneMuoEvt.py();
249  Pz_ug[0] = OneMuoEvt.pz();
250  E_ug[0] = OneMuoEvt.e();
251  //M_ug[0] = OneMuoEvt.m();
252  Vx_ug[0] = OneMuoEvt.vx();
253  Vy_ug[0] = OneMuoEvt.vy();
254  Vz_ug[0] = OneMuoEvt.vz();
255  T0_ug[0] = OneMuoEvt.t0();
256 
257  // plot variables of selected events
258 #if ROOT_INTERACTIVE
259  ene->Fill(OneMuoEvt.e());
260  the->Fill((OneMuoEvt.theta() * Rad2Deg));
261  phi->Fill((OneMuoEvt.phi() * Rad2Deg));
262  ver->Fill((OneMuoEvt.vz() / 1000.), (OneMuoEvt.vx() / 1000.), (OneMuoEvt.vy() / 1000.));
263 #endif
264  if (Debug) {
265  std::cout << "new event" << std::endl;
266  std::cout << " Px,Py,Pz,E,m = " << OneMuoEvt.px() << ", " << OneMuoEvt.py() << ", " << OneMuoEvt.pz() << ", "
267  << OneMuoEvt.e() << ", " << OneMuoEvt.m() << " GeV" << std::endl;
268  std::cout << " Vx,Vy,Vz,t0 = " << OneMuoEvt.vx() << ", " << OneMuoEvt.vy() << ", " << OneMuoEvt.vz() << ", "
269  << OneMuoEvt.t0() << " mm" << std::endl;
270  }
271  if (EventDisplay)
272  displayEv();
273 }
274 
276  if (Debug)
277  std::cout << "\nEntered CosmicMuonGenerator::nextMultiEvent()" << std::endl;
278  bool EvtRejected = true;
279  bool MuInMaxDist = false;
280  double MinDist; //[mm]
281 
282  while (EvtRejected) {
283  //read in event from SimTree
284  //ULong64_t ientry = SimTree->LoadTree(SimTree_jentry);
285  Long64_t ientry = SimTree->GetEntry(SimTree_jentry);
286  std::cout << "CosmicMuonGenerator::nextMultiEvent(): SimTree_jentry="
287  << SimTree_jentry
288  //<< " ientry=" << ientry
289  << " SimTreeEntries=" << SimTreeEntries << std::endl;
290  if (ientry < 0)
291  return false; //stop run
293  SimTree_jentry++;
294  } else {
295  std::cout << "CosmicMuonGenerator.cc::nextMultiEvent: No more events in file!" << std::endl;
296  return false; //stop run
297  }
298 
300  if (nmuons < MultiMuonNmin) {
301  std::cout << "CosmicMuonGenerator.cc: Warning! Less than " << MultiMuonNmin << " muons in event!" << std::endl;
302  std::cout << "trying next event from file" << std::endl;
304  continue; //EvtRejected while loop: get next event from file
305  }
306 
307  Px_mu.resize(nmuons);
308  Py_mu.resize(nmuons);
309  Pz_mu.resize(nmuons);
310  P_mu.resize(nmuons);
311 
312  MinDist = 99999.e9; //[mm]
313  double MuMuDist;
314  MuInMaxDist = false;
315  //check if at least one muon pair closer than 30m at surface
316  int NmuPmin = 0;
317  for (int imu = 0; imu < nmuons; ++imu) {
318  Px_mu[imu] =
320  Pz_mu[imu] =
322  Py_mu[imu] = -SimTree->particle__Pz[imu]; //Corsika down going particles defined in -z direction!
323  P_mu[imu] = sqrt(Px_mu[imu] * Px_mu[imu] + Py_mu[imu] * Py_mu[imu] + Pz_mu[imu] * Pz_mu[imu]);
324 
325  if (P_mu[imu] < MinP_CMS && AcptAllMu == false)
326  continue;
327  else if (SimTree->particle__ParticleID[imu] != 5 && SimTree->particle__ParticleID[imu] != 6)
328  continue;
329  else
330  NmuPmin++;
331 
332  for (int jmu = 0; jmu < imu; ++jmu) {
333  if (P_mu[jmu] < MinP_CMS && AcptAllMu == false)
334  continue;
335  if (SimTree->particle__ParticleID[imu] != 5 && SimTree->particle__ParticleID[imu] != 6)
336  continue;
337  MuMuDist = sqrt((SimTree->particle__x[imu] - SimTree->particle__x[jmu]) *
338  (SimTree->particle__x[imu] - SimTree->particle__x[jmu]) +
339  (SimTree->particle__y[imu] - SimTree->particle__y[jmu]) *
340  (SimTree->particle__y[imu] - SimTree->particle__y[jmu])) *
341  10.; //CORSIKA [cm] to CMSCGEN [mm]
342  if (MuMuDist < MinDist)
343  MinDist = MuMuDist;
344  if (MuMuDist < 2. * Target3dRadius)
345  MuInMaxDist = true;
346  }
347  }
348  if (MultiMuonNmin >= 2) {
349  if (MuInMaxDist) {
351  } else {
352  std::cout << "CosmicMuonGenerator.cc: Warning! No muon pair closer than " << 2. * Target3dRadius / 1000.
353  << "m MinDist=" << MinDist / 1000. << "m at surface" << std::endl;
354  std::cout << "Fraction of too wide opening angle multi muon events: "
355  << 1 - double(NcloseMultiMuonEvents) / SimTree_jentry << std::endl;
356  std::cout << "NcloseMultiMuonEvents=" << NcloseMultiMuonEvents << std::endl;
357  std::cout << "trying next event from file" << std::endl;
359  continue; //EvtRejected while loop: get next event from file
360  }
361  }
362 
363  if (NmuPmin < MultiMuonNmin && AcptAllMu == false) { //take single muon events consistently into account
365  continue; //EvtRejected while loop: get next event from file
366  }
367 
368  if (Debug)
369  if (MultiMuonNmin >= 2)
370  std::cout << "start trial do loop: MuMuDist=" << MinDist / 1000. << "[m] Nmuons=" << nmuons
371  << " NcloseMultiMuonEvents=" << NcloseMultiMuonEvents
372  << " NskippedMultiMuonEvents=" << NskippedMultiMuonEvents << std::endl;
373 
374  //int primary_id = SimTree->run_ParticleID;
376 
377  double M_at = 0.;
378  //if (Id_at == 13) {
379  Id_at = 2212; //convert from Corsika to HepPDT
380  M_at = 938.272e-3; //[GeV] mass
381  //}
382 
385  double phi_at = SimTree->shower_Phi - NorthCMSzDeltaPhi; //rotate by almost 90 degrees
386  if (phi_at < -Pi)
387  phi_at += TwoPi; //bring into interval (-Pi,Pi]
388  else if (phi_at > Pi)
389  phi_at -= TwoPi;
390  double P_at = sqrt(E_at * E_at - M_at * M_at);
391  //need to rotate about 90degrees around x->N axis => y<=>z,
392  //then rotate new x-z-plane from x->North to x->LHC centre
393  Px_at = P_at * sin(Theta_at) * sin(phi_at);
394  Py_at = -P_at * cos(Theta_at);
395  Pz_at = P_at * sin(Theta_at) * cos(phi_at);
396 
397  //compute maximal theta of secondary muons
398  double theta_mu_max = Theta_at;
399  double theta_mu_min = Theta_at;
400 
401  double phi_rel_min = 0.; //phi_mu_min - phi_at
402  double phi_rel_max = 0.; //phi_mu_max - phi_at
403 
404  //std::cout << "SimTree->shower_Energy=" << SimTree->shower_Energy <<std::endl;
405 
406  Theta_mu.resize(nmuons);
407  for (int imu = 0; imu < nmuons; ++imu) {
408  Theta_mu[imu] = acos(-Py_mu[imu] / P_mu[imu]);
409  if (Theta_mu[imu] > theta_mu_max)
410  theta_mu_max = Theta_mu[imu];
411  if (Theta_mu[imu] < theta_mu_min)
412  theta_mu_min = Theta_mu[imu];
413 
414  double phi_mu = atan2(Px_mu[imu], Pz_mu[imu]); // in (-Pi,Pi]
415  double phi_rel = phi_mu - phi_at;
416  if (phi_rel < -Pi)
417  phi_rel += TwoPi; //bring into interval (-Pi,Pi]
418  else if (phi_rel > Pi)
419  phi_rel -= TwoPi;
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;
424  }
425 
426  double h_sf = SurfaceOfEarth + PlugWidth; //[mm]
427 
428  double R_at = h_sf * tan(Theta_at);
429 
430  double JdRxzV_dR_trans = 1.;
431  double JdPhiV_dPhi_trans = 1.;
432  double JdR_trans_sqrt = 1.;
433 
434  //chose random vertex Phi and Rxz weighted to speed up and smoothen
435  double R_mu_max = (h_sf + Target3dRadius) * tan(theta_mu_max);
436  double R_max = std::min(SurfaceRadius, R_mu_max);
437  double R_mu_min = (h_sf - Target3dRadius) * tan(theta_mu_min);
438  double R_min = std::max(0., R_mu_min);
439 
440  if (R_at > SurfaceRadius) {
441  std::cout << "CosmicMuonGenerator.cc: Warning! R_at=" << R_at << " > SurfaceRadius=" << SurfaceRadius
442  << std::endl;
443  }
444 
445  //do phase space transformation for horizontal radius R
446 
447  //determine first phase space limits
448 
449  double psR1min = R_min + 0.25 * (R_max - R_min);
450  double psR1max = std::min(SurfaceRadius, R_max - 0.25 * (R_max - R_min)); //no R's beyond R_max
451  double psR1 = psR1max - psR1min;
452 
453  double psR2min = R_min;
454  double psR2max = R_max;
455  double psR2 = psR2max - psR2min;
456 
457  double psR3min = 0.;
458  double psR3max = SurfaceRadius;
459  double psR3 = psR3max - psR3min;
460 
461  //double psall = psR1+psR2+psR3;
462  double psRall = psR3;
463 
464  double fR1 = psR1 / psRall, fR2 = psR2 / psRall, fR3 = psR3 / psRall; //f1+f2+f3=130%
465  double pR1 = 0.25, pR2 = 0.7, pR3 = 0.05;
466 
467  //do phase space transformation for azimuthal angle phi
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;
472 
473  double fPh1 = psPh1 / psPhall, fPh2 = psPh2 / psPhall,
474  fPh3 = psPh3 / psPhall; //(f1+f2+f3=TwoPi+f1+f2)/(TwoPi+f1+f2)
475 
476  double pPh1 = 0.25, pPh2 = 0.7, pPh3 = 0.05;
477 
478  Trials = 0; //global int trials
479  double trials = 0.; //local weighted trials
480  Vx_mu.resize(nmuons);
481  Vy_mu.resize(nmuons);
482  Vz_mu.resize(nmuons);
483  int NmuHitTarget = 0;
484  while (NmuHitTarget < MultiMuonNmin) {
485  NmuHitTarget = 0; //re-initialize every loop iteration
486  double Nver = 0.;
487 
488  //chose phase space class
489  double RxzV;
490  double which_R_class = RanGen->flat();
491  if (which_R_class < pR1) { //pR1% in psR1
492  RxzV = psR1min + psR1 * RanGen->flat();
493  JdRxzV_dR_trans = fR1 / pR1 * SurfaceRadius / psR1;
494  } else if (which_R_class < pR1 + pR2) { //further pR2% in psR2
495  RxzV = psR2min + psR2 * RanGen->flat();
496  JdRxzV_dR_trans = fR2 / pR2 * SurfaceRadius / psR2;
497  } else { //remaining pR3% in psR3=[0., R_max]
498  RxzV = psR3min + psR3 * RanGen->flat();
499  JdRxzV_dR_trans = fR3 / pR3 * SurfaceRadius / psR3;
500  }
501 
502  JdR_trans_sqrt = 2. * RxzV / SurfaceRadius; //flat in sqrt(r) space
503 
504  //chose phase space class
505  double PhiV;
506  double which_phi_class = RanGen->flat();
507  if (which_phi_class < pPh1) { //pPh1% in psPh1
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) { //further pPh2% in psPh2
511  PhiV = phi_at + phi_rel_min + psPh2 * RanGen->flat();
512  JdPhiV_dPhi_trans = fPh2 / pPh2 * TwoPi / psPh2;
513  } else { //remaining pPh3% in psPh3=[-Pi,Pi]
514  PhiV = phi_at + phi_rel_min + psPh3 * RanGen->flat();
515  JdPhiV_dPhi_trans = fPh3 / pPh3 * TwoPi / psPh3;
516  }
517 
518  //shuffle PhiV into [-Pi,+Pi] interval
519  if (PhiV < -Pi)
520  PhiV += TwoPi;
521  else if (PhiV > Pi)
522  PhiV -= TwoPi;
523 
524  Nver++;
525  trials += JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans;
526  Trials++;
527  if (trials > max_Trials)
528  break; //while (Id_sf.size() < 2) loop
529  Ngen += (Nver - 1.); //add number of generated vertices to initial cosmic events
530 
531  Vx_at = RxzV * sin(PhiV); // [mm]
532 
533  Vy_at = h_sf; // [mm] (SurfaceOfEarth + PlugWidth; Determine primary particle height below)
534  //Vy_at = SimTree->shower_StartingAltitude*10. + h_sf; // [mm]
535  //std::cout << "SimTree->shower_StartingAltitude*10=" << SimTree->shower_StartingAltitude*10 <<std::endl;
536  Vz_at = RxzV * cos(PhiV); // [mm]
537 
538  int NmuHitTargetSphere = 0;
539  for (int imu = 0; imu < nmuons; ++imu) {
540  Vx_mu[imu] = Vx_at + (-SimTree->particle__x[imu] * sin(NorthCMSzDeltaPhi) +
542  10; //[mm] (Corsika cm to CMSCGEN mm)
543  Vy_mu[imu] = h_sf; //[mm] fixed at surface + PlugWidth
544  Vz_mu[imu] = Vz_at + (SimTree->particle__x[imu] * cos(NorthCMSzDeltaPhi) +
546  10; //[mm] (Corsika cm to CMSCGEN mm)
547 
548  //add atmospheric height to primary particle (default SimTree->shower_StartingAltitude = 0.)
549  double pt_sec = sqrt(Px_mu[imu] * Px_mu[imu] + Pz_mu[imu] * Pz_mu[imu]);
550  double theta_sec = atan2(std::fabs(Py_mu[imu]), pt_sec);
551  double r_sec = sqrt((Vx_mu[imu] - Vx_at) * (Vx_mu[imu] - Vx_at) + (Vz_mu[imu] - Vz_at) * (Vz_mu[imu] - Vz_at));
552  double h_prod = r_sec * tan(theta_sec);
553  if (h_prod + h_sf > Vy_at)
554  Vy_at = h_prod + h_sf;
555 
556  //only muons
557  if (SimTree->particle__ParticleID[imu] != 5 && SimTree->particle__ParticleID[imu] != 6)
558  continue;
559 
560  if (P_mu[imu] < MinP_CMS && AcptAllMu == false)
561  continue;
562 
563  //check here if at least 2 muons make it to the target sphere
564  double Vxz_mu = sqrt(Vx_mu[imu] * Vx_mu[imu] + Vz_mu[imu] * Vz_mu[imu]);
565  theta_mu_max = atan((Vxz_mu + Target3dRadius) / (h_sf - Target3dRadius));
566  theta_mu_min = atan((Vxz_mu - Target3dRadius) / (h_sf + Target3dRadius));
567  if (Theta_mu[imu] > theta_mu_min && Theta_mu[imu] < theta_mu_max) {
568  // check phi range (for a sphere with Target3dRadius around the target)
569  double dPhi = Pi;
570  if (Vxz_mu > Target3dRadius)
571  dPhi = asin(Target3dRadius / Vxz_mu);
572  double PhiPmu = atan2(Px_mu[imu], Pz_mu[imu]); //muon phi
573  double PhiVmu = atan2(Vx_mu[imu], Vz_mu[imu]); //muon phi
574  double rotPhi = PhiVmu + Pi;
575  if (rotPhi > Pi)
576  rotPhi -= TwoPi;
577  double disPhi = std::fabs(rotPhi - PhiPmu);
578  if (disPhi > Pi)
579  disPhi = TwoPi - disPhi;
580  if (disPhi < dPhi) {
581  NmuHitTargetSphere++;
582  }
583  }
584 
585  } //end imu for loop
586 
587  if (NmuHitTargetSphere < MultiMuonNmin)
588  continue; //while (Id_sf.size() < 2) loop
589 
590  //nmuons outgoing particle at SurFace
591  Id_sf.clear();
592  Px_sf.clear();
593  Py_sf.clear();
594  Pz_sf.clear();
595  E_sf.clear();
596  //M_sf_out.clear();
597  Vx_sf.clear();
598  Vy_sf.clear();
599  Vz_sf.clear();
600  T0_sf.clear();
601 
602  //nmuons particles at UnderGround
603  Id_ug.clear();
604  Px_ug.clear();
605  Py_ug.clear();
606  Pz_ug.clear();
607  E_ug.clear();
608  //M_ug.clear();
609  Vx_ug.clear();
610  Vy_ug.clear();
611  Vz_ug.clear();
612  T0_ug.clear();
613 
614  int Id_sf_this = 0;
615  double Px_sf_this = 0., Py_sf_this = 0., Pz_sf_this = 0.;
616  double E_sf_this = 0.;
617  //double M_sf_this=0.;
618  double Vx_sf_this = 0., Vy_sf_this = 0., Vz_sf_this = 0.;
619  double T0_sf_this = 0.;
620 
621  T0_at = SimTree->shower_GH_t0 * SpeedOfLight; // [mm]
622 
623  for (int imu = 0; imu < nmuons; ++imu) {
624  if (P_mu[imu] < MinP_CMS && AcptAllMu == false)
625  continue;
626  //for the time being only muons
627  if (SimTree->particle__ParticleID[imu] != 5 && SimTree->particle__ParticleID[imu] != 6)
628  continue;
629 
630  Id_sf_this = SimTree->particle__ParticleID[imu];
631  if (Id_sf_this == 5)
632  Id_sf_this = -13; //mu+
633  else if (Id_sf_this == 6)
634  Id_sf_this = 13; //mu-
635 
636  else if (Id_sf_this == 1)
637  Id_sf_this = 22; //gamma
638  else if (Id_sf_this == 2)
639  Id_sf_this = -11; //e+
640  else if (Id_sf_this == 3)
641  Id_sf_this = 11; //e-
642  else if (Id_sf_this == 7)
643  Id_sf_this = 111; //pi0
644  else if (Id_sf_this == 8)
645  Id_sf_this = 211; //pi+
646  else if (Id_sf_this == 9)
647  Id_sf_this = -211; //pi-
648  else if (Id_sf_this == 10)
649  Id_sf_this = 130; //KL0
650  else if (Id_sf_this == 11)
651  Id_sf_this = 321; //K+
652  else if (Id_sf_this == 12)
653  Id_sf_this = -321; //K-
654  else if (Id_sf_this == 13)
655  Id_sf_this = 2112; //n
656  else if (Id_sf_this == 14)
657  Id_sf_this = 2212; //p
658  else if (Id_sf_this == 15)
659  Id_sf_this = -2212; //pbar
660  else if (Id_sf_this == 16)
661  Id_sf_this = 310; //Ks0
662  else if (Id_sf_this == 17)
663  Id_sf_this = 221; //eta
664  else if (Id_sf_this == 18)
665  Id_sf_this = 3122; //Lambda
666 
667  else {
668  std::cout << "CosmicMuonGenerator.cc: Warning! Muon Id=" << Id_sf_this << " from file read in" << std::endl;
669  Id_sf_this = 99999; //trouble
670  }
671 
672  T0_sf_this = SimTree->particle__Time[imu] * SpeedOfLight; //in [mm]
673 
674  Px_sf_this = Px_mu[imu];
675  Py_sf_this = Py_mu[imu]; //Corsika down going particles defined in -z direction!
676  Pz_sf_this = Pz_mu[imu];
677  E_sf_this = sqrt(P_mu[imu] * P_mu[imu] + MuonMass * MuonMass);
678  Vx_sf_this = Vx_mu[imu];
679  Vy_sf_this = Vy_mu[imu]; //[mm] fixed at surface + PlugWidth
680  Vz_sf_this = Vz_mu[imu];
681 
682  OneMuoEvt.create(Id_sf_this,
683  Px_sf_this,
684  Py_sf_this,
685  Pz_sf_this,
686  E_sf_this,
687  MuonMass,
688  Vx_sf_this,
689  Vy_sf_this,
690  Vz_sf_this,
691  T0_sf_this);
692  // if angles are ok, propagate to target
693  if (goodOrientation()) {
695  }
696 
698  AcptAllMu == true) {
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);
704  //M_sf.push_back(M_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);
709  //T0_sf.push_back(0.); //synchronised arrival for 100% efficient full simulation tests
710 
711  Id_ug.push_back(OneMuoEvt.id());
712  Px_ug.push_back(OneMuoEvt.px());
713  Py_ug.push_back(OneMuoEvt.py());
714  Pz_ug.push_back(OneMuoEvt.pz());
715  E_ug.push_back(OneMuoEvt.e());
716  //M_sf.push_back(OneMuoEvt.m());
717  Vx_ug.push_back(OneMuoEvt.vx());
718  Vy_ug.push_back(OneMuoEvt.vy());
719  Vz_ug.push_back(OneMuoEvt.vz());
720  T0_ug.push_back(OneMuoEvt.t0());
721 
722  NmuHitTarget++;
723  }
724  }
725 
726  } // while (Id_sf.size() < 2); //end of do loop
727 
728  if (trials > max_Trials) {
729  std::cout << "CosmicMuonGenerator.cc: Warning! trials reach max_trials=" << max_Trials
730  << " without accepting event!" << std::endl;
731  if (Debug) {
732  std::cout << " N(mu)=" << Id_ug.size();
733  if (!Id_ug.empty())
734  std::cout << " E[0]=" << E_ug[0] << " theta="
735  << acos(-Py_ug[0] / sqrt(Px_ug[0] * Px_ug[0] + Py_ug[0] * Py_ug[0] + Pz_ug[0] * Pz_ug[0]))
736  << " R_xz=" << sqrt(Vx_sf[0] * Vx_sf[0] + Vy_sf[0] * Vy_sf[0]);
737  std::cout << " Theta_at=" << Theta_at << std::endl;
738  }
739  std::cout << "Unweighted int num of Trials = " << Trials << std::endl;
740  std::cout << "trying next event (" << SimTree_jentry << ") from file" << std::endl;
742  continue; //EvtRejected while loop: get next event from file
743  } else {
744  if (NmuHitTarget < MultiMuonNmin) {
745  std::cout << "CosmicMuonGenerator.cc: Warning! less than " << MultiMuonNmin << " muons hit target: N(mu=)"
746  << NmuHitTarget << std::endl;
747  std::cout << "trying next event (" << SimTree_jentry << ") from file" << std::endl;
749  continue; //EvtRejected while loop: get next event from file
750  } else { //if (MuInMaxDist) {
751 
752  //re-adjust T0's of surviving muons shifted to trigger time box
753  //(possible T0 increase due to propagation (loss of energy/speed + travelled distance))
754  double T0_ug_min, T0_ug_max;
755  T0_ug_min = T0_ug_max = T0_ug[0];
756  double Tbox = (MaxT0 - MinT0) * SpeedOfLight; // [mm]
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)
761  T0_ug_min = T0_this;
762  if (T0_this > T0_ug_max)
763  T0_ug_max = T0_this;
764  if (Debug)
765  std::cout << "imu=" << imu << " T0_this=" << T0_this
766  << " P=" << sqrt(pow(Px_ug[imu], 2) + pow(Py_ug[imu], 2) + pow(Pz_ug[imu], 2)) << std::endl;
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]);
770  }
771  }
772 
773  if (int(Id_ug.size()) >= MultiMuonNmin && MultiMuonNmin >= 2 && minDeltaT0 > Tbox)
774  continue; //EvtRejected while loop: get next event from file
775 
776  double T0_min = T0_ug_min + MinT0 * SpeedOfLight; //-12.5ns * c [mm]
777  double T0_max = T0_ug_max + MaxT0 * SpeedOfLight; //+12.5ns * c [mm]
778 
779  //ckeck if >= NmuMin in time box, else throw new random number + augment evt weight
780  int TboxTrials = 0;
781  int NmuInTbox;
782  double T0_offset, T0diff;
783  for (int tboxtrial = 0; tboxtrial < 1000; ++tboxtrial) { //max 1000 trials
784  T0_offset = RanGen->flat() * (T0_max - T0_min); // [mm]
785  TboxTrials++;
786  T0diff = T0_offset - T0_max; // [mm]
787  NmuInTbox = 0;
788  for (unsigned int imu = 0; imu < Id_ug.size(); ++imu) {
789  if (T0_ug[imu] + T0diff > MinT0 * SpeedOfLight && T0_ug[imu] + T0diff < MaxT0 * SpeedOfLight)
790  NmuInTbox++;
791  }
792  if (NmuInTbox >= MultiMuonNmin)
793  break;
794  }
795  if (NmuInTbox < MultiMuonNmin)
796  continue; //EvtRejected while loop: get next event from file
797 
798  if (Debug)
799  std::cout << "initial T0_at=" << T0_at << " T0_min=" << T0_min << " T0_max=" << T0_max
800  << " T0_offset=" << T0_offset;
801  T0_at += T0diff; //[mm]
802  if (Debug)
803  std::cout << " T0diff=" << T0diff << std::endl;
804  for (unsigned int imu = 0; imu < Id_ug.size(); ++imu) { //adjust @ surface + underground
805  if (Debug)
806  std::cout << "before: T0_sf[" << imu << "]=" << T0_sf[imu] << " T0_ug=" << T0_ug[imu];
807  T0_sf[imu] += T0diff;
808  T0_ug[imu] += T0diff;
809  if (Debug)
810  std::cout << " after: T0_sf[" << imu << "]=" << T0_sf[imu] << " T0_ug=" << T0_ug[imu] << std::endl;
811  }
812  if (Debug)
813  std::cout << "T0diff=" << T0diff << " T0_at=" << T0_at << std::endl;
814 
815  Nsel += 1;
816  EventWeight = JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans / (trials * TboxTrials);
817  EvtRejected = false;
818  if (Debug)
819  std::cout << "CosmicMuonGenerator.cc: Theta_at=" << Theta_at << " phi_at=" << phi_at << " Px_at=" << Px_at
820  << " Py_at=" << Py_at << " Pz_at=" << Pz_at << " T0_at=" << T0_at << " Vx_at=" << Vx_at
821  << " Vy_at=" << Vy_at << " Vz_at=" << Vz_at << " EventWeight=" << EventWeight
822  << " Nmuons=" << Id_sf.size() << std::endl;
823  }
824  }
825 
826  } //while loop EvtRejected
827 
828  return true; //write event to HepMC;
829 }
830 
832  if (NumberOfEvents > 0) {
833  std::cout << std::endl;
834  std::cout << "*********************************************************" << std::endl;
835  std::cout << "*********************************************************" << std::endl;
836  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;
838  std::cout << "*** ***" << std::endl;
839  std::cout << "*********************************************************" << std::endl;
840  std::cout << "*********************************************************" << std::endl;
841  std::cout << std::endl;
842  std::cout << " number of initial cosmic events: " << int(Ngen) << std::endl;
843  std::cout << " number of actually diced events: " << int(Ndiced) << std::endl;
844  std::cout << " number of generated and accepted events: " << int(Nsel) << std::endl;
845  double selEff = Nsel / Ngen; // selection efficiency
846  std::cout << " event selection efficiency: " << selEff * 100. << "%" << std::endl;
847  int n100cos =
848  Norm->events_n100cos(0., 0.); //get final amount of cosmics in defined range for normalisation of flux
849  std::cout << " events with ~100 GeV and 1 - cos(theta) < 1/2pi: " << n100cos << std::endl;
850  std::cout << std::endl;
851  std::cout << " momentum range: " << MinP << " ... " << MaxP << " GeV" << std::endl;
852  std::cout << " theta range: " << MinTheta * Rad2Deg << " ... " << MaxTheta * Rad2Deg << " deg"
853  << std::endl;
854  std::cout << " phi range: " << MinPhi * Rad2Deg << " ... " << MaxPhi * Rad2Deg << " deg" << std::endl;
855  std::cout << " time range: " << MinT0 << " ... " << MaxT0 << " ns" << std::endl;
856  std::cout << " energy loss: " << ElossScaleFactor * 100. << "%" << std::endl;
857  std::cout << std::endl;
858  double area = 1.e-6 * Pi * SurfaceRadius * SurfaceRadius; // area on surface [m^2]
859  if (MinTheta > 90. * Deg2Rad) //upgoing muons from neutrinos)
860  std::cout << " area of initial cosmics at CMS detector bottom surface: " << area << " m^2" << std::endl;
861  else
862  std::cout << " area of initial cosmics on Surface + PlugWidth: " << area << " m^2" << std::endl;
863  std::cout << " depth of CMS detector (from Surface): " << SurfaceOfEarth / 1000 << " m" << std::endl;
864 
865  //at least 100 evts., and
866  //downgoing inside theta parametersisation range
867  //or upgoing neutrino muons
868  if ((n100cos > 0 && MaxTheta < 84.26 * Deg2Rad) || MinTheta > 90. * Deg2Rad) {
869  // rate: corrected for area and selection-Eff. and normalized to known flux, integration over solid angle (dOmega) is implicit
870  // flux is normalised with respect to known flux of vertical 100GeV muons in area at suface level
871  // rate seen by detector is lower than rate at surface area, so has to be corrected for selection-Eff.
872  // normalisation factor has unit [1/s/m^2]
873  // rate = N/time --> normalization factor gives 1/runtime/area
874  // normalization with respect to number of actually diced events (Ndiced)
875 
876  if (MinTheta > 90. * Deg2Rad) { //upgoing muons from neutrinos)
877  double Omega = (cos(MinTheta) - cos(MaxTheta)) * (MaxPhi - MinPhi);
878  //EventRate = (Ndiced * 3.e-13) * Omega * area*1.e4 * selEff;//area in cm, flux=3.e-13cm^-2s^-1sr^-1
879  EventRate = (Ndiced * 3.e-13) * Omega * 4. * RadiusOfTarget * ZDistOfTarget * 1.e-2 *
880  selEff; //area in cm, flux=3.e-13cm^-2s^-1sr^-1
881  rateErr_stat = EventRate / sqrt((double)Ndiced); // stat. rate error
882  rateErr_syst = EventRate / 3.e-13 * 1.0e-13; // syst. rate error, from error of known flux
883  } else {
884  EventRate = (Ndiced * Norm->norm(n100cos)) * area * selEff;
885  rateErr_stat = EventRate / sqrt((double)n100cos); // stat. rate error
886  rateErr_syst = EventRate / 2.63e-3 * 0.06e-3; // syst. rate error, from error of known flux
887  }
888 
889  // normalisation in region 1.-cos(theta) < 1./(2.*Pi), if MaxTheta even lower correct for this
890  if (MaxTheta < 0.572) {
891  double spacean = 2. * Pi * (1. - cos(MaxTheta));
892  EventRate = (Ndiced * Norm->norm(n100cos)) * area * selEff * spacean;
893  rateErr_stat = EventRate / sqrt((double)n100cos); // rate error
894  rateErr_syst = EventRate / 2.63e-3 * 0.06e-3; // syst. rate error, from error of known flux
895  }
896 
897  } else {
898  EventRate = Nsel; //no info as no muons at 100 GeV
899  rateErr_stat = Nsel;
900  rateErr_syst = Nsel;
901  std::cout << std::endl;
902  if (MinP > 100.)
903  std::cout << " !!! MinP > 100 GeV. Cannot apply normalisation!" << std::endl;
904  else if (MaxTheta > 84.26 * Deg2Rad)
905  std::cout << " !!! Note: generated cosmics exceed parameterisation. No flux calculated!" << std::endl;
906 
907  else
908  std::cout << " !!! Not enough statistics to apply normalisation (rate=1 +- 1) !!!" << std::endl;
909  }
910 
911  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
912  std::cout << " rate is " << EventRate << " +-" << rateErr_stat << " (stat) "
913  << "+-" << rateErr_syst << " (syst) "
914  << " muons per second" << std::endl;
915  if (EventRate != 0)
916  std::cout << " number of events corresponds to " << Nsel / EventRate << " s"
917  << std::endl; //runtime at CMS = Nsel/rate
918  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
919  std::cout << std::endl;
920  std::cout << "*********************************************************" << std::endl;
921  std::cout << "*********************************************************" << std::endl;
922  }
923 }
924 
926  if (MinP < 0.) {
927  NumberOfEvents = 0;
928  std::cout << " CMG-ERR: min.energy is out of range (0 GeV ... inf]" << std::endl << std::endl;
929  }
930  if (MaxP < 0.) {
931  NumberOfEvents = 0;
932  std::cout << " CMG-ERR: max.energy is out of range (0 GeV ... inf]" << std::endl << std::endl;
933  }
934  if (MaxP <= MinP) {
935  NumberOfEvents = 0;
936  std::cout << " CMG-ERR: max.energy is not greater than min.energy" << std::endl << std::endl;
937  }
938  if (MinTheta < 0.) {
939  NumberOfEvents = 0;
940  std::cout << " CMG-ERR: min.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl;
941  }
942  if (MaxTheta < 0.) {
943  NumberOfEvents = 0;
944  std::cout << " CMG-ERR: max.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl;
945  }
946  if (MaxTheta <= MinTheta) {
947  NumberOfEvents = 0;
948  std::cout << " CMG-ERR: max.theta is not greater than min.theta" << std::endl << std::endl;
949  }
950  if (MinPhi < 0.) {
951  NumberOfEvents = 0;
952  std::cout << " CMG-ERR: min.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl;
953  }
954  if (MaxPhi < 0.) {
955  NumberOfEvents = 0;
956  std::cout << " CMG-ERR: max.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl;
957  }
958  if (MaxPhi <= MinPhi) {
959  NumberOfEvents = 0;
960  std::cout << " CMG-ERR: max.phi is not greater than min.phi" << std::endl << std::endl;
961  }
962  if (MaxT0 <= MinT0) {
963  NumberOfEvents = 0;
964  std::cout << " CMG-ERR: max.t0 is not greater than min.t0" << std::endl << std::endl;
965  }
966  if (ElossScaleFactor < 0.) {
967  NumberOfEvents = 0;
968  std::cout << " CMG-ERR: E-loss scale factor is out of range [0 ... inf)" << std::endl << std::endl;
969  }
970  if (MinEnu < 0.) {
971  NumberOfEvents = 0;
972  std::cout << " CMG-ERR: min.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl;
973  }
974  if (MaxEnu < 0.) {
975  NumberOfEvents = 0;
976  std::cout << " CMG-ERR: max.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl;
977  }
978  if (MaxEnu <= MinEnu) {
979  NumberOfEvents = 0;
980  std::cout << " CMG-ERR: max.Enu is not greater than min.Enu" << std::endl << std::endl;
981  }
982 }
983 
985  // check angular range (for a sphere with Target3dRadius around the target)
986  bool goodAngles = false;
987  bool phiaccepted = false;
988  bool thetaaccepted = false;
989  double RxzV = sqrt(OneMuoEvt.vx() * OneMuoEvt.vx() + OneMuoEvt.vz() * OneMuoEvt.vz());
990 
991  double rVY;
992  if (MinTheta > 90. * Deg2Rad) //upgoing muons from neutrinos
993  rVY = -sqrt(RxzV * RxzV + RadiusCMS * RadiusCMS);
994  else
995  rVY = sqrt(RxzV * RxzV + (SurfaceOfEarth + PlugWidth) * (SurfaceOfEarth + PlugWidth));
996 
997  double Phi = OneMuoEvt.phi();
998  double PhiV = atan2(OneMuoEvt.vx(), OneMuoEvt.vz()) + Pi;
999  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  //std::cout << " dPhi = " << dPhi << " (" << Phi << " <p|V> " << PhiV << ")" << std::endl;
1015  //std::cout << " dTheta = " << dTheta << " (" << Theta << " <p|V> " << ThetaV << ")" << std::endl;
1016 
1017  if (!phiaccepted && RxzV < Target3dRadius)
1018  //if (RxzV < Target3dRadius)
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 
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 
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.); // [m]
1088  float verY = float(OneMuoEvt.vy() / 1000.); // [m]
1089  float verZ = float(OneMuoEvt.vz() / 1000.); // [m]
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.; // [mm]
1100  float absZ = std::fabs(verZ) * 1000.; // [mm]
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 
1119  if (NotInitialized)
1120  NumberOfEvents = N;
1121 }
1122 
1124  if (NotInitialized)
1125  RanSeed = N;
1126 }
1127 
1129  if (NotInitialized)
1130  MinP = P;
1131 }
1132 
1134  if (NotInitialized)
1135  MinP_CMS = P;
1136 }
1137 
1139  if (NotInitialized)
1140  MaxP = P;
1141 }
1142 
1144  if (NotInitialized)
1145  MinTheta = Theta * Deg2Rad;
1146 }
1147 
1149  if (NotInitialized)
1150  MaxTheta = Theta * Deg2Rad;
1151 }
1152 
1154  if (NotInitialized)
1155  MinPhi = Phi * Deg2Rad;
1156 }
1157 
1159  if (NotInitialized)
1160  MaxPhi = Phi * Deg2Rad;
1161 }
1162 
1164  if (NotInitialized)
1165  MinT0 = T0;
1166 }
1167 
1169  if (NotInitialized)
1170  MaxT0 = T0;
1171 }
1172 
1173 void CosmicMuonGenerator::setElossScaleFactor(double ElossScaleFact) {
1174  if (NotInitialized)
1175  ElossScaleFactor = ElossScaleFact;
1176 }
1177 
1179  if (NotInitialized)
1180  RadiusOfTarget = R;
1181 }
1182 
1184  if (NotInitialized)
1185  ZDistOfTarget = Z;
1186 }
1187 
1189  if (NotInitialized)
1190  ZCentrOfTarget = Z;
1191 }
1192 
1194  if (NotInitialized)
1195  TrackerOnly = Tracker;
1196 }
1197 
1199  if (NotInitialized)
1200  MultiMuon = MultiMu;
1201 }
1203  if (NotInitialized)
1204  MultiMuonFileName = MultiMuFile;
1205 }
1207  if (NotInitialized)
1208  MultiMuonFileFirstEvent = MultiMuFile1stEvt;
1209 }
1211  if (NotInitialized)
1212  MultiMuonNmin = MultiMuNmin;
1213 }
1214 
1216  if (NotInitialized)
1217  TIFOnly_constant = TIF;
1218 }
1219 
1221  if (NotInitialized)
1222  TIFOnly_linear = TIF;
1223 }
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 
1263  if (NotInitialized)
1264  MinEnu = MinEn;
1265 }
1267  if (NotInitialized)
1268  MaxEnu = MaxEn;
1269 }
1270 void CosmicMuonGenerator::setNuProdAlt(double NuPrdAlt) {
1271  if (NotInitialized)
1272  NuProdAlt = NuPrdAlt;
1273 }
1274 
1276 
1278  if (NotInitialized)
1279  AcptAllMu = AllMu;
1280 }
Int_t shower_nParticlesWritten
Definition: sim.h:41
void setZDistOfTarget(double Z)
const double TwoPi
const double Pi
void propagate(double ElossScaleFac, double RadiusTarget, double Z_DistTarget, double Z_CentrTarget, bool TrackerOnly, bool MTCCHalf)
float norm(int n100cos)
Definition: CMSCGENnorm.cc:27
std::vector< double > Theta_mu
int events_n100cos(double energy, double theta)
Definition: CMSCGENnorm.cc:13
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_]
Definition: sim.h:62
std::vector< double > Py_mu
void setRhoAir(double VarRhoAir)
const double Z_DistTracker
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::vector< int > Id_ug
TTree * fChain
Definition: sim.h:21
int initialize(double, double, double, double, CLHEP::HepRandomEngine *, bool, bool)
Definition: CMSCGEN.cc:27
void setRandomEngine(CLHEP::HepRandomEngine *v)
Definition: CMSCGEN.cc:20
void setRadiusOfTarget(double R)
Double_t particle__Time[kMaxparticle_]
Definition: sim.h:66
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_]
Definition: sim.h:61
Int_t shower_EventID
Definition: sim.h:28
std::vector< double > Pz_mu
const double SpeedOfLight
void initialize(CLHEP::HepRandomEngine *rng=nullptr)
Double_t particle__Pz[kMaxparticle_]
Definition: sim.h:63
std::vector< double > E_ug
void setMultiMuonFileFirstEvent(int MultiMuFile1stEvt)
const double NorthCMSzDeltaPhi
virtual Int_t GetEntry(Long64_t entry)
void setRhoPlug(double VarRhoPlug)
double cos_theta()
Definition: CMSCGEN.cc:431
const double PlugWidth
Float_t shower_Energy
Definition: sim.h:29
std::vector< double > Vz_ug
int generate()
Definition: CMSCGEN.cc:249
void setMinPhi(double Phi)
const double RadiusCMS
const double SurfaceOfEarth
void setMaxPhi(double Phi)
void setRandomEngine(CLHEP::HepRandomEngine *v)
std::vector< double > Vz_sf
std::vector< double > Vz_mu
const double Deg2Rad
T sqrt(T t)
Definition: SSEVec.h:23
std::vector< double > Vx_sf
const double Z_DistCMS
void setMultiMuonNmin(int MultiMuNmin)
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Float_t shower_Phi
Definition: sim.h:34
std::vector< double > Vy_sf
void setMinTheta(double Theta)
CLHEP::HepRandomEngine * RanGen
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
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
#define N
Definition: blowfish.cc:9
std::vector< double > Px_sf
std::vector< double > Vy_ug
void setAcptAllMu(bool AllMu)
Int_t particle__ParticleID[kMaxparticle_]
Definition: sim.h:58
double momentum_times_charge()
Definition: CMSCGEN.cc:422
void setPlugVz(double PlugVtz)
void setTIFOnly_linear(bool TIF)
std::pair< OmniClusterRef, TrackingParticleRef > P
std::vector< double > Pz_ug
const double MinStepSize
Double_t particle__x[kMaxparticle_]
Definition: sim.h:64
void create(int id, double px, double py, double pz, double e, double m, double vx, double vy, double vz, double t0)
const bool EventDisplay
const int max_Trials
const double MuonMass
void setMultiMuonFileName(std::string MultiMuonFileName)
const double RadiusTracker
void setElossScaleFactor(double ElossScaleFact)
void setMaxTheta(double Theta)
ALPAKA_FN_ACC ALPAKA_FN_INLINE uint32_t iy(uint32_t id)
void setRhoWall(double VarRhoSWall)
const double Rad2Deg
void setPlugVx(double PlugVtx)
Float_t shower_GH_t0
Definition: sim.h:47
void setRhoRock(double VarRhoRock)
int initializeNuMu(double, double, double, double, double, double, double, double, double, CLHEP::HepRandomEngine *)
Definition: CMSCGEN.cc:450
std::vector< double > Vx_ug
int generateNuMu()
Definition: CMSCGEN.cc:556
std::vector< double > P_mu
std::vector< double > Py_ug
const bool Debug
Float_t shower_Theta
Definition: sim.h:33
std::vector< double > T0_ug
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
std::vector< int > Id_sf
void setRhoClay(double VarRhoClay)
std::vector< double > Pz_sf
Double_t particle__y[kMaxparticle_]
Definition: sim.h:65
void setTrackerOnly(bool Tracker)