CosmicMuonGenerator.cc

Go to the documentation of this file.

// modified by P. Biallass 29.03.2006 to implement new cosmic generator (CMSCGEN.cc) and new normalization of flux (CMSCGENnorm.cc)
// 04.12.2008 sonne: replaced Min/MaxE by Min/MaxP to get cos_sf/ug scripts working again
// 20.04.2009 sonne: Implemented mechanism to read in multi muon events and propagate each muon
#define sim_cxx


#include "GeneratorInterface/CosmicMuonGenerator/interface/CosmicMuonGenerator.h"


void CosmicMuonGenerator::runCMG(){
  initialize();
  for (unsigned int iGen=0; iGen<NumberOfEvents; ++iGen){ nextEvent(); }
  terminate();
}

void CosmicMuonGenerator::setRandomEngine(CLHEP::HepRandomEngine* v) {
  if (delRanGen)
    delete RanGen;
  RanGen = v;
  delRanGen = false;
  Cosmics->setRandomEngine(v);
}

void CosmicMuonGenerator::initialize(CLHEP::HepRandomEngine *rng){
  if (delRanGen)
    delete RanGen;
  if (!rng) {
    RanGen = new CLHEP::HepJamesRandom;
    RanGen->setSeed(RanSeed, 0); //set seed for Random Generator (seed can be controled by config-file)
    delRanGen = true;
  } else {
    RanGen = rng;
    delRanGen = false;
  }
  checkIn();
  if (NumberOfEvents > 0){
    // set up "surface geometry" dimensions
    double RadiusTargetEff = RadiusOfTarget; //get this from cfg-file
    double Z_DistTargetEff = ZDistOfTarget;  //get this from cfg-file
    //double Z_CentrTargetEff = ZCentrOfTarget;  //get this from cfg-file
    if(TrackerOnly==true){
    RadiusTargetEff = RadiusTracker;
    Z_DistTargetEff = Z_DistTracker;
    }
    Target3dRadius = sqrt(RadiusTargetEff*RadiusTargetEff + Z_DistTargetEff*Z_DistTargetEff) + MinStepSize;
    if (Debug) std::cout << "  radius of sphere  around  target = " << Target3dRadius << " mm" << std::endl;

    if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
      SurfaceRadius = (RadiusCMS)*(-tan(MinTheta)) + MinStepSize;
    else
      SurfaceRadius = (SurfaceOfEarth+PlugWidth+RadiusTargetEff)*tan(MaxTheta) + Target3dRadius;  
    if (Debug) std::cout << "  starting point radius at Surface + PlugWidth = " << SurfaceRadius << " mm" << std::endl;

    OneMuoEvt.PlugVx = PlugVx;
    OneMuoEvt.PlugVz = PlugVz;
    OneMuoEvt.RhoAir = RhoAir;
    OneMuoEvt.RhoWall = RhoWall;
    OneMuoEvt.RhoRock = RhoRock;
    OneMuoEvt.RhoClay = RhoClay;
    OneMuoEvt.RhoPlug = RhoPlug;
    //set energy and angle limits for CMSCGEN, give same seed as above 
    if (MinTheta >= 90.*Deg2Rad) //upgoing muons from neutrinos
      Cosmics->initializeNuMu(MinP, MaxP, MinTheta, MaxTheta, MinEnu, MaxEnu, 
                  MinPhi, MaxPhi, NuProdAlt, RanGen);
    else 
      Cosmics->initialize(MinP, MaxP, MinTheta, MaxTheta, RanGen, TIFOnly_constant, TIFOnly_linear);
   
#if ROOT_INTERACTIVE
  // book histos
  TH1D* ene = new TH1D("ene","generated energy",210,0.,1050.);
  TH1D* the = new TH1D("the","generated theta",90,0.,90.);
  TH1D* phi = new TH1D("phi","generated phi",120,0.,360.);
  TH3F* ver = new TH3F("ver","Z-X-Y coordinates",100,-25.,25.,20,-10.,10.,40,-10.,10.);
#endif
    if (EventDisplay) initEvDis();
    std::cout << std::endl;
    
    if (MultiMuon) {
      MultiIn = nullptr;
      
      std::cout << "MultiMuonFileName.c_str()=" << MultiMuonFileName.c_str() << std::endl;
      MultiIn = new TFile( MultiMuonFileName.c_str() );
      
      if (!MultiIn) std::cout << "MultiMuon=True: MultiMuonFileName='"
                              << MultiMuonFileName.c_str() << "' does not exist" << std::endl;
      else std::cout << "MultiMuonFile: " << MultiMuonFileName.c_str() << " opened!" << std::endl;
      //MultiTree = (TTree*) gDirectory->Get("sim");
      MultiTree = (TTree*) MultiIn->Get("sim");
      SimTree = new sim(MultiTree);
      SimTree->Init(MultiTree);
      SimTreeEntries = SimTree->fChain->GetEntriesFast();
      std::cout << "SimTreeEntries=" << SimTreeEntries << std::endl;

      if (MultiMuonFileFirstEvent <= 0)
    SimTree_jentry = 0;
      else 
    SimTree_jentry = MultiMuonFileFirstEvent - 1; //1=1st evt (SimTree_jentry=0)
      
      NcloseMultiMuonEvents = 0;
      NskippedMultiMuonEvents = 0;
    }
    
    if (!MultiMuon || (MultiMuon && MultiIn)) NotInitialized = false;
     
  }
}

