Merge branch 'main' into trig

fcl/TrkAnaExtracted.fcl

@@ -19,15 +19,14 @@ physics :
 physics.TrkAnaTrigPath : [ ]
 physics.TrkAnaEndPath : [ TrkAnaExt ]
 
-# Include more information (MC, full TrkQual and TrkPID branches)
 physics.analyzers.TrkAnaExt.candidate.options : @local::AllOpt
-physics.analyzers.TrkAnaExt.candidate.options.fillTrkQual : false
 physics.analyzers.TrkAnaExt.candidate.options.fillTrkPID : false
 physics.analyzers.TrkAnaExt.FillTriggerInfo : false
 physics.analyzers.TrkAnaExt.FitType : KinematicLine
+physics.analyzers.TrkAnaExt.candidate.branch : "trk"
 
 # for hit level diagnostics, set diagLevel to 2
-physics.analyzers.TrkAnaExt.diagLevel : 2
+physics.analyzers.TrkAnaExt.diagLevel : 1
 physics.analyzers.TrkAnaExt.FillMCInfo : true
 physics.analyzers.TrkAnaExt.FillCRVHits : true
 

inc/EventInfo.hh

@@ -8,9 +8,9 @@
 namespace mu2e
 {
   struct EventInfo {
-    int eventid = 0;
-    int runid = 0;
-    int subrunid = 0; // run/event identification
+    int event = 0;
+    int run = 0;
+    int subrun = 0; // run/event identification
     int nprotons = 0; // estimated # of protons on target for this microbunch
     float pbtime = 0.0;
     float  pbterr = 0.0; // estimated proton bunch time (and error)

inc/InfoMCStructHelper.hh

@@ -21,6 +21,10 @@
 #include "Offline/MCDataProducts/inc/KalSeedMC.hh"
 #include "BTrk/BbrGeom/HepPoint.h"
 #include "Offline/MCDataProducts/inc/PrimaryParticle.hh"
+#include "Offline/BFieldGeom/inc/BFieldManager.hh"
+#include "Offline/GeometryService/inc/GeomHandle.hh"
+#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh"
+#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh"
 
 #include <vector>
 #include <functional>

inc/SimInfo.hh

@@ -18,6 +18,7 @@ namespace mu2e
     Int_t rank = -1; // rank of this particle, by # of associated hits
     Int_t pdg = -1; // true PDG code
     Int_t proc = -1; // process code
+    Int_t stopCode = -1; // stop process code
     Int_t gen = -1; // generator code
     Float_t time = -1.0;  // Origin time of the SimParticle
     Int_t index; // index into the SimInfo vector

inc/TrkInfo.hh

@@ -89,6 +89,8 @@ namespace mu2e
     int  ndigigood =0; // number of true straw digitizations where the primary particle (not a delta-ray) crossed threshold
     int nhits =0, nactive =0; // number of hits on the track and active on the track generated by the primary particle
     int nambig =0; // number of true hits where the reconstruction assigned the correct left-right ambiguity
+    float maxr = 0.; // max radius calculated using the LoopHelix constructor
+    float rad=0,lam=0,cx=0,cy=0,phi0=0,t0=0; // truth track parameters
     void reset() { *this = TrkInfoMC(); }
   };
 }

src/InfoMCStructHelper.cc

@@ -17,6 +17,7 @@
 #include "Offline/GlobalConstantsService/inc/PhysicsParams.hh"
 #include "Offline/GeometryService/inc/GeomHandle.hh"
 #include "Offline/GeometryService/inc/DetectorSystem.hh"
+#include "Offline/TrackerGeom/inc/Tracker.hh"
 #include "art/Framework/Principal/Handle.h"
 #include "art/Framework/Principal/Event.h"
 
