Initial commit

README.md

+Detector Benchmarks
+===================
+
+First Step
+----------
+
+Create symbolic link to detectors and data directories
+
+```bash
+# Starting from the directory of this readme
+cd .. &&  git clone git@gitlab.com:Argonne_EIC/detectors.git
+cd TOF_analysis
+ln -s ../detectors
+ln -s ../data
+```
+Here we have assumed local data is organized in the directory `../data`.
+
+

bin/run_detector_bechmarks

+#!/bin/bash
+
+DIR="$(cd "$(dirname "$0")" && pwd)"
+cd $DIR/..
+pwd
+
+mkdir -p results
+
+root_script="TOF_analysis"
+echo  "Running ${root_script} ..."
+mkdir -p results/${root_script}
+root -b -q scripts/${root_script}.cxx 1> results/${root_script}/${root_script}.out 2> results/${root_script}/${root_script}.err
+
+root_script="compact_geo_test"
+echo  "Running ${root_script} ..."
+mkdir -p results/${root_script}
+root -b -q scripts/${root_script}.cxx 1> results/${root_script}/${root_script}.out 2> results/${root_script}/${root_script}.err
+

scripts/TOF_analysis.cxx

+#include "UTIL/LCTOOLS.h"
+#include "UTIL/Operators.h"
+#include <UTIL/CellIDDecoder.h>
+//#include "DD4hep/LCDD.h"
+//#include "DD4hep/DD4hepUnits.h"
+// check directory
+#include <sys/stat.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include "TMath.h"
+#include"time.h"
+#include <dirent.h>
+#include <string>
+#include <vector>
+#include <map>
+struct stat sb;
+
+std::vector<std::string> get_files_in_directory(std::string path = "data/.") {
+  std::vector<std::string> files;
+  if (stat(path.c_str(), &sb) == 0 && S_ISDIR(sb.st_mode)) { 
+    DIR*    dir;
+    dirent* pdir;
+    dir = opendir(path.c_str());
+    while( pdir = readdir(dir) ){
+        string rfile=path + "/" + pdir->d_name;
+        if (rfile.find(".slcio") != std::string::npos) 
+             files.push_back(rfile); 
+    }
+    } else {
+        cout << "Directory " << path << " does not exist!!" << endl;
+        std::exit(0); 
+        }
+    return files;
+}
+//______________________________________________________________________________
+
+
+void TOF_analysis(
+    const char* DIRNAME            = "data/electrons",
+    const char* result_file_prefix = "SID_1T",
+    const char* compact_file       = "detectors/SiD/compact/sid_working/sidloi3_v00.xml")
+{
+
+  DD4hep::Geometry::LCDD& lcdd = DD4hep::Geometry::LCDD::getInstance();
+  lcdd.fromCompact(compact_file);
+  const double position[3]={0,0,0}; // position to calculate magnetic field at (the origin in this case)
+  double bField[3]={0,0,0}; 
+  lcdd.field().magneticField(position,bField); 
+  _Bz = bField[2]/dd4hep::tesla; 
+
+  std::cout << " Magnetic Field Bz = " << _Bz << std::endl;
+  return;
+
+  //lcdd.field();
+  //CellIDDecoder<SimTrackerHit> cellid_decoder( hitCol );
+  //for(int i=0; i<nhits; i++){ 
+  //  SimTrackerHit *hit =dynamic_cast <SimTrackerHit*>( hitCol->getElementAt(i) );
+  //  int subdetectorID = cellid_decoder( hit )["subdet"];
+  //  int layerID = cellid_decoder( hit )["layer"];
+  //  int moduleID = cellid_decoder( hit )["module"];
+  //}
+
+  // Put the equivalent lines in your .rootlogon.C
+  //gSystem->AddIncludePath(" -I$EIC_SOFTWARE/local/include");
+  //gSystem->Load("liblcio.so");
+  //gSystem->Load("liblcioDict.so");
+  
+  using namespace EVENT ;
+  using namespace UTIL;
+  using namespace IMPL ;
+
+  TH1F * hTime= new TH1F("hTime","hit time ;time [ns];",100,0,100);
+  TH1F * hEDep= new TH1F("hEDep","EDep ;E [GeV];",500,0,0.001);
+  TH1F * hdEdx= new TH1F("hdEdx","dEdx ;E [GeV];",500,0,0.01);
+  TH2F * hEDepVSdEdx= new TH2F("hEDepVSdEdx","EDep Vs dEdx ;dEdx [GeV];EDep [GeV]",500,0,0.01,500,0,0.001);
+  TH2F * hEDepVSnhit= new TH2F("hEDepVSnhit","dEdx Vs nhit ;nhit ;dEdx [GeV]",20,0,20,500,0,0.01);
+  TProfile * hdEdx_vs_p= new TProfile("hdEdx_vs_p","dEdx Vs P ;P [GeV/c] ;dEdx [GeV]",20,0,20,0,0.005);
+  TH1F * h1   = new TH1F("hh1","P_{T};P_{T} [GeV];",80,0,40);
+  TH1F * hP   = new TH1F("hP"  , "P;P [GeV/c];"        , 80 , 0   , 40);
+  TH1F * hPT  = new TH1F("hPT" , "P_{T};P_{T} [GeV/c];", 80 , 0   , 40);
+  TH1F * hPDG = new TH1F("hPDG", "PDG Code; PDG code;" , 40 , -20 , 20);
+
+  auto eff_vs_Pt2 = new TEfficiency("eff_vs_Pt2","Track timing2;P_{T} [GeV];#epsilon",50,0,10);
+  auto eff_vs_Pt  = new TEfficiency("eff_vs_Pt","Track timing;P_{T} [GeV];#epsilon",50,0,10);
+  auto eff_vs_P   = new TEfficiency("eff_vs_P","Track timing;P [GeV];#epsilon",50,0,10);
+
+  auto eff_vs_Pt3 = new TEfficiency("eff_vs_Pt2","Track timing2;P_{T} [GeV];#epsilon",50,0,10);
+  auto eff_vs_Pt4 = new TEfficiency("eff_vs_Pt2","Track timing2;P_{T} [GeV];#epsilon",50,0,10);
+
+  auto chi2_vs_Pt = new TProfile("chi2_vs_Pt","#chi^{2} vs P_{T};P_{T} [GeV];#chi^{2}",20,0,10,0,100);
+  auto chi2_vs_P  = new TProfile("chi2_vs_P","#chi^{2} vs P;P [GeV];#chi^{2}",20,0,10,0,100);
+
+  int  max_files    = 10;
+  int  nEvents      = 0  ;
+  int  n_files_read = 0;
+  auto files        = get_files_in_directory(DIRNAME);
+
+  IO::LCReader* lcReader = IOIMPL::LCFactory::getInstance()->createLCReader() ;
+  EVENT::LCEvent* evt = 0 ;
+  std::string mcpType("std::vector<IMPL::MCParticleImpl*>") ;
+  std::string mcpName("MCParticle") ;
+
+  std::sort(files.begin(),files.end());
+
+  for(auto aFile : files) {
+
+    std::cout << "Opening " << aFile << '\n';
+    lcReader->open( aFile.c_str() ) ;
+    n_files_read++;
+
+    //----------- the event loop -----------
+    while( (evt = lcReader->readNextEvent()) != 0 ) {
+
+      if(nEvents%1000==0) {
+        if(nEvents==0)UTIL::LCTOOLS::dumpEvent( evt ) ;
+        std::cout << "Event " << nEvents << '\n';
+      }
+
+      IMPL::LCCollectionVec* thrown_col = (IMPL::LCCollectionVec*) evt->getCollection( "MCParticle"  ) ;
+      IMPL::LCCollectionVec* tracks_col = (IMPL::LCCollectionVec*) evt->getCollection( "Tracks"  ) ;
+      IMPL::LCCollectionVec* simtrackhits_col = (IMPL::LCCollectionVec*) evt->getCollection("SiTrackerEndcapHits") ;
+
+      //UTIL::LCTOOLS::printMCParticles( col ) ;
+
+      int nMCP    = thrown_col->getNumberOfElements() ;
+      int nTracks = tracks_col->getNumberOfElements() ;
+      bool hasOneTrack  = (nTracks==1);
+      bool hasTwoTracks = (nTracks==2);
+      int nSimTrackHits = simtrackhits_col->getNumberOfElements() ;
+
+      for(int i_thrown=0 ; i_thrown<nMCP ; ++i_thrown) {
+
+        auto  part = (EVENT::MCParticle*) thrown_col->getElementAt(i_thrown);
+
+        if( (part->getGeneratorStatus() == 1) ) {
+          // (part->getPDG() == 11) 
+          //std::cout << int(part->getGeneratorStatus()) << std::endl;;
+          //std::cout << *part << std::endl;;
+          //std::cout << (UTIL::lcio_short<EVENT::MCParticle>(*part));
+          //<< std::endl;;
+          TVector3 p(part->getMomentum());
+          double eta = p.PseudoRapidity();
+
+          hPT->Fill(p.Pt());
+          hP->Fill(p.Mag());
+          hPDG->Fill(part->getPDG());
+
+          eff_vs_P  ->Fill(hasOneTrack, p.Mag());
+          eff_vs_Pt ->Fill(hasOneTrack, p.Pt());
+          eff_vs_Pt2->Fill(hasOneTrack||hasTwoTracks, p.Pt());
+          if( TMath::Abs(eta) < 1 ) {
+            eff_vs_Pt3 ->Fill(hasOneTrack, p.Pt());
+          }
+          if( TMath::Abs(TMath::Abs(eta)-1.5) < 0.5) {
+            eff_vs_Pt4 ->Fill(hasOneTrack, p.Pt());
+          }
+
+
+          //std::cout << "Tracks: " << nTracks << "\n";
+          for(int i_track=0 ; i_track<nTracks ; ++i_track) {
+            auto  track = (EVENT::Track*) tracks_col->getElementAt(i_track);
+            chi2_vs_Pt->Fill(p.Pt(), track->getChi2(), 1);
+            chi2_vs_P->Fill( p.Mag(),track->getChi2(), 1);
+
+            double dEdx = track->getdEdx();
+            hdEdx->Fill(dEdx);
+            auto hits = track->getTrackerHits();
+            hEDepVSnhit->Fill( hits.size(),dEdx);
+            hdEdx_vs_p->Fill(p.Mag(), dEdx);
+
+            for(auto ahit : hits) {
+              double hit_time  = ahit->getTime();
+              double EDep  = ahit->getEDep();
+              if(hit_time != 0.0) {
+                std::cout << " time " << hit_time << '\n';
+              }
+              //if(EDep != 0.0) {
+              //  std::cout << " EDep " << EDep << '\n';
+              //}
+              hEDep->Fill(EDep); 
+              hEDepVSdEdx->Fill(dEdx,EDep);
+            }
+
+            //if(nTracks>1) {
+            //  std::cout <<  (*track) << std::endl;
+            //}
+          }
+
+          //std::cout << "Tracks: " << nTracks << "\n";
+          for(int i_simtrack=0 ; i_simtrack<nSimTrackHits ; ++i_simtrack) {
+            auto  sim_track_hit = (EVENT::SimTrackerHit*) simtrackhits_col->getElementAt(i_simtrack);
+
+            double ht = sim_track_hit->getTime();
+            //if(ht != 0.0) {
+            //  std::cout << " time " << ht << '\n';
+            //}
+            hTime->Fill(ht); 
+          }
+
+        }
+      }
+
+      nEvents++;
+    }
+    // -------- end of event loop -----------
+
+
+
+    lcReader->close() ;
+
+    if(n_files_read >= max_files) break;
+  }
+
+  std::cout << std::endl 
+    <<  "  " <<  nEvents 
+    << " events read from file: " 
+    << DIRNAME << std::endl  ;
+
+  delete lcReader ;
+
+  gSystem->mkdir("results/timing",true);
+
+  TCanvas * c = new TCanvas();
+  eff_vs_Pt ->SetLineWidth(2);
+  eff_vs_Pt2->SetLineWidth(2);
+  eff_vs_Pt3->SetLineWidth(2);
+  eff_vs_Pt4->SetLineWidth(2);
+  eff_vs_Pt ->SetLineColor(1);
+  eff_vs_Pt2->SetLineColor(2);
+  eff_vs_Pt3->SetLineColor(4);
+  eff_vs_Pt4->SetLineColor(8);
+  eff_vs_Pt->Draw();
+  eff_vs_Pt2->Draw("same");
+  eff_vs_Pt3->Draw("same");
+  eff_vs_Pt4->Draw("same");
+  c->SaveAs(Form("results/timing/%s_effVsPT.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  eff_vs_P->Draw();
+  c->SaveAs(Form("results/timing/%s_effVsP.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  chi2_vs_Pt->Draw();
+  c->SaveAs(Form("results/timing/%s_Chi2VsPT.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  chi2_vs_P->Draw();
+  c->SaveAs(Form("results/timing/%s_Chi2VsP.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  hTime->Draw();
+  c->SaveAs(Form("results/timing/%s_Time.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  hEDep->Draw();
+  c->SaveAs(Form("results/timing/%s_EDep.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  hdEdx->Draw();
+  c->SaveAs(Form("results/timing/%s_dEdx.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  hEDepVSdEdx->Draw("colz");
+  c->SaveAs(Form("results/timing/%s_EDepVSdEdx.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  hEDepVSnhit->Draw("colz");
+  c->SaveAs(Form("results/timing/%s_EDepVSnhit.pdf",result_file_prefix));
+
+  c = new TCanvas();
+  hdEdx_vs_p->Draw("colz");
+  c->SaveAs(Form("results/timing/%s_dEdx_vs_p.pdf",result_file_prefix));
+
+  //c = new TCanvas();
+  //hPDG->Draw();
+  //c->SaveAs(Form("results/resolution/%s_dPTOverPT.pdf",result_file_prefix));
+
+}

scripts/compact_geo_test.cxx

+#include "UTIL/LCTOOLS.h"
+#include "UTIL/Operators.h"
+#include <UTIL/CellIDDecoder.h>
+//#include "DD4hep/LCDD.h"
+//#include "DD4hep/DD4hepUnits.h"
+// check directory
+#include <sys/stat.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include "TMath.h"
+#include"time.h"
+#include <dirent.h>
+#include <string>
+#include <vector>
+#include <map>
+struct stat sb;
+
+std::vector<std::string> get_files_in_directory(std::string path = "data/.") {
+  std::vector<std::string> files;
+  if (stat(path.c_str(), &sb) == 0 && S_ISDIR(sb.st_mode)) { 
+    DIR*    dir;
+    dirent* pdir;
+    dir = opendir(path.c_str());
+    while( pdir = readdir(dir) ){
+        string rfile=path + "/" + pdir->d_name;
+        if (rfile.find(".slcio") != std::string::npos) 
+             files.push_back(rfile); 
+    }
+    } else {
+        cout << "Directory " << path << " does not exist!!" << endl;
+        std::exit(0); 
+        }
+    return files;
+}
+//______________________________________________________________________________
+
+
+void compact_geo_test(
+    const char* DIRNAME            = "data/electrons",
+    const char* result_file_prefix = "SID_1T",
+    const char* compact_file       = "detectors/SiD/compact/sid_working/sidloi3_v00.xml")
+{
+
+  DD4hep::Geometry::LCDD& lcdd = DD4hep::Geometry::LCDD::getInstance();
+  lcdd.fromCompact(compact_file);
+  const double position[3]={0,0,0}; // position to calculate magnetic field at (the origin in this case)
+  double bField[3]={0,0,0}; 
+  lcdd.field().magneticField(position,bField); 
+  _Bz = bField[2]/dd4hep::tesla; 
+
+  std::cout << " Magnetic Field Bz = " << _Bz << std::endl;
+  return;
+
+}