Brittany 1pi analysis #158

jtenavidal · 2024-02-21T13:08:31Z

References of interest:
We want to have a look at what is done in neutrino physics in the case of pion production.

https://arxiv.org/pdf/1011.3572.pdf
https://cds.cern.ch/record/2623949?ln=ca Page 118 - 131 ( good overview of variables of interest )
Some of these variables are not available in the code yet. We want to add the missing variables here: https://github.com/e4nu/e4nuanalysiscode/blob/master/src/utils/KinematicUtils.cxx and https://github.com/e4nu/e4nuanalysiscode/blob/master/src/utils/KinematicUtils.h
We will have to store them in the analysed ROOT file
Gide Brittany on the plotting scripts. Get cross-section measurement for different energies and targets
Plotting scripts. As of now, we plot RES (Delta + Above Delta contribution). We want to have the option to turn off Non-RES contribution as it is very off.

Quantify Background subtraction systematics with fiducial variations

Discuss with Adi and Larry what approach they used. We will have to run it from the data. Need to understand limitations

jtenavidal · 2024-02-22T08:14:49Z

Exercise with ROOT file
So following our chat it will be great if you could take a look at a few kinematic variables, save their distribution and think a bit about the reason they look this way:

For events with 1 negative pion take a look at:

The momentum of the outgoing electron
The scattering angle of the electron
The azimuthal angle (phi) of the electron

The momentum of the outgoing pion
The scattering angle of the outgoing pion
The azimuthal angel (phi) of the outgoing pion

W
Q2

How many events are QE/MEC/RES/DIS

The momentum of the outgoing electron vs. its scattering angle
Q2 vs. W

For the 2d histograms it’s advisable to draw them with the "colz” option

A reminder: The cross section can be calculated as the number of events divided by the luminosity or take the events in your histogram and normalize it to the cross section like Julia showed you from the spline file

inclusive cross section you take from the spline root file by taking Eval of the relevant TGraph in the value of the incoming energy.

histogram->Scale(inclusive cross section / total generated number of events)

jtenavidal · 2024-02-22T08:25:32Z

To plot the gst W branch for 1pim only:
gst->Draw("W","nfpim==1&&nfpip==0&&nfem==0")

jtenavidal · 2024-02-22T09:09:54Z

jtenavidal · 2024-02-22T09:17:34Z

//////////////////////////////////////////////////////////////////////////////////

