Project

General

Profile

AnalysisTree variables » History » Version 243

« Previous - Version 243/246 (diff) - Next » - Current version
En-Chuan Huang, 03/27/2017 12:07 PM


What is AnalysisTree?

In LArSoft jargon, it is an art::EDAnalyzer module that resides in the “uboonecode” repository. In simple words, it makes a root TTree with tons of useful information that enables the user to perform wide spectrum of studies: Physics analysis (main goal for its existence!), Calibration studies, reconstruction validation studies, truth-level studies, Cosmic ray studies etc.

Starting MCC2 Monte Carlo production, AnalysisTree is being run with every MC production.

The AnalysisTree code resides in the MicroBooNE maintained uboonecode repository: https://cdcvs.fnal.gov/redmine/projects/uboonecode/repository/revisions/develop/entry/uboone/AnalysisTree/AnalysisTree_module.cc

The corresponding fhicl file to run this module all resides in the same repository:
https://cdcvs.fnal.gov/redmine/projects/uboonecode/repository/revisions/develop/changes/uboone/AnalysisTree/anatree_uboone.fcl

Maintainers: Tingjun, Sowjanya and Gianluca
Contributors: Tingjun, Sowjanya, Gianluca (and others)

More technical details on how to check-out analysis tree, or develop etc. can be found in the Tutorials section of this wiki page.

Configure what gets saved in the AnalysisTree

User can choose what information to dump in the ntuples by using fcl parameters. Analysistree can save following types of information:

1. MC truth (generator-level)
        -> GENIE info
        -> CRY info
2. MC particle (GEANT4 level)
        -> MC Shower information
        -> MC Track information
        -> Auxiliary detector info
3. Hit information
4. Hits from Raw Digits
5. Hits from Calibrated Waveforms
6. Hits from Sim Channels
7. Raw Waveform information
8. Calibrated Waveform information
9. Optical flash information
10. Cluster information
11. Track information (multiple trackers as provided by the user in the fcl file)
12. Vertex information
13. Shower information

Except for CRY and Auxiliary detector information, all other information is saved by default (see analysistreemodule.fcl in the srcs/uboonecode/uboone/AnalysisTree directory). User can turn on/off the information they like by setting the following fcl parameters to true or false:

SaveAuxDetInfo            //save auxiliary information?
SaveCryInfo                  //save CRY information?
SaveGenieInfo              //save GENIE information? 
SaveGeantInfo              //save GEANT information?
SaveMCShowerInfo       //save MC Shower information?
SaveMCTrackInfo       //save MC Shower information?
SaveHitInfo                  //save Hit information? 
SaveRawDigitInfo        //hit extraction at the raw-digit level (Sowjanya and Tingjun's recipe, if you would like to use it, contact Sowjanya). By default off
SaveCalWireInfo        //hit extraction at the CalWire level (Sowjanya's recipe, if you would like to use it, contact Sowjanya). By default off
SaveSimChannelInfo        //hit extraction at the Sim Channel level (Sowjanya's recipe, if you would like to use it, contact Sowjanya). By default off
SaveTrackInfo              //save Track information?
SaveShowerInfo           //save reconstructed Shower information?
SaveClusterInfo           //save reconstructed Cluster information?
SaveVertexInfo            //save Vertex information?
SaveFlashInfo              //save Flash information?
SaveCaloCosmics        //save calorimetry for cosmics?
SaveRawWaveFormInfo   //save raw waveform information? (by default off)
SaveCalibWaveFormInfo  //save calibrated waveform information? (by default off)

The following example shows how to set true or false values to the above fcl parameters in your anatree_uboone.fcl file:

physics.analyzers.analysistree.SaveCryInfo: true
physics.analyzers.analysistree.SaveGenieInfo: false
physics.analyzers.analysistree.SaveHitInfo: false
physics.analyzers.analysistree.SaveTrackInfo: true
physics.analyzers.analysistree.SaveVertexInfo: false
physics.analyzers.analysistree.SaveClusterInfo: false
physics.analyzers.analysistree.SaveShowerInfo: false
physics.analyzers.analysistree.SaveGeantInfo: true
physics.analyzers.analysistree.SaveMCShowerInfo: true
physics.analyzers.analysistree.SaveMCTrackInfo: true
physics.analyzers.analysistree.SaveAuxDetInfo: false
...
...

Important note: In the case of GEANT information (doesn't apply to MC Shower information), an energy cut of 10 MeV is imposed on the secondary particles. All primary particles are saved with no energy cuts. One can always configure the energy cut by accessing the fcl parameter G4MinE as follows:

physics.analyzers.analysistree.G4MinE: 0.01

AnalysisTree variables

Run information

Analysis tree holds the following run information, these are saved on a per-event basis:

    Int_t      run;                  //run number
    Int_t      subrun;               //subrun number
    Int_t      event;                //event number
    Double_t   evttime;              //event time in sec
    Double_t   beamtime;             //beam time in sec
    Double_t   taulife;              //electron lifetime put in simulation (3ms)
    Char_t     isdata;               //flag, 0=MC 1=data

Information from the sub-run

Although, you will see that the pot information is entered per event in the tree, the pot information in the analysis tree is actually per sub-run not per event.

    Double_t pot; //protons on target

For example, in the plot below, which shows the pot information from Analysistree for a GENIE BNB run, the sub-run has 100 events and you will see the "same pot value (2.353e+17)" is entered 100 times in the histogram. This one value (2.353e+17) corresponds to the total pot for 100 events. One DOES NOT have to multiply the mean 2.353e+17 by 100 to get the total pot. This detail is especially relevant when you run jobs on grid, in which case, each grid job/submission is considered a sub-run. And the pot information you see per job corresponds to the total number of events in that grid job.

Software trigger information

The software trigger results are saved in AnalysisTree if the flag SaveSWTriggerInfo is set to true and a valid software trigger producer is specified in SWTriggerModuleLabel. The software trigger module in uboonecode may produce several outputs as it may run several algorithms and check is the event is triggered based on these algorithm logic. For this reason, the trigger output is stored in std::vector objects. You can retrieve the name of the trigger algorithm and a flag. The flag is a boolean and tells is the event is triggered or not. The AnalysisTree variables are:

  std::vector<std::string> swtrigger_name;       // the name of the trigger algorithm
  std::vector<bool>        swtrigger_triggered;  // true = event is triggered; false = event is not triggered  based on the relative algorithm logic

An example on how to access the software trigger information:

    //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        int nSWTriggerAlgos = swtrigger_name.size();
        for (int trigger = 0; trigger < nSWTriggerAlgos; trigger++) {
          std::cout << "SW Trigger name: " << swtrigger_name[trigger] << std::endl;
          std::cout << "  -> was triggered? " << swtrigger_triggered[trigger] << std::endl;
        }
        . . .
    }

Hit information

Note that only hit information from the hit module specified by HitsModuleLabel in the fcl file used to produce this analysis tree will be stored here. (Information about the fcl files used to produce Monte Carlo Challenge (MCC) files can be found in pages available from the MCC index page)

The value of kMaxHits is set at compile time. The data stored is non-resizeable, e.g., 45x kMaxHits = 900k bytes worth for kMaxHits=20000 (tag v02_05_01, used for MCC5), proportionally higher for kMaxHits=25000 (v03_03_00).

The actual number of hits is given by the no_hits variable. (Here "no" is an abbreviation for "number", not the word "no".)

Note: hit times (hit_peakT, hit_startT, hit_endT) are all in time ticks not in real time. Using ticks as units is less ambiguous because, there are multiple ways to interpret what you mean by time, is it w.r.t. beam trigger or something else. One way to correspond the ticks to real time is: each tick corresponds to 500ns sampling time, so if you multiply the hit times by 500ns, this will give you the generation time (not the time w.r.t. the trigger, but the actual generation time).

    Int_t    no_hits;                  //number of hits
    Short_t   hit_plane[kMaxHits];      //plane number corresponding to this hit
    Short_t  hit_wire[kMaxHits];        //index of wire number corresponding to this hit
    Short_t  hit_channel[kMaxHits];   //index of channel ID (raw digit) corresponding to this hit
    Float_t  hit_peakT[kMaxHits];    //time corresponding to the peak charge deposition (in units of time ticks)
    Float_t  hit_charge[kMaxHits];     //area (total charge (ADC) deposited for hit in the tdc range) 
    Float_t  hit_ph[kMaxHits];         //amplitude (maximum ADC value in hit window)
    Float_t  hit_startT[kMaxHits];     //hit start time (initial tdc) in units of time ticks 
    Float_t  hit_endT[kMaxHits];       //hit end time (final tdc) in units of time ticks
    Float_t  hit_goodnessOfFit[kMaxHits]; //chi2/dof goodness of fit 
    Short_t  hit_multiplicity[kMaxHits];  //multiplicity of the given hit                                         
    Float_t  hit_trueX[kMaxHits];      // hit true X (cm)
    Float_t  hit_nelec[kMaxHits];     //hit number of electrons
    Float_t  hit_energy[kMaxHits];       //hit energy
    Short_t  hit_trkid[kMaxHits];      //is this hit associated with a reco track?
    Short_t  hit_trkKey[kMaxHits];      //is this hit associated with a reco track,  if so associate a unique track key ID?
    Short_t  hit_clusterid[kMaxHits];  //is this hit associated with a reco cluster?
    Short_t  hit_clusterKey[kMaxHits];  //is this hit associated with a reco cluster, if so associate a unique cluster key ID?

The hit_trkid[] variable as defined above (which tells whether a given hit is associated to a track or not) is only stored for the first tracker currently (usually trackkalmanhit). For example, if you ran the Analysistree on two trackers, say, trackkalmanhit and trackkalmansps, you will see the hit_trkid entries listed in the ntuple only for trackkalmanhit tracker. Since the hit finder is not uniform from tracker to tracker, it is not easy to save this information for all trackers. Also, trkID or ClusterID are saved purely based on author's preference and choice and in some cases they may not be unique from event to event, so, to deal with such cases, we save hit_trkKey and hit_clusterKey which are always unique objects to map track associated hits or cluster associated hits to hit collection.

An example template below shows how to access the hit information in a root macro,

    //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the total number of hits
        for(Int_t j=0;j<no_hits;j++) {
             //for e.g., if you want to get the hit width
             float width = hit_startT[j] - hit_endT[j];
             //look at hit charge information,
             std::cout<<"\n Total hit charge ="<<hit_charge[j]<<"\t Hit_peakcharge="<<hit_ph[j];
             . . .
        }
    }

