Skip to content

Commit

Permalink
Merge pull request #51 from soleti/livetime_dataproduct
Browse files Browse the repository at this point in the history
Added CosmicLivetime data product; Clean-up of CORSIKA code
  • Loading branch information
kutschke authored Nov 21, 2019
2 parents cc2da20 + 193a1bf commit a6d18b2
Show file tree
Hide file tree
Showing 10 changed files with 331 additions and 104 deletions.
31 changes: 31 additions & 0 deletions Analyses/src/CORSIKAGenPlots_module.cc
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include "GlobalConstantsService/inc/PhysicsParams.hh"
#include "GlobalConstantsService/inc/ParticleDataTable.hh"
#include "MCDataProducts/inc/GenParticleCollection.hh"
#include "Mu2eUtilities/inc/TwoLinePCA.hh"

#include "TH1F.h"
#include "TH2F.h"
Expand Down Expand Up @@ -74,6 +75,7 @@ class mu2e::CORSIKAGenPlots : public art::EDAnalyzer {
TH1F *_hNSecondaries = nullptr;

TTree *_cosmicTree = nullptr;
TTree *_eventTree = nullptr;

float _x = std::numeric_limits<float>::lowest();
float _y = std::numeric_limits<float>::lowest();
Expand All @@ -86,6 +88,7 @@ class mu2e::CORSIKAGenPlots : public art::EDAnalyzer {
float _KE = std::numeric_limits<float>::lowest();
float _p = std::numeric_limits<float>::lowest();
float _t = std::numeric_limits<float>::lowest();
std::vector<float> _dca;
PDGCode::type _pdgId;

void bookHists(art::ServiceHandle<art::TFileService> &);
Expand Down Expand Up @@ -118,6 +121,8 @@ mu2e::CORSIKAGenPlots::CORSIKAGenPlots(fhicl::ParameterSet const &p)
_cosmicTree->Branch("t", &_t, "t/F");
_cosmicTree->Branch("pdgId", &_pdgId, "pdgId/I");

_eventTree = tfs->make<TTree>("eventTree", "TTree with cosmic ray info per shower");
_eventTree->Branch("dca", &_dca);
}


Expand All @@ -137,8 +142,34 @@ void mu2e::CORSIKAGenPlots::analyze(art::Event const & e)

const auto & particles = *gpHandle;

// Store the point of closest approach between target box and roof for
// events with more than one particle
_dca.clear();
for (GenParticleCollection::const_iterator i = particles.begin(); i != particles.end(); ++i)
{
for (GenParticleCollection::const_iterator j = i+1; j != particles.end(); ++j)
{
GenParticle const &particle1 = *i;
GenParticle const &particle2 = *j;

TwoLinePCA twoLine(particle1.position(), -particle1.momentum().vect(), particle2.position(), -particle2.momentum().vect());
const CLHEP::Hep3Vector point1 = twoLine.point1();
const CLHEP::Hep3Vector point2 = twoLine.point2();


if (point1.y() > 5000 &&
point1.y() < 15365.4 &&
point2.y() > 5000 &&
point2.y() < 15365.4)
{
_dca.push_back(twoLine.dca());
}
}
}
_eventTree->Fill();

_hNSecondaries->Fill(particles.size());

for(const auto & p : particles)
{
_hXZ->Fill(p.position().x(), p.position().z());
Expand Down
130 changes: 99 additions & 31 deletions EventGenerator/src/CORSIKAEventGenerator_module.cc
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
#include "GeometryService/inc/DetectorSystem.hh"
#include "GeometryService/inc/Mu2eEnvelope.hh"
#include "MCDataProducts/inc/GenParticle.hh"
#include "MCDataProducts/inc/CosmicLivetime.hh"
#include "MCDataProducts/inc/GenParticleCollection.hh"
#include "CalorimeterGeom/inc/Calorimeter.hh"
#include "ExtinctionMonitorFNAL/Geometry/inc/ExtMonFNAL.hh"
Expand All @@ -41,7 +42,7 @@
#include "CLHEP/Random/RandFlat.h"

#include "Mu2eUtilities/inc/VectorVolume.hh"

#include "Mu2eUtilities/inc/TwoLinePCA.hh"

using CLHEP::Hep3Vector;
using CLHEP::HepLorentzVector;
Expand All @@ -60,23 +61,19 @@ namespace mu2e {
fhicl::Atom<std::string> corsikaModuleLabel{Name("corsikaModuleLabel"), Comment("Reference point in the coordinate system"), "FromCorsikaBinary"};
fhicl::Atom<std::string> refPointChoice{Name("refPointChoice"), Comment("Reference point in the coordinate system"), "UNDEFINED"};
fhicl::Atom<bool> projectToTargetBox{Name("projectToTargetBox"), Comment("Store only events that cross the target box"), false};
fhicl::Atom<float> targetBoxXmin{Name("targetBoxXmin"), Comment("Extension of the generation plane"), -5000};
fhicl::Atom<float> targetBoxXmax{Name("targetBoxXmax"), Comment("Extension of the generation plane"), 5000};
fhicl::Atom<float> targetBoxYmin{Name("targetBoxYmin"), Comment("Extension of the generation plane"), -5000};
fhicl::Atom<float> targetBoxYmax{Name("targetBoxYmax"), Comment("Extension of the generation plane"), 5000};
fhicl::Atom<float> targetBoxZmin{Name("targetBoxZmin"), Comment("Extension of the generation plane"), -5000};
fhicl::Atom<float> targetBoxZmax{Name("targetBoxZmax"), Comment("Extension of the generation plane"), 5000};
fhicl::Atom<float> targetBoxYmax{Name("targetBoxYmax"), Comment("Top coordinate of the target box on the Y axis")};
fhicl::Atom<float> intDist{Name("intDist"), Comment("Maximum distance of closest approach for backtracked trajectories")};
};
typedef art::EDProducer::Table<Config> Parameters;

explicit CorsikaEventGenerator(const Parameters &conf);
// Accept compiler written d'tor. Modules are never moved or copied.
virtual void produce (art::Event& e);
virtual void endSubRun(art::SubRun &sr);
virtual void beginSubRun(art::SubRun &sr);

private:

std::vector<CLHEP::Hep3Vector> _targetBoxIntersections;
std::vector<CLHEP::Hep3Vector> _worldIntersections;

bool _geomInfoObtained = false;
Expand All @@ -89,36 +86,45 @@ namespace mu2e {
float _worldYmax = 0;
float _worldZmin = 0;
float _worldZmax = 0;
float _targetBoxXmin = 0;
float _targetBoxXmax = 0;
float _targetBoxYmin = 0;
float _targetBoxYmax = 0;
float _targetBoxZmin = 0;
float _targetBoxZmax = 0;

float _targetBoxYmax;

bool _projectToTargetBox = false;
Hep3Vector _cosmicReferencePointInMu2e;
std::string _refPointChoice;

float _intDist = 0;

unsigned int _primaries = 0;
float _area = 0;
float _lowE = 0;
float _highE = 0;
float _fluxConstant = 0;

};

CorsikaEventGenerator::CorsikaEventGenerator(const Parameters &conf) : EDProducer{conf},
_conf(conf()),
_corsikaModuleLabel(_conf.corsikaModuleLabel()),
_targetBoxXmin(_conf.targetBoxXmin()),
_targetBoxXmax(_conf.targetBoxXmax()),
_targetBoxYmin(_conf.targetBoxYmin()),
_targetBoxYmax(_conf.targetBoxYmax()),
_targetBoxZmin(_conf.targetBoxZmin()),
_targetBoxZmax(_conf.targetBoxZmax()),
_projectToTargetBox(_conf.projectToTargetBox()),
_refPointChoice(_conf.refPointChoice())
_refPointChoice(_conf.refPointChoice()),
_intDist(_conf.intDist())
{
produces<GenParticleCollection>();
produces<CosmicLivetime,art::InSubRun>();
}

void CorsikaEventGenerator::beginSubRun(art::SubRun &subrun)
{
_primaries = 0;
}

void CorsikaEventGenerator::beginSubRun( art::SubRun &subrun){

void CorsikaEventGenerator::endSubRun(art::SubRun &subrun)
{
std::unique_ptr<CosmicLivetime> livetime(new CosmicLivetime(_primaries, _area, _lowE, _highE, _fluxConstant));
std::cout << *livetime << std::endl;
subrun.put(std::move(livetime));
}

void CorsikaEventGenerator::produce(art::Event &evt)
Expand Down Expand Up @@ -165,30 +171,92 @@ namespace mu2e {
const GenParticleCollection &particles(*corsikaParticles);
std::unique_ptr<GenParticleCollection> genParticles(new GenParticleCollection);

art::Handle<mu2e::CosmicLivetime> livetime;
evt.getByLabel(_corsikaModuleLabel, livetime);
_area = (*livetime).area();
_lowE = (*livetime).lowE();
_highE = (*livetime).highE();
_fluxConstant = (*livetime).fluxConstant();
_primaries += (*livetime).primaries();

/*
* In this block we check if there are two or more particles that intersect
* before the roof of the Mu2e world. If this happens, the particles are assumed to come
* from a common mother particle. This is to avoid simulating unphysical processes
* that can happen between the roof and the interaction point.
*/
std::vector<CLHEP::Hep3Vector> endPoints(particles.size());
std::vector<bool> doNotBacktrack(particles.size(), false);

if (_intDist >= 0) {
for (GenParticleCollection::const_iterator i = particles.begin(); i != particles.end(); ++i)
{
for (GenParticleCollection::const_iterator j = i+1; j != particles.end(); ++j)
{

GenParticle const &particle1 = *i;
GenParticle const &particle2 = *j;

TwoLinePCA twoLine(particle1.position(), -particle1.momentum().vect(), particle2.position(), -particle2.momentum().vect());
const CLHEP::Hep3Vector point1 = twoLine.point1();
const CLHEP::Hep3Vector point2 = twoLine.point2();

/*
* If their distance of closest approach is smaller than _intDist
* and the point of closest approach is between the roof and the
* top of the target box then we stop the backtracking at interaction
* point
*/
if (twoLine.dca() < _intDist &&
point1.y() > _targetBoxYmax &&
point1.y() < _worldYmax &&
point2.y() > _targetBoxYmax &&
point2.y() < _worldYmax)
{
doNotBacktrack[i - particles.begin()] = true;
doNotBacktrack[j - particles.begin()] = true;
endPoints[i - particles.begin()] = point1;
endPoints[j - particles.begin()] = point2;
}
}
}
}

for (GenParticleCollection::const_iterator i = particles.begin(); i != particles.end(); ++i)
{
GenParticle const &particle = *i;
const HepLorentzVector mom4 = particle.momentum();

const Hep3Vector position(
particle.position().x() + _cosmicReferencePointInMu2e.x(),
particle.position().y(),
particle.position().z() + _cosmicReferencePointInMu2e.z());

if (_projectToTargetBox)
{
_worldIntersections.clear();
VectorVolume particleWorld(particle.position(), particle.momentum().vect(),
_worldXmin, _worldXmax, _worldYmin, _worldYmax, _worldZmin, _worldZmax);
particleWorld.calIntersections(_worldIntersections);
// Check if the particle needs to be backtracked to the roof or not
if (doNotBacktrack[i - particles.begin()]) {
Hep3Vector projectedPos = endPoints[i - particles.begin()];

projectedPos.setX(projectedPos.x() + _cosmicReferencePointInMu2e.x());
projectedPos.setZ(projectedPos.z() + _cosmicReferencePointInMu2e.z());

if (_worldIntersections.size() > 0)
{
// Being inside the world volume, the intersection can be only one.
const Hep3Vector projectedPos = _worldIntersections.at(0);
genParticles->push_back(GenParticle(static_cast<PDGCode::type>(particle.pdgId()),
GenId::cosmicCORSIKA, projectedPos, mom4,
particle.time()));
} else {
_worldIntersections.clear();
VectorVolume particleWorld(position, particle.momentum().vect(),
_worldXmin, _worldXmax, _worldYmin, _worldYmax, _worldZmin, _worldZmax);
particleWorld.calIntersections(_worldIntersections);

if (_worldIntersections.size() > 0)
{
// Being inside the world volume, the intersection can be only one.
const Hep3Vector projectedPos = _worldIntersections.at(0);
genParticles->push_back(GenParticle(static_cast<PDGCode::type>(particle.pdgId()),
GenId::cosmicCORSIKA, projectedPos, mom4,
particle.time()));
}
}
} else {
genParticles->push_back(GenParticle(static_cast<PDGCode::type>(particle.pdgId()),
Expand Down
12 changes: 9 additions & 3 deletions EventGenerator/test/corsikaTest.fcl
Original file line number Diff line number Diff line change
Expand Up @@ -17,14 +17,18 @@ targetParams: {

source: {
module_type: FromCorsikaBinary
fileNames: ["/mu2e/app/users/srsoleti/corsika-77100/run/DAT110001"]
fileNames: ["/pnfs/mu2e/persistent/users/srsoleti/corsika/DAT110001"]
runNumber : 1
firstSubRunNumber : 0
firstEventNumber : 1
showerAreaExtension : 10000
maxEvents: -1
@table::targetParams
seed: 1
resample: false
fluxConstant: 1.8e4
lowE: 1.3
highE: 1e6
}

customCutPhoton: {
Expand Down Expand Up @@ -70,7 +74,9 @@ physics : {
module_type : CORSIKAEventGenerator
corsikaModuleLabel: "FromCorsikaBinary"
refPointChoice: "UNDEFINED"
@table::targetParams
projectToTargetBox : true
targetBoxYmax : 5000
intDist: -1
}

g4run :
Expand Down Expand Up @@ -128,7 +134,7 @@ physics : {

}

trigFilter: [corsikaGen, g4run, g4run, randomsaver]
trigFilter: [corsikaGen, g4run, randomsaver]
e1 : [outfile]
ana : [corsikaGenPlots]

Expand Down
78 changes: 78 additions & 0 deletions MCDataProducts/inc/CosmicLivetime.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
#ifndef MCDataProducts_CosmicLivetime_hh
#define MCDataProducts_CosmicLivetime_hh

// Original author Stefano Roberto Soleti

// C++ includes
#include <iostream>
#include <vector>
#include <math.h>

// Mu2e includes

namespace mu2e {

struct CosmicLivetime{

public:

CosmicLivetime():
_primaries(0),
_area(0),
_lowE(0),
_highE(0),
_fluxConstant(0) {
}

CosmicLivetime( unsigned int primaries,
float area,
float lowE,
float highE,
float fluxConstant ):
_primaries(primaries),
_area(area),
_lowE(lowE),
_highE(highE),
_fluxConstant(fluxConstant) {
const float eslope = -2.7;

const float EiToOneMinusGamma = pow(_lowE, 1 + eslope);
const float EfToOneMinusGamma = pow(_highE, 1 + eslope);
// http://pdg.lbl.gov/2018/reviews/rpp2018-rev-cosmic-rays.pdf eq. 29.2
_livetime = _primaries / (M_PI * _area * _fluxConstant * (EfToOneMinusGamma - EiToOneMinusGamma) / (1. + eslope));
}

// Accessors
unsigned int primaries() const { return _primaries; }
float area() const { return _area; }
float lowE() const { return _lowE;}
float highE() const { return _highE; }
float fluxConstant() const { return _fluxConstant; }
float liveTime() const { return _livetime; }


// Accept compiler generated versions of d'tor, copy c'tor, assignment operator.

// Print contents of the object.
void print( std::ostream& ost = std::cout, bool doEndl = true ) const;

private:

unsigned int _primaries;
float _area; // m2
float _lowE; // (GeV)
float _highE; // (GeV)
float _fluxConstant;
float _livetime; // s

};

inline std::ostream& operator<<( std::ostream& ost,
CosmicLivetime const& lt){
lt.print(ost,false);
return ost;
}

} // namespace mu2e

#endif /* MCDataProducts_CosmicLivetime_hh */
Loading

0 comments on commit a6d18b2

Please sign in to comment.