CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 
8 
10 
11 
13  initialize();
14  for (unsigned int iGen=0; iGen<NumberOfEvents; ++iGen){ nextEvent(); }
15  terminate();
16 }
17 
18 void CosmicMuonGenerator::initialize(CLHEP::HepRandomEngine *rng){
19  if (delRanGen)
20  delete RanGen;
21  if (!rng) {
22  RanGen = new CLHEP::HepJamesRandom;
23  RanGen->setSeed(RanSeed, 0); //set seed for Random Generator (seed can be controled by config-file)
24  delRanGen = true;
25  } else {
26  RanGen = rng;
27  delRanGen = false;
28  }
29  checkIn();
30  if (NumberOfEvents > 0){
31  // set up "surface geometry" dimensions
32  double RadiusTargetEff = RadiusOfTarget; //get this from cfg-file
33  double Z_DistTargetEff = ZDistOfTarget; //get this from cfg-file
34  //double Z_CentrTargetEff = ZCentrOfTarget; //get this from cfg-file
35  if(TrackerOnly==true){
36  RadiusTargetEff = RadiusTracker;
37  Z_DistTargetEff = Z_DistTracker;
38  }
39  Target3dRadius = sqrt(RadiusTargetEff*RadiusTargetEff + Z_DistTargetEff*Z_DistTargetEff) + MinStepSize;
40  if (Debug) std::cout << " radius of sphere around target = " << Target3dRadius << " mm" << std::endl;
41 
42  if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
44  else
46  if (Debug) std::cout << " starting point radius at Surface + PlugWidth = " << SurfaceRadius << " mm" << std::endl;
47 
55  //set energy and angle limits for CMSCGEN, give same seed as above
56  if (MinTheta >= 90.*Deg2Rad) //upgoing muons from neutrinos
59  else
61 
62 #if ROOT_INTERACTIVE
63  // book histos
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.);
68 #endif
69  if (EventDisplay) initEvDis();
70  std::cout << std::endl;
71 
72  if (MultiMuon) {
73  MultiIn = 0;
74 
75  std::cout << "MultiMuonFileName.c_str()=" << MultiMuonFileName.c_str() << std::endl;
76  MultiIn = new TFile( MultiMuonFileName.c_str() );
77 
78  if (!MultiIn) std::cout << "MultiMuon=True: MultiMuonFileName='"
79  << MultiMuonFileName.c_str() << "' does not exist" << std::endl;
80  else std::cout << "MultiMuonFile: " << MultiMuonFileName.c_str() << " opened!" << std::endl;
81  //MultiTree = (TTree*) gDirectory->Get("sim");
82  MultiTree = (TTree*) MultiIn->Get("sim");
83  SimTree = new sim(MultiTree);
85  SimTreeEntries = SimTree->fChain->GetEntriesFast();
86  std::cout << "SimTreeEntries=" << SimTreeEntries << std::endl;
87 
88  if (MultiMuonFileFirstEvent <= 0)
89  SimTree_jentry = 0;
90  else
91  SimTree_jentry = MultiMuonFileFirstEvent - 1; //1=1st evt (SimTree_jentry=0)
92 
95  }
96 
97  if (!MultiMuon || (MultiMuon && MultiIn)) NotInitialized = false;
98 
99  }
100 }
101 
103 
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;
106  // generate cosmic (E,theta,phi)
107  bool notSelected = true;
108  while (notSelected){
109  bool badMomentumGenerated = true;
110  while (badMomentumGenerated){
111 
112  if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
114  else
115  Cosmics->generate(); //dice one event now
116 
118  Theta = TMath::ACos( Cosmics->cos_theta() ) ; //angle has to be in RAD here
119  Ngen+=1.; //count number of initial cosmic events (in surface area), vertices will be added later
120  badMomentumGenerated = false;
121  Phi = RanGen->flat()*(MaxPhi-MinPhi) + MinPhi;
122  }
123  Norm->events_n100cos(E, Theta); //test if this muon is in normalization range
124  Ndiced += 1; //one more cosmic is diced
125 
126  // generate vertex
127  double Nver = 0.;
128  bool badVertexGenerated = true;
129  while (badVertexGenerated){
130  RxzV = sqrt(RanGen->flat())*SurfaceRadius;
131  PhiV = RanGen->flat()*TwoPi;
132  // check phi range (for a sphere with Target3dRadius around the target)
133  double dPhi = Pi; if (RxzV > Target3dRadius) dPhi = asin(Target3dRadius/RxzV);
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;
137  Nver+=1.;
138  }
139  Ngen += (Nver-1.); //add number of generated vertices to initial cosmic events
140 
141  // complete event at surface
142  int id = 13; // mu-
143  if (Cosmics->momentum_times_charge() >0.) id = -13; // mu+
144  double absMom = sqrt(E*E - MuonMass*MuonMass);
145  double verMom = absMom*cos(Theta);
146  double horMom = absMom*sin(Theta);
147  double Px = horMom*sin(Phi); // [GeV/c]
148  double Py = -verMom; // [GeV/c]
149  double Pz = horMom*cos(Phi); // [GeV/c]
150  double Vx = RxzV*sin(PhiV); // [mm]
151 
152  double Vy;
153  if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
154  Vy = -RadiusCMS;
155  else
156  Vy = SurfaceOfEarth + PlugWidth; // [mm]
157 
158  double Vz = RxzV*cos(PhiV); // [mm]
159  double T0 = (RanGen->flat()*(MaxT0-MinT0) + MinT0)*SpeedOfLight; // [mm/c];
160 
161  Id_at = id;
162  Px_at = Px; Py_at = Py; Pz_at = Pz; E_at = E; //M_at = MuonMass;
163  Vx_at = Vx; Vy_at = Vy; Vz_at = Vz; T0_at = T0;
164 
165  OneMuoEvt.create(id, Px, Py, Pz, E, MuonMass, Vx, Vy, Vz, T0);
166  // if angles are ok, propagate to target
167  if (goodOrientation()) {
168  if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
170  else
172  }
173 
174  if ( (OneMuoEvt.hitTarget() && sqrt(OneMuoEvt.e()*OneMuoEvt.e() - MuonMass*MuonMass) > MinP_CMS)
175  || AcptAllMu==true){
176  Nsel+=1.; //count number of generated and accepted events
177  notSelected = false;
178  }
179  }
180 
181  EventWeight = 1.;
182 
183  //just one outgoing particle at SurFace
184  Id_sf.resize(1);
185  Px_sf.resize(1);
186  Py_sf.resize(1);
187  Pz_sf.resize(1);
188  E_sf.resize(1);
189  //M_sf.resize(1);
190  Vx_sf.resize(1);
191  Vy_sf.resize(1);
192  Vz_sf.resize(1);
193  T0_sf.resize(1);
194 
195  Id_sf[0] = Id_at;
196  Px_sf[0] = Px_at; Py_sf[0] = Py_at; Pz_sf[0] = Pz_at; E_sf[0] = E_at; //M_fs[0] = MuonMass;
197  Vx_sf[0] = Vx_at; Vy_sf[0] = Vy_at; Vz_sf[0] = Vz_at; T0_sf[0] = T0_at;
198 
199 
200  //just one particle at UnderGround
201  Id_ug.resize(1);
202  Px_ug.resize(1);
203  Py_ug.resize(1);
204  Pz_ug.resize(1);
205  E_ug.resize(1);
206  //M_ug.resize(1);
207  Vx_ug.resize(1);
208  Vy_ug.resize(1);
209  Vz_ug.resize(1);
210  T0_ug.resize(1);
211 
212  Id_ug[0] = OneMuoEvt.id();
213  Px_ug[0] = OneMuoEvt.px();
214  Py_ug[0] = OneMuoEvt.py();
215  Pz_ug[0] = OneMuoEvt.pz();
216  E_ug[0] = OneMuoEvt.e();
217  //M_ug[0] = OneMuoEvt.m();
218  Vx_ug[0] = OneMuoEvt.vx();
219  Vy_ug[0] = OneMuoEvt.vy();
220  Vz_ug[0] = OneMuoEvt.vz();
221  T0_ug[0] = OneMuoEvt.t0();
222 
223  // plot variables of selected events
224 #if ROOT_INTERACTIVE
225  ene->Fill(OneMuoEvt.e());
226  the->Fill((OneMuoEvt.theta()*Rad2Deg));
227  phi->Fill((OneMuoEvt.phi()*Rad2Deg));
228  ver->Fill((OneMuoEvt.vz()/1000.),(OneMuoEvt.vx()/1000.),(OneMuoEvt.vy()/1000.));
229 #endif
230  if (Debug){
231  std::cout << "new event" << std::endl;
232  std::cout << " Px,Py,Pz,E,m = " << OneMuoEvt.px() << ", " << OneMuoEvt.py() << ", "
233  << OneMuoEvt.pz() << ", " << OneMuoEvt.e() << ", " << OneMuoEvt.m() << " GeV" << std::endl;
234  std::cout << " Vx,Vy,Vz,t0 = " << OneMuoEvt.vx() << ", " << OneMuoEvt.vy() << ", "
235  << OneMuoEvt.vz() << ", " << OneMuoEvt.t0() << " mm" << std::endl;
236  }
237  if (EventDisplay) displayEv();
238 
239 }
240 
241 
242 
244 
245  if (Debug) std::cout << "\nEntered CosmicMuonGenerator::nextMultiEvent()" << std::endl;
246  bool EvtRejected = true;
247  bool MuInMaxDist = false;
248  double MinDist; //[mm]
249 
250  while (EvtRejected) {
251 
252  //read in event from SimTree
253  //ULong64_t ientry = SimTree->LoadTree(SimTree_jentry);
254  Long64_t ientry = SimTree->GetEntry(SimTree_jentry);
255  std::cout << "CosmicMuonGenerator::nextMultiEvent(): SimTree_jentry=" << SimTree_jentry
256  //<< " ientry=" << ientry
257  << " SimTreeEntries=" << SimTreeEntries << std::endl;
258  if (ientry < 0) return false; //stop run
260  SimTree_jentry++;
261  }
262  else {
263  std::cout << "CosmicMuonGenerator.cc::nextMultiEvent: No more events in file!" << std::endl;
264  return false; //stop run
265  }
266 
267 
268 
269  int nmuons = SimTree->shower_nParticlesWritten;
270  if (nmuons<MultiMuonNmin) {
271  std::cout << "CosmicMuonGenerator.cc: Warning! Less than " << MultiMuonNmin <<
272  " muons in event!" << std::endl;
273  std::cout << "trying next event from file" << std::endl;
275  continue; //EvtRejected while loop: get next event from file
276  }
277 
278 
279 
280  Px_mu.resize(nmuons); Py_mu.resize(nmuons); Pz_mu.resize(nmuons);
281  P_mu.resize(nmuons);
282 
283  MinDist = 99999.e9; //[mm]
284  double MuMuDist;
285  MuInMaxDist = false;
286  //check if at least one muon pair closer than 30m at surface
287  int NmuPmin = 0;
288  for (int imu=0; imu<nmuons; ++imu) {
289 
294  Py_mu[imu] = -SimTree->particle__Pz[imu]; //Corsika down going particles defined in -z direction!
295  P_mu[imu] = sqrt(Px_mu[imu]*Px_mu[imu] + Py_mu[imu]*Py_mu[imu] + Pz_mu[imu]*Pz_mu[imu]);
296 
297  if (P_mu[imu] < MinP_CMS && AcptAllMu==false) continue;
298  else if (SimTree->particle__ParticleID[imu] != 5 &&
299  SimTree->particle__ParticleID[imu] != 6) continue;
300  else NmuPmin++;
301 
302  for (int jmu=0; jmu<imu; ++jmu) {
303  if (P_mu[jmu] < MinP_CMS && AcptAllMu==false) continue;
304  if (SimTree->particle__ParticleID[imu] != 5 &&
305  SimTree->particle__ParticleID[imu] != 6) continue;
306  MuMuDist = sqrt( (SimTree->particle__x[imu]-SimTree->particle__x[jmu])*
307  (SimTree->particle__x[imu]-SimTree->particle__x[jmu])
308  +(SimTree->particle__y[imu]-SimTree->particle__y[jmu])*
309  (SimTree->particle__y[imu]-SimTree->particle__y[jmu])
310  )*10.; //CORSIKA [cm] to CMSCGEN [mm]
311  if (MuMuDist < MinDist) MinDist = MuMuDist;
312  if (MuMuDist < 2.*Target3dRadius) MuInMaxDist = true;
313  }
314  }
315  if (MultiMuonNmin>=2) {
316  if (MuInMaxDist) {
318  }
319  else {
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: "
323  << 1 - double(NcloseMultiMuonEvents)/SimTree_jentry << std::endl;
324  std::cout << "NcloseMultiMuonEvents=" << NcloseMultiMuonEvents << std::endl;
325  std::cout << "trying next event from file" << std::endl;
327  continue; //EvtRejected while loop: get next event from file
328  }
329  }
330 
331  if (NmuPmin < MultiMuonNmin && AcptAllMu==false) { //take single muon events consistently into account
333  continue; //EvtRejected while loop: get next event from file
334  }
335 
336  if (Debug)
337  if (MultiMuonNmin>=2)
338  std::cout << "start trial do loop: MuMuDist=" << MinDist/1000. << "[m] Nmuons="
339  << nmuons << " NcloseMultiMuonEvents=" << NcloseMultiMuonEvents
340  << " NskippedMultiMuonEvents=" << NskippedMultiMuonEvents << std::endl;
341 
342 
343  //int primary_id = SimTree->run_ParticleID;
345 
346  double M_at = 0.;
347  //if (Id_at == 13) {
348  Id_at = 2212; //convert from Corsika to HepPDT
349  M_at = 938.272e-3; //[GeV] mass
350  //}
351 
354  double phi_at = SimTree->shower_Phi - NorthCMSzDeltaPhi; //rotate by almost 90 degrees
355  if (phi_at < -Pi) phi_at +=TwoPi; //bring into interval (-Pi,Pi]
356  else if (phi_at > Pi) phi_at -= TwoPi;
357  double P_at = sqrt(E_at*E_at - M_at*M_at);
358  //need to rotate about 90degrees around x->N axis => y<=>z,
359  //then rotate new x-z-plane from x->North to x->LHC centre
360  Px_at = P_at*sin(Theta_at)*sin(phi_at);
361  Py_at = -P_at*cos(Theta_at);
362  Pz_at = P_at*sin(Theta_at)*cos(phi_at);
363 
364  //compute maximal theta of secondary muons
365  double theta_mu_max = Theta_at;
366  double theta_mu_min = Theta_at;
367 
368  double phi_rel_min = 0.; //phi_mu_min - phi_at
369  double phi_rel_max = 0.; //phi_mu_max - phi_at
370 
371  //std::cout << "SimTree->shower_Energy=" << SimTree->shower_Energy <<std::endl;
372 
373  Theta_mu.resize(nmuons);
374  for (int imu=0; imu<nmuons; ++imu) {
375  Theta_mu[imu] = acos(-Py_mu[imu]/P_mu[imu]);
376  if (Theta_mu[imu]>theta_mu_max) theta_mu_max = Theta_mu[imu];
377  if (Theta_mu[imu]<theta_mu_min) theta_mu_min = Theta_mu[imu];
378 
379  double phi_mu = atan2(Px_mu[imu],Pz_mu[imu]); // in (-Pi,Pi]
380  double phi_rel = phi_mu - phi_at;
381  if (phi_rel < -Pi) phi_rel += TwoPi; //bring into interval (-Pi,Pi]
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;
385 
386 
387  }
388 
389 
390  double h_sf = SurfaceOfEarth + PlugWidth; //[mm]
391 
392  double R_at = h_sf*tan(Theta_at);
393 
394  double JdRxzV_dR_trans = 1.;
395  double JdPhiV_dPhi_trans = 1.;
396  double JdR_trans_sqrt = 1.;
397 
398  //chose random vertex Phi and Rxz weighted to speed up and smoothen
399  double R_mu_max = (h_sf+Target3dRadius)*tan(theta_mu_max);
400  double R_max = std::min(SurfaceRadius, R_mu_max);
401  double R_mu_min = (h_sf-Target3dRadius)*tan(theta_mu_min);
402  double R_min = std::max(0., R_mu_min);
403 
404  if (R_at>SurfaceRadius) {
405  std::cout << "CosmicMuonGenerator.cc: Warning! R_at=" << R_at
406  << " > SurfaceRadius=" << SurfaceRadius <<std::endl;
407  }
408 
409  //do phase space transformation for horizontal radius R
410 
411  //determine first phase space limits
412 
413  double psR1min = R_min + 0.25*(R_max-R_min);
414  double psR1max = std::min(SurfaceRadius,R_max-0.25*(R_max-R_min)); //no R's beyond R_max
415  double psR1 = psR1max - psR1min;
416 
417  double psR2min = R_min;
418  double psR2max = R_max;
419  double psR2 = psR2max - psR2min;
420 
421  double psR3min = 0.;
422  double psR3max = SurfaceRadius;
423  double psR3 = psR3max - psR3min;
424 
425  //double psall = psR1+psR2+psR3;
426  double psRall = psR3;
427 
428  double fR1=psR1/psRall, fR2=psR2/psRall, fR3=psR3/psRall; //f1+f2+f3=130%
429  double pR1=0.25, pR2=0.7, pR3=0.05;
430 
431 
432  //do phase space transformation for azimuthal angle phi
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;
437 
438  double fPh1=psPh1/psPhall, fPh2=psPh2/psPhall, fPh3=psPh3/psPhall; //(f1+f2+f3=TwoPi+f1+f2)/(TwoPi+f1+f2)
439 
440  double pPh1=0.25, pPh2=0.7, pPh3=0.05;
441 
442  int temp = 0;
443  bool btemp = false;
444 
445  Trials = 0; //global int trials
446  double trials = 0.; //local weighted trials
447  Vx_mu.resize(nmuons); Vy_mu.resize(nmuons); Vz_mu.resize(nmuons);
448  int NmuHitTarget = 0;
449  while (NmuHitTarget < MultiMuonNmin) {
450 
451  NmuHitTarget = 0; //re-initialize every loop iteration
452  double Nver = 0.;
453 
454 
455  //chose phase space class
456  double RxzV;
457  double which_R_class = RanGen->flat();
458  if (which_R_class < pR1) { //pR1% in psR1
459  RxzV = psR1min + psR1 * RanGen->flat();
460  JdRxzV_dR_trans = fR1/pR1 * SurfaceRadius/psR1;
461  }
462  else if (which_R_class < pR1+pR2) { //further pR2% in psR2
463  RxzV = psR2min + psR2 * RanGen->flat();
464  JdRxzV_dR_trans = fR2/pR2 * SurfaceRadius/psR2;
465  }
466  else { //remaining pR3% in psR3=[0., R_max]
467  RxzV = psR3min + psR3 * RanGen->flat();
468  JdRxzV_dR_trans = fR3/pR3 * SurfaceRadius/psR3;
469  }
470 
471  JdR_trans_sqrt = 2.*RxzV/SurfaceRadius; //flat in sqrt(r) space
472 
473  //chose phase space class
474  double PhiV;
475  double which_phi_class = RanGen->flat();
476  if (which_phi_class < pPh1) { //pPh1% in psPh1
477  PhiV = phi_at + phi_rel_min + psPh1 * RanGen->flat();
478  JdPhiV_dPhi_trans = fPh1/pPh1 * TwoPi/psPh1;
479  }
480  else if (which_phi_class < pPh1+pPh2) { //further pPh2% in psPh2
481  PhiV = phi_at + phi_rel_min + psPh2 * RanGen->flat();
482  JdPhiV_dPhi_trans = fPh2/pPh2 * TwoPi/psPh2;
483  }
484  else { //remaining pPh3% in psPh3=[-Pi,Pi]
485  PhiV = phi_at + phi_rel_min + psPh3 * RanGen->flat();
486  JdPhiV_dPhi_trans = fPh3/pPh3 * TwoPi/psPh3;
487  }
488 
489  //shuffle PhiV into [-Pi,+Pi] interval
490  if (PhiV < -Pi) PhiV+=TwoPi;
491  else if (PhiV > Pi) PhiV-=TwoPi;
492 
493 
494  Nver++;
495  trials += JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans;
496  Trials++;
497  if (trials > max_Trials) break; //while (Id_sf.size() < 2) loop
498  Ngen += (Nver-1.); //add number of generated vertices to initial cosmic events
499 
500 
501  Vx_at = RxzV*sin(PhiV); // [mm]
502 
503  Vy_at = h_sf; // [mm] (SurfaceOfEarth + PlugWidth; Determine primary particle height below)
504  //Vy_at = SimTree->shower_StartingAltitude*10. + h_sf; // [mm]
505  //std::cout << "SimTree->shower_StartingAltitude*10=" << SimTree->shower_StartingAltitude*10 <<std::endl;
506  Vz_at = RxzV*cos(PhiV); // [mm]
507 
508  int NmuHitTargetSphere = 0;
509  for (int imu=0; imu<nmuons; ++imu) {
510 
512  +SimTree->particle__y[imu]*cos(NorthCMSzDeltaPhi) )*10; //[mm] (Corsika cm to CMSCGEN mm)
513  Vy_mu[imu] = h_sf; //[mm] fixed at surface + PlugWidth
515  +SimTree->particle__y[imu]*sin(NorthCMSzDeltaPhi) )*10; //[mm] (Corsika cm to CMSCGEN mm)
516 
517 
518  //add atmospheric height to primary particle (default SimTree->shower_StartingAltitude = 0.)
519  double pt_sec = sqrt(Px_mu[imu]*Px_mu[imu]+Pz_mu[imu]*Pz_mu[imu]);
520  double theta_sec = atan2(std::fabs(Py_mu[imu]),pt_sec);
521  double r_sec = sqrt((Vx_mu[imu]-Vx_at)*(Vx_mu[imu]-Vx_at)
522  +(Vz_mu[imu]-Vz_at)*(Vz_mu[imu]-Vz_at));
523  double h_prod = r_sec * tan(theta_sec);
524  if (h_prod + h_sf > Vy_at) Vy_at = h_prod + h_sf;
525 
526  //only muons
527  if (SimTree->particle__ParticleID[imu] != 5 &&
528  SimTree->particle__ParticleID[imu] != 6) continue;
529 
530  if (P_mu[imu] < MinP_CMS && AcptAllMu==false) continue;
531 
532  //check here if at least 2 muons make it to the target sphere
533  double Vxz_mu = sqrt(Vx_mu[imu]*Vx_mu[imu] + Vz_mu[imu]*Vz_mu[imu]);
534  theta_mu_max = atan((Vxz_mu+Target3dRadius)/(h_sf-Target3dRadius));
535  theta_mu_min = atan((Vxz_mu-Target3dRadius)/(h_sf+Target3dRadius));
536  if (Theta_mu[imu] > theta_mu_min && Theta_mu[imu] < theta_mu_max) {
537 
538  // check phi range (for a sphere with Target3dRadius around the target)
539  double dPhi = Pi; if (Vxz_mu > Target3dRadius) dPhi = asin(Target3dRadius/Vxz_mu);
540  double PhiPmu = atan2(Px_mu[imu],Pz_mu[imu]); //muon phi
541  double PhiVmu = atan2(Vx_mu[imu],Vz_mu[imu]); //muon phi
542  double rotPhi = PhiVmu + Pi; if (rotPhi > Pi) rotPhi -= TwoPi;
543  double disPhi = std::fabs(rotPhi - PhiPmu); if (disPhi > Pi) disPhi = TwoPi - disPhi;
544  if (disPhi < dPhi) {
545  NmuHitTargetSphere++;
546  }
547 
548  }
549 
550  } //end imu for loop
551 
552 
553  if (temp > 0) btemp = true;
554 
555 
556  if (NmuHitTargetSphere < MultiMuonNmin) continue; //while (Id_sf.size() < 2) loop
557 
558  //nmuons outgoing particle at SurFace
559  Id_sf.clear();
560  Px_sf.clear();
561  Py_sf.clear();
562  Pz_sf.clear();
563  E_sf.clear();
564  //M_sf_out.clear();
565  Vx_sf.clear();
566  Vy_sf.clear();
567  Vz_sf.clear();
568  T0_sf.clear();
569 
570  //nmuons particles at UnderGround
571  Id_ug.clear();
572  Px_ug.clear();
573  Py_ug.clear();
574  Pz_ug.clear();
575  E_ug.clear();
576  //M_ug.clear();
577  Vx_ug.clear();
578  Vy_ug.clear();
579  Vz_ug.clear();
580  T0_ug.clear();
581 
582  int Id_sf_this =0;
583  double Px_sf_this =0., Py_sf_this=0., Pz_sf_this=0.;
584  double E_sf_this=0.;
585  //double M_sf_this=0.;
586  double Vx_sf_this=0., Vy_sf_this=0., Vz_sf_this=0.;
587  double T0_sf_this=0.;
588 
589  T0_at = SimTree->shower_GH_t0 * SpeedOfLight; // [mm]
590 
591  for (int imu=0; imu<nmuons; ++imu) {
592 
593  if (P_mu[imu] < MinP_CMS && AcptAllMu==false) continue;
594  //for the time being only muons
595  if (SimTree->particle__ParticleID[imu] != 5 &&
596  SimTree->particle__ParticleID[imu] != 6) continue;
597 
598  Id_sf_this = SimTree->particle__ParticleID[imu];
599  if (Id_sf_this == 5) Id_sf_this = -13; //mu+
600  else if (Id_sf_this == 6) Id_sf_this = 13; //mu-
601 
602  else if (Id_sf_this == 1) Id_sf_this = 22; //gamma
603  else if (Id_sf_this == 2) Id_sf_this = -11; //e+
604  else if (Id_sf_this == 3) Id_sf_this = 11; //e-
605  else if (Id_sf_this == 7) Id_sf_this = 111; //pi0
606  else if (Id_sf_this == 8) Id_sf_this = 211; //pi+
607  else if (Id_sf_this == 9) Id_sf_this = -211; //pi-
608  else if (Id_sf_this == 10) Id_sf_this = 130; //KL0
609  else if (Id_sf_this == 11) Id_sf_this = 321; //K+
610  else if (Id_sf_this == 12) Id_sf_this = -321; //K-
611  else if (Id_sf_this == 13) Id_sf_this = 2112; //n
612  else if (Id_sf_this == 14) Id_sf_this = 2212; //p
613  else if (Id_sf_this == 15) Id_sf_this = -2212; //pbar
614  else if (Id_sf_this == 16) Id_sf_this = 310; //Ks0
615  else if (Id_sf_this == 17) Id_sf_this = 221; //eta
616  else if (Id_sf_this == 18) Id_sf_this = 3122; //Lambda
617 
618  else {
619  std::cout << "CosmicMuonGenerator.cc: Warning! Muon Id=" << Id_sf_this
620  << " from file read in" <<std::endl;
621  Id_sf_this = 99999; //trouble
622  }
623 
624  T0_sf_this = SimTree->particle__Time[imu] * SpeedOfLight; //in [mm]
625 
626  Px_sf_this = Px_mu[imu];
627  Py_sf_this = Py_mu[imu]; //Corsika down going particles defined in -z direction!
628  Pz_sf_this = Pz_mu[imu];
629  E_sf_this = sqrt(P_mu[imu]*P_mu[imu] + MuonMass*MuonMass);
630  Vx_sf_this = Vx_mu[imu];
631  Vy_sf_this = Vy_mu[imu]; //[mm] fixed at surface + PlugWidth
632  Vz_sf_this = Vz_mu[imu];
633 
634 
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);
636  // if angles are ok, propagate to target
637  if (goodOrientation()) {
639  }
640 
641  if ( (OneMuoEvt.hitTarget()
642  && sqrt(OneMuoEvt.e()*OneMuoEvt.e() - MuonMass*MuonMass) > MinP_CMS)
643  || AcptAllMu==true ) {
644 
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);
650  //M_sf.push_back(M_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);
655  //T0_sf.push_back(0.); //synchronised arrival for 100% efficient full simulation tests
656 
657  Id_ug.push_back(OneMuoEvt.id());
658  Px_ug.push_back(OneMuoEvt.px());
659  Py_ug.push_back(OneMuoEvt.py());
660  Pz_ug.push_back(OneMuoEvt.pz());
661  E_ug.push_back(OneMuoEvt.e());
662  //M_sf.push_back(OneMuoEvt.m());
663  Vx_ug.push_back(OneMuoEvt.vx());
664  Vy_ug.push_back(OneMuoEvt.vy());
665  Vz_ug.push_back(OneMuoEvt.vz());
666  T0_ug.push_back(OneMuoEvt.t0());
667 
668  NmuHitTarget++;
669  }
670  }
671 
672 
673  } // while (Id_sf.size() < 2); //end of do loop
674 
675 
676  if (trials > max_Trials) {
677  std::cout << "CosmicMuonGenerator.cc: Warning! trials reach max_trials=" << max_Trials
678  << " without accepting event!" << std::endl;
679  if (Debug) {
680  std::cout << " N(mu)=" << Id_ug.size();
681  if (Id_ug.size()>=1)
682  std::cout << " E[0]=" << E_ug[0] << " theta="
683  << acos(-Py_ug[0]/sqrt(Px_ug[0]*Px_ug[0]+Py_ug[0]*Py_ug[0]+Pz_ug[0]*Pz_ug[0]))
684  << " R_xz=" << sqrt(Vx_sf[0]*Vx_sf[0]+Vy_sf[0]*Vy_sf[0]);
685  std::cout << " Theta_at=" << Theta_at << std::endl;
686  }
687  std::cout << "Unweighted int num of Trials = " << Trials << std::endl;
688  std::cout << "trying next event (" << SimTree_jentry << ") from file" << std::endl;
690  continue; //EvtRejected while loop: get next event from file
691  }
692  else {
693  if (NmuHitTarget < MultiMuonNmin) {
694  std::cout << "CosmicMuonGenerator.cc: Warning! less than " << MultiMuonNmin <<
695  " muons hit target: N(mu=)" << NmuHitTarget << std::endl;
696  std::cout << "trying next event (" << SimTree_jentry << ") from file" << std::endl;
698  continue; //EvtRejected while loop: get next event from file
699  }
700  else { //if (MuInMaxDist) {
701 
702  //re-adjust T0's of surviving muons shifted to trigger time box
703  //(possible T0 increase due to propagation (loss of energy/speed + travelled distance))
704  double T0_ug_min, T0_ug_max;
705  T0_ug_min = T0_ug_max = T0_ug[0];
706  double Tbox = (MaxT0 - MinT0) * SpeedOfLight; // [mm]
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;
712  if (Debug) std::cout << "imu=" << imu << " T0_this=" << T0_this
713  << " P=" << sqrt(pow(Px_ug[imu],2) + pow(Py_ug[imu],2) + pow(Pz_ug[imu],2))
714  << std::endl;
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]);
717  }
718  }
719 
720  if (int(Id_ug.size()) >= MultiMuonNmin && MultiMuonNmin>=2 && minDeltaT0 > Tbox)
721  continue; //EvtRejected while loop: get next event from file
722 
723  double T0_min = T0_ug_min +MinT0*SpeedOfLight; //-12.5ns * c [mm]
724  double T0_max = T0_ug_max +MaxT0*SpeedOfLight; //+12.5ns * c [mm]
725 
726  //ckeck if >= NmuMin in time box, else throw new random number + augment evt weight
727  int TboxTrials = 0;
728  int NmuInTbox;
729  double T0_offset, T0diff;
730  for (int tboxtrial=0; tboxtrial<1000; ++tboxtrial) { //max 1000 trials
731  T0_offset = RanGen->flat()*(T0_max -T0_min); // [mm]
732  TboxTrials++;
733  T0diff = T0_offset - T0_max; // [mm]
734  NmuInTbox = 0;
735  for (unsigned int imu=0; imu<Id_ug.size(); ++imu) {
736  if (T0_ug[imu]+T0diff > MinT0*SpeedOfLight && T0_ug[imu]+T0diff < MaxT0*SpeedOfLight)
737  NmuInTbox++;
738  }
739  if (NmuInTbox >= MultiMuonNmin) break;
740 
741  }
742  if (NmuInTbox < MultiMuonNmin) continue; //EvtRejected while loop: get next event from file
743 
744 
745  if (Debug) std::cout << "initial T0_at=" << T0_at << " T0_min=" << T0_min << " T0_max=" << T0_max
746  << " T0_offset=" << T0_offset;
747  T0_at += T0diff; //[mm]
748  if (Debug) std::cout << " T0diff=" << T0diff << std::endl;
749  for (unsigned int imu=0; imu<Id_ug.size(); ++imu) { //adjust @ surface + underground
750  if (Debug) std::cout << "before: T0_sf[" << imu << "]=" << T0_sf[imu] << " T0_ug=" << T0_ug[imu];
751  T0_sf[imu] += T0diff;
752  T0_ug[imu] += T0diff;
753  if (Debug)
754  std::cout << " after: T0_sf[" << imu << "]=" << T0_sf[imu] << " T0_ug=" << T0_ug[imu] << std::endl;
755  }
756  if (Debug) std::cout << "T0diff=" << T0diff << " T0_at=" << T0_at << std::endl;
757 
758 
759 
760  Nsel += 1;
761  EventWeight = JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans
762  / (trials * TboxTrials);
763  EvtRejected = false;
764  if (Debug) std::cout << "CosmicMuonGenerator.cc: Theta_at=" << Theta_at << " phi_at=" << phi_at
765  << " Px_at=" << Px_at << " Py_at=" << Py_at << " Pz_at=" << Pz_at
766  << " T0_at=" << T0_at
767  << " Vx_at=" << Vx_at << " Vy_at=" << Vy_at << " Vz_at=" << Vz_at
768  << " EventWeight=" << EventWeight << " Nmuons=" << Id_sf.size() << std::endl;
769  }
770  }
771 
772 
773  } //while loop EvtRejected
774 
775 
776  return true; //write event to HepMC;
777 
778 }
779 
780 
781 
783  if (NumberOfEvents > 0){
784  std::cout << std::endl;
785  std::cout << "*********************************************************" << std::endl;
786  std::cout << "*********************************************************" << std::endl;
787  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;
789  std::cout << "*** ***" << std::endl;
790  std::cout << "*********************************************************" << std::endl;
791  std::cout << "*********************************************************" << std::endl;
792  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;
796  double selEff = Nsel/Ngen; // selection efficiency
797  std::cout << " event selection efficiency: " << selEff*100. << "%" << std::endl;
798  int n100cos = Norm->events_n100cos(0., 0.); //get final amount of cosmics in defined range for normalisation of flux
799  std::cout << " events with ~100 GeV and 1 - cos(theta) < 1/2pi: " << n100cos << std::endl;
800  std::cout << std::endl;
801  std::cout << " momentum range: " << MinP << " ... " << MaxP << " GeV" << std::endl;
802  std::cout << " theta range: " << MinTheta*Rad2Deg << " ... " << MaxTheta*Rad2Deg << " deg" << std::endl;
803  std::cout << " phi range: " << MinPhi*Rad2Deg << " ... " << MaxPhi*Rad2Deg << " deg" << std::endl;
804  std::cout << " time range: " << MinT0 << " ... " << MaxT0 << " ns" << std::endl;
805  std::cout << " energy loss: " << ElossScaleFactor*100. << "%" << std::endl;
806  std::cout << std::endl;
807  double area = 1.e-6*Pi*SurfaceRadius*SurfaceRadius; // area on surface [m^2]
808  if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos)
809  std::cout << " area of initial cosmics at CMS detector bottom surface: " << area << " m^2" << std::endl;
810  else
811  std::cout << " area of initial cosmics on Surface + PlugWidth: " << area << " m^2" << std::endl;
812  std::cout << " depth of CMS detector (from Surface): " << SurfaceOfEarth/1000 << " m" << std::endl;
813 
814  //at least 100 evts., and
815  //downgoing inside theta parametersisation range
816  //or upgoing neutrino muons
817  if( (n100cos>0 && MaxTheta<84.26*Deg2Rad)
818  || MinTheta>90.*Deg2Rad) {
819  // rate: corrected for area and selection-Eff. and normalized to known flux, integration over solid angle (dOmega) is implicit
820  // flux is normalised with respect to known flux of vertical 100GeV muons in area at suface level
821  // rate seen by detector is lower than rate at surface area, so has to be corrected for selection-Eff.
822  // normalisation factor has unit [1/s/m^2]
823  // rate = N/time --> normalization factor gives 1/runtime/area
824  // normalization with respect to number of actually diced events (Ndiced)
825 
826  if (MinTheta > 90.*Deg2Rad) {//upgoing muons from neutrinos)
827  double Omega = (cos(MinTheta)-cos(MaxTheta)) * (MaxPhi-MinPhi);
828  //EventRate = (Ndiced * 3.e-13) * Omega * area*1.e4 * selEff;//area in cm, flux=3.e-13cm^-2s^-1sr^-1
829  EventRate = (Ndiced * 3.e-13) * Omega * 4.*RadiusOfTarget*ZDistOfTarget*1.e-2 * selEff;//area in cm, flux=3.e-13cm^-2s^-1sr^-1
830  rateErr_stat = EventRate/sqrt( (double) Ndiced); // stat. rate error
831  rateErr_syst = EventRate/3.e-13 * 1.0e-13; // syst. rate error, from error of known flux
832  }
833  else {
834  EventRate= (Ndiced * Norm->norm(n100cos)) * area * selEff;
835  rateErr_stat = EventRate/sqrt( (double) n100cos); // stat. rate error
836  rateErr_syst = EventRate/2.63e-3 * 0.06e-3; // syst. rate error, from error of known flux
837  }
838 
839  // normalisation in region 1.-cos(theta) < 1./(2.*Pi), if MaxTheta even lower correct for this
840  if(MaxTheta<0.572){
841  double spacean = 2.*Pi*(1.-cos(MaxTheta));
842  EventRate= (Ndiced * Norm->norm(n100cos)) * area * selEff * spacean;
843  rateErr_stat = EventRate/sqrt( (double) n100cos); // rate error
844  rateErr_syst = EventRate/2.63e-3 * 0.06e-3; // syst. rate error, from error of known flux
845  }
846 
847  }else{
848  EventRate=Nsel; //no info as no muons at 100 GeV
851  std::cout << std::endl;
852  if (MinP > 100.)
853  std::cout << " !!! MinP > 100 GeV. Cannot apply normalisation!" << std::endl;
854  else if (MaxTheta > 84.26*Deg2Rad)
855  std::cout << " !!! Note: generated cosmics exceed parameterisation. No flux calculated!" << std::endl;
856 
857  else
858  std::cout << " !!! Not enough statistics to apply normalisation (rate=1 +- 1) !!!" << std::endl;
859  }
860 
861  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
862  std::cout << " rate is " << EventRate << " +-" << rateErr_stat <<" (stat) " << "+-" <<
863  rateErr_syst << " (syst) " <<" muons per second" << std::endl;
864  if(EventRate!=0) std::cout << " number of events corresponds to " << Nsel/EventRate << " s" << std::endl; //runtime at CMS = Nsel/rate
865  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
866  std::cout << std::endl;
867  std::cout << "*********************************************************" << std::endl;
868  std::cout << "*********************************************************" << std::endl;
869  }
870 }
871 
873  if (MinP < 0.){ NumberOfEvents = 0;
874  std::cout << " CMG-ERR: min.energy is out of range (0 GeV ... inf]" << std::endl << std::endl; }
875  if (MaxP < 0.){ NumberOfEvents = 0;
876  std::cout << " CMG-ERR: max.energy is out of range (0 GeV ... inf]" << std::endl << std::endl; }
877  if (MaxP <= MinP){ NumberOfEvents = 0;
878  std::cout << " CMG-ERR: max.energy is not greater than min.energy" << std::endl << std::endl; }
879  if (MinTheta < 0.){ NumberOfEvents = 0;
880  std::cout << " CMG-ERR: min.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl; }
881  if (MaxTheta < 0.){ NumberOfEvents = 0;
882  std::cout << " CMG-ERR: max.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl; }
883  if (MaxTheta <= MinTheta){ NumberOfEvents = 0;
884  std::cout << " CMG-ERR: max.theta is not greater than min.theta" << std::endl << std::endl; }
885  if (MinPhi < 0.){ NumberOfEvents = 0;
886  std::cout << " CMG-ERR: min.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl; }
887  if (MaxPhi < 0.){ NumberOfEvents = 0;
888  std::cout << " CMG-ERR: max.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl; }
889  if (MaxPhi <= MinPhi){ NumberOfEvents = 0;
890  std::cout << " CMG-ERR: max.phi is not greater than min.phi" << std::endl << std::endl; }
891  if (MaxT0 <= MinT0){ NumberOfEvents = 0;
892  std::cout << " CMG-ERR: max.t0 is not greater than min.t0" << std::endl << std::endl; }
893  if (ElossScaleFactor < 0.){ NumberOfEvents = 0;
894  std::cout << " CMG-ERR: E-loss scale factor is out of range [0 ... inf)" << std::endl << std::endl; }
895  if (MinEnu < 0.){ NumberOfEvents = 0;
896  std::cout << " CMG-ERR: min.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl; }
897  if (MaxEnu < 0.){ NumberOfEvents = 0;
898  std::cout << " CMG-ERR: max.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl; }
899  if (MaxEnu <= MinEnu){ NumberOfEvents = 0;
900  std::cout << " CMG-ERR: max.Enu is not greater than min.Enu" << std::endl << std::endl; }
901 
902 }
903 
905  // check angular range (for a sphere with Target3dRadius around the target)
906  bool goodAngles = false;
907  bool phiaccepted = false;
908  bool thetaaccepted = false;
909  double RxzV = sqrt(OneMuoEvt.vx()*OneMuoEvt.vx() + OneMuoEvt.vz()*OneMuoEvt.vz());
910 
911  double rVY;
912  if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
913  rVY = -sqrt(RxzV*RxzV + RadiusCMS*RadiusCMS);
914  else
915  rVY = sqrt(RxzV*RxzV + (SurfaceOfEarth+PlugWidth)*(SurfaceOfEarth+PlugWidth));
916 
917  double Phi = OneMuoEvt.phi();
918  double PhiV = atan2(OneMuoEvt.vx(),OneMuoEvt.vz()) + Pi; if (PhiV > TwoPi) PhiV -= TwoPi;
919  double disPhi = std::fabs(PhiV - Phi); if (disPhi > Pi) disPhi = TwoPi - disPhi;
920  double dPhi = Pi; if (RxzV > Target3dRadius) dPhi = asin(Target3dRadius/RxzV);
921  if (disPhi < dPhi) phiaccepted = true;
922  double Theta = OneMuoEvt.theta();
923  double ThetaV = asin(RxzV/rVY);
924  double dTheta = Pi; if (std::fabs(rVY) > Target3dRadius) dTheta = asin(Target3dRadius/std::fabs(rVY));
925  //std::cout << " dPhi = " << dPhi << " (" << Phi << " <p|V> " << PhiV << ")" << std::endl;
926  //std::cout << " dTheta = " << dTheta << " (" << Theta << " <p|V> " << ThetaV << ")" << std::endl;
927 
928  if (!phiaccepted && RxzV < Target3dRadius)
929  //if (RxzV < Target3dRadius)
930  std::cout << "Rejected phi=" << Phi << " PhiV=" << PhiV
931  << " dPhi=" << dPhi << " disPhi=" << disPhi
932  << " RxzV=" << RxzV << " Target3dRadius=" << Target3dRadius
933  << " Theta=" << Theta << std::endl;
934 
935  if (std::fabs(Theta-ThetaV) < dTheta) thetaaccepted = true;
936  if (phiaccepted && thetaaccepted) goodAngles = true;
937  return goodAngles;
938 }
939 
941 #if ROOT_INTERACTIVE
942  float rCMS = RadiusCMS/1000.;
943  float zCMS = Z_DistCMS/1000.;
944  if(TrackerOnly==true){
945  rCMS = RadiusTracker/1000.;
946  zCMS = Z_DistTracker/1000.;
947 }
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);
954  disC->Divide(2,1);
955  disC->cd(1);
956  gPad->SetTicks(1,1);
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]");
967  disC->cd(2);
968  gPad->SetGrid(1,1);
969  gPad->SetTicks(1,1);
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]");
980 #endif
981 }
982 
984 #if ROOT_INTERACTIVE
985  double RadiusDet=RadiusCMS;
986  double Z_DistDet=Z_DistCMS;
987  if(TrackerOnly==true){
988  RadiusDet = RadiusTracker;
989  Z_DistDet = Z_DistTracker;
990  }
991  disXY->Reset();
992  disZY->Reset();
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);
996  float energy = float(OneMuoEvt.e());
997  float verX = float(OneMuoEvt.vx()/1000.); // [m]
998  float verY = float(OneMuoEvt.vy()/1000.); // [m]
999  float verZ = float(OneMuoEvt.vz()/1000.); // [m]
1000  float dirX = float(OneMuoEvt.px())/std::fabs(OneMuoEvt.py());
1001  float dirY = float(OneMuoEvt.py())/std::fabs(OneMuoEvt.py());
1002  float dirZ = float(OneMuoEvt.pz())/std::fabs(OneMuoEvt.py());
1003  float yStep = disXY->GetYaxis()->GetBinWidth(1);
1004  int NbinY = disXY->GetYaxis()->GetNbins();
1005  for (int iy=0; iy<NbinY; ++iy){
1006  verX += dirX*yStep;
1007  verY += dirY*yStep;
1008  verZ += dirZ*yStep;
1009  float rXY = sqrt(verX*verX + verY*verY)*1000.; // [mm]
1010  float absZ = std::fabs(verZ)*1000.; // [mm]
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");
1017  gPad->Update();
1018  }
1019  }
1020 #endif
1021 }
1022 
1024 
1026 
1028 
1030 
1032 
1034 
1036 
1038 
1040 
1042 
1044 
1045 void CosmicMuonGenerator::setElossScaleFactor(double ElossScaleFact){ if (NotInitialized) ElossScaleFactor = ElossScaleFact; }
1046 
1048 
1050 
1052 
1054 
1055 void CosmicMuonGenerator::setMultiMuon(bool MultiMu){ if (NotInitialized) MultiMuon = MultiMu; }
1056 void CosmicMuonGenerator::setMultiMuonFileName(std::string MultiMuFile){ if (NotInitialized) MultiMuonFileName = MultiMuFile; }
1057 void CosmicMuonGenerator::setMultiMuonFileFirstEvent(int MultiMuFile1stEvt){ if (NotInitialized) MultiMuonFileFirstEvent = MultiMuFile1stEvt; }
1058 void CosmicMuonGenerator::setMultiMuonNmin(int MultiMuNmin){ if (NotInitialized) MultiMuonNmin = MultiMuNmin; }
1059 
1061 
1064 
1065 void CosmicMuonGenerator::setPlugVx(double PlugVtx){ if (NotInitialized) PlugVx = PlugVtx; }
1066 void CosmicMuonGenerator::setPlugVz(double PlugVtz){ if (NotInitialized) PlugVz = PlugVtz; }
1067 void CosmicMuonGenerator::setRhoAir(double VarRhoAir){ if (NotInitialized) RhoAir = VarRhoAir; }
1068 void CosmicMuonGenerator::setRhoWall(double VarRhoWall){ if (NotInitialized) RhoWall = VarRhoWall; }
1069 void CosmicMuonGenerator::setRhoRock(double VarRhoRock){ if (NotInitialized) RhoRock = VarRhoRock; }
1070 void CosmicMuonGenerator::setRhoClay(double VarRhoClay){ if (NotInitialized) RhoClay = VarRhoClay; }
1071 void CosmicMuonGenerator::setRhoPlug(double VarRhoPlug){ if (NotInitialized) RhoPlug = VarRhoPlug; }
1072 void CosmicMuonGenerator::setClayWidth(double ClayLayerWidth){ if (NotInitialized) ClayWidth = ClayLayerWidth; }
1073 
1074 void CosmicMuonGenerator::setMinEnu(double MinEn){ if (NotInitialized) MinEnu = MinEn; }
1075 void CosmicMuonGenerator::setMaxEnu(double MaxEn){ if (NotInitialized) MaxEnu = MaxEn; }
1076 void CosmicMuonGenerator::setNuProdAlt(double NuPrdAlt){ if (NotInitialized) NuProdAlt = NuPrdAlt; }
1077 
1079 
const double Z[kNumberCalorimeter]
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)
void initialize(CLHEP::HepRandomEngine *rng=0)
float norm(int n100cos)
Definition: CMSCGENnorm.cc:31
std::vector< double > Theta_mu
int events_n100cos(double energy, double theta)
Definition: CMSCGENnorm.cc:15
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:23
void setRadiusOfTarget(double R)
Double_t particle__Time[kMaxparticle_]
Definition: sim.h:66
void setNumberOfEvents(unsigned int N)
#define P
std::vector< double > Py_sf
#define min(a, b)
Definition: mlp_lapack.h:161
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
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:462
const double PlugWidth
Float_t shower_Energy
Definition: sim.h:29
std::vector< double > Vz_ug
int generate()
Definition: CMSCGEN.cc:256
void setMinPhi(double Phi)
const double RadiusCMS
const double SurfaceOfEarth
void setMaxPhi(double Phi)
std::vector< double > Vz_sf
std::vector< double > Vz_mu
const double Deg2Rad
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
const T & max(const T &a, const T &b)
T sqrt(T t)
Definition: SSEVec.h:28
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
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:447
void setPlugVz(double PlugVtz)
void setTIFOnly_linear(bool TIF)
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)
tuple cout
Definition: gather_cfg.py:41
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:495
std::vector< double > Vx_ug
int generateNuMu()
Definition: CMSCGEN.cc:588
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:40
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)
Definition: DDAxes.h:10