Skip to content
Snippets Groups Projects
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));

}