Hits from Raw Digits

This is not for everyone, take caution when using this and contact Sowjanya if you would like to use this.

The recipe used to extract raw digits is described, for example in this talk on Slide 16:
http://microboone-docdb.fnal.gov:8080/cgi-bin/RetrieveFile?docid=5798&filename=DiffusionUpdateMay12.pdf&version=1

The raw digit ROI is currently taken to be 3 times the reconstructed hit RMS on both sides of the reconstructed hit peak time.
One can set the "n" sigma RMS ROI using the fcl parameter "SetRawDigitROI" (for example if you set this to 4, the window for raw digits will go from -4*sigma to +4*sigma)
One can also set the ADC threshold above which to integrate raw digits to extract hit area and hit RMS using the fcl parameter "SetThreshold"
(this threshold is set to a certain percent of the pulse height, for example if it is set to 20, all digits with an ADC above 20% of the pulse height are used)

By default,
SetRawDigitROI is set to 3
SetRawDigitThresh is set to 10 (which means 10% of the raw digit pulse height)

//RawDigit information
    Float_t  rawD_ph[kMaxHits];  
    Int_t    rawD_peakT[kMaxHits];  
    Float_t  rawD_charge[kMaxHits];  
    Float_t  rawD_fwhh[kMaxHits];  
    Double_t rawD_rms[kMaxHits]; 

Hits from Calibrated signals

This is not for everyone, take caution when using this and contact Sowjanya if you would like to use this.
Same description as raw digits hold here.

By default,
SetCalWireROI is set to 3
SetCalWireThresh is set to 10 (which means 10% of the raw digit pulse height)

    //CalWire information
    Float_t  calwire_ph[kMaxHits];  
    Int_t    calwire_peakT[kMaxHits];  
    Float_t  calwire_charge[kMaxHits];  
    Float_t  calwire_fwhh[kMaxHits];  
    Double_t calwire_rms[kMaxHits]; 

Hits from Sim Channels

This is not for everyone, take caution when using this and contact Sowjanya if you would like to use this.
This is basically representing hits in number of electrons. Consider this as truth information and can be used as a cross-check. See the analysis tree code to understand how this is done.

   //SimChannel information
    Float_t  sim_ph[kMaxHits];  
    Int_t    sim_tdc[kMaxHits];  
    Float_t  sim_charge[kMaxHits];  
    Float_t  sim_fwhh[kMaxHits];  
    Double_t sim_rms[kMaxHits]; 

Raw and Calibrated Wave form information

One can save the raw waveform and calibrated waveform information in AnalysisTree by setting the "SaveRawWaveFormInfo" and "SaveCalibWaveFormInfo" to true in the analysis tree fcl file.
(see section "Configure what gets saved in AnalysisTree" for details on how to do this).
The number of samples for which the waveform is saved can be set using the variable "kMaxTicks". This is currently set to 9600.
The raw_channelId (calib_channelId) are 1D vectors with the size "raw_nchannels" (calib_nchannels). The raw_wf (calib_wf) is 2D vectors with the first index the size of raw_nchannels (calib_nchnanels) and the second index is currently set to have a max size of 9600. The data structures for both forms of information is implemented as boxed arrays, the same structure we used for tracks.

To extract raw waveforms, the code looks for "DigitModuleLabel".

Short_t                         raw_nchannels;             //number of raw channels
RawData_t<Short_t>              raw_channelId;
RawWaveFormData_t<Short_t>      raw_wf; 

To extract calibrated waveforms, the code looks for "CalDataModuleLabel"

Short_t                        calib_nchannels;        //number of calibration channels
CalibData_t<Short_t>           calib_channelId;
CalibWaveFormData_t<Short_t>     calib_wf; 

Cluster information

Note that only cluster information from the cluster module specified by ClusterModuleLabel in the fcl file used to produce this analysis tree will be stored here. Currently can only store 1000 clusters per event (as specified by the kMaxClusters variable in analysis tree). The actual number of clusters is given by the nclusters variable.

    Short_t nclusters;     //number of clusters in a  given event
    Short_t clusterId[kMaxClusters];   //ID of this cluster      
    Short_t clusterView[kMaxClusters];   //which plane this cluster belongs to          
    Int_t   cluster_isValid[kMaxClusters];  //is this cluster valid? will have a value of -1 if it is not valid
    Float_t cluster_StartCharge[kMaxClusters];  //charge on the first wire of the cluster in ADC
    Float_t cluster_StartAngle[kMaxClusters];    //starting angle of the cluster
    Float_t cluster_EndCharge[kMaxClusters];   //charge on the last wire of the cluster in ADC
    Float_t cluster_EndAngle[kMaxClusters];    //ending angle of the cluster
    Float_t cluster_Integral[kMaxClusters];      //returns the total charge of the cluster from hit shape in ADC
    Float_t cluster_IntegralAverage[kMaxClusters]; //average charge of the cluster hits in ADC
    Float_t cluster_SummedADC[kMaxClusters];  //total charge of the cluster from signal ADC counts
    Float_t cluster_SummedADCaverage[kMaxClusters]; //average signal ADC counts of the cluster hits.
    Float_t cluster_MultipleHitDensity[kMaxClusters]; //Density of wires in the cluster with more than one hit. 
    Float_t cluster_Width[kMaxClusters];  //cluster width in ? units
    Short_t cluster_NHits[kMaxClusters]; //Number of hits in the cluster
    Short_t cluster_StartWire[kMaxClusters]; //wire coordinate of the start of the cluster 
    Short_t cluster_StartTick[kMaxClusters];  //tick coordinate of the start of the cluster in time ticks
    Short_t cluster_EndWire[kMaxClusters];  //wire coordinate of the end of the cluster
    Short_t cluster_EndTick[kMaxClusters];  //tick coordinate of the end of the cluster in time ticks    
Cosmic Cluster tagging information in analysis tree:
  • In addition to the above information, Analysis tree also saves one type of cosmic tagging information for clusters: Geometry_based cluster tagger (variables will have the word "tagger" as part of their name). Unlike track information, cluster cosmic tagging information is only saved for one cluster tagging algorithm. One can specify the cosmic cluster tagging label using CosmicClusterTaggerAssocLabel.
    For example, the following line needs to be added in the analysis tree fhicl file:

physics.analyzers.analysistree.CosmicClusterTaggerAssocLabel: "ccclustertag"

  • For each cluster, the following cosmic tag information is saved:
    1. number of cosmic tags associated to a given cluster (almost always there is only one tag per cluster, saved for sanity check)
    2. cosmic score (1=definite cosmic; 0=non-cosmic ( or unknown); 0.5= unsure; -1=very likely neutrino).
    3. Type of cosmic tag (-1=unknown type; 0=untagged; 1=crossed 2 Y boundaries; 21=crossed one Y boundary; 100=partially outside drift;)

Additional details of cosmic scores and tags can be found here:
https://cdcvs.fnal.gov/redmine/projects/lardata/repository/revisions/develop/entry/AnalysisBase/CosmicTag.h

    Short_t cluncosmictags_tagger[kMaxClusters];      //No. of cosmic tags associated to this cluster    
    Float_t clucosmicscore_tagger[kMaxClusters];      //Cosmic score associated to this cluster. In the case of more than one tag, the first one is associated.    
    Short_t clucosmictype_tagger[kMaxClusters];       //Cosmic tag type for this cluster.    

Vertex information

Note that only vertex information from the vertex module specified by VertexModuleLabel in the fcl file used to produce this analysis tree will be stored here. (Information about the fcl files used to produce Monte Carlo Challenge (MCC) files can be found in pages available from the MCC index page) The results of several different vertex reconstruction algorithms can be stored. The naming convention is (variablename)_(vertexRecoAlgName).

    Short_t  nvtx_vertexRecoAlgName;                     //number of vertices
    Float_t  vtxx_vertexRecoAlgName[kMaxVertices];     //the X location (in cm) for a given vertex
    Float_t  vtxy_vertexRecoAlgName[kMaxVertices];     //the Y location (in cm) for a given vertex
    Float_t  vtxz_vertexRecoAlgName[kMaxVertices];     //the Z location (in cm) for a given vertex

An example template below shows how to access the vertex nformation in a root macro, assuming the cluster reco algorithm was "cccluster".

    //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the total number of vertices
        for(Int_t j=0;j<nvtx_cccluster;j++) {
             //for e.g., if you want to print the vertex position information, you would do,                
             std::cout<<"\n Vertex X ="<<vtxx_cccluster[j]<<"\t Vertex Y="<<vtxy_cccluster[j] <<"\t Vertex Z="<<vtxz_cccluster[j];
             . . .
        }
    }

Neutrino vertex information

The following variables store only the reconstructed vertices that are thought to be due to a neutrino interaction. They are produced by the pandoraNu algorithm. They are saved selecting only the neutrino PFParticles and retrieving the vertices associated to those. One may have several instances of pandoraNu, for this reason it is possible to save neutrino reconstructed vertices for different pandoraNu algorithms. They will be saved as (variablename)_(pandoraNuInstanceName).

Note that only neutrino vertex information from the pandoraNu module specified by PandoraNuVertexModuleLabel in the fcl file used to produce this analysis tree will be stored here. Also ensure that the fcl parameter SavePandoraNuVertexInfo is set to true.

    Short_t  nnuvtx_ pandoraNuInstanceName;                              //number of neutrino recon vertices
    Float_t  nuvtxId_ pandoraNuInstanceName[kMaxVertices];               //the ID of the vertex
    Float_t  nuvtxx_ pandoraNuInstanceName[kMaxVertices];                //the X location (in cm) for a given vertex
    Float_t  nuvtxy_ pandoraNuInstanceName[kMaxVertices];                //the Y location (in cm) for a given vertex
    Float_t  nuvtxz_ pandoraNuInstanceName[kMaxVertices];                //the Z location (in cm) for a given vertex
    Short_t  nuvtxpdg_ pandoraNuInstanceName[kMaxVertices];              //the pdf of the PFP associated to this vertex (a neutrino PFP)
    Short_t  nuvtxhasPFParticle_ pandoraNuInstanceName[kMaxVertices];    //if this vertex has a PFP (will always be true) (0=false; 1=true)
    Short_t  nuvtxPFParticleID_ pandoraNuInstanceName[kMaxVertices];     //the ID of the associated PFP (a neutrino PFP)