@@ -50,11 +51,33 @@ namespace mu2e {
   void InfoMCStructHelper::fillTrkInfoMC(const KalSeed& kseed, const KalSeedMC& kseedmc, TrkInfoMC& trkinfomc) {
     // use the primary match of the track
     // primary associated SimParticle
+    GeomHandle<DetectorSystem> det;
     if(kseedmc.simParticles().size() > 0){
       auto const& simp = kseedmc.simParticles().front();
       trkinfomc.valid = true;
       trkinfomc.nhits = simp._nhits; // number of hits from the primary particle
       trkinfomc.nactive = simp._nactive; // number of active hits from the primary particle
+
+      static GlobalConstantsHandle<ParticleDataList> pdt;
+      auto charge = pdt->particle(simp._pdg).charge();
+
+      XYZTVectorF mom = XYZTVectorF(simp._mom);
+      ROOT::Math::XYZTVector pos0(simp._pos.x(), simp._pos.y(), simp._pos.z(), simp._pos.t());
+      ROOT::Math::PxPyPzMVector mom0(mom.x(), mom.y(), mom.z(),  pdt->particle(simp._pdg).mass());
+
+      GeomHandle<BFieldManager> bfmgr;
+      mu2e::GeomHandle<mu2e::Tracker> tracker;
+      auto tracker_origin = det->toMu2e(tracker->origin());
+      ROOT::Math::XYZVector bnom(bfmgr->getBField(tracker_origin));
+
+      KinKal::LoopHelix lh(pos0, mom0, charge, bnom);
+      trkinfomc.maxr =sqrt(lh.cx()*lh.cx()+lh.cy()*lh.cy())+fabs(lh.rad());
+      trkinfomc.rad = lh.rad();
+      trkinfomc.lam = lh.lam();
+      trkinfomc.cx = lh.cx();
+      trkinfomc.cy = lh.cy();
+      trkinfomc.phi0= lh.phi0();
+      trkinfomc.t0 = lh.t0();
     }
 
     fillTrkInfoMCDigis(kseed, kseedmc, trkinfomc);
@@ -209,12 +232,13 @@ namespace mu2e {
     siminfo.valid = true;
     if(sp.genParticle().isNonnull())siminfo.gen = sp.genParticle()->generatorId().id();
     siminfo.proc = sp.creationCode();
+    siminfo.stopCode = sp.stoppingCode();
     siminfo.pdg = sp.pdgId();
     siminfo.time = sp.startGlobalTime();
     siminfo.mom = XYZVectorF(sp.startMomentum());
     siminfo.pos = XYZVectorF(det->toDetector(sp.startPosition()));
-    siminfo.endmom = XYZVectorF(sp.endMomentum());
     siminfo.endpos = XYZVectorF(det->toDetector(sp.endPosition()));
+    siminfo.endmom = XYZVectorF(sp.endMomentum());
   }
 
   void InfoMCStructHelper::fillVDInfo(const KalSeed& kseed, const KalSeedMC& kseedmc, std::vector<MCStepInfo>& vdinfos) {

src/TrkAnaTreeMaker_module.cc

@@ -465,8 +465,8 @@ namespace mu2e {
             auto& mcssi = mcssis.at(ixtra);
             auto const& tag = _extraMCStepTags[ixtra];
             auto inst = tag.instance();
-            std::string  mcsiname = branch +"mcsic_" + inst;
-            std::string  mcssiname = branch + "mcssi_" + inst;
+            std::string  mcsiname = branch +"mcsic_" + inst + ".";
+            std::string  mcssiname = branch + "mcssi_" + inst + ".";
             _trkana->Branch(mcsiname.c_str(),mcsic,_buffsize,_splitlevel);
             _trkana->Branch(mcssiname.c_str(),mcssi,_buffsize,_splitlevel);
           }
@@ -701,9 +701,9 @@ namespace mu2e {
 
   void TrkAnaTreeMaker::fillEventInfo( const art::Event& event) {
     // fill basic event information
-    _einfo.eventid = event.event();
-    _einfo.runid = event.run();
-    _einfo.subrunid = event.subRun();
+    _einfo.event = event.event();
+    _einfo.run = event.run();
+    _einfo.subrun = event.subRun();
     // currently no reco nproton estimate TODO
 
     auto PBThandle = event.getValidHandle<mu2e::ProtonBunchTime>(_PBTTag);

tutorial/pages/genealogy.md

@@ -47,7 +47,8 @@ c1->SetGridy(true);
 Now we can draw the X-Z starting position of the conversion electrons like this:
 
 ```
-  trkana->Draw("demmcsim.pos.x():demmcsim.pos.z()>>hist(1000,-20000,1000, 1000,-700,9000)", "demmcsim.prirel._rem==1 && demmcsim.prirel._rel==2", "COLZ");
+  TH2D* hist = new TH2D("hist", "", 1000,-20000,1000, 1000,-700,9000);
+  trkana->Draw("demmcsim.pos.x():demmcsim.pos.z()>>hist", "demmcsim.prirel._rem==1 && demmcsim.prirel._rel==2", "COLZ");
 ```
 where ```prirel._rel==2``` is for parent. You should see the S-bend of the transport solenoid and also the production solenoid.
 

tutorial/pages/mom-res.md

@@ -47,7 +47,8 @@ c1->SetGridy(true);
 We can first plot the reconstructed momentum and the MC truth information separately:
 
 ```
-trkana->Draw("demfit.mom.R()>>hist(100,100,110)", "demfit.sid==0", "HIST");
+TH1D* hist = new TH1D("hist", "", 100,100,110));
+trkana->Draw("demfit.mom.R()>>hist", "demfit.sid==0", "HIST");
 trkana->Draw("demmcvd.mom.R()>>hist2", "demmcvd.sid==0", "HIST SAMES");
 ```
 

tutorial/pages/n-hits.md

@@ -58,7 +58,8 @@ c1->SetGridy(true);
 We can then plot a histogram of the number of hits on each track like so:
 
 ```
-trkana->Draw("dem.nhits>>hist(100,0,100)", "", "goff");
+TH1D* hist = new TH1D("hist", "", 100,0,100)
+trkana->Draw("dem.nhits>>hist", "", "goff");
 ```
 
 where ```goff``` is the drawing option to now draw since we are still missing some important information: axis labels

tutorial/pages/reco-mom.md

@@ -89,15 +89,16 @@ We want to plot the momentum of the tracks. If you ```Print()``` the ```demfit``
 We want to plot the magnitude of the momentum vector. We could calculate it ourselves from the components but in ROOT, we can use member functions of [XYZVectorF](https://root.cern.ch/doc/v628/namespaceROOT_1_1Math.html#a767e8c52a85dc9538fe00603961eab98). We will use ```R()``` function:
 
 ```
-trkana->Draw("demfit.mom.R()>>hist(100,100,110)", "", "HIST");
+TH1D* hist = new TH1D("hist", "", 100,100,110")
+trkana->Draw("demfit.mom.R()>>hist, "", "HIST");
 ```
 
 You may notice that the peak is rather broad... That's because we are plotting the reconstructed track momentum at all the intersections - the entrance, middle, and exit of the tracker. 
 
 To see just the track momentum at one of the intersections, we need to apply a cut using the second argument to the ```Draw()``` function:
 
 ```
-trkana->Draw("demfit.mom.R()>>hist(100,100,110)", "demfit.sid==0", "HIST")
+trkana->Draw("demfit.mom.R()>>hist", "demfit.sid==0", "HIST")
 ```
 
 where ```sid``` is the surface id of the intersection we want to plot (```sid=0``` is the entrance to the tracker).

tutorial/pages/start-pos.md

@@ -54,7 +54,8 @@ c1->SetGridy(true);
 Now we can draw the X-Y starting position of the conversion electrons like this:
 
 ```
-  trkana->Draw("demmcsim.pos.y():demmcsim.pos.x()>>hist(200,-100,100, 200,-100,100)", "demmcsim.prirel._rem==0", "COLZ");
+  TH2D* hist = new TH2D("hist", "", 100-100,100, 200,-100,100)
+  trkana->Draw("demmcsim.pos.y():demmcsim.pos.x()>>hist", "demmcsim.prirel._rem==0", "COLZ");
 ```
 
 where ```COLZ``` tells ROOT to draw with a color map. Note that the coordinate system used by TrkAna is the detector coordinate system (origin = center of tracker) and not the global Mu2e coordinate system (origin = center of TS3).

validation/create_val_file.C

@@ -0,0 +1,74 @@
+void create_val_file() {
+
+  TFile* trkana_file = new TFile("/pnfs/mu2e/tape/phy-nts/nts/mu2e/CeEndpointMix1BBSignal/MDC2020z1_best_v1_1_std_v04_01_00/tka/85/7b/nts.mu2e.CeEndpointMix1BBSignal.MDC2020z1_best_v1_1_std_v04_01_00.001210_00000289.tka", "READ");
+  TTree* trkana = (TTree*) trkana_file->Get("TrkAnaNeg/trkana");
+
+  TFile* file = new TFile("val-trkana-v4.root", "RECREATE");
+
+  // evtinfo histograms
+  trkana->Draw("evtinfo.eventid>>h_evtinfo_eventid", "", "");
+  trkana->Draw("evtinfo.subrunid>>h_evtinfo_subrunid", "", "");
+  trkana->Draw("evtinfo.runid>>h_evtinfo_runid", "", "");
+  trkana->Draw("evtinfo.nprotons>>h_evtinfo_nprotons", "", "");
+  trkana->Draw("evtinfo.pbtime>>h_evtinfo_pbtime", "", "");
+  trkana->Draw("evtinfo.pbterr>>h_evtinfo_pbterr", "", "");
+
+  // evtinfomc histograms
+  trkana->Draw("evtinfomc.nprotons>>h_evtinfomc_nprotons", "", "");
+  trkana->Draw("evtinfomc.pbtime>>h_evtinfomc_pbtime", "", "");
+
+  // demfit histograms
+  trkana->Draw("demfit.mom.R()>>h_demfit_mom_all", "", "goff");
+  trkana->Draw("demfit.mom.R()>>h_demfit_mom_ent", "demfit.sid==0", "goff");
+  trkana->Draw("demfit.mom.R()>>h_demfit_mom_mid", "demfit.sid==1", "goff");
+  trkana->Draw("demfit.mom.R()>>h_demfit_mom_xit", "demfit.sid==2", "goff");
+  trkana->Draw("demfit.mom.R()>>h_demfit_mom_ent_t0cut", "demfit.sid==0 && demlh.t0>=700", "goff");
+  trkana->Draw("demfit.mom.R()>>h_demfit_mom_mid_t0cut", "demfit.sid==1 && demlh.t0>=700", "goff");
+  trkana->Draw("demfit.mom.R()>>h_demfit_mom_xit_t0cut", "demfit.sid==2 && demlh.t0>=700", "goff");
+
+  // resolution histograms
+  trkana->Draw("(demfit[demmcvd.iinter].mom.R() - demmcvd.mom.R())>>h_demfit_momres_all", "", "goff");
+  trkana->Draw("(demfit[demmcvd.iinter].mom.R() - demmcvd.mom.R())>>h_demfit_momres_ent", "demmcvd.sid==0", "goff");
+  trkana->Draw("(demfit[demmcvd.iinter].mom.R() - demmcvd.mom.R())>>h_demfit_momres_mid", "demmcvd.sid==1", "goff");
+  trkana->Draw("(demfit[demmcvd.iinter].mom.R() - demmcvd.mom.R())>>h_demfit_momres_xit", "demmcvd.sid==2", "goff");
+
+  // trkcalohit histograms
+  trkana->Draw("demtch.ctime>>h_demtch_ctime_all", "", "goff");
+  trkana->Draw("demtch.ctime>>h_demtch_ctime_active", "demtch.active==1", "goff");
+  trkana->Draw("demtch.ctime>>h_demtch_ctime_inactive", "demtch.active==0", "goff");
+  trkana->Draw("demtch.ctime-demlh.t0>>h_demtch_demfit_dt_ent", "demfit.sid==0 && demtch.active==1", "goff");
+
+  // crv histograms
+  trkana->Draw("crvsummary.totalPEs>>h_crvsummary_totalPEs_all", "", "goff");
+  trkana->Draw("crvhit.pos.fCoordinates.fX>>h_crvhit_pos_x", "", "goff");
+  trkana->Draw("crvhit.pos.fCoordinates.fY>>h_crvhit_pos_y", "", "goff");
+  trkana->Draw("crvhit.pos.fCoordinates.fZ>>h_crvhit_pos_z", "", "goff");
+  trkana->Draw("crvhitmc.primary.fCoordinates.fX>>h_crvhitmc_primary_x", "", "goff");
+  trkana->Draw("crvhitmc.primary.fCoordinates.fY>>h_crvhitmc_primary_y", "", "goff");
+  trkana->Draw("crvhitmc.primary.fCoordinates.fZ>>h_crvhitmc_primary_z", "", "goff");
+  trkana->Draw("crvhitmc.depositedEnergy>>h_crvhitmc_depostedEnergy", "", "goff");
+
+  // demmcsim histograms
+  trkana->Draw("demmcsim.pos.x()>>h_demmcsim_pos_x_all", "", "goff");
+  trkana->Draw("demmcsim.pos.y()>>h_demmcsim_pos_y_all", "", "goff");
+  trkana->Draw("demmcsim.pos.z()>>h_demmcsim_pos_z_all", "", "goff");
+  trkana->Draw("demmcsim.pos.x()>>h_demmcsim_pos_x_evtprimary", "demmcsim.prirel._rel==0 && demmcsim.prirel._rem==0", "goff");
+  trkana->Draw("demmcsim.pos.y()>>h_demmcsim_pos_y_evtprimary", "demmcsim.prirel._rel==0 && demmcsim.prirel._rem==0", "goff");
+  trkana->Draw("demmcsim.pos.z()>>h_demmcsim_pos_z_evtprimary", "demmcsim.prirel._rel==0 && demmcsim.prirel._rem==0", "goff");
+  trkana->Draw("demmcsim.pos.x()>>h_demmcsim_pos_x_trkprimary", "demmcsim.trkrel._rel==0 && demmcsim.trkrel._rem==0", "goff");
+  trkana->Draw("demmcsim.pos.y()>>h_demmcsim_pos_y_trkprimary", "demmcsim.trkrel._rel==0 && demmcsim.trkrel._rem==0", "goff");
+  trkana->Draw("demmcsim.pos.z()>>h_demmcsim_pos_z_trkprimary", "demmcsim.trkrel._rel==0 && demmcsim.trkrel._rem==0", "goff");
+  trkana->Draw("demmcsim.pos.x()>>h_demmcsim_pos_x_trkparent", "demmcsim.trkrel._rel==2 && demmcsim.trkrel._rem==1", "goff");
+  trkana->Draw("demmcsim.pos.y()>>h_demmcsim_pos_y_trkparent", "demmcsim.trkrel._rel==2 && demmcsim.trkrel._rem==1", "goff");
+  trkana->Draw("demmcsim.pos.z()>>h_demmcsim_pos_z_trkparent", "demmcsim.trkrel._rel==2 && demmcsim.trkrel._rem==1", "goff");
+  trkana->Draw("demmcsim.pos.x()>>h_demmcsim_pos_x_trkgparent", "demmcsim.trkrel._rel==5 && demmcsim.trkrel._rem==2", "goff");
+  trkana->Draw("demmcsim.pos.y()>>h_demmcsim_pos_y_trkgparent", "demmcsim.trkrel._rel==5 && demmcsim.trkrel._rem==2", "goff");
+  trkana->Draw("demmcsim.pos.z()>>h_demmcsim_pos_z_trkgparent", "demmcsim.trkrel._rel==5 && demmcsim.trkrel._rem==2", "goff");
+  trkana->Draw("demmcsim.pos.x()>>h_demmcsim_pos_x_trkggparent", "demmcsim.trkrel._rel==5 && demmcsim.trkrel._rem==3", "goff");
+  trkana->Draw("demmcsim.pos.y()>>h_demmcsim_pos_y_trkggparent", "demmcsim.trkrel._rel==5 && demmcsim.trkrel._rem==3", "goff");
+  trkana->Draw("demmcsim.pos.z()>>h_demmcsim_pos_z_trkggparent", "demmcsim.trkrel._rel==5 && demmcsim.trkrel._rem==3", "goff");
+
+
+  file->Write();
+  file->Close();
+}