Added cosmics example with analysis code

# 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. Because the time output of the CRY simulation code needs to be known for the analysis, we store the log in a *.txt* file:

```bash

allpix -c cosmic_flux.conf >> cosmic_flux_log.txt

```

To analyse the tracks of the *MCParticles*, we convert them into a text list. Open the root framework and issue

```

.L PATH_TO_ALLPIX_INSTALL/lib/libAllpixObjects.so

```

Next, load the conversion macro via

```

.L analysis/rootMacro.cc

```

Finally, execute it with

```

rootMacro()

```

and close root with

```

.q

```

Now you can analyse the result by issuing

```bash

python analysis/analysis.py

```

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

from track import Track

from parse import compile

from tqdm import tqdm

from simulatedTime import SimulatedTime

from histogram import Histogram

"""

Analyse the MCTracks.txt file and plot the observed muon flux

"""

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

trackParser = compile("({:g},{:g},{:g}) ({:g},{:g},{:g})")

def getTracks():

    time = SimulatedTime("cosmic_flux_log.txt").time

    tracks = []

    input = "MCTracks.txt"

    for l in open(input, "r").readlines():

        l = l.strip()

        res = trackParser.parse(l)

        tracks.append(Track(res[0:3], res[3:6], 0, 0, 0, 0))

    return tracks, time

tracks, time = getTracks()

histogram.addTracks(tracks, time)

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)

from bin import Bin

import numpy as np

import matplotlib.pyplot as plt

class Histogram(object):

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

        self.azimuthalLinspace = np.linspace(

            azimuthal[0], azimuthal[1], num=granularity[0] + 1, endpoint=True)

        self.zenithLinspace = np.linspace(

            zenith[0], zenith[1], num=granularity[1] + 1, endpoint=True)

        self.bins = [[None for j in range(granularity[1])]

                     for i in range(granularity[0])]

        for i in range(granularity[0]):

            for j in range(granularity[1]):

                azimuthalRange = [self.azimuthalLinspace[i],

                                  self.azimuthalLinspace[i + 1]]

                zenithRange = [self.zenithLinspace[j],

                               self.zenithLinspace[j + 1]]

                b = Bin(azimuthal=azimuthalRange, zenith=zenithRange, deg=deg)

                self.bins[i][j] = b

        self.granularity = granularity

    def addTracks(self, tracks, time):

        for a in self.bins:

            for b in a:

                b.addTracks(tracks, time)

    def plot(self):

        for a in self.bins:

            for b in a:

                b.getCOMPoints()

    def printFlux(self):

        for a in self.bins:

            for b in a:

                print("Bin", b.zenith, b.azimuthal,

                      ":", b.flux(), "Muons/s/m^2/sr")

    def plotZenith(self):

        assert self.granularity[0] == 1

        x = [(self.zenithLinspace[i] + self.zenithLinspace[i + 1]) /

             2 for i in range(len(self.zenithLinspace) - 1)]

        xerr = [(self.zenithLinspace[i + 1] - self.zenithLinspace[i]

                 ) / 2 for i in range(len(self.zenithLinspace) - 1)]

        flux = []

        for a in self.bins:

            for b in a:

                flux.append(b.flux())

        y = [f.n for f in flux]

        yerr = [f.s for f in flux]

        print(x)

        plt.errorbar(x, y, xerr=xerr, yerr=yerr, linestyle="None", marker='.')

        xcos = np.linspace(

            self.zenithLinspace[0], self.zenithLinspace[-1], 200, endpoint=True)

        ycos = y[0] * np.cos(np.deg2rad(xcos))**2

        #plt.plot(xcos, ycos)

        plt.xlabel(r"Zenith Angle [$^\circ$]")

        plt.ylabel(r"Muon Flux [Muons $s^{-1}$ $sr^{-1}$ $m^{-2}$]")

        x = np.linspace(0, 90, 100)

        y = np.cos(np.deg2rad(x))**2 * max(y)

        plt.plot(x, y)

        plt.savefig("ZenithFlux.pdf")

        plt.show()

 No newline at end of file

#include <TCanvas.h>

#include <TFile.h>

#include <TH1D.h>

#include <TH2D.h>

#include <TSystem.h>

#include <TTree.h>

#include <fstream>

#include <iostream>

#include <memory>

#include <stdlib.h>

using namespace std;

/**

 * An example of how to get iterate over data from Allpix Squared,

 * how to obtain information on individual objects and the object history.

 */

void rootMacro() {

    std::string detector = "telescope0_0";

    TFile* file = new TFile("output/cosmicsMC.root");

    // Initialise reading of the MCParticle TTrees

    TTree* mc_particle_tree = static_cast<TTree*>(file->Get("MCParticle"));

    if(!mc_particle_tree) {

        std::cout << "Could not read tree MCParticle" << std::endl;

        return;

    }

    TBranch* mc_particle_branch = mc_particle_tree->FindBranch(detector.c_str());

    if(!mc_particle_branch) {

        std::cout << "Could not find the branch on tree MCParticle for the corresponding detector, cannot continue"

                  << std::endl;

        return;

    }

    // Bind the information to a predefined vector

    std::vector<allpix::MCParticle*> input_particles;

    mc_particle_branch->SetObject(&input_particles);

    ofstream myfile;

    myfile.open("MCTracks.txt");

    // Iterate over all events

    for(int i = 0; i < mc_particle_tree->GetEntries(); ++i) {

        if(i % 100 == 0) {

            std::cout << "Processing event " << i << std::endl;

        }

        // Access next event. Pushes information into input_*

        mc_particle_tree->GetEntry(i);

        // Loop over all particles in the event

        for(auto& mc_part : input_particles) {

            // Only look at (anti-)muons (PDG-ID: (-)13)

            if(abs(mc_part->getParticleID()) == 13) {

                // Get track info

                auto startPoint = mc_part->getLocalStartPoint();

                auto endPoint = mc_part->getLocalEndPoint();

                auto direction = startPoint - endPoint;

                std::cout << startPoint << " " << direction / direction.z() << std::endl;

                // Write to file

                myfile << startPoint << " " << direction / direction.z() << std::endl;

            }

        }

    }

    myfile.close();

}