Flash information

Flash information is saved only if the fSaveFlashInfo is set to true in the fhicl file.
The results of several different flash reconstruction algorithms can be stored. The naming convention is (variablename)_(flashRecoAlgName).
The following information is saved for optical flashes:

    Int_t    nfls_flashAlgName;                                    // number of flashes in the event
    Float_t  flsTime_flashAlgName[kMaxFlashes];                    // flash time (in microseconds)
    Float_t  flsPe_flashAlgName[kMaxFlashes];                      // total number of photoelectrons corresponding to the flash
    Float_t  flsPePerOpDet_flashAlgName[kMaxFlashes][kNOpDets];    // photoelectrons per each optical detector (PMT ID, not channel)
    Float_t  flsYcenter_flashAlgName[kMaxFlashes];                 // y center of flash (in cm)
    Float_t  flsZcenter_flashAlgName[kMaxFlashes];                 // z center of flash (in cm)
    Float_t  flsYwidth_flashAlgName[kMaxFlashes];                  // y width of flash (in cm)
    Float_t  flsZwidth_flashAlgName[kMaxFlashes];                  // z width of flash (in cm)

Note that if you are using uboonecode v06_14_00 or below, only one OpFlash data product is saved (the one generated with the OpFlash algorithm).
SimpleFlash and OpFlash are two different flash reconstruction algorithms. For MCC7, OpFlash was the default algorithm, while by MCC8, SimpleFlash was recently developed and has current support. Furthermore, "Beam" and "Cosmic" in the flashRecoAlgName refers to the time window the flashes were reconstructed. The Beam window is ~24 us around the beam spill, while the Cosmic window is much larger before and after the Beam window. For the Cosmic waveforms, we only record 40 samples of them, while the Beam waveforms are recorded entirely.

Shower information

Shower information is saved only if the fSaveShowerInfo is set to true in the fhicl file (by default fSaveShowerInfo is set to true). User is also required to supply a valid ShowerModuleLabel(s) in the anatree_uboone.fcl file to extract the showers. Variables listed below appear in tree with the naming convention (variablename)_(algorithm name). E.g., the nshowers variable appears as nshowers_beziertracker for the "showerrecofuzzy" shower algorithm, nshowers_showerrecopandora for "showerrecopandora" algorithm, etc.

The shower variables are arrays of length nshowers. These are STL vectors in the original TTree, so the number of showers is in principle unlimited; (note: however, be aware that code made by TTree::MakeClass() will have a fixed maximum array length in the retrieved data based on the specific tree on which TTree::MakeClass() was run. See warning in ROOT documentation of TTree::MakeClass() about this effect if MakeClass() is used on just one TTree of a TChain.)

      Short_t  nshowers;                      ///< number of showers
      ShowerData_t<Int_t>    showerID;        ///< Shower ID
      ShowerData_t<Short_t>  shwr_bestplane;  ///< Shower best plane
      ShowerData_t<Float_t>  shwr_length;     ///< Shower length
      ShowerData_t<Float_t>  shwr_startdcosx; ///< X directional cosine at start of shower
      ShowerData_t<Float_t>  shwr_startdcosy; ///< Y directional cosine at start of shower
      ShowerData_t<Float_t>  shwr_startdcosz; ///< Z directional cosine at start of shower
      ShowerData_t<Float_t>  shwr_startx;     ///< startx of shower
      ShowerData_t<Float_t>  shwr_starty;     ///< starty of shower
      ShowerData_t<Float_t>  shwr_startz;     ///< startz of shower
      PlaneData_t<Float_t>   shwr_totEng;     ///< Total energy of the shower per plane
      PlaneData_t<Float_t>   shwr_dedx;       ///< dE/dx of the shower per plane
      PlaneData_t<Float_t>   shwr_mipEng;     ///< Total MIP energy of the shower per plane

Important note: One will notice that, often for some mysteriously failing reason, there is an input event that has no shower branch (expert should debug why this is happening!). In cases like this, the code currently throws a "ProductNotFound" exception. Below is the strategy currently implemented in analysis tree to deal with cases like this:
The analysis tree code now throws an informative exception (instead of just ProductNotFound exception) when a shower is missing (An example of the screen output is shown below) and there will be no shower entry in the output ntuple which means that the user would never know that a shower was missing unless they carefully see the log messages.

Module failed due to an exception
---- ProductNotFound BEGIN
  Showers with input tag 'showerrecofuzzy' were not found in the event. If you really know what you are doing, set AnalysisTree's configuration parameter IgnoreMissingShowers to "true"; the lack of any shower set will be tolerated, and the shower list in the corresponding event will be set to a list of one shower, with an invalid ID.
---- ProductNotFound END

We wanted to provide a provision in analysis tree for the users to know when such things happen and still be able to count the correct number of showers. This can be done by using the fhicl parameter "IgnoreMissingShowers". If the parameter "physics.analyzers.analysistree.IgnoreMissingShowers: true" is set, an error message is printed on the log, and the processing
continues and the event with the missing shower information will have nshowers =1 and some "magic" values for other shower variables:

showerid[0]: -9999
(all the rest of the information is also -999xxx)

The screen output will look like this when the "IgnoreMissingShowers" parameter is set to true,

%MSG-e AnalysisTree:  AnalysisTree:analysistree 09-Aug-2015 14:04:56 CDT  run: 1 subRun: 23 event: 2208
No showers found for input tag 'showerrecofuzzy' ; FILLING WITH FAKE DATA AS FOR USER'S REQUEST
%MSG
%MSG-e AnalysisTree:  AnalysisTree:analysistree 09-Aug-2015 14:04:56 CDT  run: 1 subRun: 23 event: 2208
No showers found for input tag 'showerrecopandora' ; FILLING WITH FAKE DATA AS FOR USER'S REQUEST
%MSG

By default, the fhicl parameter "IgnoreMissingShowers" is set to False, it is recommended to set this to true to be able to keep the correct count of showers and hopefully provide input to the shower algorithm authors or relevant experts that this is happening.

Track information

Variables listed below appear in tree with the naming convention (variablename)_(algorithmname).
E.g., the ntracks variable appears as ntracks_beziertracker for the "beziertracker" tracking algorithm, ntracks_trackkalmanhit for "trackkalmanhit", etc.

The image below shows the naming convention for track variables when the analysis tree is run on multiple trackers:

There are four types of track information: TrackData, PlaneData, HitData, and HitCoordData, as described in the subsections below.

    Short_t  ntracks;             //number of reconstructed tracks

Track "TrackData"

These are arrays of length ntracks. These are STL vectors in the original TTree, so the number of tracks is in principle unlimited; however, be aware that code made by TTree::MakeClass() will have a fixed maximum array length in the retrieved data based on the specific tree on which TTree::MakeClass() was run. (See warning in ROOT documentation of TTree::MakeClass() about this effect if MakeClass() is used on just one TTree of a TChain.)

      TrackData_t<Short_t> trkId;          //Geant track ID of the track
      TrackData_t<Float_t> trkstartx;     // starting x position of the track in cm
      TrackData_t<Float_t> trkstarty;     // starting y position of the track in cm
      TrackData_t<Float_t> trkstartz;     // starting z position of the track in cm
      TrackData_t<Float_t> trkstartd;     // starting distance to boundary in cm
      TrackData_t<Float_t> trkendx;       // ending x position of the track in cm
      TrackData_t<Float_t> trkendy;       // ending y position of the track in cm
      TrackData_t<Float_t> trkendz;       // ending z position of the track in cm
      TrackData_t<Float_t> trkendd;       // ending distance to boundary in cm
      TrackData_t<Float_t> trkACpierceT0;   // t0 per track from anode or cathode piercing tracks (in micro seconds) 
      TrackData_t<Float_t> trkflashT0;   // t0 per track from matching tracks to flashes (in ns)
                                          // using the PhotonCounterT0Matching_module.cc module in larana/T0Finder     
      TrackData_t<Float_t> trktrueT0;    // t0 per track from truth information (in ns) 
                                          // using the MCTruthT0Matching_module.cc  module in larana/T0Finder
      TrackData_t<Float_t> trktheta;      // 3D angle theta w.r.t. the vertex direction in radians
      TrackData_t<Float_t> trkphi;        //  3D angle phi w.r.t. the vertex direction in radians
      TrackData_t<Float_t> trkstartdcosx; //the X direction cosine for track Start
      TrackData_t<Float_t> trkstartdcosy; //the Y direction cosine for track Start
      TrackData_t<Float_t> trkstartdcosz; //the Z direction cosine  for track Start
      TrackData_t<Float_t> trkenddcosx;  //the X direction cosine for track End
      TrackData_t<Float_t> trkenddcosy;  //the Y direction cosine for track End
      TrackData_t<Float_t> trkenddcosz;  //the Z direction cosine for track End
      TrackData_t<Float_t> trkthetaxz;    // Orthogonal projection of theta on to the XZ axis (units in radians)
      TrackData_t<Float_t> trkthetayz;    // Orthogonal projection of theta on to the YZ axis (units in radians)
      TrackData_t<Float_t> trkmom;        // vertex momentum (or start momentum) of the track in GeV
      TrackData_t<Float_t> trklen;        // length of the track in cm (length calculated going trajectory by trajectory)
      TrackData_t<Float_t> trkmomrange;    // track momentum from range using CSDA tables in MeV
      TrackData_t<Float_t> trkmommschi2;   // track momentum from multiple scattering (in GeV)
      TrackData_t<Float_t> trkmommsllhd;    // track momentum from the MCS LLHD method (in GeV) 
      TrackData_t<Short_t> trksvtxid;     // Vertex ID associated with the track start
      TrackData_t<Short_t> trkevtxid;     // Vertex ID associated with the track end
Additional notes on some of the above variables:
  • trkthetaxz variable is defined as trkthetaxz= atan2(dir_start.X(), dir_start.Z()), where dir_start is the vertex direction. Similarly, trkthetayz is defined as trkthetayz=atan2(dir_start.Y(), dir_start.Z()).
  • The trkmomrange variable calculated using sowjanya's track momentum calculator is currently only valid for muons and protons.
  • It is important to note that the FeatureVertex finder (which is ran in MCC5) makes vertices at the cluster-level, so the variables trksvtxid and trkevtxid which show the vertexID associated to a track end or track start, do not come directly from the vertexing or tracking algorithm. The decision of associating a vertex to a track start or track end comes from a simple method written in analysis tree. The distance between Vertex position and track end points is calculated, then the vertex ID is associated to the end point whichever lies closest to the vertex.