int plot_GENIEgst() {
  std::string MC_file_name = "/Users/juliatenavidal/Downloads/GENIE_MC_File_example.gst.root" ;
  TCanvas * c1 = new TCanvas("c1","c1",800,600);

  TFile * file_true_MC = new TFile((MC_file_name).c_str(),"ROOT");
  TTree * my_gst_tree = (TTree*)file_true_MC->Get("gst");

  if( !file_true_MC ) { std::cout << "ERROR: the "<< file_true_MC << " does not exist." <<std::endl; return 0;}
  if( !my_gst_tree) std::cout << " Tree does not exist " << std::endl;

//  TH1D * myh = new TH1D( "myh", "title", 100, 0, 50) ;
  TH1D * myh_1p0pi = new TH1D( "myh_1p0pi", "myh_1p0pi", 100, 0, 5) ;
  TH1D * myh_1p0pi_qel = new TH1D( "myh_1p0pi_qel", "myh_1p0pi_qel", 100, 0, 5) ;
  TH2D * myh = new TH2D( "myh", "title", 100, 0, 4, 100,0,4) ;

   // OBSERVABLE DEFINITION:
  int nf ;
  double Ev, W, Q2 ; // Electron energy
  Int_t pdgf[120] ;
  double pxf[120],pyf[120],pzf[120],Ef[120];
  double pl ;
  bool res, dis, qel;
  my_gst_tree->SetBranchAddress("Ev", &Ev);
  my_gst_tree->SetBranchAddress("nf", &nf);
  my_gst_tree->SetBranchAddress("pl", &pl);
  my_gst_tree->SetBranchAddress("pdgf", &pdgf);
  my_gst_tree->SetBranchAddress("pxf", &pxf);
  my_gst_tree->SetBranchAddress("pyf", &pyf);
  my_gst_tree->SetBranchAddress("pzf", &pzf);
  my_gst_tree->SetBranchAddress("Ef", &Ef);
  my_gst_tree->SetBranchAddress("qel", &qel);
  my_gst_tree->SetBranchAddress("res", &res);
  my_gst_tree->SetBranchAddress("dis", &dis);
  my_gst_tree->SetBranchAddress("W", &W);
  my_gst_tree->SetBranchAddress("Q2", &Q2);
  for( int j = 0 ; j < NEntries ; ++j ) {
    my_gst_tree->GetEntry(j) ;
    int number_pip = 0 ;
    int number_p = 0 ;
    int number_pim = 0 ;
    int number_em = 0 ;

    double mom = 0 ;
    double phi = 0 ;
    double theta = 0 ;
    for( unsigned int k = 0 ; k < nf ; ++k ){

      TLorentzVector my4mom;
      my4mom.SetPxPyPzE (pxf[k], pyf[k], pzf[k], Ef[k]);

      if(  pdgf[k] == -211 ) { mom = my4mom.P(); theta = my4mom.Theta(); phi = my4mom.Phi() ; }

      if( pdgf[k] == 211 ) ++number_pip; // number_pip = number_pip + 1;
      if( pdgf[k] == -211 ) ++number_pim;
      if( pdgf[k] == 2212 ) ++number_p;
      if( pdgf[k] == 22 ) ++number_em; // photon
    }
    if( number_pip == 0 && number_pim ==1 && number_em == 0 ) myh->Fill(W,Q2);
  }

  myh->SetLineColor(kBlack);
  myh_1p0pi->SetLineColor(kBlue);
  myh_1p0pi_qel->SetLineColor(kRed);

  myh->Draw("hist colz");
  myh_1p0pi->Draw("hist same");
  myh_1p0pi_qel->Draw("hist same");

  c1->SaveAs("test.pdf");
  return 0;

jtenavidal · 2024-02-22T11:02:12Z

jtenavidal · 2024-03-06T16:21:08Z

Add Reconstructed Ev from pion and outgoing electron momentum

Add new function in Kinematics.h (definition in Kinematics.cxx): https://github.com/e4nu/e4nuanalysiscode/blob/master/src/utils/KinematicUtils.h
We want to store the results in the analysis output root file. To do so, we need to modify the analysis code for the MC and the data. Both have a Store Tree function which will define the branches and add the event content in the branch. See:
- e4nuanalysiscode/src/analysis/CLAS6AnalysisI.cxx
  
  Line 132 in afa1ec3
  
  bool CLAS6AnalysisI::StoreTree(Event event){
- e4nuanalysiscode/src/analysis/MCCLAS6AnalysisI.cxx
  
  Line 300 in afa1ec3
  
  bool MCCLAS6AnalysisI::StoreTree(Event event){
The next step is to run the analysis with the MC and CLAS data files. As you don't have access, I will run it once the implementation is complete. For testing, you can use a smaller root file ().

./e4nuanalysis --conf-file ConfFiles/mc_conf/clas6mc_1pimanalysis_eC12_1GeV.txt --root-file example_root_file.gst.root --output-file test_RecoEvPion --analysis-type ComputeTrueAccCorr --xsec-file example_xsec.root

jtenavidal · 2024-03-10T09:22:45Z

Interesting paper
https://arxiv.org/pdf/2102.03346.pdf

jtenavidal · 2024-03-11T14:38:01Z

For the plots,

1D Plot:

Stack plot of distribution for true W and Reconstructed W
For reconstructed W distribution, get breakdown into processes

2D Plot:

For Julia:

QEL events can produce pions if the proton created at the interaction vertex has enough energy to excite a nucleon from the target.
Give Brittany Monte Carlo files for 1GeV, 2GeV and 4 GeV initial electron energy. She will repeat the plots above for each case. For now she has the 1GeV file.

jtenavidal · 2024-03-11T14:53:39Z

Neutrino summer school recordings
https://npc.fnal.gov/inss/
Day 1 day 2, are the most relevant: https://indico.cern.ch/event/1011452/timetable/

jtenavidal · 2024-03-13T13:27:34Z

Here you can find my latest slides:
https://indico.jlab.org/event/829/contributions/14080/attachments/10742/16277/CLAS6_introduction_analysis_1p1pi_JTena_March2024.pdf

It is a good summary of what you will do, you will just look at 1pi events instead.

There's also work done by collaborators of our on similar topologies using different data. This is the latest update of their work: https://indico.jlab.org/event/829/contributions/14082/attachments/10739/16269/2024_cfogler_NPWG_v1.pdf

jtenavidal · 2024-06-30T08:39:04Z

Include New observables in data format
We want to add the reconstructed Ev for pion production, and the reconstructed W. To do so, we need two new functions here:

https://github.com/e4nu/e4nuanalysiscode/blob/master/src/utils/KinematicUtils.cxx

and the definition in the .h file.

double utils::GetRecoEvPionProduction( const TVector3 p1 /*out electron*/, const TVector3 p2/*pion*/ );

double utils::GetRecoWPionProduction( const TVector3 p1 /*out electron*/, const TVector3 p2/*pion*/ );

and maybe also
double utils::GetRecoQ2PionProduction( const TVector3 p1 /*out electron*/, const TVector3 p2/*pion*/ );

Once this is done, we just need to compute these for all events that pass our selection criteria. This can be done by computing the variables (calling the functions above) in this function as well as storing the new variables in the TTree with a new name:

e4nuanalysiscode/src/analysis/MCCLAS6AnalysisI.cxx

Line 412 in 027d824

if( n == true ) {

The above is for the MC, and for the data,

e4nuanalysiscode/src/analysis/CLAS6AnalysisI.cxx

Line 294 in 027d824

kAnalysisTree -> Branch( "InputROOTFile", &InputROOTFile, "InputROOTFile/C");

They must have the same name

Once this is done, I will run the code and send you the files

jtenavidal · 2024-06-30T11:23:40Z

To get the Costh:

brCosthl = TMath::Cos( k2.Vect().Angle(k1.Vect()) );

jtenavidal · 2024-06-30T12:00:00Z

TLorentzVector pi_mom(0,0,0,0) ;
  if( topology_has_pip ) {
    double max_mom = 0 ;
    for( unsigned int i = 0 ; i < hadron_map[conf::kPdgPiP].size() ; ++i ) {
      if( hadron_map[conf::kPdgPiP][i].P() > max_mom ) {
	max_mom = hadron_map[conf::kPdgPiP][i].P() ;
	pi_max = hadron_map[conf::kPdgPiP][i] ;
      }
    }
  } else  if( topology_has_pim ) {
    double max_mom = 0 ;
    for( unsigned int i = 0 ; i < hadron_map[conf::kPdgPiM].size() ; ++i ) {
      if( hadron_map[conf::kPdgPiM][i].P() > max_mom ) {
	max_mom = hadron_map[conf::kPdgPiM][i].P() ;
	pi_max = hadron_map[conf::kPdgPiM][i] ;
      }
    }

jtenavidal · 2024-07-04T09:06:52Z

To plot the MC predictions out of the analysis code you need

src/plotting_apps/plot_e4nuanalysis.cxx

Example:

MCLOC=/storage/t3_data/tvjulia/e4nuanalysis_files/1pi_analysis/analised_MC/
OUTLOC=/storage/t3_data/brittany/e4nuanalysis_files/1pi_analysis/Plots

./plot_e4nuanalysis --mc_location $MCLOC --output_location $OUTLOC --output_name clas6analysis_1pim_1GeV --input_mc_files e4nuanalysis_1pim_G18_10a_Dipole_LFG_Q2_01_1GeV_eCarbon_NoRad --model_names G18_10a --title "CLAS6 C^{12} 1#pi^{-} 1GeV" --observable_list "pim_mom"  --analysis_id "1pim"

scp -o "ProxyJump jtena@login.jlab.org " -r jtena@ifarm1802.jlab.org:/storage/t3_data/tvjulia/e4nuanalysis_files/1pi_analysis/Plots/TotalXSec/clas6analysis_1pim_1GeV_with_breakdown_dxsec_dpim_mom.pdf .

jtenavidal · 2024-07-04T09:18:29Z

Add Name of variable in the Tree as new observable

Change GetAxisName ( add new entry ) :

e4nuanalysiscode/src/plotting/PlottingUtils.cxx

Line 65 in 027d824

    
           else if ( observable == "Angleqvshad") { x_axis = "#theta_{#vec{q}#dot#vec{p}_{had}}[deg]"; y_axis = "d#sigma/d#theta_{#vec{q}#dot#vec{p}_{had}} #left[#mub (deg)^{-1}#right]"; }

Add new analysis key "pim" and within that if statement, add range for the new variables for each energy:

e4nuanalysiscode/src/plotting/PlottingUtils.cxx

Line 246 in 027d824

else if( analysis_key == "1p1pip" ) {

Add definition of new variables and read from TBranch:

e4nuanalysiscode/src/plotting/XSecUtils.cxx

Line 235 in 027d824

double AdlerAngleThetaP, AdlerAnglePhiP, AdlerAngleThetaPi, AdlerAnglePhiPi ;

and also here

e4nuanalysiscode/src/plotting/AcceptanceUtils.cxx

Line 103 in 027d824

double pim_mom,pim_theta,pim_phi;

jtenavidal · 2024-07-21T10:26:59Z

Location of data files:

/storage/t3_data/tvjulia/e4nuanalysis_files/1pi_analysis/analised_data

jtenavidal · 2024-07-25T11:45:48Z

jtenavidal · 2024-07-25T14:25:10Z

To do for Julia

Send NoFSI files
Review 4 GeV normalization
2 GeV G18 empty re-run

jtenavidal · 2024-07-31T11:08:52Z

Small update

The plotting script should now be able to run without data when using your branch
I am re-running the G18_10a model at 2GeV, we should have the updated results in an hour from now
The Rarita MC files are missing. I have to re-generate them and then run them through the e4nu analysis code. Once this is finished, we should be able to see the results in the final root files
I still don't understand what is wrong for the 4GeV files. Need to investigate it further

jtenavidal · 2024-07-31T13:00:16Z

New G18_10a files at 2GeV available in
/storage/t3_data/tvjulia/e4nuanalysis_files/1pi_analysis/analised_MC/e4nuanalysis_1pim_G18_10a_Dipole_LFG_Q2_04_2GeV_eCarbon_NoRad_truereco.root
and
/storage/t3_data/tvjulia/e4nuanalysis_files/1pi_analysis/analised_MC/e4nuanalysis_1pim_G18_10a_Dipole_LFG_Q2_04_2GeV_eCarbon_NoRad_true.root

jtenavidal · 2024-07-31T13:57:35Z

Ok I double checked. We do not have ElectronPT or pimom_T. It would be great if you could add them. To get the transverse momentum, you have to consider that the beam goes in the z direction, as you said. Here is an example function you can use:

e4nuanalysiscode/src/utils/KinematicUtils.cxx

Line 80 in 9a84573

TVector3 utils::GetPT( const TVector3 p ) {

For instance, for the electron case it could look like this, after

e4nuanalysiscode/src/analysis/MCCLAS6AnalysisI.cxx

Line 259 in 9a84573

TLorentzVector out_mom = event.GetOutLepton4Mom();

double ElectronPT = utils::GetPT( out_mom.Vect() ).Mag() ; # in line 265
and then just add it in the branches like we did for the other variables (for data and MC)

The same should be done for pip (L369) and pim (L389)

jtenavidal · 2024-07-31T16:36:57Z

There are new files available at the analysis_MC folder:

e4nuanalysis_1pim_G18_10a_Rarita_LFG_Q2_01_1GeV_eCarbon_NoRad_true.root
e4nuanalysis_1pim_G18_10a_Rarita_LFG_Q2_04_2GeV_eCarbon_NoRad_true.root
e4nuanalysis_1pim_G18_10a_Rarita_LFG_Q2_08_4GeV_eCarbon_NoRad_true.root
e4nuanalysis_1pim_G18_10a_Dipole_LFG_Q2_04_2GeV_eCarbon_NoRad_true.root

and the corresponding reconstructed ones

jtenavidal · 2024-07-31T16:41:13Z

I will update the NOFSI ones as soon as I can, not sure I will have time tomorrow but you have plenty for now

Let me know if the new 4GeV ones make sense or look also wrong

Other things to do are adding the new axes in the plotting (like Adi mentioned)

jtenavidal changed the title ~~Jessica 1pi analysis~~ Britany 1pi analysis Feb 21, 2024

jtenavidal assigned BrittanyKCohen96 Feb 22, 2024

jtenavidal changed the title ~~Britany 1pi analysis~~ Brittany 1pi analysis Mar 11, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Brittany 1pi analysis #158

Brittany 1pi analysis #158

jtenavidal commented Feb 21, 2024 •

edited by BrittanyKCohen96

Loading

jtenavidal commented Feb 22, 2024 •

edited

Loading

jtenavidal commented Feb 22, 2024

jtenavidal commented Feb 22, 2024

jtenavidal commented Feb 22, 2024

jtenavidal commented Feb 22, 2024

jtenavidal commented Mar 6, 2024 •

edited

Loading

jtenavidal commented Mar 10, 2024 •

edited

Loading

jtenavidal commented Mar 11, 2024 •

edited

Loading

jtenavidal commented Mar 11, 2024

jtenavidal commented Mar 13, 2024 •

edited

Loading

jtenavidal commented Jun 30, 2024 •

edited

Loading

jtenavidal commented Jun 30, 2024

jtenavidal commented Jun 30, 2024 •

edited

Loading

jtenavidal commented Jul 4, 2024 •

edited by BrittanyKCohen96

Loading

jtenavidal commented Jul 4, 2024

jtenavidal commented Jul 21, 2024

jtenavidal commented Jul 25, 2024 •

edited by BrittanyKCohen96

Loading

jtenavidal commented Jul 25, 2024 •

edited

Loading

jtenavidal commented Jul 31, 2024 •

edited

Loading

jtenavidal commented Jul 31, 2024

jtenavidal commented Jul 31, 2024 •

edited

Loading

jtenavidal commented Jul 31, 2024

jtenavidal commented Jul 31, 2024

Brittany 1pi analysis #158

Brittany 1pi analysis #158

Comments

jtenavidal commented Feb 21, 2024 • edited by BrittanyKCohen96 Loading

jtenavidal commented Feb 22, 2024 • edited Loading

jtenavidal commented Feb 22, 2024

jtenavidal commented Feb 22, 2024

jtenavidal commented Feb 22, 2024

jtenavidal commented Feb 22, 2024

jtenavidal commented Mar 6, 2024 • edited Loading

jtenavidal commented Mar 10, 2024 • edited Loading

jtenavidal commented Mar 11, 2024 • edited Loading

jtenavidal commented Mar 11, 2024

jtenavidal commented Mar 13, 2024 • edited Loading

jtenavidal commented Jun 30, 2024 • edited Loading

jtenavidal commented Jun 30, 2024

jtenavidal commented Jun 30, 2024 • edited Loading

jtenavidal commented Jul 4, 2024 • edited by BrittanyKCohen96 Loading

jtenavidal commented Jul 4, 2024

jtenavidal commented Jul 21, 2024

jtenavidal commented Jul 25, 2024 • edited by BrittanyKCohen96 Loading

jtenavidal commented Jul 25, 2024 • edited Loading

jtenavidal commented Jul 31, 2024 • edited Loading

jtenavidal commented Jul 31, 2024

jtenavidal commented Jul 31, 2024 • edited Loading

jtenavidal commented Jul 31, 2024

jtenavidal commented Jul 31, 2024

jtenavidal commented Feb 21, 2024 •

edited by BrittanyKCohen96

Loading

jtenavidal commented Feb 22, 2024 •

edited

Loading

jtenavidal commented Mar 6, 2024 •

edited

Loading

jtenavidal commented Mar 10, 2024 •

edited

Loading

jtenavidal commented Mar 11, 2024 •

edited

Loading

jtenavidal commented Mar 13, 2024 •

edited

Loading

jtenavidal commented Jun 30, 2024 •

edited

Loading

jtenavidal commented Jun 30, 2024 •

edited

Loading

jtenavidal commented Jul 4, 2024 •

edited by BrittanyKCohen96

Loading

jtenavidal commented Jul 25, 2024 •

edited by BrittanyKCohen96

Loading

jtenavidal commented Jul 25, 2024 •

edited

Loading

jtenavidal commented Jul 31, 2024 •

edited

Loading

jtenavidal commented Jul 31, 2024 •

edited

Loading