sPHENIX tracking discussion

[Hugo Pereira Da Costa <hugo.pereira-da-costa AT cea.fr>]

Hello all, especially Christof   

Following yesterday's discussion, please find attached

- a trimmed down Fun4All_sPHENIX steering macro to generate flat pions (2 per events) and the TrackEvaluation module, suited for studying pulls in our Tracking detectors

- a post-processing macro, for inspiration that generate the rphi pulls distributions in one of the TPC layers (there are several sets of possible variables to draw in there that you can comment/uncomment at will, depending on the pulls you are interested in)

You would need to update your local code, if any, to make sure it includes https://github.com/sPHENIX-Collaboration/coresoftware/pull/1383

In principle you could then run the two macros in a row and get the pulls in layer 23 of the TPC

Hope this helps,

Hugo

On 1/9/22 19:24, Anthony Frawley via sPHENIX-tracking-l wrote:

       Hi All,
       
      
       We will have our usual
        tracking software meeting at 10 am ET tomorrow, January 10. The
        Zoom link is:
       https://fsu.zoom.us/j/98126117840
      
       
      
       There is an Indico page
        at:

https://indico.bnl.gov/event/14115/

Cheers

Tony

_______________________________________________
sPHENIX-tracking-l mailing list
sPHENIX-tracking-l AT lists.bnl.gov
https://lists.bnl.gov/mailman/listinfo/sphenix-tracking-l
    
  
#include <fun4all/SubsysReco.h>
#include <fun4all/Fun4AllServer.h>
#include <fun4all/Fun4AllDummyInputManager.h>
#include <fun4all/Fun4AllDstOutputManager.h>
#include <g4main/PHG4SimpleEventGenerator.h>
#include <phool/PHRandomSeed.h>
#include <phool/recoConsts.h>

#include <g4eval/TrackEvaluation.h>

// local macros
#include "G4Setup_sPHENIX.C"
#include "G4_Bbc.C"
#include "G4_Global.C"
#include "G4_Tracking.C"

R__LOAD_LIBRARY(libfun4all.so)
R__LOAD_LIBRARY(libqa_modules.so)

//____________________________________________________________________
int Fun4All_G4_sPHENIX_cr(
  const int nEvents = 2000,
  const char *outputFile = "dst_eval_cr.root"
  )
{

  // options
  Enable::PIPE = true;
  Enable::BBC = true;
  Enable::MAGNET = true;
  Enable::PLUGDOOR = false;

  // enable all absorbers
  // this is equivalent to the old "absorberactive" flag
  Enable::ABSORBER = true;

  // central tracking
  Enable::MVTX = true;
  Enable::INTT = true;
  Enable::TPC = true;
  Enable::MICROMEGAS = true;
  Enable::BLACKHOLE = true;

  // TPC
  // space charge distortions
  G4TPC::ENABLE_STATIC_DISTORTIONS = false;
  G4TPC::static_distortion_filename = "/star/u/rcorliss/sphenix/trackingStudySampleNov2021/static_only.distortion_map.hist.root";
    
  // space charge corrections
  G4TPC::ENABLE_CORRECTIONS = false;
  G4TPC::correction_filename = "distortion_maps-new/static_only_inverted_10-new.root";
  
  // micromegas configuration
  G4MICROMEGAS::CONFIG = G4MICROMEGAS::CONFIG_BASELINE;

  // tracking configuration
  G4TRACKING::use_full_truth_track_seeding = false;
  G4TRACKING::use_truth_tpc_seeding = false;
  G4TRACKING::SC_CALIBMODE = false;

  // server
  auto se = Fun4AllServer::instance();
  se->Verbosity(1);

  // make sure to printout random seeds for reproducibility
  PHRandomSeed::Verbosity(1);

  // reco const
  auto rc = recoConsts::instance();

  {
    // event generation
    auto gen = new PHG4SimpleEventGenerator;
    gen->add_particles("pi+",1);
    gen->add_particles("pi-",1);

    gen->set_eta_range(-1.0, 1.0);
    gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());

    // flat pt distribution
    gen->set_pt_range(0.2, 20.0);
    
    // vertex
    gen->set_vertex_distribution_function(
      PHG4SimpleEventGenerator::Uniform,
      PHG4SimpleEventGenerator::Uniform,
      PHG4SimpleEventGenerator::Uniform);
    gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
    gen->set_vertex_distribution_width(0.0, 0.0, 5.0);
    gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
    gen->set_vertex_size_parameters(0.0, 0.0);

    gen->Embed(2);
    se->registerSubsystem(gen);
  }

  // Geant4 initialization
  G4Init();
  G4Setup();

  // BBC
  BbcInit();
  Bbc_Reco();

  // cells
  Mvtx_Cells();
  Intt_Cells();
  TPC_Cells();
  Micromegas_Cells();

  // digitizer and clustering
  Mvtx_Clustering();
  Intt_Clustering();
  TPC_Clustering();
  Micromegas_Clustering();

  TrackingInit();
  Tracking_Reco();

  // track evaluation
  se->registerSubsystem( new TrackEvaluation );
  
  // for single particle generators we just need something which drives
  // the event loop, the Dummy Input Mgr does just that
  se->registerInputManager( new Fun4AllDummyInputManager("JADE"));

  // output manager
  se->registerOutputManager( new Fun4AllDstOutputManager("DSTOUT", outputFile) );

  // process events
  se->run(nEvents);
  // terminate
  se->End();
  se->PrintTimer();

  std::cout << "All done" << std::endl;
  delete se;
  gSystem->Exit(0);
  return 0;
}