Cosmic track tagging information in analysis tree:
  • Analysis tree implements 2 types of cosmic taggers (both at track level): Geometry_based tagger (Tagger) & Light-based tagger (Flash track tagger)
  • For each track, the following cosmic tag information is saved:
    1. number of cosmic tags associated to a given track (almost always there is only one tag per track)
    2. cosmic score (1=definite cosmic; 0=non-cosmic ( or unknown); 0.5= unsure; -1=very likely neutrino).
    3. Type of cosmic tag (-1=unknown type; 0=untagged; 1=crossed 2 Y boundaries; 21=crossed one Y boundary; 100=partially outside drift;
      200=flash incompatible to beam;300=flash matched to beam)

Additional details of cosmic scores and tags can be found here:
https://cdcvs.fnal.gov/redmine/projects/lardata/repository/revisions/develop/entry/AnalysisBase/CosmicTag.h

      TrackData_t<Short_t> trkncosmictags_tagger;        // number of cosmic tags associated to a given track using Geometry based tagger
      TrackData_t<Float_t> trkcosmicscore_tagger;         // assigned cosmic score by Geometry based tagger
      TrackData_t<Short_t> trkcosmictype_tagger;          // Type of cosmic tag assigned by the Geometry-based tagger
      TrackData_t<Short_t> trkncosmictags_flashmatch; // number of cosmic tags associated to a given track using Flash tagger
      TrackData_t<Float_t> trkcosmicscore_flashmatch;  // cosmic score assigned by Flash tagger
      TrackData_t<Short_t> trkcosmictype_flashmatch;   // cosmic tag type assigned by Flash tagger

The cosmic tag information is saved per tracker in the analysis tree. For example, if the analysis tree is run on trackkalmanhit and beziertracker, then the taggers are saved for both trackers following the convention like other tracking variables (variablename)_(algorithmname). In the tree viewer image below, the cosmic tag variables in the variable list are high lighted to show how these are stored:

The example template below shows how to access track information inside a ROOT macro:

  //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the number of tracks (depends on the tracking algorithm, ntracks_trackkalmanhit or ntracks_trackkalsps or ntracks_beziertracker etc.)
        // For now, lets assume here we want to look at trackkalmanhit tracker
        for(Int_t j=0;j<ntracks_trackkalmanhit;j++) {
             ...
            //for example, print trackID and track length of the track
            std::cout<<"\n"<<trackId_trackkalmanhit[j]<<"\t"<<trklen_trackkalmanhit[j];
           //get the difference between momentum calculated from range and multiple scattering
           float mom_diff = trkmomrange[j] - trkmommschi2[j];
           ....
           //check if the track is a cosmic
           if (trkcosmicscore_tagger_trackkalmanhit[j] == 1)
                         std::cout<<"\n Geometry tagged this track as cosmic";
           else
                         std::cout<<"\n most likey not a cosmic!";
           if (trkcosmicscore_fashmatch_trackkalmanhit[j] == 1)
                         std::cout<<"\n Flash matcher tagged this track as cosmic";
           .....

       } 
    }

Track "PlaneData"

These appear as 2-index arrays of dimension [ntracks][3]. The first index selects the track, and the second index is for the plane. The plane number 0 corresponds to U, 1 corresponds to V and 2 corresponds to W (or Z). See note in "TrackData" section above regarding maximum number of tracks. The Track "Plane Data" is also saved for each tracker, following the (variable name)_(algorithm) convention as mentioned above.

      PlaneData_t<Float_t>    trkke;       //kinetic energy of the track in each plane from calorimetry in MeV
      PlaneData_t<Float_t>    trkrange;    //total range of the track in each plane from calorimetry in cm
      PlaneData_t<Int_t>      trkidtruth;  //true geant trackid assigned to the reconstructed track in a given plane
                                           //(more in the additional information section below)
      PlaneData_t<Short_t>    trkorigin;   //mc origin type: 0=unknown, 1=beam neutrino, 2=cosmic, 3=supernova neutrino, 4=single particle 
                                           //(more in the additional information section below)
      PlaneData_t<Int_t>      trkpdgtruth; //true pdg code associated to the particle with GEANT trackID  @trkidtruth@
      PlaneData_t<Float_t>    trkefftruth; //completeness of the track in a given plane (defined in the additional information section below)
      PlaneData_t<Float_t>    trkpurtruth; //purity of track in a given plane (defined in the additional information section below)
      PlaneData_t<Float_t>    trkpitchc;   //track pitch in the collection plane from Calorimetry (in cm)
      PlaneData_t<Short_t>    ntrkhits;    //number of hits associated to the track in that plane
      //PID information
      PlaneData_t<Int_t> trkpidpdg;         //pdg code assigned by the @Chi2PIDAlg@ PID algorithm
      PlaneData_t<Float_t> trkpidchi;       //Minimum reduced Chi2 assigned by the @Chi2PIDAlg@ PID algorithm
      PlaneData_t<Float_t> trkpidchipr;     //Chi-square for proton assigned by the @Chi2PIDAlg@ PID algorithm
      PlaneData_t<Float_t> trkpidchika;     //Chi-square for kaon assigned by the @Chi2PIDAlg@ PID algorithm
      PlaneData_t<Float_t> trkpidchipi;     //Chi-square for pion assigned by the @Chi2PIDAlg@ PID algorithm
      PlaneData_t<Float_t> trkpidchimu;     //Chi-square for muon assigned by the @Chi2PIDAlg@ PID algorithm
      PlaneData_t<Float_t> trkpidpida;      //PIDA value assigned to the track by Bruce's @PIDA@ algorithm
      TrackData_t<Short_t> trkpidbestplane; //Currently, this is defined as the plane with most hits     
Additional information for some of the above variables:
  • trkidtruth is assigned in the following way: extract all the hits corresponding to the track in a given plane, associate hits to true particles using the simIDE object and find which particle's signal (hit) deposited the most energy for this track in that plane. The GEANT trackID associated to this particle is then stored in the trkidtruth variable for that plane. More information on the simIDE object can be found here: http://nusoft.fnal.gov/larsoft/doxsvn/html/SimChannel_8h_source.html#l00030
  • The trkefftruth variable is defined as the ratio of the energy of the particle that contributed most to this track in a given plane to the total true energy associated to this trackID in that plane
  • The trkpurtruth variable is defined as the ratio of the energy of the particle that contributed most to this track in a given plane to the total energy coming from all particles that contribute to this track in that plane

An example ROOT macro below shows how to access the Track Plane data information:

  //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the number of tracks (depends on the tracking algorithm, ntracks_trackkalmanhit or ntracks_trackkalsps or ntracks_beziertracker etc.)
        // For now, lets assume we want to look at trackkalmanhit tracker
        for(Int_t j=0;j<ntracks_trackkalmanhit;j++) {
             ...
            //loop over the number of wire planes (which is always 3: 0 is U plane, 1 is V plane, and 2 is W (or Z) plane
            for(Int_t k=0;k<3;k++) {
                  //lets look at number of hits associated to the track in a given plane
                  std::cout<<"\n"<<ntrkhits[j][k];
                  //lets pick tracks with kinetic energy greater than 0.5 GeV
                  if (trkke[j][k]> 500) {
                         std::cout<<"\n Range of this track is="<<trkrange[j][k];
                         //now lets see if it is a muon using the PID information
                         if (trkpidpdg[j][k]==13)
                               std::cout<<"\n Found a muon with energy greater than 20 GeV with a track length of"<<trkrange[j][k];
                  }
             }//end of wire-plane loop                
           .....
        }//end of ntracks 
    }//end of ntuple entries

Track "HitData" and "HitCoordData"

The HitData variables appear as 2-index arrays of dimension [ntracks][kMaxTrackHits], where kMaxTrackHits is a fixed number set at compile time by a constant in AnalysisTree_module.cc -- kMaxTrackHits=2000 since v01_01_00 through at least v03_03_00. (Replacing the fixed-length array with a vector is on a to-do list, so this may change.)

The HitCoordData variables appear as 3-index arrays of dimension [ntracks][kMaxTrackHits][3]. Currently the only HitCoordData variable is trkxyz. The first index selects the track, and the second index is for the hit on the track. The third index for trkxyz selects the coordinate (X, Y or Z).

The actual number of hits on a given track is given by ntrkhits for that track. (See TrackData section above. See also the note in "TrackData" section above regarding maximum number of tracks.)

      HitData_t<Float_t>      trkdedx;  // dedx of the track hit in a given plane (in MeV/cm)
      HitData_t<Float_t>      trkdqdx;  // dqdx of the track hit in a given plane (in ADC/cm)
      HitData_t<Float_t>      trkresrg; // range from end of track (in cm)
      HitCoordData_t<Float_t> trkxyz;   // (X,Y,Z) position of the track hit in a given plane (in cm)

An example root macro template below shows how to access "Hit" and "HitCoord" data:

  //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the number of tracks (depends on the tracking algorithm, ntracks_trackkalmanhit or ntracks_trackkalsps or ntracks_beziertracker etc.)
        // For now, lets assume we want to look at trackkalmanhit tracker
        for(Int_t j=0;j<ntracks_trackkalmanhit;j++) {
                   ......
                  //loop ONLY over the number of hits associated to kalmanhit tracker in plane 2 (or Z plane). 
                  //This doesn't require running a loop over wire planes after the ntracks loop
                  //require at least 100 hits in the track in the Z plane, otherwise continue
                  if (ntrkhits_trackkalmanhit[j][2] < 100)      continue;      
                  for(Int_t k=0;k<ntrkhits_trackkalmanhit[j][2];k++){                                     
                        //for tracks with > 100 hits in Z plane, lets print "dedx" of hit as a function of hit X position
                        std::cout<<"\n Hit X ="<< trkxyz_trackkalmanhit[j][2][k][0]<<"\t Hit dedx ="<<trkdedx_trackkalmanhit[j][2][k];
                            //One can convert the Hit X to time by dividing it by drift velocity, 160 cm/ms
                            float hit_time = trkxyz_trackkalmanhit[j][2][k][0]/160;
                  }              
           .....
        }//end of ntracks 
    }//end of ntuple entries

