Commit e3a3d28f authored by Maximilian Caspar's avatar Maximilian Caspar
Browse files

Switched from ROOT macro conversion to pure python analysis

parent f1423bd7
Loading
Loading
Loading
Loading
+4 −19
Original line number Diff line number Diff line
@@ -6,24 +6,9 @@ First of all, change into the example directory. Because the time output of the
```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
To analyse the tracks of the *MCParticles*, issue
```bash
python analysis/analysis.py
python analysis/analysis.py -l PATH_TO_ALLPIX_INSTALL/lib/libAllpixObjects.so
```
Notice that the analysis script relies on the python modules *uncertainties*, *parse*, *matplotlib* and *numpy.*
The library flag is required when your Allpix installation is not in a standard path.
Notice that the analysis script relies on the python modules *uncertainties*, *tqdm*, *matplotlib* and *numpy.*
+51 −10
Original line number Diff line number Diff line
from track import Track
from parse import compile
from tqdm import tqdm
from simulatedTime import SimulatedTime
from histogram import Histogram
import argparse
import os
import os.path as path
import ROOT
from ROOT import gSystem

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

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

trackParser = compile("({:g},{:g},{:g}) ({:g},{:g},{:g})")
# 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))
    if (not os.path.isfile(lib_file_name)):
        print("WARNING: ", lib_file_name, " does not exist, exiting")
        exit(1)

elif os.path.isfile(path.abspath(path.join(__file__, "..", "..", "opt", "allpix-squared", "lib", "libAllpixObjects.so"))):  # For native installs
    lib_file_name = path.abspath(path.join(
        __file__, "..", "..", "opt", "allpix-squared", "lib", "libAllpixObjects.so"))

elif os.path.isfile(path.join(path.sep, "opt", "allpix-squared", "lib", "libAllpixObjects.so")):  # For Docker installs
    lib_file_name = path.join(
        path.sep, "opt", "allpix-squared", "lib", "libAllpixObjects.so")


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 = 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

    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())

tracks, time = getTracks()
        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()
+0 −59
Original line number Diff line number Diff line
#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();
}