// This function is only used to test if we can load this as root6 macro
// without running into unresolved libraries and include files
void RunFFALoadTest() {}

#include <TCanvas.h>
#include <TChain.h>
#include <TCut.h>
#include <TGraphErrors.h>
#include <TH1.h>
#include <TH2.h>
#include <TStyle.h>

//____________________________________________________________________________
Float_t delta_phi( Float_t phi )
{
  if( phi >= M_PI ) return phi - 2*M_PI;
  else if( phi < -M_PI ) return phi + 2*M_PI; 
  else return phi;
}

//____________________________________________________________________________
void Pulls_cr()
{
  
  const TString inputFile = "dst_eval_cr.root";
  const TString pdfFile = "Pulls_cr.pdf";
  std::cout << "Pulls_cr - inputFile: " << inputFile << std::endl;
  std::cout << "Pulls_cr - pdfFile: " << pdfFile << std::endl;
  
  auto tfile = TFile::Open( inputFile );
  if( !tfile ) {
    std::cout << "Pulls_cr - invalid input file: " << inputFile << std::endl;
    return;
  }
  
  auto tree = dynamic_cast<TTree*>( tfile->Get("T") );
  if( !tree )
  {
    std::cout << "Pulls_cr - invalid tre \"T\"" << std::endl;
    return;
  }
  
  // define variables and cuts
  /* this is cluster - truth pulls, for clusters on tracks */
  const TString var( "delta_phi(m_tracks.clusters.phi - m_tracks.clusters.truth_phi)/m_tracks.clusters.phi_error" );

  /* same thing with z */
  // const TString var( "(m_tracks.clusters.z - clusters.truth_z)/m_tracks.clusters.z_error" );
  
  /* 
   * this is cluster - track pulls, for clusters on tracks  
   * in principle one should also add the track error quadratically to the denominator.
   * It is expected to be small and its validity should be checked anyway
   */
  // const TString var( "delta_phi(m_tracks.clusters.phi - m_tracks.clusters.trk_phi)/m_tracks.clusters.phi_error" );

  /* same thing with z */  
  // const TString var( "m_tracks.clusters.z - m_tracks.clusters.trk_z)/m_tracks.clusters.z_error" );

  // layer selection
  /* 23 is the first layer of the R2 GEMS in TPC */
  const int ilayer = 23;
  const TCut layer_cut = Form( "m_tracks.clusters.layer == %i", ilayer );
  
  // cluster size cut
  /* no cut */
  const TCut csize_cut = ("m_tracks.clusters.phi_size > 0 && m_tracks.clusters.z_size > 0" );
  
  // create histogram
  const TString hname = "h_pullsrph";
  auto h = new TH1F( hname, "", 100, -5, 5 );
  h->GetXaxis()->SetTitle( "#Delta#phi_{clus-truth}/#sigma#phi" );
  // h->GetXaxis()->SetTitle( "#Deltaz_{clus-truth}/#sigmaz" );
  
  // project tree
  tree->Project( hname, var, layer_cut && csize_cut );

  // fit
  h->Fit( "gaus", "0" );
  
  // draw
  auto cv = new TCanvas( "cv", "cv", 900, 900 );
  h->Draw();
  
  auto f = h->GetFunction("gaus");
  f->SetLineColor(2);
  f->Draw("same");
  
  cv->SaveAs( pdfFile );
}