As an example, the graph below shows the Hit dedx drawn as a function of Hit time. The graph is generated using the analysis tree ntuples ( by roughly following the above template). The plot also shows a polynomial curve fitted to the dedx distribution:

MC Truth information

This set of truth information specifically holds neutrino-interaction specific information which would be very useful, for e.g., if you are doing a cross-section analysis. These variables are arrays with up to kMaxTruth(=10) entries per event. This array structure allows for recording multiple neutrinos to interacting per spill. Multiple neutrinos are very common for studies which include the dirt geometry because of the large mass and volume, but less common for simulations only considering the detector volume. In the very rare case, where there are more than 10 neutrinos interacting per spill, we just save the interaction information for the first 10 neutrinos.

Below are the list of variables saved in analysis tree for neutrino interaction information:


    Int_t    mcevts_truth;               //number of neutrino interactions in the spill
    Int_t    nuPDG_truth[kMaxTruth];     //neutrino PDG code (nue=12; anti-nue=-12; numu=14; anti-numu=-14; nutau=16; anti-nutau=-16)
    Int_t    ccnc_truth[kMaxTruth];      //neutrino interaction type: 0=Charged current (CC), 1=Neutral current (NC)
    Int_t    mode_truth[kMaxTruth];      //neutrino nucleus 0=Quasi-elastic or Elastic, 1=Resonant (RES), 2=DIS, 3=Coherent production
    Float_t  enu_truth[kMaxTruth];       //true neutrino energy in GeV
    Float_t  Q2_truth[kMaxTruth];        //Momentum transfer squared in GeV^2
    Float_t  W_truth[kMaxTruth];         //hadronic invariant mass in GeV
    Float_t  X_truth[kMaxTruth];         //Bjorken X
    Float_t  Y_truth[kMaxTruth];         //Inelasticity
    Int_t    hitnuc_truth[kMaxTruth];    //neutrino scattering off of which nucleon (proton or neutron); holds the pdg of the nucleon
    Float_t  nuvtxx_truth[kMaxTruth];    //neutrino vertex x in cm
    Float_t  nuvtxy_truth[kMaxTruth];    //neutrino vertex y in cm
    Float_t  nuvtxz_truth[kMaxTruth];    //neutrino vertex z in cm
    Float_t  nu_dcosx_truth[kMaxTruth];  //neutrino x directional cosine for neutrino track Start
    Float_t  nu_dcosy_truth[kMaxTruth];  //neutrino y directional cosine for neutrino track Start
    Float_t  nu_dcosz_truth[kMaxTruth];  //neutrino z directional cosine for neutrino track Start
    Float_t  lep_mom_truth[kMaxTruth];   //lepton momentum in GeV
    Float_t  lep_dcosx_truth[kMaxTruth]; //lepton x directional cosine for lepton track Start
    Float_t  lep_dcosy_truth[kMaxTruth]; //lepton y directional cosine for lepton track Start
    Float_t  lep_dcosz_truth[kMaxTruth]; //lepton z directional cosine for lepton track Start

The example template below shows how to access this truth information in a ROOT macro:

//loop over the ntuple entries
    for(Int_t i = 0; i < ntuple_entries; i++) {
        . . .
        //loop over neutrino interactions
        for(Int_t n = 0; n < mcevts_truth && n < 10; n++) {      
          //lets find charged current numu events
          if (ccnc_truth[n] == 0 && nuPDG_truth[n] == 14)
                std::cout<<"\n found charged current nu mu events";
              .....
        }//end of neutrino interaction loop   
    }//end of ntuple entries

Flux information

The Flux variables hold information per each neutrino interaction. Accordingly, each entry in the MC Truth (_truth) variable collection has a corresponding _flux variable collection. A brief description of variables stored in the analysis tree are listed below:

    Float_t  vx_flux[kMaxTruth];          //X position of hadron/muon decay (cm)
    Float_t  vy_flux[kMaxTruth];          //Y position of hadron/muon decay (cm)
    Float_t  vz_flux[kMaxTruth];          //Z position of hadron/muon decay (cm)
    Float_t  pdpx_flux[kMaxTruth];        //Parent X momentum at decay point (GeV)
    Float_t  pdpy_flux[kMaxTruth];        //Parent Y momentum at decay point (GeV)
    Float_t  pdpz_flux[kMaxTruth];        //Parent Z momentum at decay point (GeV)
    Float_t  ppdxdz_flux[kMaxTruth];      //Parent dxdz direction at production
    Float_t  ppdydz_flux[kMaxTruth];      //Parent dydz direction at production
    Float_t  pppz_flux[kMaxTruth];        //Parent Z momentum at production (GeV)

    Int_t    ptype_flux[kMaxTruth];        //Parent GEANT code particle ID
    Float_t  ppvx_flux[kMaxTruth];        //Parent production vertex X (cm)
    Float_t  ppvy_flux[kMaxTruth];        //Parent production vertex Y (cm)
    Float_t  ppvz_flux[kMaxTruth];        //Parent production vertex Z (cm)
    Float_t  muparpx_flux[kMaxTruth];     //Muon neutrino parent production vertex X (cm)
    Float_t  muparpy_flux[kMaxTruth];     //Muon neutrino parent production vertex Y (cm)
    Float_t  muparpz_flux[kMaxTruth];     //Muon neutrino parent production vertex Z (cm)
    Float_t  mupare_flux[kMaxTruth];      //Muon neutrino parent energy (GeV)

    Int_t    tgen_flux[kMaxTruth];        //Parent generation in cascade. (1 = primary proton, 
                                         //2=particles produced by proton interaction, 3 = particles 
                                         //produced by interactions of the 2's, ...
    Int_t    tgptype_flux[kMaxTruth];     //Type of particle that created a particle flying of the target  
    Float_t  tgppx_flux[kMaxTruth];       //X Momentum of a particle, that created a particle that flies 
                                            //off the target, at the interaction point. (GeV)
    Float_t  tgppy_flux[kMaxTruth];       //Y Momentum of a particle, that created a particle that flies  
                                          //off the target, at the interaction point. (GeV)
    Float_t  tgppz_flux[kMaxTruth];       //Z Momentum of a particle, that created a particle that flies
                                          //off the target, at the interaction point. (GeV)
    Float_t  tprivx_flux[kMaxTruth];      //Primary particle interaction vertex X (cm)  
    Float_t  tprivy_flux[kMaxTruth];      //Primary particle interaction vertex Y (cm)  
    Float_t  tprivz_flux[kMaxTruth];      //Primary particle interaction vertex Z (cm) 
    Float_t  dk2gen_flux[kMaxTruth];      //distance from decay to ray origin (cm) 
    Float_t  gen2vtx_flux[kMaxTruth];     //distance from ray origin to event vtx (cm)

    Float_t  tpx_flux[kMaxTruth];        //Px of parent particle leaving BNB/NuMI target (GeV)
    Float_t  tpy_flux[kMaxTruth];        //Py of parent particle leaving BNB/NuMI target (GeV)
    Float_t  tpz_flux[kMaxTruth];        //Pz of parent particle leaving BNB/NuMI target (GeV)
    Int_t    tptype_flux[kMaxTruth];     //Type of parent particle leaving BNB/NuMI target   

The template below shows how to access the flux information in a root macro:

    unsigned int count =0;
    //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        //loop over neutrino interactions
        for(Int_t n = 0; n < mcevts_truth && n < 10; n++) {
            . . .
            //as an example, lets count how many particles came from K meson decays in the beam and look at the Kaon Pz distribution
            if (abs(tptype_flux[n] ==321){
                  count++;  
                  std::cout<<"\n"<<tpx_flux;
                  .....
            }
        }//end loop over neutrino interactions
    }//end of ntuple entries

GENIE information

MicroBooNE uses GENIE as the neutrino MC generator. The genie variables saved in the analysis tree are 1-index variables of dimension [genie_no_primaries]. Only the genie primary particles from the first neutrino recorded in MC Truth is recorded here. To get the primary particle information for any remaining neutrino interaction recorded in MC Truth AnalysisTree variables, please see the section on GEANT Information about the MCTruthIndex variable.

Below are the list of variables that are saved in the analysis tree for GENIE:

    Int_t     genie_no_primaries;              //number of primary particles generated by GENIE 
    std::vector<Int_t>    genie_primaries_pdg;    //particle type (pdg) of the GENIE particle
    std::vector<Float_t>  genie_Eng;           //Energy of the GENIE particle in GeV
    std::vector<Float_t>  genie_Px;            //Px of the GENIE particle in GeV
    std::vector<Float_t>  genie_Py;            //Py of the GENIE particle in GeV
    std::vector<Float_t>  genie_Pz;            //Pz of the GENIE particle in GeV
    std::vector<Float_t>  genie_P;             //Magnitude of the momentum vector of the GENIE particle in GeV
    std::vector<Int_t>    genie_status_code;   //particle status code of the GENIE particle 
                                               //  only genie_status_code==1 particles are to be tracked in detector MC,
                                               //  see genie.hepforge.org/doxygen/html/GHepStatus_8h_source.html
    std::vector<Float_t>  genie_mass;          //mass of the GENIE particle in GeV
    std::vector<Int_t>    genie_trackID;       //trackID of the GENIE particle (different from the GEANT-assigned track ID)
    std::vector<Int_t>    genie_ND;            //number of daughters of the GENIE particle 
    std::vector<Int_t>    genie_mother;        //mother trackID of the GENIE particle

An example template below shows how to access GENIE information in a ROOT macro.
Be aware that for general studies you will more likely want to use only GENIE particles flagged as "stable final state", i.e. those particles that have genie_status_code == 1.

    //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the total number of @GENIE@ particles
        for(Int_t j=0;j<genie_no_primaries; j++) {
             ...
            // We want stable final states. Those are the particles used for detector simulation
            if (genie_status_code != 1)  continue;

            //Lets select muons with an energy greater than 500 MeV  
            if  (abs(genie_pdg[j])==13 && genie_Eng[j]>=0.5)  
                    std::cout<<"\n Found a muon with energy greater than 500 MeV!";
            .......       
        }//end of genie loop 
    }//end of ntuple entries

Cosmic CRY information

