search

e4nu / e4nuanalysiscode

Placeholder for all e4nu analyses
1 stars 2 forks source link

Brittany 1pi analysis #158

Open jtenavidal opened 7 months ago

jtenavidal commented 7 months ago

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

Quantify Background subtraction systematics with fiducial variations

jtenavidal commented 7 months ago

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 commented 7 months ago

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

jtenavidal commented 7 months ago
image
jtenavidal commented 7 months ago 
//////////////////////////////////////////////////////////////////////////////////

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 commented 7 months ago
image
jtenavidal commented 7 months ago

Add Reconstructed Ev from pion and outgoing electron momentum

jtenavidal commented 6 months ago

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

jtenavidal commented 6 months ago

For the plots,

1D Plot:

2D Plot:

For Julia:

jtenavidal commented 6 months ago

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 commented 6 months ago

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 commented 3 months ago

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: https://github.com/e4nu/e4nuanalysiscode/blob/027d82485321c83859499beb00bdb700771fbf43/src/analysis/MCCLAS6AnalysisI.cxx#L412

The above is for the MC, and for the data, https://github.com/e4nu/e4nuanalysiscode/blob/027d82485321c83859499beb00bdb700771fbf43/src/analysis/CLAS6AnalysisI.cxx#L294

They must have the same name

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

jtenavidal commented 3 months ago

To get the Costh:

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

jtenavidal commented 3 months ago 
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 commented 3 months ago

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 commented 3 months ago

Add Name of variable in the Tree as new observable

Change GetAxisName ( add new entry ) : https://github.com/e4nu/e4nuanalysiscode/blob/027d82485321c83859499beb00bdb700771fbf43/src/plotting/PlottingUtils.cxx#L65

Add new analysis key "pim" and within that if statement, add range for the new variables for each energy: https://github.com/e4nu/e4nuanalysiscode/blob/027d82485321c83859499beb00bdb700771fbf43/src/plotting/PlottingUtils.cxx#L246

Add definition of new variables and read from TBranch: https://github.com/e4nu/e4nuanalysiscode/blob/027d82485321c83859499beb00bdb700771fbf43/src/plotting/XSecUtils.cxx#L235 and also here https://github.com/e4nu/e4nuanalysiscode/blob/027d82485321c83859499beb00bdb700771fbf43/src/plotting/AcceptanceUtils.cxx#L103

jtenavidal commented 2 months ago

Location of data files:

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

jtenavidal commented 2 months ago
jtenavidal commented 2 months ago

To do for Julia

jtenavidal commented 2 months ago

Small update

jtenavidal commented 2 months ago

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 commented 2 months ago

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: https://github.com/e4nu/e4nuanalysiscode/blob/9a845737f4d4743d054f17fac6488fa8c622a6a7/src/utils/KinematicUtils.cxx#L80

For instance, for the electron case it could look like this, after https://github.com/e4nu/e4nuanalysiscode/blob/9a845737f4d4743d054f17fac6488fa8c622a6a7/src/analysis/MCCLAS6AnalysisI.cxx#L259

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 commented 2 months ago

There are new files available at the analysis_MC folder:

and the corresponding reconstructed ones

jtenavidal commented 2 months ago

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)

BrittanyKCohen96 commented 3 weeks ago

TDL after meeting 12 September 2024:

ElectronPT:

image

PionPT:

image

Total (HadDeltaPT):

image
jtenavidal commented 3 weeks ago

The data files are now available on the tau cluster:

Location: /storage/t3_data/tvjulia/e4nuanalysis_files/1pi_analysis/analised_data/ Files:

BrittanyKCohen96 commented 3 weeks ago

PionPT 1GeV Dipole breakdown w/ data:

image

ElectronPT 1GeV Dipole breakdown w/ data:

image
BrittanyKCohen96 commented 3 weeks ago

I opened a Drive containing all of the current PDF images and ROOT files where applicable: https://drive.google.com/drive/u/0/folders/1JqKBIaOI7pBHejtzVa3y7xcdjmH23OsW

BrittanyKCohen96 commented 3 weeks ago

Current progress and problems:

2GeV Dipole comparison (Angleqvshad):

image

4GeV Dipole comparison (Angleqvshad)

image

Example when plotting using the Rarita model on its own: ./plot_e4nuanalysis --mc_location $MCLOC --output_location $OUTLOC --output_name clas6analysis_1pim_1GeV --input_mc_files e4nuanalysis_1pim_G18_10a_Rarita_LFG_Q2_01_1GeV_eCarbon_NoRad --model_names G18_10a --title "CLAS6 C^{12} 1#pi^{-} 1GeV" --observable_list ECal --analysis_id "1pim" --data_location $DATALOC --input_data_file e4nuanalysis_1pimanalysis_e_on_1000060120_1161MeV_4MaxBkgMult_clas6data This command prints Warning in <TH1::TH1>: nbins is <=0 - set to nbins = 1 to the screen and then we get a segmentation fault. Other variables plot with very big cross sections instead (like ElectronPT).

Example when trying to plot Rarita and Dipole models on the same axis: ./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,e4nuanalysis_1pim_G18_10a_Rarita_LFG_Q2_01_1GeV_eCarbon_NoRad --model_names G18_10a,G18_10a_Rarita --title "CLAS6 C^{12} 1#pi^{-} 1GeV" --observable_list Angleqvshad --analysis_id "1pim" --data_location $DATALOC --input_data_file e4nuanalysis_1pimanalysis_e_on_1000060120_1161MeV_4MaxBkgMult_clas6data --store_root true In this case we get no errors and the graph looks like this:

image

I'm not sure what to make of this so I only did it for 1GeV for Angleqvshad as of now

jtenavidal commented 2 weeks ago

No FSI files location