void CosmicMuonGenerator::nextEvent(){

  double E = 0.; double Theta = 0.; double Phi = 0.; double RxzV = 0.; double PhiV = 0.;
  if (int(Nsel)%100 == 0) std::cout << "    generated " << int(Nsel) << " events" << std::endl;
  // generate cosmic (E,theta,phi)
  bool   notSelected = true;
  while (notSelected){
    bool   badMomentumGenerated = true;
    while (badMomentumGenerated){
      
      if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
    Cosmics->generateNuMu();
      else
    Cosmics->generate(); //dice one event now
      
      E = sqrt(Cosmics->momentum_times_charge()*Cosmics->momentum_times_charge() + MuonMass*MuonMass);
      Theta = TMath::ACos( Cosmics->cos_theta() ) ; //angle has to be in RAD here
      Ngen+=1.;   //count number of initial cosmic events (in surface area), vertices will be added later
      badMomentumGenerated = false;
      Phi = RanGen->flat()*(MaxPhi-MinPhi) + MinPhi;
    }
    Norm->events_n100cos(E, Theta); //test if this muon is in normalization range
    Ndiced += 1; //one more cosmic is diced

    // generate vertex
    double Nver = 0.;
    bool   badVertexGenerated = true;
    while (badVertexGenerated){
      RxzV = sqrt(RanGen->flat())*SurfaceRadius;
      PhiV = RanGen->flat()*TwoPi;
      // check phi range (for a sphere with Target3dRadius around the target)
      double dPhi = Pi; if (RxzV > Target3dRadius) dPhi = asin(Target3dRadius/RxzV);
      double rotPhi = PhiV + Pi; if (rotPhi > TwoPi) rotPhi -= TwoPi;
      double disPhi = std::fabs(rotPhi - Phi); if (disPhi > Pi) disPhi = TwoPi - disPhi;
      if (disPhi < dPhi || AcptAllMu) badVertexGenerated = false;
      Nver+=1.;
    }
    Ngen += (Nver-1.); //add number of generated vertices to initial cosmic events
    
    // complete event at surface
    int id =  13; // mu-
    if (Cosmics->momentum_times_charge() >0.) id = -13; // mu+
    double absMom = sqrt(E*E - MuonMass*MuonMass);
    double verMom = absMom*cos(Theta);
    double horMom = absMom*sin(Theta);
    double Px = horMom*sin(Phi); // [GeV/c]
    double Py = -verMom;         // [GeV/c]
    double Pz = horMom*cos(Phi); // [GeV/c]
    double Vx = RxzV*sin(PhiV);  // [mm]

    double Vy;
    if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
      Vy = -RadiusCMS;
    else
      Vy = SurfaceOfEarth + PlugWidth;  // [mm]

    double Vz = RxzV*cos(PhiV);  // [mm]
    double T0 = (RanGen->flat()*(MaxT0-MinT0) + MinT0)*SpeedOfLight; // [mm/c];

    Id_at = id;
    Px_at = Px; Py_at = Py; Pz_at = Pz; E_at = E; //M_at = MuonMass;
    Vx_at = Vx; Vy_at = Vy; Vz_at = Vz; T0_at = T0;

    OneMuoEvt.create(id, Px, Py, Pz, E, MuonMass, Vx, Vy, Vz, T0); 
    // if angles are ok, propagate to target
    if (goodOrientation()) { 
      if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
    OneMuoEvt.propagate(0., RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
      else
    OneMuoEvt.propagate(ElossScaleFactor, RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
    }

    if ( (OneMuoEvt.hitTarget() && sqrt(OneMuoEvt.e()*OneMuoEvt.e() - MuonMass*MuonMass) > MinP_CMS)
    || AcptAllMu==true){
      Nsel+=1.; //count number of generated and accepted events  
      notSelected = false;
    }
  }

  EventWeight = 1.;

  //just one outgoing particle at SurFace
  Id_sf.resize(1); 
  Px_sf.resize(1); 
  Py_sf.resize(1); 
  Pz_sf.resize(1);
  E_sf.resize(1); 
  //M_sf.resize(1);
  Vx_sf.resize(1);
  Vy_sf.resize(1);
  Vz_sf.resize(1);
  T0_sf.resize(1);

  Id_sf[0] = Id_at;
  Px_sf[0] = Px_at; Py_sf[0] = Py_at; Pz_sf[0] = Pz_at; E_sf[0] = E_at; //M_fs[0] = MuonMass;
  Vx_sf[0] = Vx_at; Vy_sf[0] = Vy_at; Vz_sf[0] = Vz_at; T0_sf[0] = T0_at;
  

  //just one particle at UnderGround  
  Id_ug.resize(1); 
  Px_ug.resize(1); 
  Py_ug.resize(1); 
  Pz_ug.resize(1);
  E_ug.resize(1); 
  //M_ug.resize(1);
  Vx_ug.resize(1);
  Vy_ug.resize(1);
  Vz_ug.resize(1);
  T0_ug.resize(1);

  Id_ug[0] = OneMuoEvt.id();
  Px_ug[0] = OneMuoEvt.px(); 
  Py_ug[0] = OneMuoEvt.py(); 
  Pz_ug[0] = OneMuoEvt.pz();
  E_ug[0] = OneMuoEvt.e(); 
  //M_ug[0] = OneMuoEvt.m();
  Vx_ug[0] = OneMuoEvt.vx(); 
  Vy_ug[0] = OneMuoEvt.vy(); 
  Vz_ug[0] = OneMuoEvt.vz();
  T0_ug[0] = OneMuoEvt.t0();

  // plot variables of selected events
#if ROOT_INTERACTIVE
  ene->Fill(OneMuoEvt.e());
  the->Fill((OneMuoEvt.theta()*Rad2Deg));
  phi->Fill((OneMuoEvt.phi()*Rad2Deg));
  ver->Fill((OneMuoEvt.vz()/1000.),(OneMuoEvt.vx()/1000.),(OneMuoEvt.vy()/1000.));
#endif
  if (Debug){
    std::cout << "new event" << std::endl;
    std::cout << "  Px,Py,Pz,E,m = " << OneMuoEvt.px() << ", " << OneMuoEvt.py() << ", "
         << OneMuoEvt.pz() << ", " << OneMuoEvt.e() << ", " << OneMuoEvt.m() << " GeV" << std::endl;
    std::cout << "  Vx,Vy,Vz,t0  = " << OneMuoEvt.vx() << ", " << OneMuoEvt.vy() << ", " 
         << OneMuoEvt.vz() << ", " << OneMuoEvt.t0() << " mm" << std::endl;
  }
  if (EventDisplay) displayEv();
  
}



bool CosmicMuonGenerator::nextMultiEvent() {

  if (Debug) std::cout << "\nEntered CosmicMuonGenerator::nextMultiEvent()" << std::endl;
  bool EvtRejected = true;
  bool MuInMaxDist = false;
  double MinDist; //[mm] 

  while (EvtRejected) {

    //read in event from SimTree
    //ULong64_t ientry = SimTree->LoadTree(SimTree_jentry);
    Long64_t ientry = SimTree->GetEntry(SimTree_jentry);
    std::cout << "CosmicMuonGenerator::nextMultiEvent(): SimTree_jentry=" << SimTree_jentry
      //<< " ientry=" << ientry 
          << " SimTreeEntries=" << SimTreeEntries << std::endl;
    if (ientry < 0) return false; //stop run
    if (SimTree_jentry < SimTreeEntries) {
      SimTree_jentry++;
    }
    else {
      std::cout << "CosmicMuonGenerator.cc::nextMultiEvent: No more events in file!" << std::endl;
      return false; //stop run
    }
    


    int nmuons = SimTree->shower_nParticlesWritten;
    if (nmuons<MultiMuonNmin) {
      std::cout << "CosmicMuonGenerator.cc: Warning!  Less than " << MultiMuonNmin << 
    " muons in event!" << std::endl;
      std::cout << "trying next event from file" << std::endl;
      NskippedMultiMuonEvents++;
      continue; //EvtRejected while loop: get next event from file
    }



    Px_mu.resize(nmuons); Py_mu.resize(nmuons); Pz_mu.resize(nmuons);
    P_mu.resize(nmuons);

    MinDist = 99999.e9; //[mm] 
    double MuMuDist;
    MuInMaxDist = false;
    //check if at least one muon pair closer than 30m at surface
    int NmuPmin = 0;
    for (int imu=0; imu<nmuons; ++imu) {

      Px_mu[imu] = -SimTree->particle__Px[imu]*sin(NorthCMSzDeltaPhi)
    + SimTree->particle__Py[imu]*cos(NorthCMSzDeltaPhi);
      Pz_mu[imu] = SimTree->particle__Px[imu]*cos(NorthCMSzDeltaPhi)
    + SimTree->particle__Py[imu]*sin(NorthCMSzDeltaPhi);
      Py_mu[imu] = -SimTree->particle__Pz[imu]; //Corsika down going particles defined in -z direction!
      P_mu[imu] = sqrt(Px_mu[imu]*Px_mu[imu] + Py_mu[imu]*Py_mu[imu] + Pz_mu[imu]*Pz_mu[imu]);
      
      if (P_mu[imu] < MinP_CMS  && AcptAllMu==false) continue;
      else if (SimTree->particle__ParticleID[imu] != 5 &&
           SimTree->particle__ParticleID[imu] != 6) continue;
      else NmuPmin++;
      
      for (int jmu=0; jmu<imu; ++jmu) {
    if (P_mu[jmu] < MinP_CMS && AcptAllMu==false) continue;
    if (SimTree->particle__ParticleID[imu] != 5 &&
        SimTree->particle__ParticleID[imu] != 6) continue;
    MuMuDist = sqrt( (SimTree->particle__x[imu]-SimTree->particle__x[jmu])*
             (SimTree->particle__x[imu]-SimTree->particle__x[jmu]) 
             +(SimTree->particle__y[imu]-SimTree->particle__y[jmu])*
             (SimTree->particle__y[imu]-SimTree->particle__y[jmu])
             )*10.; //CORSIKA [cm] to CMSCGEN [mm]
    if (MuMuDist < MinDist) MinDist = MuMuDist;
    if (MuMuDist < 2.*Target3dRadius) MuInMaxDist = true;
      }
    }
    if (MultiMuonNmin>=2) {
      if (MuInMaxDist) { 
    NcloseMultiMuonEvents++;
      }
      else {
    std::cout << "CosmicMuonGenerator.cc: Warning! No muon pair closer than " 
          << 2.*Target3dRadius/1000. << "m   MinDist=" << MinDist/1000. << "m at surface" << std::endl;
    std::cout << "Fraction of too wide opening angle multi muon events: "
          << 1 - double(NcloseMultiMuonEvents)/SimTree_jentry << std::endl;
    std::cout << "NcloseMultiMuonEvents=" << NcloseMultiMuonEvents << std::endl;
    std::cout << "trying next event from file" << std::endl;
    NskippedMultiMuonEvents++;
    continue; //EvtRejected while loop: get next event from file
      }
    }
  
    if (NmuPmin < MultiMuonNmin && AcptAllMu==false) { //take single muon events consistently into account
      NskippedMultiMuonEvents++;
      continue; //EvtRejected while loop: get next event from file
    }

    if (Debug) 
      if (MultiMuonNmin>=2)
    std::cout << "start trial do loop: MuMuDist=" << MinDist/1000. << "[m]   Nmuons=" 
          << nmuons << "  NcloseMultiMuonEvents=" << NcloseMultiMuonEvents 
          << "  NskippedMultiMuonEvents=" << NskippedMultiMuonEvents << std::endl;  
    
    
    //int primary_id = SimTree->run_ParticleID;
    Id_at = SimTree->shower_EventID;
    
    double M_at = 0.;
    //if (Id_at == 13) {
    Id_at = 2212; //convert from Corsika to HepPDT
    M_at = 938.272e-3; //[GeV] mass
    //}
    
    E_at = SimTree->shower_Energy;
    Theta_at = SimTree->shower_Theta;
    double phi_at = SimTree->shower_Phi - NorthCMSzDeltaPhi; //rotate by almost 90 degrees
    if (phi_at < -Pi) phi_at +=TwoPi; //bring into interval (-Pi,Pi]
    else if (phi_at > Pi) phi_at -= TwoPi; 
    double P_at = sqrt(E_at*E_at - M_at*M_at);
    //need to rotate about 90degrees around x->N axis => y<=>z, 
    //then rotate new x-z-plane from x->North to x->LHC centre
    Px_at = P_at*sin(Theta_at)*sin(phi_at);
    Py_at = -P_at*cos(Theta_at);
    Pz_at = P_at*sin(Theta_at)*cos(phi_at);

    //compute maximal theta of secondary muons
    double theta_mu_max = Theta_at;
    double theta_mu_min = Theta_at;

    double phi_rel_min = 0.; //phi_mu_min - phi_at
    double phi_rel_max = 0.; //phi_mu_max - phi_at

    //std::cout << "SimTree->shower_Energy=" << SimTree->shower_Energy <<std::endl;

    Theta_mu.resize(nmuons);
    for (int imu=0; imu<nmuons; ++imu) {
      Theta_mu[imu] = acos(-Py_mu[imu]/P_mu[imu]);
      if (Theta_mu[imu]>theta_mu_max) theta_mu_max = Theta_mu[imu];
      if (Theta_mu[imu]<theta_mu_min) theta_mu_min = Theta_mu[imu];

      double phi_mu = atan2(Px_mu[imu],Pz_mu[imu]); // in (-Pi,Pi]
      double phi_rel = phi_mu - phi_at;
      if (phi_rel < -Pi) phi_rel += TwoPi; //bring into interval (-Pi,Pi] 
      else if (phi_rel > Pi) phi_rel -= TwoPi;
      if (phi_rel < phi_rel_min) phi_rel_min = phi_rel;
      else if (phi_rel > phi_rel_max) phi_rel_max =phi_rel;


    }
    
    
    double h_sf = SurfaceOfEarth + PlugWidth; //[mm]

    double R_at = h_sf*tan(Theta_at);
    
    double JdRxzV_dR_trans = 1.; 
    double JdPhiV_dPhi_trans = 1.; 
    double JdR_trans_sqrt = 1.;
    
    //chose random vertex Phi and Rxz weighted to speed up and smoothen
    double R_mu_max = (h_sf+Target3dRadius)*tan(theta_mu_max);
    double R_max = std::min(SurfaceRadius, R_mu_max);
    double R_mu_min = (h_sf-Target3dRadius)*tan(theta_mu_min);
    double R_min = std::max(0., R_mu_min);
    
    if (R_at>SurfaceRadius) {
      std::cout << "CosmicMuonGenerator.cc: Warning! R_at=" << R_at 
       << " > SurfaceRadius=" << SurfaceRadius <<std::endl;
    }
    
    //do phase space transformation for horizontal radius R
    
    //determine first phase space limits
    
    double psR1min = R_min + 0.25*(R_max-R_min); 
    double psR1max = std::min(SurfaceRadius,R_max-0.25*(R_max-R_min)); //no R's beyond R_max
    double psR1 = psR1max - psR1min;

    double psR2min = R_min;
    double psR2max = R_max;
    double psR2 = psR2max - psR2min;

    double psR3min = 0.;
    double psR3max = SurfaceRadius;
    double psR3 = psR3max - psR3min;

    //double psall = psR1+psR2+psR3;
    double psRall = psR3;

    double fR1=psR1/psRall, fR2=psR2/psRall, fR3=psR3/psRall; //f1+f2+f3=130%
    double pR1=0.25, pR2=0.7, pR3=0.05;


    //do phase space transformation for azimuthal angle phi
    double psPh1 = 0.5*(phi_rel_max - phi_rel_min);
    double psPh2 = phi_rel_max - phi_rel_min;
    double psPh3 = TwoPi;
    double psPhall = psPh3;
    
    double fPh1=psPh1/psPhall, fPh2=psPh2/psPhall, fPh3=psPh3/psPhall; //(f1+f2+f3=TwoPi+f1+f2)/(TwoPi+f1+f2) 

    double pPh1=0.25, pPh2=0.7, pPh3=0.05;

    Trials = 0; //global int trials
    double trials = 0.; //local weighted trials
    Vx_mu.resize(nmuons); Vy_mu.resize(nmuons); Vz_mu.resize(nmuons);
    int NmuHitTarget = 0;
    while (NmuHitTarget < MultiMuonNmin) { 
   
      NmuHitTarget = 0; //re-initialize every loop iteration
      double Nver = 0.;

      
      //chose phase space class
      double RxzV;
      double which_R_class = RanGen->flat();
      if (which_R_class < pR1) { //pR1% in psR1
    RxzV = psR1min + psR1 * RanGen->flat();
    JdRxzV_dR_trans = fR1/pR1 * SurfaceRadius/psR1;
      }
      else if (which_R_class < pR1+pR2) { //further pR2% in psR2 
    RxzV = psR2min + psR2 * RanGen->flat();   
    JdRxzV_dR_trans = fR2/pR2 * SurfaceRadius/psR2;
      }
      else { //remaining pR3% in psR3=[0., R_max]
    RxzV = psR3min + psR3 * RanGen->flat();
    JdRxzV_dR_trans = fR3/pR3 * SurfaceRadius/psR3;
      }
      
      JdR_trans_sqrt = 2.*RxzV/SurfaceRadius; //flat in sqrt(r) space

      //chose phase space class
      double PhiV;
      double which_phi_class = RanGen->flat();
      if (which_phi_class < pPh1) { //pPh1% in psPh1 
    PhiV = phi_at + phi_rel_min + psPh1 * RanGen->flat();
    JdPhiV_dPhi_trans = fPh1/pPh1 * TwoPi/psPh1;
      }
      else if (which_phi_class < pPh1+pPh2) { //further pPh2% in psPh2
    PhiV = phi_at + phi_rel_min + psPh2 * RanGen->flat();     
    JdPhiV_dPhi_trans = fPh2/pPh2 * TwoPi/psPh2;
      }
      else { //remaining pPh3% in psPh3=[-Pi,Pi]
    PhiV = phi_at + phi_rel_min + psPh3 * RanGen->flat();
    JdPhiV_dPhi_trans = fPh3/pPh3 * TwoPi/psPh3;
      }

      //shuffle PhiV into [-Pi,+Pi] interval
      if (PhiV < -Pi) PhiV+=TwoPi;
      else if (PhiV > Pi) PhiV-=TwoPi;
      
      
      Nver++;
      trials += JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans;
      Trials++;
      if (trials > max_Trials) break; //while (Id_sf.size() < 2) loop 
      Ngen += (Nver-1.); //add number of generated vertices to initial cosmic events
         
   
      Vx_at = RxzV*sin(PhiV); // [mm]

      Vy_at = h_sf; // [mm] (SurfaceOfEarth + PlugWidth; Determine primary particle height below)
      //Vy_at = SimTree->shower_StartingAltitude*10. + h_sf; // [mm]
      //std::cout << "SimTree->shower_StartingAltitude*10=" << SimTree->shower_StartingAltitude*10 <<std::endl;
      Vz_at = RxzV*cos(PhiV); // [mm]
      
      int NmuHitTargetSphere = 0;
      for (int imu=0; imu<nmuons; ++imu) {
    
    Vx_mu[imu] = Vx_at + (-SimTree->particle__x[imu]*sin(NorthCMSzDeltaPhi)
                  +SimTree->particle__y[imu]*cos(NorthCMSzDeltaPhi) )*10; //[mm] (Corsika cm to CMSCGEN mm)
    Vy_mu[imu] = h_sf; //[mm] fixed at surface + PlugWidth
    Vz_mu[imu] = Vz_at + ( SimTree->particle__x[imu]*cos(NorthCMSzDeltaPhi)
                   +SimTree->particle__y[imu]*sin(NorthCMSzDeltaPhi) )*10; //[mm] (Corsika cm to CMSCGEN mm)

    
    //add atmospheric height to primary particle (default SimTree->shower_StartingAltitude = 0.)
    double pt_sec = sqrt(Px_mu[imu]*Px_mu[imu]+Pz_mu[imu]*Pz_mu[imu]);
    double theta_sec = atan2(std::fabs(Py_mu[imu]),pt_sec);
    double r_sec = sqrt((Vx_mu[imu]-Vx_at)*(Vx_mu[imu]-Vx_at)
                +(Vz_mu[imu]-Vz_at)*(Vz_mu[imu]-Vz_at));
    double h_prod =  r_sec * tan(theta_sec);
    if (h_prod + h_sf > Vy_at) Vy_at = h_prod + h_sf;

    //only muons
    if (SimTree->particle__ParticleID[imu] != 5 &&
        SimTree->particle__ParticleID[imu] != 6) continue;

    if (P_mu[imu] < MinP_CMS && AcptAllMu==false) continue; 

    //check here if at least 2 muons make it to the target sphere
    double Vxz_mu = sqrt(Vx_mu[imu]*Vx_mu[imu] + Vz_mu[imu]*Vz_mu[imu]);
    theta_mu_max = atan((Vxz_mu+Target3dRadius)/(h_sf-Target3dRadius));
    theta_mu_min = atan((Vxz_mu-Target3dRadius)/(h_sf+Target3dRadius));
    if (Theta_mu[imu] > theta_mu_min && Theta_mu[imu] < theta_mu_max) {
      
      // check phi range (for a sphere with Target3dRadius around the target)
      double dPhi = Pi; if (Vxz_mu > Target3dRadius) dPhi = asin(Target3dRadius/Vxz_mu);
      double PhiPmu = atan2(Px_mu[imu],Pz_mu[imu]); //muon phi
      double PhiVmu = atan2(Vx_mu[imu],Vz_mu[imu]); //muon phi
      double rotPhi = PhiVmu + Pi; if (rotPhi > Pi) rotPhi -= TwoPi;
      double disPhi = std::fabs(rotPhi - PhiPmu); if (disPhi > Pi) disPhi = TwoPi - disPhi;
      if (disPhi < dPhi) {
        NmuHitTargetSphere++;
      }

    }

      } //end imu for loop

      


      if (NmuHitTargetSphere < MultiMuonNmin) continue; //while (Id_sf.size() < 2) loop
    
      //nmuons outgoing particle at SurFace
      Id_sf.clear(); 
      Px_sf.clear(); 
      Py_sf.clear(); 
      Pz_sf.clear();
      E_sf.clear(); 
      //M_sf_out.clear();
      Vx_sf.clear(); 
      Vy_sf.clear(); 
      Vz_sf.clear();
      T0_sf.clear();
      
      //nmuons particles at UnderGround  
      Id_ug.clear(); 
      Px_ug.clear(); 
      Py_ug.clear(); 
      Pz_ug.clear();
      E_ug.clear(); 
      //M_ug.clear();
      Vx_ug.clear();
      Vy_ug.clear();
      Vz_ug.clear();
      T0_ug.clear();
      
      int Id_sf_this =0;
      double Px_sf_this =0., Py_sf_this=0., Pz_sf_this=0.;
      double E_sf_this=0.;
      //double M_sf_this=0.;
      double Vx_sf_this=0., Vy_sf_this=0., Vz_sf_this=0.;
      double T0_sf_this=0.;
      
      T0_at = SimTree->shower_GH_t0 * SpeedOfLight; // [mm]

      for (int imu=0; imu<nmuons; ++imu) {

    if (P_mu[imu] < MinP_CMS && AcptAllMu==false) continue;
    //for the time being only muons
    if (SimTree->particle__ParticleID[imu] != 5 &&
        SimTree->particle__ParticleID[imu] != 6) continue;
    
    Id_sf_this = SimTree->particle__ParticleID[imu];
    if (Id_sf_this == 5) Id_sf_this = -13; //mu+
    else if (Id_sf_this == 6) Id_sf_this = 13; //mu-

    else if (Id_sf_this == 1) Id_sf_this = 22; //gamma
    else if (Id_sf_this == 2) Id_sf_this = -11; //e+
    else if (Id_sf_this == 3) Id_sf_this = 11; //e-
    else if (Id_sf_this == 7) Id_sf_this = 111; //pi0
    else if (Id_sf_this == 8) Id_sf_this = 211; //pi+
    else if (Id_sf_this == 9) Id_sf_this = -211; //pi-
    else if (Id_sf_this == 10) Id_sf_this = 130; //KL0
    else if (Id_sf_this == 11) Id_sf_this = 321; //K+
    else if (Id_sf_this == 12) Id_sf_this = -321; //K-
    else if (Id_sf_this == 13) Id_sf_this = 2112; //n
    else if (Id_sf_this == 14) Id_sf_this = 2212; //p
    else if (Id_sf_this == 15) Id_sf_this = -2212; //pbar
    else if (Id_sf_this == 16) Id_sf_this = 310; //Ks0
    else if (Id_sf_this == 17) Id_sf_this = 221; //eta
    else if (Id_sf_this == 18) Id_sf_this = 3122; //Lambda
    
    else {
      std::cout << "CosmicMuonGenerator.cc: Warning! Muon Id=" << Id_sf_this 
           << " from file read in" <<std::endl;
      Id_sf_this = 99999; //trouble
    }
    
    T0_sf_this = SimTree->particle__Time[imu] * SpeedOfLight; //in [mm]

    Px_sf_this = Px_mu[imu];
    Py_sf_this = Py_mu[imu]; //Corsika down going particles defined in -z direction!
    Pz_sf_this = Pz_mu[imu];
    E_sf_this = sqrt(P_mu[imu]*P_mu[imu] + MuonMass*MuonMass);
    Vx_sf_this = Vx_mu[imu];
    Vy_sf_this = Vy_mu[imu]; //[mm] fixed at surface + PlugWidth
    Vz_sf_this = Vz_mu[imu];


    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); 
    // if angles are ok, propagate to target
    if (goodOrientation()) {
      OneMuoEvt.propagate(ElossScaleFactor, RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
    }
    
    if ( (OneMuoEvt.hitTarget() 
          && sqrt(OneMuoEvt.e()*OneMuoEvt.e() - MuonMass*MuonMass) > MinP_CMS)
         || AcptAllMu==true ) {
      
      Id_sf.push_back(Id_sf_this);
      Px_sf.push_back(Px_sf_this);
      Py_sf.push_back(Py_sf_this);
      Pz_sf.push_back(Pz_sf_this);
      E_sf.push_back(E_sf_this);
      //M_sf.push_back(M_sf_this);
      Vx_sf.push_back(Vx_sf_this);
      Vy_sf.push_back(Vy_sf_this);
      Vz_sf.push_back(Vz_sf_this);
      T0_sf.push_back(T0_sf_this);
      //T0_sf.push_back(0.); //synchronised arrival for 100% efficient full simulation tests
      
      Id_ug.push_back(OneMuoEvt.id());
      Px_ug.push_back(OneMuoEvt.px());
      Py_ug.push_back(OneMuoEvt.py());
      Pz_ug.push_back(OneMuoEvt.pz());
      E_ug.push_back(OneMuoEvt.e());
      //M_sf.push_back(OneMuoEvt.m());
      Vx_ug.push_back(OneMuoEvt.vx());
      Vy_ug.push_back(OneMuoEvt.vy());
      Vz_ug.push_back(OneMuoEvt.vz());
      T0_ug.push_back(OneMuoEvt.t0());

      NmuHitTarget++;
    }
      }
      

    } // while (Id_sf.size() < 2); //end of do loop
    

    if (trials > max_Trials) {
      std::cout << "CosmicMuonGenerator.cc: Warning! trials reach max_trials=" << max_Trials
        << " without accepting event!" << std::endl;
      if (Debug) {
    std::cout << " N(mu)=" << Id_ug.size();
    if (!Id_ug.empty()) 
      std::cout << " E[0]=" << E_ug[0] << " theta=" 
            << acos(-Py_ug[0]/sqrt(Px_ug[0]*Px_ug[0]+Py_ug[0]*Py_ug[0]+Pz_ug[0]*Pz_ug[0]))
            << " R_xz=" << sqrt(Vx_sf[0]*Vx_sf[0]+Vy_sf[0]*Vy_sf[0]);
    std::cout << " Theta_at=" << Theta_at << std::endl; 
      }
      std::cout << "Unweighted int num of Trials = " << Trials << std::endl;
      std::cout << "trying next event (" << SimTree_jentry << ") from file" << std::endl;
      NskippedMultiMuonEvents++;
      continue; //EvtRejected while loop: get next event from file
    }
    else {
      if (NmuHitTarget < MultiMuonNmin) {
    std::cout << "CosmicMuonGenerator.cc: Warning! less than " << MultiMuonNmin << 
      " muons hit target: N(mu=)" << NmuHitTarget << std::endl;
    std::cout << "trying next event (" << SimTree_jentry << ") from file" << std::endl;
    NskippedMultiMuonEvents++;
    continue; //EvtRejected while loop: get next event from file
      }
      else { //if (MuInMaxDist) {

    //re-adjust T0's of surviving muons shifted to trigger time box
    //(possible T0 increase due to propagation (loss of energy/speed + travelled distance))
    double T0_ug_min, T0_ug_max;
    T0_ug_min = T0_ug_max = T0_ug[0];
    double Tbox = (MaxT0 - MinT0) * SpeedOfLight; // [mm]
    double minDeltaT0 = 2*Tbox;
    for (unsigned int imu=0; imu<Id_ug.size(); ++imu) {
      double T0_this = T0_ug[imu];
      if (T0_this < T0_ug_min) T0_ug_min = T0_this;
      if (T0_this > T0_ug_max) T0_ug_max = T0_this;
      if (Debug) std::cout << "imu=" << imu << " T0_this=" << T0_this 
                   << " P=" << sqrt(pow(Px_ug[imu],2) + pow(Py_ug[imu],2) + pow(Pz_ug[imu],2)) 
                   << std::endl;
      for (unsigned int jmu=0; jmu<imu; ++jmu) {
        if (std::fabs(T0_ug[imu]-T0_ug[jmu]) < minDeltaT0) minDeltaT0 = std::fabs(T0_ug[imu]-T0_ug[jmu]);
      }
    }
    
    if (int(Id_ug.size()) >= MultiMuonNmin && MultiMuonNmin>=2 && minDeltaT0 > Tbox)
      continue; //EvtRejected while loop: get next event from file

    double T0_min = T0_ug_min +MinT0*SpeedOfLight; //-12.5ns * c [mm]
    double T0_max = T0_ug_max +MaxT0*SpeedOfLight; //+12.5ns * c [mm]

    //ckeck if >= NmuMin in time box, else throw new random number + augment evt weight
    int TboxTrials = 0;
    int NmuInTbox;
    double T0_offset, T0diff;
    for (int tboxtrial=0; tboxtrial<1000; ++tboxtrial) { //max 1000 trials
      T0_offset = RanGen->flat()*(T0_max -T0_min); // [mm]
      TboxTrials++;
      T0diff = T0_offset - T0_max; // [mm] 
      NmuInTbox = 0;
      for (unsigned int imu=0; imu<Id_ug.size(); ++imu) {
        if (T0_ug[imu]+T0diff > MinT0*SpeedOfLight && T0_ug[imu]+T0diff < MaxT0*SpeedOfLight)
          NmuInTbox++;
      }
      if (NmuInTbox >= MultiMuonNmin) break;

    }
    if (NmuInTbox < MultiMuonNmin) continue; //EvtRejected while loop: get next event from file


    if (Debug) std::cout << "initial T0_at=" << T0_at << " T0_min=" << T0_min << " T0_max=" << T0_max 
            << " T0_offset=" << T0_offset;
    T0_at += T0diff; //[mm]
    if (Debug) std::cout << " T0diff=" << T0diff << std::endl;
    for (unsigned int imu=0; imu<Id_ug.size(); ++imu) { //adjust @ surface + underground
      if (Debug) std::cout << "before: T0_sf[" << imu << "]=" << T0_sf[imu] << "  T0_ug=" <<  T0_ug[imu]; 
      T0_sf[imu] += T0diff;
      T0_ug[imu] += T0diff;
      if (Debug) 
        std::cout << "   after: T0_sf[" << imu << "]=" << T0_sf[imu] << "  T0_ug=" <<  T0_ug[imu] << std::endl;
    }   
    if (Debug) std::cout << "T0diff=" << T0diff << " T0_at=" << T0_at << std::endl;



    Nsel += 1;
    EventWeight = JdR_trans_sqrt * JdRxzV_dR_trans * JdPhiV_dPhi_trans 
      / (trials * TboxTrials);
    EvtRejected = false;
    if (Debug) std::cout << "CosmicMuonGenerator.cc: Theta_at=" << Theta_at << " phi_at=" << phi_at 
                 << " Px_at=" << Px_at << " Py_at=" << Py_at << " Pz_at=" << Pz_at 
                 << " T0_at=" << T0_at
                 << " Vx_at=" << Vx_at << " Vy_at=" << Vy_at << " Vz_at=" << Vz_at 
                 << " EventWeight=" << EventWeight << " Nmuons=" << Id_sf.size() << std::endl;
      }
    }
    
    
  } //while loop EvtRejected
  

  return true; //write event to HepMC;

}



void CosmicMuonGenerator::terminate(){
  if (NumberOfEvents > 0){
    std::cout << std::endl;
    std::cout << "*********************************************************" << std::endl;
    std::cout << "*********************************************************" << std::endl;
    std::cout << "***                                                   ***" << std::endl;
    std::cout << "***    C O S M I C   M U O N   S T A T I S T I C S    ***" << std::endl;
    std::cout << "***                                                   ***" << std::endl;
    std::cout << "*********************************************************" << std::endl;
    std::cout << "*********************************************************" << std::endl;
    std::cout << std::endl;  
    std::cout << "       number of initial cosmic events:  " << int(Ngen) << std::endl;
    std::cout << "       number of actually diced events:  " << int(Ndiced) << std::endl;
    std::cout << "       number of generated and accepted events:  " << int(Nsel) << std::endl;
    double selEff = Nsel/Ngen; // selection efficiency
    std::cout << "       event selection efficiency:  " << selEff*100. << "%" << std::endl;
    int n100cos =  Norm->events_n100cos(0., 0.); //get final amount of cosmics in defined range for normalisation of flux
    std::cout << "       events with ~100 GeV and 1 - cos(theta) < 1/2pi: " << n100cos << std::endl;
    std::cout << std::endl;
    std::cout << "       momentum range: " << MinP             << " ... " << MaxP << " GeV" << std::endl;
    std::cout << "       theta  range:   " << MinTheta*Rad2Deg << " ... " << MaxTheta*Rad2Deg << " deg" << std::endl; 
    std::cout << "       phi    range:   " << MinPhi*Rad2Deg   << " ... " << MaxPhi*Rad2Deg << " deg" << std::endl;
    std::cout << "       time   range:   " << MinT0            << " ... " << MaxT0 << " ns" << std::endl;
    std::cout << "       energy  loss:   " << ElossScaleFactor*100. << "%" << std::endl;
    std::cout << std::endl;
    double area = 1.e-6*Pi*SurfaceRadius*SurfaceRadius; // area on surface [m^2] 
    if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos)
      std::cout << "       area of initial cosmics at CMS detector bottom surface:   " << area << " m^2" << std::endl;
    else
      std::cout << "       area of initial cosmics on Surface + PlugWidth:   " << area << " m^2" << std::endl;
    std::cout << "       depth of CMS detector (from Surface):   " << SurfaceOfEarth/1000 << " m" << std::endl;
       
    //at least 100 evts., and
    //downgoing inside theta parametersisation range
    //or upgoing neutrino muons 
    if( (n100cos>0 && MaxTheta<84.26*Deg2Rad) 
       || MinTheta>90.*Deg2Rad) {
      // rate: corrected for area and selection-Eff. and normalized to known flux, integration over solid angle (dOmega) is implicit
      // flux is normalised with respect to known flux of vertical 100GeV muons in area at suface level 
      // rate seen by detector is lower than rate at surface area, so has to be corrected for selection-Eff.
      // normalisation factor has unit [1/s/m^2] 
      // rate = N/time --> normalization factor gives 1/runtime/area 
      // normalization with respect to number of actually diced events (Ndiced)

      if (MinTheta > 90.*Deg2Rad) {//upgoing muons from neutrinos) 
    double Omega = (cos(MinTheta)-cos(MaxTheta)) * (MaxPhi-MinPhi);
    //EventRate = (Ndiced * 3.e-13) * Omega * area*1.e4 * selEff;//area in cm, flux=3.e-13cm^-2s^-1sr^-1
    EventRate = (Ndiced * 3.e-13) * Omega * 4.*RadiusOfTarget*ZDistOfTarget*1.e-2 * selEff;//area in cm, flux=3.e-13cm^-2s^-1sr^-1
    rateErr_stat = EventRate/sqrt( (double) Ndiced);  // stat. rate error 
    rateErr_syst = EventRate/3.e-13 * 1.0e-13;  // syst. rate error, from error of known flux 
      }
      else {
    EventRate= (Ndiced * Norm->norm(n100cos)) * area * selEff;
    rateErr_stat = EventRate/sqrt( (double) n100cos);  // stat. rate error 
    rateErr_syst = EventRate/2.63e-3 * 0.06e-3;  // syst. rate error, from error of known flux 
      }

      // normalisation in region 1.-cos(theta) < 1./(2.*Pi), if MaxTheta even lower correct for this
      if(MaxTheta<0.572){
    double spacean = 2.*Pi*(1.-cos(MaxTheta));
    EventRate= (Ndiced * Norm->norm(n100cos)) * area * selEff * spacean;
    rateErr_stat = EventRate/sqrt( (double) n100cos);  // rate error 
    rateErr_syst = EventRate/2.63e-3 * 0.06e-3;  // syst. rate error, from error of known flux 
      }

    }else{
      EventRate=Nsel; //no info as no muons at 100 GeV
      rateErr_stat =Nsel;
      rateErr_syst =Nsel;
      std::cout << std::endl;
      if (MinP > 100.)
    std::cout << " !!! MinP > 100 GeV. Cannot apply normalisation!" << std::endl;
      else if (MaxTheta > 84.26*Deg2Rad)
    std::cout << " !!! Note: generated cosmics exceed parameterisation. No flux calculated!" << std::endl;

      else
    std::cout << " !!! Not enough statistics to apply normalisation (rate=1 +- 1) !!!" << std::endl;
    } 
    
    std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
    std::cout << "       rate is " << EventRate << " +-" << rateErr_stat <<" (stat) " << "+-" << 
      rateErr_syst << " (syst) " <<" muons per second" << std::endl;
    if(EventRate!=0) std::cout << "       number of events corresponds to " << Nsel/EventRate << " s" << std::endl;  //runtime at CMS = Nsel/rate
    std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
    std::cout << std::endl;
    std::cout << "*********************************************************" << std::endl;
    std::cout << "*********************************************************" << std::endl;
  }
}

void CosmicMuonGenerator::checkIn(){
  if (MinP < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: min.energy is out of range (0 GeV ... inf]" << std::endl << std::endl; }
  if (MaxP < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.energy is out of range (0 GeV ... inf]" << std::endl << std::endl; }
  if (MaxP <= MinP){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.energy is not greater than min.energy" << std::endl << std::endl; }
  if (MinTheta < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: min.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl; }
  if (MaxTheta < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.theta is out of range [0 deg ... 90 deg)" << std::endl << std::endl; }
  if (MaxTheta <= MinTheta){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.theta is not greater than min.theta" << std::endl << std::endl; }
  if (MinPhi < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: min.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl; }
  if (MaxPhi < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.phi is out of range [0 deg ... 360 deg]" << std::endl << std::endl; }
  if (MaxPhi <= MinPhi){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.phi is not greater than min.phi" << std::endl << std::endl; }
  if (MaxT0 <= MinT0){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.t0 is not greater than min.t0" << std::endl << std::endl; }
  if (ElossScaleFactor < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: E-loss scale factor is out of range [0 ... inf)" << std::endl << std::endl; }
  if (MinEnu < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: min.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl; }
  if (MaxEnu < 0.){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.Enu is out of range [0 GeV ... inf]" << std::endl << std::endl; }
  if (MaxEnu <= MinEnu){ NumberOfEvents = 0;
    std::cout << "  CMG-ERR: max.Enu is not greater than min.Enu" << std::endl << std::endl; }

}

bool CosmicMuonGenerator::goodOrientation(){
  // check angular range (for a sphere with Target3dRadius around the target)
  bool goodAngles = false;
  bool phiaccepted = false;
  bool thetaaccepted = false;
  double RxzV = sqrt(OneMuoEvt.vx()*OneMuoEvt.vx() + OneMuoEvt.vz()*OneMuoEvt.vz());

  double rVY;
  if (MinTheta > 90.*Deg2Rad) //upgoing muons from neutrinos
    rVY = -sqrt(RxzV*RxzV + RadiusCMS*RadiusCMS);
  else 
    rVY = sqrt(RxzV*RxzV + (SurfaceOfEarth+PlugWidth)*(SurfaceOfEarth+PlugWidth));

  double Phi = OneMuoEvt.phi();
  double PhiV = atan2(OneMuoEvt.vx(),OneMuoEvt.vz()) + Pi; if (PhiV > TwoPi) PhiV -= TwoPi;
  double disPhi = std::fabs(PhiV - Phi); if (disPhi > Pi) disPhi = TwoPi - disPhi;
  double dPhi = Pi; if (RxzV > Target3dRadius) dPhi = asin(Target3dRadius/RxzV);
  if (disPhi < dPhi) phiaccepted = true;
  double Theta = OneMuoEvt.theta();
  double ThetaV = asin(RxzV/rVY);
  double dTheta = Pi; if (std::fabs(rVY) > Target3dRadius) dTheta = asin(Target3dRadius/std::fabs(rVY));
  //std::cout << "    dPhi = " <<   dPhi << "  (" <<   Phi << " <p|V> " <<   PhiV << ")" << std::endl;
  //std::cout << "  dTheta = " << dTheta << "  (" << Theta << " <p|V> " << ThetaV << ")" << std::endl;

  if (!phiaccepted && RxzV < Target3dRadius)
  //if (RxzV < Target3dRadius)
    std::cout << "Rejected phi=" << Phi << "  PhiV="  << PhiV 
     << "  dPhi=" << dPhi << "  disPhi=" << disPhi
     << "  RxzV=" << RxzV << "  Target3dRadius=" << Target3dRadius 
     << "  Theta=" << Theta << std::endl;

  if (std::fabs(Theta-ThetaV) < dTheta) thetaaccepted = true;
  if (phiaccepted && thetaaccepted) goodAngles = true;
  return goodAngles;
}

void CosmicMuonGenerator::initEvDis(){
#if ROOT_INTERACTIVE
  float rCMS = RadiusCMS/1000.;
  float zCMS = Z_DistCMS/1000.;
  if(TrackerOnly==true){
    rCMS = RadiusTracker/1000.;
    zCMS = Z_DistTracker/1000.;
}
  TH2F* disXY = new TH2F("disXY","X-Y view",160,-rCMS,rCMS,160,-rCMS,rCMS);
  TH2F* disZY = new TH2F("disZY","Z-Y view",150,-zCMS,zCMS,160,-rCMS,rCMS);
  gStyle->SetPalette(1,0);
  gStyle->SetMarkerColor(1);
  gStyle->SetMarkerSize(1.5);
  TCanvas *disC = new TCanvas("disC","Cosmic Muon Event Display",0,0,800,410);
  disC->Divide(2,1);
  disC->cd(1);
  gPad->SetTicks(1,1);
  disXY->SetMinimum(log10(MinP));
  disXY->SetMaximum(log10(MaxP));
  disXY->GetXaxis()->SetLabelSize(0.05);
  disXY->GetXaxis()->SetTitleSize(0.05);
  disXY->GetXaxis()->SetTitleOffset(1.0);
  disXY->GetXaxis()->SetTitle("X [m]");
  disXY->GetYaxis()->SetLabelSize(0.05);
  disXY->GetYaxis()->SetTitleSize(0.05);
  disXY->GetYaxis()->SetTitleOffset(0.8);
  disXY->GetYaxis()->SetTitle("Y [m]");
  disC->cd(2);
  gPad->SetGrid(1,1);
  gPad->SetTicks(1,1);
  disZY->SetMinimum(log10(MinP));
  disZY->SetMaximum(log10(MaxP));
  disZY->GetXaxis()->SetLabelSize(0.05);
  disZY->GetXaxis()->SetTitleSize(0.05);
  disZY->GetXaxis()->SetTitleOffset(1.0);
  disZY->GetXaxis()->SetTitle("Z [m]");
  disZY->GetYaxis()->SetLabelSize(0.05);
  disZY->GetYaxis()->SetTitleSize(0.05);
  disZY->GetYaxis()->SetTitleOffset(0.8);
  disZY->GetYaxis()->SetTitle("Y [m]");
#endif
}

void CosmicMuonGenerator::displayEv(){
#if ROOT_INTERACTIVE
  double RadiusDet=RadiusCMS;
  double Z_DistDet=Z_DistCMS;
  if(TrackerOnly==true){
    RadiusDet = RadiusTracker;
    Z_DistDet = Z_DistTracker;
  }
  disXY->Reset();
  disZY->Reset();
  TMarker* InteractionPoint = new TMarker(0.,0.,2);
  TArc* r8m = new TArc(0.,0.,(RadiusDet/1000.));
  TLatex* logEaxis = new TLatex(); logEaxis->SetTextSize(0.05);
  float energy = float(OneMuoEvt.e());
  float verX = float(OneMuoEvt.vx()/1000.); // [m]
  float verY = float(OneMuoEvt.vy()/1000.); // [m]
  float verZ = float(OneMuoEvt.vz()/1000.); // [m]
  float dirX = float(OneMuoEvt.px())/std::fabs(OneMuoEvt.py());
  float dirY = float(OneMuoEvt.py())/std::fabs(OneMuoEvt.py());
  float dirZ = float(OneMuoEvt.pz())/std::fabs(OneMuoEvt.py());
  float yStep = disXY->GetYaxis()->GetBinWidth(1);
  int   NbinY = disXY->GetYaxis()->GetNbins();
  for (int iy=0; iy<NbinY; ++iy){
    verX += dirX*yStep;
    verY += dirY*yStep;
    verZ += dirZ*yStep;
    float rXY = sqrt(verX*verX + verY*verY)*1000.; // [mm]
    float absZ = std::fabs(verZ)*1000.;                 // [mm]
    if (rXY < RadiusDet && absZ < Z_DistDet){
      disXY->Fill(verX,verY,log10(energy));
      disZY->Fill(verZ,verY,log10(energy));
      disC->cd(1); disXY->Draw("COLZ"); InteractionPoint->Draw("SAME"); r8m->Draw("SAME");
      logEaxis->DrawLatex((0.65*RadiusDet/1000.),(1.08*RadiusDet/1000.),"log_{10}E(#mu^{#pm})");
      disC->cd(2); disZY->Draw("COL"); InteractionPoint->Draw("SAME");
      gPad->Update();
    }
  }
#endif
}

void CosmicMuonGenerator::setNumberOfEvents(unsigned int N){ if (NotInitialized) NumberOfEvents = N; }

void CosmicMuonGenerator::setRanSeed(int N){ if (NotInitialized) RanSeed = N; }

void CosmicMuonGenerator::setMinP(double P){ if (NotInitialized) MinP = P; }

void CosmicMuonGenerator::setMinP_CMS(double P){ if (NotInitialized) MinP_CMS = P; }

void CosmicMuonGenerator::setMaxP(double P){ if (NotInitialized) MaxP = P; }

void CosmicMuonGenerator::setMinTheta(double Theta){ if (NotInitialized) MinTheta = Theta*Deg2Rad; }

void CosmicMuonGenerator::setMaxTheta(double Theta){ if (NotInitialized) MaxTheta = Theta*Deg2Rad; } 

void CosmicMuonGenerator::setMinPhi(double Phi){ if (NotInitialized) MinPhi = Phi*Deg2Rad; }

void CosmicMuonGenerator::setMaxPhi(double Phi){ if (NotInitialized) MaxPhi = Phi*Deg2Rad; }

void CosmicMuonGenerator::setMinT0(double T0){ if (NotInitialized) MinT0 = T0; }

void CosmicMuonGenerator::setMaxT0(double T0){ if (NotInitialized) MaxT0 = T0; }

void CosmicMuonGenerator::setElossScaleFactor(double ElossScaleFact){ if (NotInitialized) ElossScaleFactor = ElossScaleFact; }

void CosmicMuonGenerator::setRadiusOfTarget(double R){ if (NotInitialized) RadiusOfTarget = R; }

void CosmicMuonGenerator::setZDistOfTarget(double Z){ if (NotInitialized) ZDistOfTarget = Z; }

void CosmicMuonGenerator::setZCentrOfTarget(double Z){ if (NotInitialized) ZCentrOfTarget = Z; }

void CosmicMuonGenerator::setTrackerOnly(bool Tracker){ if (NotInitialized) TrackerOnly = Tracker; }

void CosmicMuonGenerator::setMultiMuon(bool MultiMu){ if (NotInitialized) MultiMuon = MultiMu; }
void CosmicMuonGenerator::setMultiMuonFileName(std::string MultiMuFile){ if (NotInitialized) MultiMuonFileName = MultiMuFile; }
void CosmicMuonGenerator::setMultiMuonFileFirstEvent(int MultiMuFile1stEvt){ if (NotInitialized) MultiMuonFileFirstEvent = MultiMuFile1stEvt; }
void CosmicMuonGenerator::setMultiMuonNmin(int MultiMuNmin){ if (NotInitialized) MultiMuonNmin = MultiMuNmin; }

void CosmicMuonGenerator::setTIFOnly_constant(bool TIF){ if (NotInitialized) TIFOnly_constant = TIF; }

void CosmicMuonGenerator::setTIFOnly_linear(bool TIF){ if (NotInitialized) TIFOnly_linear = TIF; }
void CosmicMuonGenerator::setMTCCHalf(bool MTCC){ if (NotInitialized) MTCCHalf = MTCC; }

void CosmicMuonGenerator::setPlugVx(double PlugVtx){ if (NotInitialized) PlugVx = PlugVtx; }
void CosmicMuonGenerator::setPlugVz(double PlugVtz){ if (NotInitialized) PlugVz = PlugVtz; }
void CosmicMuonGenerator::setRhoAir(double VarRhoAir){ if (NotInitialized) RhoAir = VarRhoAir; }
void CosmicMuonGenerator::setRhoWall(double VarRhoWall){ if (NotInitialized) RhoWall = VarRhoWall; }
void CosmicMuonGenerator::setRhoRock(double VarRhoRock){ if (NotInitialized) RhoRock = VarRhoRock; }
void CosmicMuonGenerator::setRhoClay(double VarRhoClay){ if (NotInitialized) RhoClay = VarRhoClay; }
void CosmicMuonGenerator::setRhoPlug(double VarRhoPlug){ if (NotInitialized) RhoPlug = VarRhoPlug; }
void CosmicMuonGenerator::setClayWidth(double ClayLayerWidth){ if (NotInitialized) ClayWidth = ClayLayerWidth; }

void CosmicMuonGenerator::setMinEnu(double MinEn){ if (NotInitialized) MinEnu = MinEn; }
void CosmicMuonGenerator::setMaxEnu(double MaxEn){ if (NotInitialized) MaxEnu = MaxEn; }
void CosmicMuonGenerator::setNuProdAlt(double NuPrdAlt){ if (NotInitialized) NuProdAlt = NuPrdAlt; }

double CosmicMuonGenerator::getRate(){ return EventRate; }

void CosmicMuonGenerator::setAcptAllMu(bool AllMu){ if (NotInitialized) AcptAllMu = AllMu; }