MicroBooNE uses two geneators, CRY and CORSIKA, to generate cosmic ray fluxes. CRY variables saved in the analysis tree are 1-index arrays of dimension [cry_no_primaries]. Since it is not very common that one would go down to the CRY-level information, this information is only extracted and saved in the tree if the flag SaveCryInfo is set to true in the fhicl file. Note that even though the labels begin with "cry_" they do in fact apply to whichever generator was run.

Below are the list of CRY variables saved in the analysis tree along with a brief description of what those are:

    Int_t     cry_no_primaries;                  //number of primary particles generated by CRY 
    std::vector<Int_t>    cry_primaries_pdg;            //pdg of the CRY particle
    std::vector<Float_t>  cry_Eng;           //Energy of the CRY particle in GeV
    std::vector<Float_t>  cry_Px;             //Momentum in X direction of the CRY particle in GeV
    std::vector<Float_t>  cry_Py;             //Momentum in Y direction of the CRY particle in GeV
    std::vector<Float_t>  cry_Pz;             //Momentum in Z direction of the CRY particle in GeV
    std::vector<Float_t>  cry_P;               //Magnitude of the Momentum vector of the CRY particle in GeV
    std::vector<Float_t>  cry_StartPointx;    //X StartPoint of the CRY particle in cm
    std::vector<Float_t>  cry_StartPointy;    //Y StartPoint of the CRY particle in cm
    std::vector<Float_t>  cry_StartPointz;    //Z StartPoint of the CRY particle in cm
    std::vector<Float_t>  cry_StartPointt;    //Start time of the CRY particle in ns
    std::vector<Int_t>    cry_status_code;    //particle status code of the CRY particle
    std::vector<Float_t>  cry_mass;           //Mass of the CRY particle in GeV  
    std::vector<Int_t>    cry_trackID;         //trackID of the CRY particle (different from the GEANT assigned trackID) 

There seem to be some confusion where one understands that cry_no_primaries refer to the primary particle spectra that CRY/CORSIKA uses to generate air-showers. Thats not the case, these are all the air-shower particles that CRY generates and are input to GEANT for further simulation. And like in the case of any other generator, in the context of GEANT, these particles are referred to as primary particles.

An example template below shows how to access CRY/CORSIKA information in a ROOT macro:

    //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the total number of CRY particles
        for(Int_t j=0;j<cry_no_primaries; j++) {
             ...
            //see how many anti-muons CRY generated with an energy greater than 10 GeV  
            if  (cry_pdg[j]==-13 && cry_Eng[j]>=10)  
                    std::cout<<"\n Found an anti-muon with energy greater than 10 GeV!";
            .......       
        }//end of cry loop 
    }//end of ntuple entries

MC Shower information

The MC Shower variables are 1-index arrays of dimension [no_mcshowers]. One can turn on/off saving MC Shower information in the analysis tree tuples using the switch "SaveMCShowerInfo" from the fhicl file. User also needs to pass the appropriate mc shower module label "MCShowerModuleLabel" (currently set to "mcreco") to successfully extract and store mc shower information. All showers irrespective of whether they deposit energy inside the detector or not are saved to give a proper count of the total number of showers in the full event window. The variable "mcshwr_isEngDeposited" tells the user whether a given shower deposited energy inside the detector or not. A value of one means yes and zero means no. The "Combined Energy Deposition" variables (see below) will have valid values only if a shower deposits energy inside the detector. Otherwise, the usual place holder values "-99999." will populate these particular variables.

 //MC Shower information
    Int_t     no_mcshowers;                         //number of MC Showers in this event.
    //MC Shower particle information
    std::vector<Int_t>       mcshwr_origin;        //MC Shower origin information ( //mc origin type: 0=unknown, 1=beam neutrino, 2=cosmic, 3=supernova neutrino, 4=single particle)
    std::vector<Int_t>       mcshwr_pdg;        //MC Shower particle PDG code.   
    std::vector<Int_t>       mcshwr_TrackId;        //MC Shower particle G4 track ID.
    std::vector<std::string> mcshwr_Process;        //MC Shower particle's creation process. 
    std::vector<Float_t>     mcshwr_startX;        //MC Shower particle G4 startX (cm)
    std::vector<Float_t>     mcshwr_startY;        //MC Shower particle G4 startY (cm)
    std::vector<Float_t>     mcshwr_startZ;        //MC Shower particle G4 startZ (cm)
    std::vector<Float_t>     mcshwr_endX;        //MC Shower particle G4 endX (cm)
    std::vector<Float_t>     mcshwr_endY;        //MC Shower particle G4 endY (cm)
    std::vector<Float_t>     mcshwr_endZ;        //MC Shower particle G4 endZ (cm)
    std::vector<Float_t>    mcshwr_CombEngX;        //MC Shower Combined energy deposition information, Start Point X Position (cm)
    std::vector<Float_t>    mcshwr_CombEngY;        //MC Shower Combined energy deposition information, Start Point Y Position (cm)
    std::vector<Float_t>    mcshwr_CombEngZ;        //MC Shower Combined energy deposition information, Start Point Z Position (cm)
    std::vector<Float_t>     mcshwr_CombEngPx;        //MC Shower Combined energy deposition information, Momentum X direction (MeV/c)
    std::vector<Float_t>     mcshwr_CombEngPy;        //MC Shower Combined energy deposition information, Momentum X direction (MeV/c)
    std::vector<Float_t>     mcshwr_CombEngPz;        //MC Shower Combined energy deposition information, Momentum X direction (MeV/c)
    std::vector<Float_t>     mcshwr_CombEngE;        //MC Shower Combined energy deposition information, Energy (MeV/c)
    std::vector<Int_t>       mcshwr_isEngDeposited;  //tells whether if this shower deposited energy in the detector or not.
                                                            //yes = 1; no =0;    
    //MC Shower mother information
    std::vector<Int_t>       mcshwr_Motherpdg;       //MC Shower's mother PDG code. 
    std::vector<Int_t>       mcshwr_MotherTrkId;     //MC Shower's mother G4 track ID.
    std::vector<std::string> mcshwr_MotherProcess;   //MC Shower's mother creation process. 
    std::vector<Float_t>     mcshwr_MotherstartX;    //MC Shower's mother  G4 startX  (cm)
    std::vector<Float_t>     mcshwr_MotherstartY;    //MC Shower's mother  G4 startY  (cm)
    std::vector<Float_t>     mcshwr_MotherstartZ;    //MC Shower's mother  G4 startZ  (cm)
    std::vector<Float_t>     mcshwr_MotherendX;         //MC Shower's mother  G4 endX (cm)
    std::vector<Float_t>     mcshwr_MotherendY;         //MC Shower's mother  G4 endY (cm)  
    std::vector<Float_t>     mcshwr_MotherendZ;         //MC Shower's mother  G4 endZ (cm)
    //MC Shower ancestor information
    std::vector<Int_t>       mcshwr_Ancestorpdg;       //MC Shower's ancestor PDG code. 
    std::vector<Int_t>       mcshwr_AncestorTrkId;     //MC Shower's ancestor G4 track ID.
    std::vector<std::string> mcshwr_AncestorProcess;   //MC Shower's ancestor creation process. 
    std::vector<Float_t>     mcshwr_AncestorstartX;    //MC Shower's ancestor  G4 startX (cm)
    std::vector<Float_t>     mcshwr_AncestorstartY;    //MC Shower's ancestor  G4 startY (cm)
    std::vector<Float_t>     mcshwr_AncestorstartZ;    //MC Shower's ancestor  G4 startZ (cm)
    std::vector<Float_t>     mcshwr_AncestorendX;      //MC Shower's ancestor  G4 endX (cm)
    std::vector<Float_t>     mcshwr_AncestorendY;      //MC Shower's ancestor  G4 endY (cm) 
    std::vector<Float_t>     mcshwr_AncestorendZ;      //MC Shower's ancestor  G4 endZ (cm) 

MC Track information

The MC Track variables are 1-index arrays of dimension [no_mctracks]. One can turn on/off saving MC Track information in the analysis tree tuples using the switch "SaveMCTrackInfo" from the fhicl file. User also needs to pass the appropriate mc track module label "MCTrackModuleLabel" (currently set to "mcreco") to successfully extract and store mc track information. The meaning of the variables with "_drifted" suffix is explained in the next section.

 //MC Track information
    Int_t     no_mctracks;                         //number of MC Tracks in this event.
    //MC Track particle information
    std::vector<Int_t>       mctrk_origin;        //MC Track origin information ( //mc origin type: 0=unknown, 1=beam neutrino, 2=cosmic, 3=supernova neutrino, 4=single particle)
    std::vector<Int_t>       mctrk_pdg;        //MC Track particle PDG code.   
    std::vector<Int_t>       mctrk_TrackId;        //MC Track particle G4 track ID.
    std::vector<std::string> mctrk_Process;        //MC Track particle's creation process. 
    std::vector<Float_t>     mctrk_startX;        //MC Track particle G4 startX (cm)
    std::vector<Float_t>     mctrk_startY;        //MC Track particle G4 startY (cm)
    std::vector<Float_t>     mctrk_startZ;        //MC Track particle G4 startZ (cm)
    std::vector<Float_t>     mctrk_endX;        //MC Track particle G4 endX (cm)
    std::vector<Float_t>     mctrk_endY;        //MC Track particle G4 endY (cm)
    std::vector<Float_t>     mctrk_endZ;        //MC Track particle G4 endZ (cm)
    std::vector<Float_t>     mctrk_startX_drifted;        //MC track particle first step in drifted TPC volume x (cm)
    std::vector<Float_t>     mctrk_startY_drifted;        //MC track particle first step in drifted TPC volume y (cm)
    std::vector<Float_t>     mctrk_startZ_drifted;        //MC track particle first step in drifted TPC volume z (cm)
    std::vector<Float_t>     mctrk_endX_drifted;          //MC track particle last step in drifted TPC volume x (cm)
    std::vector<Float_t>     mctrk_endY_drifted;          //MC track particle last step in drifted TPC volume y (cm)
    std::vector<Float_t>     mctrk_endZ_drifted;          //MC track particle last step in drifted TPC volume z (cm)
    std::vector<Float_t>     mctrk_len_drifted;           //MC track length within drifted TPC volume (cm)
    std::vector<Float_t>     mctrk_p_drifted;             //MC track momentum at start point in drifted TPC volume (MeV/c)
    std::vector<Float_t>     mctrk_px_drifted;            //MC track x momentum at start point in drifted TPC volume (MeV/c)
    std::vector<Float_t>     mctrk_py_drifted;            //MC track y momentum at start point in drifted TPC volume (MeV/c)
    std::vector<Float_t>     mctrk_pz_drifted;            //MC track z momentum at start point in drifted TPC volume (MeV/c)

    //MC Track mother information
    std::vector<Int_t>       mctrk_Motherpdg;       //MC Track's mother PDG code. 
    std::vector<Int_t>       mctrk_MotherTrkId;     //MC Track's mother G4 track ID.
    std::vector<std::string> mctrk_MotherProcess;   //MC Track's mother creation process. 
    std::vector<Float_t>     mctrk_MotherstartX;    //MC Track's mother  G4 startX  (cm)
    std::vector<Float_t>     mctrk_MotherstartY;    //MC Track's mother  G4 startY  (cm)
    std::vector<Float_t>     mctrk_MotherstartZ;    //MC Track's mother  G4 startZ  (cm)
    std::vector<Float_t>     mctrk_MotherendX;         //MC Track's mother  G4 endX (cm)
    std::vector<Float_t>     mctrk_MotherendY;         //MC Track's mother  G4 endY (cm)  
    std::vector<Float_t>     mctrk_MotherendZ;         //MC Track's mother  G4 endZ (cm)
    //MC Track ancestor information
    std::vector<Int_t>       mctrk_Ancestorpdg;       //MC Track's ancestor PDG code. 
    std::vector<Int_t>       mctrk_AncestorTrkId;     //MC Track's ancestor G4 track ID.
    std::vector<std::string> mctrk_AncestorProcess;   //MC Track's ancestor creation process. 
    std::vector<Float_t>     mctrk_AncestorstartX;    //MC Track's ancestor  G4 startX (cm)
    std::vector<Float_t>     mctrk_AncestorstartY;    //MC Track's ancestor  G4 startY (cm)
    std::vector<Float_t>     mctrk_AncestorstartZ;    //MC Track's ancestor  G4 startZ (cm)
    std::vector<Float_t>     mctrk_AncestorendX;      //MC Track's ancestor  G4 endX (cm)
    std::vector<Float_t>     mctrk_AncestorendY;      //MC Track's ancestor  G4 endY (cm) 
    std::vector<Float_t>     mctrk_AncestorendZ;      //MC Track's ancestor  G4 endZ (cm) 

