#include <root.h>#include "amschain.h"#include <TFile.h>#include <TH1F.h>#include <TH2.h>#include <map>#include <vector>#include <iostream>#include <unordered_map>#include "Point.h"#include "TRDConstants.h"#include "Utilities.h"#include "AMSTrdMCTrack.h"#include "Intersection3D.h"#include "ClusterFitter.h"#include "LinearAlgebra.h"#include "AMSTracker.h"
Go to the source code of this file.
Functions | |
| std::pair< int, double > | GetTRDVertexClassifierOutput (AMSEventR *pev, TH2D *h_TrackerL1_XY, bool debugMode=false) |
Function Documentation
◆ GetTRDVertexClassifierOutput()
| std::pair<int, double> GetTRDVertexClassifierOutput | ( | AMSEventR * | pev, |
| TH2D * | h_TrackerL1_XY, | ||
| bool | debugMode = false |
||
| ) |
Author: Anirvan Shukla
Simple header to be included in the miniDST production code. GetTRDVertexClassifierOutput() will classify an AMS event, and return std::pair<int, int> containing:
- TRDVertexStatus flag
- PrimaryParticleL1HitDistance value.
Make sure that the Tracker L1 geometry file "geometry_TrkL1_repaired.root" is available.
--— Define variables
TRDVertexStatus: To save as a branch in the miniDST. TRDVertexStatus = -99, no counters needed TRDVertexStatus = -39, nTRDVertexStatus_catchElse (catch final Else) TRDVertexStatus = -29, nTRDVertexStatus_error (If empty/bad/problematic event: e.g., no TRD tracks found, or no fitted tracks we found.) TRDVertexStatus = -19, nTRDVertexStatus_backg (All fitted track pairs in this event come from a single vertex) TRDVertexStatus = -9, nTRDVertexStatus_noFitPrimMC (If no fitted primary MC track found; replacing the old nTRDVertexStatus_noFitProton) TRDVertexStatus = -2, nTRDVertexStatus_multiStars1 (A fitted track intersect with some fitted tracks, but NOT with others. Implies multiple stars/vertices/Cascade. TRDVertexStatus = -1; nTRDVertexStatus_multiStars2 (many scenarios are possible, including multiple stars/vertices/Cascade. Some of these could be good events. TRDVertexStatus = 0, nTRDVertexStatus_mixed (At least one track is NOT pairing up the other tracks in the event. Could be good.) TRDVertexStatus = 1, nTRDVertexStatus_signal1 (No vertex/star found. Clean single track. This is our main signal event.) TRDVertexStatus = 2, nTRDVertexStatus_signal2 (No vertex/star found. Clean multiple tracks.) TRDVertexStatus = 3, nTRDVertexStatus_signal3 (Two tracks, one of them is an electron, and it is coming from somewhere on the other track.)
default/init value
PrimaryParticleL1HitDistance is the closest distance (on the X-Y plane, at the Z position of L1) between the extrapolated TRD track and MC L1 hit (if it exists). PrimaryParticleL1HitDistance = -9, bad/error PrimaryParticleL1HitDistance = -5, no L1 hit found! PrimaryParticleL1HitDistance > 0, L1 hit found, and closest distance is provided.
init as bad/error.
define per-event counters
Analyze the global MC track clusters check if particle enters the inner tracker: —> Define the inner tracker 3D polygon, —> Then check if any MCEventgR global cluster is inside it.
this is the number of global hits/clusters in this MC event.
container to store (x, y) of L1 hits, to compare to fitted 3D TRD tracks later.
Fill only for hits in L1 z-position = 158.920 cm from Travis thesis. For X and Y, see TH2D hQA_TrackerL1_XY.
Z (158, 159), defined in AMSTracker.h
X, Y acceptance map
add point to event vector
x, y, trkID, PID, closedis
end loop over nMCEventg clusters in this event
Print the unique elements and their frequencies
Now find the unique trkIDs of MC TRD track clusters, then put MC track clusters with the same trkIDs together as a track. This needs to be done because, unfortunately, the AMSEventR does not store MC tracks as one object. Instead it stores all track clusters, and their associted trkIDs (and PIDs etc.).
Remove entries with frequency less than 2 and shrink the map, as we need at least 2 clusters to make a linear fit
store separately, to correctly account for tracks with bad fits or bad G3PIDs
Print the unique elements and their frequencies
Need at least 2 clusters to make a linear fit
Now that the number of unique trkIDs have been identified in this event we re-loop over iTRDMCCluster, and group these clusters under separate AMSTrdMCTrack::Tracks.
container to store MC cluster under identified MC TRD tracks
Need at least 2 clusters to make a linear fit
Set the G3PID for each track. Take care of edge cases or bad clusters
This can happen some times in the AMS data.
Sanity checks
This should never happen!
Print the number of identified/saved tracks and clusters:
Now perform track fits:
container to store the three "good" (XZ, YZ and XY) lin fits of each AMSTrdMCTrack::Track
if no tracks could be fitted, move to next event
if the primary MC track could be fitted inside the TRD, move to next event
all three plane projections have same trkID and g3PID. choose any one, and check if it is the primary MC track:
this should be at most 1.
if primary MC track is not properly fitted in the TRD, then set this to background
Now we find if the fitted 3D lines have a matched point on the Tracker L1
We have a 3D line, but we do not know its equation. However we do know the equations of its three projections in the XZ, YZ and XY plane: XZ projection: z = m1*x + b1 YZ projection: z = m2*y + b2 XY projection: y = m3*x + b3 Use this information to get the equation of the 3D line, in the vector/parametric form.
XZ projection: z = m1*x + b1
YZ projection: z = m2*y + b2
XY projection: y = m3*x + b3
Perform the fit
Check if the vectors are not empty and then use them
Good fit found!
The line3D object was default-constructed and its point and direction vectors were not set.
move to next fitted track
For each fitted line3D, loop over the global MC hits from MCEventgR and find the closest cluster.
Check if the returned pair is empty
Also check if the trkIDs and G3PIDs match
see if this is the primary MC particle as well:
Now we find how many of the fitted AMSTrdMCTrack::Track pairs intersect
For N fitted tracks, we have nMAX_POSSIBLE_PAIRS = N(N-1)/2 (max possible track pairs). If nIntersecting_3Dlines < N(N-1)/2, we may have one clean TRD track.
For example: 4 tracks. Expected pairs = 6. This is the nMAX_POSSIBLE_PAIRS. If pairs found = 3, i.e., nNEXTTOMAX_POSSIBLE_PAIRS = 3(3-1)/2: Then we have one track that does not intersect with the other 3 tracks. This is a clean TRD event then, but with a special flag.
But if 5 pairs were found (i.e., 5 > nNEXTTOMAX_POSSIBLE_PAIRS), then this implies multiple stars/vertices in this event. Set this with another special flag.
We could consider even more NEXTTONEXTTOMAX_POSSIBLE_PAIRS and so on, and try to reconstruct the multiple showers/stars... Not now, perhaps in future. :-)
Already ensured that nFittedTracks > 0
Use the fits to get the number of intersecting 3D line pairs The if-else statements have been checked for zero overlap :-)
Decide if signal or background event i.e., the TRDVertexStatus
No vertex/star found. This is our main signal event.
Clean single track: good!
Clean multiple tracks: good!
special case for 2 intersecting tracks one track is just a secondary e+ splitting off the main track
2 fitted tracks
at least one of them is a secondary e-
one track is primary, one is secondary
All fitted track pairs in this event come from a single vertex This is a background event.
At least one track is NOT pairing up the other tracks in the event. The other fitted track pairs in this event come from a single vertex, UNLESS if the event as only two tracks. Then both tracks are clean and independent! But this is already covered by TRDVertexStatus = 10 / nTRDVertexStatus_signal2. So set this as a mixed event (it could be good!)
A fitted track intersect with a few fitted tracks, but does NOT interact with others. This implies multiple stars/vertices in a cascade Set this with a separate background event code.
if nFittedTracks == 3, then nNEXTTOMAX_POSSIBLE_PAIRS = 1, and nIntersecting_3Dlines = 0. This is already covered by TRDVertexStatus = 10 / nTRDVertexStatus_signal2. But if nFittedTracks > 3, then many scenarios are possible, including multiple stars/vertices in a cascade. Set this to background (even though some of these could be good events).
Ideally, this should never happen!
clear per-event containers to reset before next event
