Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:28 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!