GEANT (MC Particle) information

The GEANT variables are 1-index arrays of dimension [MaxGeantParticles]. Analysis tree saves an extensive list of GEANT-level information as listed below. Information for all primary particles are saved. In the case of secondary particles, an energy cut of 10 MeV is imposed in order to reduce the very large size of arrays. One can configure the energy cut on G4 secondary particles using the "G4MinE" fcl parameter.

There are two sets of variables that describe particle kinematics related to the TPC geometry. The first set is flagged by inTPCActive and labelled with "_tpcAV" and simply states whether the particle intersected the physical TPC active volume (the geometry calls this volTPCActive). For instance the StartPointx_tpcAV variable is the x coordinate of the first trajectory point that was found inside the TPC. These tpcAV variables are useful for studies where the actual particle trajectory is needed.

The other set is flagged by inTPCDrifted and labelled with "_drifted". The flag defines whether or not the particle intersects the time projected TPC active volume. Or, put another way, this flag defines whether or not the ionization deposited by a particle would reach the anode plane during the readout window. The StartPointx_drifted variable here is the projected x coordinate (trajectory's x coordinate + time multiplied by the drift velocity). The drift velocity and the inTPCDrifted projected x range are pulled from util::DetectorProperties. These variables are useful for comparison with reconstruction.

Note: previous versions of AnalysisTree have used both definitions of inTPCActive defined above. Most recently the inTPCActive was equivalent to inTPCDrifted.

    Int_t     no_primaries;                  //number of primary geant particles
    Int_t     geant_list_size;               //list of all geant particles
    Int_t     geant_list_size_in_tpcAV;      //total number of geant particles crossing TPC boundaries (requires at least two trajectory points inside the TPC active volume)
    std::vector<Int_t>    pdg;               //pdgcode of the particle
    std::vector<Int_t>    status;            //status code of secondary particles assigned by geant
    std::vector<Float_t>  Eng;               //Energy of the particle at its start in GeV
    std::vector<Float_t>  EndE;             //Energy of the particle at its end point in GeV                   
    std::vector<Float_t>  Mass;             // Rest mass of the particle in GeV
    std::vector<Float_t>  Px;                //Particle start momentum in X direction in GeV
    std::vector<Float_t>  Py;                //Particle start momentum in X direction in GeV 
    std::vector<Float_t>  Pz;                //Particle start momentum in X direction in GeV
    std::vector<Float_t>  P;                 //Scalar momentum of the particle in GeV
    std::vector<Float_t>  StartPointx;      //starting 'X' position of the particle in cm at its origin
    std::vector<Float_t>  StartPointy;      //starting 'Y' position of the particle in cm at its origin
    std::vector<Float_t>  StartPointz;      //starting 'Z' position of the particle in cm at its origin
    std::vector<Float_t>  StartT;           //start time of the particle (time corresponding to the first trajectory point) in ns
    std::vector<Float_t>  EndT;             //end time of the particle (time corresponding to the last trajectory point) in ns       
    std::vector<Float_t>  EndPointx;        //ending 'X' position of the particle in cm
    std::vector<Float_t>  EndPointy;        //ending 'Y' position of the particle in cm
    std::vector<Float_t>  EndPointz;        //ending 'Z' position of the particle in cm
    std::vector<Float_t>  theta;           //3D angle theta w.r.t. the vertex (start) direction in radians    
    std::vector<Float_t>  phi;              //3D angle phi w.r.t. the vertex (start) direction in radians
    std::vector<Float_t>  theta_xz;           //Orthogonal projection of theta on to the XZ axis (units in radians)  
    std::vector<Float_t>  theta_yz;          //Orthogonal projection of theta on to the YZ axis (units in radians) 
    std::vector<Float_t>  pathlen;                 //Particle length in cm calculated trajectory-by-trajectory
    std::vector<Int_t>    inTPCActive;            //Did the particle pass through TPC physical boundaries? (Yes=1; No=0)
    std::vector<Float_t>  StartPointx_tpcAV;    //'X' position of particle's first step inside the TPC boundaries in cm
    std::vector<Float_t>  StartPointy_tpcAV;    //'Y' position of particle's first step inside the TPC boundaries in cm
    std::vector<Float_t>  StartPointz_tpcAV;    //'Z' position of particle's first step inside the TPC boundaries in cm
    std::vector<Float_t>  StartT_tpcAV;            //Time associated with particle's first step inside the TPC [ns]
    std::vector<Float_t>  StartE_tpcAV;            //Energy of the particle at the first step inside the TPC [GeV]
    std::vector<Float_t>  StartP_tpcAV;            //Momentum of the particle at the first step inside the TPC [GeV/c]
    std::vector<Float_t>  StartPx_tpcAV;        //X component of P at the first step inside the TPC [GeV/c]
    std::vector<Float_t>  StartPy_tpcAV;        //Y component of P at the first step inside the TPC [GeV/c]
    std::vector<Float_t>  StartPz_tpcAV;        //Z component of P at the first step inside the TPC [GeV/c]
    std::vector<Float_t>  EndPointx_tpcAV;        //'X' position of particle's last step inside the TPC boundaries in cm
    std::vector<Float_t>  EndPointy_tpcAV;        //'Y' position of particle's first step inside the TPC boundaries in cm
    std::vector<Float_t>  EndPointz_tpcAV;        //'Z' position of particle's first step inside the TPC boundaries in cm
    std::vector<Float_t>  EndT_tpcAV;            //Time associated with particle's last step inside the TPC [ns]
    std::vector<Float_t>  EndE_tpcAV;            //Energy of the particle at the last step inside the TPC [GeV]
    std::vector<Float_t>  EndP_tpcAV;            //Momentum of the particle at its last step inside the TPC [GeV/c]
    std::vector<Float_t>  EndPx_tpcAV;            //X component of P at the last step inside the TPC [GeV/c]
    std::vector<Float_t>  EndPy_tpcAV;            //Y component of P at the last step inside the TPC [GeV/c]
    std::vector<Float_t>  EndPz_tpcAV;            //Z component of P at the last step inside the TPC [GeV/c]
    std::vector<Float_t>  pathlen_drifted;        //Particle length [cm] in drifted volume calculated trajectory-by-trajectory
    std::vector<Int_t>    inTPCDrifted;            //Was the particle in the drifted volume? (Yes=1; No=0)
    std::vector<Float_t>  StartPointx_drifted;    //'X' position of particle's first step inside the drifted volume in cm
    std::vector<Float_t>  StartPointy_drifted;    //'Y' position of particle's first step inside the drifted volume in cm
    std::vector<Float_t>  StartPointz_drifted;    //'Z' position of particle's first step inside the drifted volume in cm
    std::vector<Float_t>  StartT_drifted;        //Time associated with particle's first step inside drifted volume [ns]
    std::vector<Float_t>  StartE_drifted;        //Energy of the particle at the first step inside the drifted volume [GeV]
    std::vector<Float_t>  StartP_drifted;        //Momentum of the particle at the first step inside the drifted volume [GeV/c]
    std::vector<Float_t>  StartPx_drifted;        //X component of P at the first step inside the drifted volume [GeV/c]
    std::vector<Float_t>  StartPy_drifted;        //Y component of P at the first step inside the drifted volume [GeV/c]
    std::vector<Float_t>  StartPz_drifted;        //Z component of P at the first step inside the drifted volume [GeV/c]
    std::vector<Float_t>  EndPointx_drifted;    //'X' position of particle's last step inside the drifted volume in cm
    std::vector<Float_t>  EndPointy_drifted;    //'Y' position of particle's last step inside the drifted volume in cm
    std::vector<Float_t>  EndPointz_drifted;    //'Z' position of particle's last step inside the drifted volume in cm
    std::vector<Float_t>  EndT_drifted;            //Time associated with particle's last step inside the TPC [ns]
    std::vector<Float_t>  EndE_drifted;            //Energy of the particle at the last step inside the drifted volume [GeV]
    std::vector<Float_t>  EndP_drifted;            //Momentum of the particle at its last step inside the TPC [GeV/c]
    std::vector<Float_t>  EndPx_drifted;        //X component of P at the last step inside the TPC [GeV/c]
    std::vector<Float_t>  EndPy_drifted;        //Y component of P at the last step inside the TPC [GeV/c]
    std::vector<Float_t>  EndPz_drifted;        //Z component of P at the last step inside the TPC [GeV/c]
    std::vector<Int_t>    NumberDaughters;        //Number of daughters of the particle
    std::vector<Int_t>    TrackId;                //TrackId assigned by Geant4
    std::vector<Int_t>    Mother;                //track ID of the Mother particle
    std::vector<Int_t>    process_primary;      //process_primary =1 for primary particle; process_primary=0 for secondary, tertiary,...,particles
    std::vector<std::string> processname;       //Physics process by which the particle was created
    std::vector<Int_t>    MergedId;             //geant track segments, which belong to the same particle, get the same MergedId
    std::vector<Int_t>    origin;               //What type of MC made this particle? 0: unknown, 1: neutrino, 2: cosmic, 3: supernova, 4: singles 
    std::vector<Int_t>    MCTruthIndex;         //this geant particle comes from the neutrino interaction of the _truth variables with this index
Additional notes:
  • In the case of single particle samples (or samples generated with particle gun), as mentioned in the MC Truth section above, we don't save the primary particle information separately like we do for GENIE or CRY generators. The reason is because, for single particle samples, there is only one primary parent particle and the information about this particle is available from the first index of the GEANT particles that are saved in the analysis tree. So, for e.g., in any given event, pdg[0], Eng[0], phi[0], processname[0], StartPointx[0] etc. will all have the information corresponding to the primary single particle that is fed into GEANT4.

The TTree Viewer below shows some of the GEANT variables listed in the tree:

An example root template that shows how to access GEANT information in a ROOT macro is given below:

  //loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the total number of geant particles
        for(Int_t j=0;j<geant_list_size; j++) {
             ...
            //check how many primary muons enter the TPC
            if  (inTPCActive[j]==1 && abs(pdg[j])==13 && process_primary[j]==1)  
                    std::cout<<"\n Found a primary muon!";
            //now check how many compton electrons are present
            if  (abs(pdg[j])==11  &&  processname->at(j)=="compt")
                   std::cout<<"\n Found a compton electron!"; 
            .......       
        }//end of geant_list_size 
    }//end of ntuple entries