jtenavidal commented 2 weeks ago

This paper has information on the transverse variables: https://arxiv.org/pdf/1512.05748 https://journals.aps.org/prd/pdf/10.1103/PhysRevD.108.053002 (For neutrinos) https://www.nature.com/articles/s41586-021-04046-5 (Nature paper e4nu collaboration)

image

Boosting angle - AlphaPT (HadAlphaPT) Angle between delta PT and electron

image

In the code: DeltaPT = | Vector_Electron_Transverse_Momentum + Vector_Highest_momentum_Pion_Transverse_Momentum|

HadDeltaPT = | Vector_Electron_Transverse_Momentum + Vector_Sum_Final_state_Hadrons_Transverse_Momentum| = (in your case, only 1 pion ) = | Vector_Electron_Transverse_Momentum + Vector_Pion_Transverse_Momentum|

you should get the cross section plots for:

jtenavidal commented 2 weeks ago
image image
jtenavidal commented 2 weeks ago
image
jtenavidal commented 2 weeks ago
jtenavidal commented 2 weeks ago

1 - Introduction Goals of the analysis. We want to study single pion production with electron-scattering to understand nuclear effects and how to correctly model it for neutrinos DUNE flux -

image

Near detector flux = neutrino peaks at 3 GeV Far detector - peaks at 1 and 4, shifts in MC prediction due to incorrect pion modelling will affect the results of the oscillation fit.

We fit far detector to near detector ratio and fit it to the MC prediction. We need good understanding of Ev, also cross section. Exclusive is important as it helps us determine bias in Ev!

https://arxiv.org/abs/2009.07228

2 - Signal definition and run details Introduce GENIE (event generator used widely for neutrino community. Can run electron scattering as well as neutrinos) Constraining electron models we can improve neutrinos. You should refer to this paper: https://arxiv.org/pdf/2009.07228

You are using the G2018 model which uses the Rosenbluth cross section with the Local Fermi Gas for QE scattering and the Empirical MEC model [21]. RES is modeled using the Berger- Sehgal model [23] and DIS reactions are modeled using Bodek and Yang[24]. The models are described in more detail in Ref(). This model set is formally marked as the G18 10a 02 11a configura- tion of GENIE v3.

Define your event selection

Run events on Deuterium, what energy Run events on Carbon what energy

Study Target Energy Cuts

True study D 1, 4 GeV 1pi (+or-), no photons, cut on Q2 C 1, 4 GeV 1pi (+or-), no photons, cut on Q2

Data analysis C 1,2 (4). 1pi (+or-), no photons, cut on Q2 + momentum cuts We account for gaps in the detector (it has six sectors) details on Sec. ().

Particle cuts defined here: https://github.com/e4nu/e4nuanalysiscode/blob/c2c85c429729316a11cb54706475d4c1c93eacfd/src/conf/AnalysisCutsI.cxx#L13

Angle cut

image

3 - GENIE True level studies On hydrogen (and deuterium), all events with W> m_pi have a pion in the final state. Neutrino experiments, we use heavy nuclear targets. In this case, we have significant nuclear effects which alter the final state particles we detect and their kinematics. In this case, separating pion production events by using reconstructed W is not possible. In neutrino physics we can't compute the true W (depends on neutrino kinematics) but we can only guess it from the detected final state particles (muon, pion). Equations. With electrons we can do both, we can use it to understand the bias in Wreco.

https://journals.aps.org/prd/pdf/10.1103/PhysRevD.104.072009

image

GENIE for deuterium doesnt compute FSI.

Here show GENIE true level studies.

for all energies. Mention in two lines what you got from the multiplicity study.

The method used shows that it is not possible to reconstruct Wreco accurately and cuts on Wreco should not be used to isolate pion production.

We could attempt to correct the data from the reconstructed space (Wreco) back to the true variable using Monte carlo. This assumes our model can correctly describe this bias. We will show in the next section that this is not possible because our models are not good/precise/ enough. We will quantify how bad it is using electron data because we can compute both true and reconstructed only with data. bias = ( Wreco - W ) / W (We can use electrons to understand if we can correctly quantify the bias using just neutrinos, or if the event generators are wrong.)

4 - CLAS (maybe in introduction) and data analysis

Empathise that it is a large acceptance spectrometer (angular acceptance is 2pi).

image

Put the picture of clas with the sectors

image

Here you should say again what targets and energies CLAS6 used. (He, C, Fe) and 1.161, 2.261 and 4.461 GeV. Say we focus on Carbon

data analysis: here maybe put a table saying what momentum cuts. We correct for gaps in the detector, background contamination (corrects for additional particles in detector gaps), we account for particles momentum smearing ( we account for it in MC prediction).

How we convert to cross section:

image

5 - Results with data 5a ) Simple variables - momentum and angle of particles 5b ) Beam energy reconstruction - ECal and RecoEvPion Here mention they are the same equations than used in Sec 3

5c ) Event invariant mass (W) bias Here mention they are the same equations than used in Sec 3

5d ) Deeper look into nuclear effects Transverse invariance variables and theta qpi

  1. Conclusions!! Great work!!
jtenavidal commented 2 weeks ago

No FSI files for 4GeV now also available at: /storage/t3_data/tvjulia/e4nuanalysis_files/1pi_analysis/analised_MC/e4nuanalysis_1pim_GEM21_11a_Dipole_LFG_NoFSI_Q2_08_4GeV_eCarbon_NoRad_true.root /storage/t3_data/tvjulia/e4nuanalysis_files/1pi_analysis/analised_MC/e4nuanalysis_1pim_GEM21_11a_Dipole_LFG_NoFSI_Q2_08_4GeV_eCarbon_NoRad_truereco.root

jtenavidal commented 2 weeks ago

Presentation2.pptx Slides