Merge branch 'cosmicsExample' into 'master'

* Tobias Bisanz, CERN, @tbisanz

* Marco Bomben, Université de Paris, @mbomben

* Koen van den Brandt, Nikhef, @kvandenb

* Carsten Daniel Burgard, DESY @cburgard

* Carsten Daniel Burgard, DESY, @cburgard

* Maximilian Felix Caspar, DESY, @mcaspar

* Liejian Chen, Institute of High Energy Physics Beijing, @chenlj

* Katharina Dort, University of Gie\ss en, @kdort

* Neal Gauvin, Université de Genève, @ngauvin

\item Marco Bomben, Université de Paris

\item Koen van den Brandt, Nikhef

\item Carsten Daniel Burgard, DESY

\item Maximilian Felix Caspar, DESY

\item Liejian Chen, Institute of High Energy Physics Beijing

\item Katharina Dort, University of Gie\ss en

\item Neal Gauvin, Université de Genève

# Cosmic Muon Flux Simulation

This example illustrates how the *depositionCosmics* module is used to model the flux of cosmic muons in Allpix Squared. Python analysis code is included to calculate the flux through the detector from the *MCParticle* objects

## Usage

First of all, change into the example directory. Run the simulation example from here:

```bash

allpix -c cosmic_flux.conf

```

To analyse the tracks of the *MCParticles*, issue

```bash

python analysis/analysis.py -l PATH_TO_ALLPIX_INSTALL/lib/libAllpixObjects.so

```

The library flag is only required when your *allpix-squared/lib* path is not in *LD_LIBRARY_PATH*.

Notice that the analysis script relies on the python modules *uncertainties*, *tqdm*, *matplotlib* and *numpy.*

from track import Track

from tqdm import tqdm

from histogram import Histogram

import argparse

import os

import os.path as path

import ROOT

from ROOT import gSystem

"""

Analyse the cosmicsMC.root file and plot the observed muon flux

"""

# argument parser

parser = argparse.ArgumentParser()

parser.add_argument("-l", metavar='libAllpixObjects', required=False,

                    help="specify path to the libAllpixObjects library (generally in allpix-squared/lib/)) ")

args = parser.parse_args()

if args.l is not None:  # Try to find Allpix Library

    lib_file_name = (str(args.l))

else:  # Look in LD_LIBRARY_PATH

    libraryPaths = os.environ['LD_LIBRARY_PATH'].split(':')

    for p in libraryPaths:

        if path.isfile(path.join(p, "libAllpixObjects.so")):

            lib_file_name = path.join(p, "libAllpixObjects.so")

            break

if (not os.path.isfile(lib_file_name)):

    print("WARNING: ", lib_file_name, " does not exist, exiting")

    exit(1)

histogram = Histogram(granularity=[1, 15])

# load library and rootfile

gSystem.Load(lib_file_name)

rootFile = ROOT.TFile("output/cosmicsMC.root")

McParticle = rootFile.Get('MCParticle')

def getTracks():

    "Load tracks from the ROOT file"

    # time is given in ns

    time = float(

        str(rootFile.config.DepositionCosmics.total_time_simulated)) / 1e9

    tracks = []

    for i in tqdm(range(0, McParticle.GetEntries())):

        McParticle.GetEntry(i)

        McParticle_branch = McParticle.GetBranch("telescope0_0")

        br_mc_part = getattr(McParticle, McParticle_branch.GetName())

        startPoint = None

        endPoint = None

        for mc_part in br_mc_part:

            # Only consider muons

            if abs(mc_part.getParticleID()) == 13:

                startPoint = mc_part.getLocalStartPoint()

                endPoint = mc_part.getLocalEndPoint()

                direction = startPoint - endPoint

                direction = direction / direction.z()

                tracks.append(Track([startPoint.x(), startPoint.y(), startPoint.z()], [

                              direction.x(), direction.y(), direction.z()], 0, 0, 0, 0))

    return tracks, time

# Get tracks from file

tracks, time = getTracks()

histogram.addTracks(tracks, time)

# Do the analysis and show the plot

histogram.printFlux()

histogram.plotZenith()

import numpy as np

from solidAngle import SolidAngle

from uncertainties import ufloat

import pickle

class Bin(object):

    def __init__(self, azimuthal=[0, 360], zenith=[0, 90], deg=True):

        self.solidAngle = SolidAngle(azimuthal, zenith, deg=deg)

        self.zenith = zenith

        self.azimuthal = azimuthal

        self.deg = deg

        self.numberOfEntries = 0

        self.tracks = []

        self.time = 0  # in seconds

    def _addTrack(self, t):

        if self.azimuthal[0] <= t.azimuthalAngle <= self.azimuthal[1]:

            if self.zenith[0] <= t.zenithAngle <= self.zenith[1]:

                self.tracks.append(t)

    def addTracks(self, tracks, time):

        self.time += time

        for t in tracks:

            self._addTrack(t)

        #print("Now have", len(self.tracks),"tracks in bin", self.zenith, self.azimuthal)

    def getCOMPoints(self):

        return [t.COMIntersection() for t in self.tracks]

    def plot(self):

        import matplotlib.pyplot as plt

        x = [t.COMIntersection()[0] for t in self.tracks]

        y = [t.COMIntersection()[1] for t in self.tracks]

        plt.scatter(x, y)

        plt.show()

    def realSurface(self):

        return np.cos(np.deg2rad(np.mean(self.zenith))) * 0.3455

    def flux(self):

        area = self.realSurface()

        if area <= 0:

            return ufloat(0, 0)

        n = len(self.tracks)

        # / self.efficiency

        return ufloat(n, n**0.5) / self.time / self.solidAngle.value / area

if __name__ == '__main__':

    b = Bin(azimuthal=[0, 360], zenith=[0, 20], deg=True)

    print(b)