Auxiliary detector information

The Auxiliary detector variables are 2-index arrays with dimensions [MaxGeantParticles][kMaxAuxDets]. The first index is the number of GEANT particles and the second index is the number of auxiliary detectors which the particle traversed. In general, given how the auxiliary detectors are arranged in MicroBooNE, it is highly unusual for a particle to interact with more than 4 detectors along its path. So, the maximum size for the second index kMaxAuxDets is set to 4 in the analysis tree (this may change in the future). Also, since the auxiliary detector information is not of interest to every user, this information is only saved in the tree if the flag SaveAuxdetInfo is set to true in the fhicl file.

The list of variables saved in the analysis tree for auxiliary detectors along with a brief description of what they are is given below:

    std::vector<UShort_t> NAuxDets;            // Number of Auxiliary detectors crossed by this particle
    AuxDetMCData_t<Short_t> AuxDetID;          // Which auxiliary detector this particle went through
    AuxDetMCData_t<Float_t> entryX;            // Entry X position of particle into AuxDet in cm
    AuxDetMCData_t<Float_t> entryY;            // Entry Y position of particle into AuxDet in cm
    AuxDetMCData_t<Float_t> entryZ;            // Entry Z position of particle into AuxDet in cm
    AuxDetMCData_t<Float_t> entryT;            // Entry time of particle into AuxDet in ns (from time G4 takes control)
    AuxDetMCData_t<Float_t> exitX;             // Exit X position of particle out of AuxDet in cm
    AuxDetMCData_t<Float_t> exitY;             // Exit Y position of particle out of AuxDet in cm
    AuxDetMCData_t<Float_t> exitZ;             // Exit Z position of particle out of AuxDet in cm
    AuxDetMCData_t<Float_t> exitT;             // Exit time of particle out of AuxDet in ns (from time G4 takes control)
    AuxDetMCData_t<Float_t> exitPx;            // Exit x momentum of particle out of AuxDet
    AuxDetMCData_t<Float_t> exitPy;            // Exit y momentum of particle out of AuxDet
    AuxDetMCData_t<Float_t> exitPz;            // Exit z momentum of particle out of AuxDet
    AuxDetMCData_t<Float_t> CombinedEnergyDep; //Sum energy (GeV) of all particles with this trackID (+ID or -ID) in AuxDet

Additional information for some of the above variables:

  • AuxDetID: Each auxiliary detector is given an integer 0, 1, 2, ...
  • CombinedEnergyDep: As a particle passes through the detector, energy is lost on each simulated step it takes through. This particle has a particular trackID assigned by Geant4. If a particle initiates an electromagnetic shower, Geant4 in LArSoft applies a short-cut to save memory. New particles created in the electromagnetic shower will be assigned a new trackID equal to the initial particle's trackID x (-1), these are called "untracked particles" because they do not enter into the G4 particle list. The CombinedEnergyDep variable sums up the combined energy deposited for all particles in this AuxDet with + or - trackID of the current particle.

Example for accessing AuxDet information:

//loop over the ntuple entries
    for(Int_t i=0;i<ntuple_entries;i++) {
        . . .
        //loop over the total number of geant particles
        for(Int_t j=0;j<geant_list_size; j++) {
             ...
            //loop over the number of AuxDets this particle traversed
           for(Int_t k=0; k<NAuxDets[j]; k++) {
               float energyDep=CombinedEnergyDep[j][k];
               std::cout<<"\n This particle is responsible for "<<energyDep<<"GeV of energy in AuxDet["<<k<<"]!";
           }
            .......       
        }//end of geant_list_size 
    }//end of ntuple entries

MCEventWeight Information

If the art event contains MCEventWeights, those will be automatically saved in AnalysisTree. MCEventWeights are weights that the user may apply to the event in order to study systematics. See the MCEventWeight wiki page for more information.

    Int_t evtwgt_nfunc;                                // number of functions used
    std::vector<std::string> evtwgt_funcname;          // the name of the functions used
    std::vector<int> evtwgt_nweight;                   // number of weights for each function
    std::vector<std::vector<double>> evtwgt_weight;    // the weights (a vector for each function used)

The art event may contain several sets of weights corresponding to different configurations. For example, for the GENIE case, one may have decided to vary one or more GENIE parameters separately. Each of these configurations are called functions here. The number of functions that have been used is stored in evtwgt_nfunc, the name of those functions is stored in evtwgt_funcname. The number of weights generated for each function is stored in evtwgt_nweight and finally the actual weights are stored in vectors inside evtwgt_weight.

For example, let's assume that weights have been generated varying two GENIE parameters (called "QEMA" and "NCELaxial") separately and varying them by +-1 sigma. In this case the above variables will look like this:

evtwgt_nfunc    = 2
evtwgt_funcname = { “QEMA”,   “NCELaxial” }
evtwgt_nweight  = { 2,         2          }
evtwgt_weight   = { {w1, w2}, {w'1, w'2}  }

Skeleton code for analyzing Ntuples

ROOT provides the TTree::MakeClass method to create a framework for analyzing the ntuple data. The MakeClass() utility creates a .C and .h file with the name that you supply to the MakeClass method. The .h file declares and defines classes and Tree branches, sets branch addresses for reading all tree variables etc. which otherwise one would write from scratch if one were to create their own root macro. The MakeClass method eliminates such requirement for the users and does everything for you. The .C file is kept as short as possible and contains only the code structure to loop over the entries of the ntuple using the Loop() method. User can then modify the .C file and add their analysis code inside the tree loop and do their analysis.

The example root macro below shows how to generate .C and .h files for AnalysisTree ntuples (users are encouraged to copy the below code to their area and run them on their favorite input file):

void test()
{
//lets pick a analysis tree file from the MCC5 BNB sample located on dCache
TFile *f = new TFile("/pnfs/uboone/scratch/uboonepro/v02_05_01/mergeana/prodgenie_bnb_nu_uboone/75646_0/ana_hist.root ");
f->cd("analysistree");
TTree *v = (TTree*)f->Get("anatree");
anatree->MakeClass("AnaBNB")
}

The above code generates AnaBNB.h and AnaBNB.C files. The AnaBNB.C file will look something like this:

#define AnaBNB_cxx
#include "AnaBNB.h" 

void AnaBNB::Loop()
{
   if (fChain == 0) return;
   Long64_t nentries = fChain->GetEntriesFast();
   Long64_t nbytes = 0, nb = 0;

   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   
       nbytes += nb;
   }
}

At the beginning of the file are instructions on how to compile/run it and about reading selected branches. They are not reprinted here, but please read them from your own file. Inside the nentries loop one can then add code to analyze the ntuple data. Using the template above, lets add in some code in the AnaBNB.C file inside the jentry loop that selects and counts the number of reconstructed tracks with a length 50% or higher than the truth track. Using analysis tree variables, the code that does this is shown here:

#define AnaBNB_cxx
#include "AnaBNB.h" 

void AnaBNB::Loop()
{
   if (fChain == 0) return;
   Long64_t nentries = fChain->GetEntriesFast();
   Long64_t nbytes = 0, nb = 0;

   unsigned int count=0; 
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   
      nbytes += nb;
      //loop over the list of reconstructed tracks (lets pick @kalmanhit@ tracking algorithm for this
      for(Int_t i=0;i<ntracks_trackkalmanhit;i++){       
          //now loop over the geant particle list
          for(Int_t j=0;j<geant_list_size;j++){          
             if (trklen_trackkalmanhit[i]>0.5*pathlen[j])
                               count++;
         }
      }
   }
  std::cout<<"\n Found "<<count<<" number of reconstructed tracks with a length greater than 50% of the truth track!";
}

Now, save the file and lets compile and run it and get some results:

<uboonegpvm03.fnal.gov> root -l
root [0] .L AnaBNB.C+                           (this compiles the root macro)
Info in <TUnixSystem::ACLiC>: creating shared library /uboone/app/users/sowjanya/MyCFiles/./AnaBNB_C.so
root [1] AnaBNB  t
root [2] t.Loop()

Tutorials

See MicroBooNE-doc-3794-v1 for an AnalysisTree tutorial.