-
Whitney Armstrong authoredWhitney Armstrong authored
TOF_analysis.cxx 9.36 KiB
#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));
}