Functions | Variables
doAna.cpp File Reference
#include <iostream>
#include "include.h"
#include "AnalysisTree.h"

Go to the source code of this file.

Functions

double dEdx (float dqdx, float Efield)
 
double tot_Ef (double xval, double yval, double zval)
 
int EventSelection (const char *filename, const char *outfile)
 
int main (int argc, char *argv[])
 

Variables

double Rho = 1.383
 
double betap = 0.212
 
double alpha = 0.93
 
double Wion = 23.6e-6
 
double normalisation_factor =1.0
 
double calib_factor = 5.6e-3
 
double recombination =0.6417
 
TFile * ef =new TFile("/dune/data/users/higuera/data/SCE_DataDriven_180kV_v3.root")
 
TH3F * xneg =(TH3F*)ef->Get("Reco_ElecField_X_Neg")
 
TH3F * yneg =(TH3F*)ef->Get("Reco_ElecField_Y_Neg")
 
TH3F * zneg =(TH3F*)ef->Get("Reco_ElecField_Z_Neg")
 
TH3F * xpos =(TH3F*)ef->Get("Reco_ElecField_X_Pos")
 
TH3F * ypos =(TH3F*)ef->Get("Reco_ElecField_Y_Pos")
 
TH3F * zpos =(TH3F*)ef->Get("Reco_ElecField_Z_Pos")
 

Function Documentation

double dEdx ( float  dqdx,
float  Efield 
)

Definition at line 21 of file doAna.cpp.

21  {
22  return (exp(dqdx*(betap/(Rho*Efield)*Wion))-alpha)/(betap/(Rho*Efield));
23 }
betap
double betap
Definition: doAna.cpp:14
Rho
double Rho
Definition: doAna.cpp:13
alpha
double alpha
Definition: doAna.cpp:15
Wion
double Wion
Definition: doAna.cpp:16
int EventSelection ( const char *  filename,
const char *  outfile 
)

Definition at line 49 of file doAna.cpp.

49  {
50 
51  TFile *fyz_corr=new TFile("/dune/data/users/higuera/data/YZcalo_ProtoDUNE.root");
52  TH2F *YZ_neg=(TH2F*)fyz_corr->Get("run_YZ_negativeX_2");
53  TH2F *YZ_pos=(TH2F*)fyz_corr->Get("run_YZ_positiveX_2");
54 
55  //X Corr.
56  TFile *fx_corr=new TFile("/dune/data/users/higuera/data/Xcalo_ProtoDUNE.root");
57  TH2F *X_correction= (TH2F*)fx_corr->Get("X_correction_2");
58 
59  TFile *f_data = new TFile(filename,"READ");
60  TFile *f =new TFile(outfile,"recreate");
61 
62  TH1D *h_costheta = new TH1D("h_costheta",";cos#theta",100,-1,1);
63  TH1D *h_nhits = new TH1D("h_nhits",";",500,0,2000);
64  TH1D *h_dEdx = new TH1D("h_dEdx",";",50,0,10);
65  TH1D *h_sh_length = new TH1D("h_sh_length",";",100,0,400);
66  TH1D *h_sh_length_afcut = new TH1D("h_sh_length_afcut",";",100,0,400);
67  TH1D *h_P = new TH1D("h_P","; Beamline Momentum [GeV/c]",300,0,3);
68  TH1D *h_ratio = new TH1D("h_ratio",";(Shower Energy - Beamline Momentum)/Beamline Momentum",100,-1,1);
69  TH1D *h_Ecal = new TH1D("h_Ecal",";Shower Energy [Gev]",300,0,3);
70 
71  //========================================================
72  TTree *tree =(TTree*)f_data->Get("PandoraBeam");
73  AnalysisTree *signal = new AnalysisTree(tree);
74  int n_entries = tree->GetEntries();
75 
76 
77  for(long int jentry=0; jentry<n_entries;jentry++) {
78 
79  signal->GetEntry(jentry);
80  if( signal->beamtrackMomentum == -999.0 ) continue;
81  signal->beamtrackMomentum *= 1000.0;
82  //============= electrons =========
83  if( signal->cerenkovStatus[1] == 1 ){
84  //=======shower ===============================
85  if( signal->primaryIsshower ==1 ){
86  h_nhits->Fill(signal->primaryNHits);
87  h_sh_length->Fill(signal->primaryLength);
88  //select complete showers
89  if( signal->primaryNHits < 200 ) continue;
90 
91  TVector3 zaxis(signal->beamtrackDir[0],signal->beamtrackDir[1],signal->beamtrackDir[2]);
92  TVector3 shdir(signal->primaryStartDirection[0],signal->primaryStartDirection[1],signal->primaryStartDirection[2]);
93  double costheta = (shdir.Dot(zaxis))/(shdir.Mag()*zaxis.Mag());
94  h_costheta->Fill(costheta);
95  //remove potential cosmic-ray background
96  if( costheta < 0.9 ) continue;
97 
98  //calculate shower energy
99  double ave_Y=0;
100  double ave_Z=0;
101  int n_hits_y =0;
102  int n_hits_z =0;
103  //calculate average Y and average Z for cases we could not find a space point for a given hit
104  for( size_t h=0; h<signal->primaryShower_nHits; h++){
105  if( signal->primaryShower_hit_Y[h] != -999 ){
106  ave_Y += signal->primaryShower_hit_Y[h];
107  n_hits_y ++;
108  }
109  if(signal->primaryShower_hit_Z[h] != -999 ){
110  ave_Z += signal->primaryShower_hit_Z[h];
111  n_hits_z ++;
112  }
113  }
114  ave_Y /= n_hits_y;
115  ave_Z /= n_hits_z;
116  double YZcorrection_factor =1.0;
117  double Xcorrection_factor=1.0;
118  double Ecal =0;
119  //loop over hits in the shower
120  for( size_t h=0; h<signal->primaryShower_nHits; h++){
121  Xcorrection_factor=X_correction->GetBinContent(X_correction->FindBin(5834,signal->primaryShower_hit_X[h]));
122  if( signal->primaryShower_hit_Y[h] != -999 && signal->primaryShower_hit_Z[h] != -999 ){
123  if(signal->primaryShower_hit_X[h]<0) { YZcorrection_factor=YZ_neg->GetBinContent(YZ_neg->FindBin(5834,signal->primaryShower_hit_Z[h],signal->primaryShower_hit_Y[h])); }
124  if(signal->primaryShower_hit_X[h]>=0) { YZcorrection_factor=YZ_pos->GetBinContent(YZ_pos->FindBin(5834,signal->primaryShower_hit_Z[h],signal->primaryShower_hit_Y[h])); }
125  }
126  else {
127  if(signal->primaryShower_hit_X[h]<0) { YZcorrection_factor=YZ_neg->GetBinContent(5834,YZ_neg->FindBin(ave_Z,ave_Y)); }
128  if(signal->primaryShower_hit_X[h]>=0) { YZcorrection_factor=YZ_pos->GetBinContent(5834,YZ_pos->FindBin(ave_Z,ave_Y)); }
129  }
130  double E = (Xcorrection_factor*YZcorrection_factor*signal->primaryShower_hit_q[h]*Wion)/(calib_factor*recombination);
131  E *= normalisation_factor;
132  Ecal += E;
133  }
134  //this part calculates dE/dx at the beginning of the shower
135  int idx=0;
136  double dEdx_tmp[signal->primarynCal];
137  //tempral vectors to be used as vertex
138  //SCE corrections do not modify the shower vertex, therefore we need to use the first (X,Y,Z) of the dQ/dx calorimetry point as vertex
139  vector<double> tmp_x;
140  vector<double> tmp_y;
141  vector<double> tmp_z;
142  for( int k=0; k<signal->primarynCal-2; ++k){ //skip first two hits due to SCE
143  k +=2;
144  //skip calorimetry points where calorimetry information is not available
145  if( signal->primarydQdx[k] == 0 ) continue;
146  if( signal->primary_calX[k] == 0. && signal->primary_calY[k] == 0. && signal->primary_calZ[k] ==0 ) continue;
147  if( signal->primary_calZ[k] < 0 ) continue;
148 
149  tmp_x.push_back(signal->primary_calX[k]);
150  tmp_y.push_back(signal->primary_calY[k]);
151  tmp_z.push_back(signal->primary_calZ[k]);
152  }
153  //we found a "vertex"
154  for( int k=0; k<signal->primarynCal-2; ++k){
155  k += 2;
156  if( signal->primarydQdx[k] == 0 ) continue;
157  if( signal->primary_calZ[k] == 0. && signal->primary_calY[k] == 0. && signal->primary_calZ[k] ==0 ) continue;
158  if( signal->primary_cal_pitch[k] > 0 )
159  double YZcorrection_factor =1.0;
160  double Xcorrection_factor=1.0;
161  if(signal->primary_calX[k]<0) { YZcorrection_factor=YZ_neg->GetBinContent(YZ_neg->FindBin(5834,signal->primary_calZ[k],signal->primary_calY[k])); }
162  if(signal->primary_calX[k]>=0) { YZcorrection_factor=YZ_pos->GetBinContent(YZ_pos->FindBin(5834,signal->primary_calZ[k],signal->primary_calY[k])); }
163  if( signal->primary_calZ[k] < 0 ) continue; //no YZcorrection_factor;
164 
165  Xcorrection_factor=X_correction->GetBinContent(X_correction->FindBin(5834,signal->primary_calX[k]));
166  double corrected_dQdx=signal->primarydQdx[k]*normalisation_factor*Xcorrection_factor*YZcorrection_factor;
167 
168  double scaled_dQdx=corrected_dQdx/calib_factor;
169  double Efield=tot_Ef(signal->primary_calX[k],signal->primary_calY[k],signal->primary_calZ[k]);
170  double cali_dEdx=dEdx(scaled_dQdx, Efield); //this is your final dE/dx value
171  //calculate dE/dx at the first 4 cm of the shower
172  double newvtx_x = tmp_x.at(0);
173  double newvtx_y = tmp_y.at(0);
174  double newvtx_z = tmp_z.at(0);
175  double phi = TMath::ATan2(signal->primaryStartDirection[1],signal->primaryStartDirection[0]);
176  double theta = TMath::ACos(signal->primaryStartDirection[2]);
177  double x = 4.0*sin(theta)*cos(phi)+newvtx_x;
178  double y = 4.0*sin(theta)*sin(phi)+newvtx_y;
179  double z = 4.0*cos(theta)+newvtx_z;
180 
181  TVector3 x1(newvtx_x,newvtx_y,newvtx_z);
182  TVector3 x2(x,y,z);
183  TVector3 x0(signal->primary_calX[k],signal->primary_calY[k],signal->primary_calZ[k]);
184  TVector3 diff1 = x2-x1;
185  TVector3 diff2 = x1-x0;
186  TVector3 cross = diff1.Cross(diff2);
187  double r= cross.Mag()/diff1.Mag();
188  double d = sqrt( pow(signal->primary_calX[k]-newvtx_x,2)+ pow(signal->primary_calY[k]-newvtx_y,2) + pow(signal->primary_calZ[k]-newvtx_z,2));
189  //save dE/dx for the first 4 cm of the shower
190  if( r < 1.0 && d < 4.0){
191  dEdx_tmp[idx]=cali_dEdx;
192  idx ++;
193  }
194 
195 
196  }
197  //calculate the median
198  double dEdx_m =0;
199  dEdx_m = TMath::Median(idx,dEdx_tmp);
200 
201  if( dEdx_m == 0 ) continue;
202  h_dEdx->Fill(dEdx_m);
203  h_sh_length_afcut->Fill(signal->primaryLength);
204  h_P->Fill(signal->beamtrackMomentum/1000.0);
205  h_Ecal->Fill(Ecal/1000.0);
206  h_ratio->Fill((Ecal-signal->beamtrackMomentum)/signal->beamtrackMomentum);
207 
208  }//showers
209  }//electrons
210  }
211 
212  f->Write();
213  f->Close();
214  return 0;
215 }
AnalysisTree::primaryShower_hit_Z
Double_t primaryShower_hit_Z[4847]
Definition: AnalysisTree.h:82
AnalysisTree::primaryShower_hit_X
Double_t primaryShower_hit_X[4847]
Definition: AnalysisTree.h:80
keras_to_tensorflow.f
f
Definition: keras_to_tensorflow.py:162
AnalysisTree::primaryIsshower
Int_t primaryIsshower
Definition: AnalysisTree.h:54
truth_ana.theta
theta
Definition: truth_ana.py:138
AnalysisTree::primaryShower_nHits
Int_t primaryShower_nHits
Definition: AnalysisTree.h:76
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
bump_copyright.d
string d
Definition: bump_copyright.py:70
AnalysisTree::cerenkovStatus
Int_t cerenkovStatus[2]
Definition: AnalysisTree.h:35
AnalysisTree::primary_calZ
Double_t primary_calZ[4847]
Definition: AnalysisTree.h:93
AnalysisTree::primarydQdx
Double_t primarydQdx[4847]
Definition: AnalysisTree.h:90
AnalysisTree::primary_cal_pitch
Double_t primary_cal_pitch[4847]
Definition: AnalysisTree.h:94
reco_momentum_tuples.x2
x2
Definition: reco_momentum_tuples.py:83
train.filename
string filename
Definition: train.py:213
make_norm_csv.outfile
outfile
Definition: make_norm_csv.py:55
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
dEdx
double dEdx(float dqdx, float Efield)
Definition: doAna.cpp:21
AnalysisTree::primary_calY
Double_t primary_calY[4847]
Definition: AnalysisTree.h:92
recombination
double recombination
Definition: doAna.cpp:19
AnalysisTree::primaryShower_hit_Y
Double_t primaryShower_hit_Y[4847]
Definition: AnalysisTree.h:81
AnalysisTree::primaryStartDirection
Double_t primaryStartDirection[3]
Definition: AnalysisTree.h:65
calib_factor
double calib_factor
Definition: doAna.cpp:18
AnalysisTree::beamtrackMomentum
Double_t beamtrackMomentum
Definition: AnalysisTree.h:38
AnalysisTree::beamtrackDir
Double_t beamtrackDir[3]
Definition: AnalysisTree.h:42
normalisation_factor
double normalisation_factor
Definition: doAna.cpp:17
E
E
Definition: 018_def.c:13
AnalysisTree::GetEntry
virtual Int_t GetEntry(Long64_t entry)
train.x
list x
Definition: train.py:276
tot_Ef
double tot_Ef(double xval, double yval, double zval)
Definition: doAna.cpp:33
geo::vect::cross
Vector cross(Vector const &a, Vector const &b)
Return cross product of two vectors.
Definition: geo_vectors_utils.h:985
Wion
double Wion
Definition: doAna.cpp:16
AnalysisTree::primaryNHits
Int_t primaryNHits
Definition: AnalysisTree.h:56
AnalysisTree::primarynCal
Int_t primarynCal
Definition: AnalysisTree.h:89
AnalysisTree::primaryLength
Double_t primaryLength
Definition: AnalysisTree.h:59
AnalysisTree::primaryShower_hit_q
Double_t primaryShower_hit_q[4847]
Definition: AnalysisTree.h:77
make_resolution.tree
tree
Definition: make_resolution.py:22
plot_model.r
r
Definition: plot_model.py:57
reco_momentum_tuples.x1
x1
Definition: reco_momentum_tuples.py:82
AnalysisTree::primary_calX
Double_t primary_calX[4847]
Definition: AnalysisTree.h:91
int main ( int  argc,
char *  argv[] 
)

Definition at line 217 of file doAna.cpp.

217  {
218  cout<<"**************************** "<<endl;
219  cout<<"* WELCOME TO JAMAICA * "<<endl;
220  cout<<"* HAVE A NICE DAY * "<<endl;
221  cout<<"**************************** "<<endl;
222 
223  if( argc == 3 ) return EventSelection(argv[1],argv[2]);
224  else {
225  cout<<"please provide an input and output file "<<endl;
226  return 0;
227  }
228 
229 }
EventSelection
int EventSelection(const char *filename, const char *outfile)
Definition: doAna.cpp:49
freeze_graph.argv
argv
Definition: freeze_graph.py:218
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
double tot_Ef ( double  xval,
double  yval,
double  zval 
)

Definition at line 33 of file doAna.cpp.

33  {
34  if(xval>=0){
35  double ex=0.5+0.5*xpos->GetBinContent(xpos->FindBin(xval,yval,zval));
36  double ey=0.5*ypos->GetBinContent(ypos->FindBin(xval,yval,zval));
37  double ez=0.5*zpos->GetBinContent(zpos->FindBin(xval,yval,zval));
38  return sqrt(ex*ex+ey*ey+ez*ez);
39  }
40  else{
41  double ex=0.5+0.5*xneg->GetBinContent(xneg->FindBin(xval,yval,zval));
42  double ey=0.5*yneg->GetBinContent(yneg->FindBin(xval,yval,zval));
43  double ez=0.5*zneg->GetBinContent(zneg->FindBin(xval,yval,zval));
44  return sqrt(ex*ex+ey*ey+ez*ez);
45  }
46 }
yneg
TH3F * yneg
Definition: doAna.cpp:27
zneg
TH3F * zneg
Definition: doAna.cpp:28
xneg
TH3F * xneg
Definition: doAna.cpp:26
xpos
TH3F * xpos
Definition: doAna.cpp:29
ypos
TH3F * ypos
Definition: doAna.cpp:30
zpos
TH3F * zpos
Definition: doAna.cpp:31

Variable Documentation

double alpha = 0.93

Definition at line 15 of file doAna.cpp.

double betap = 0.212

Definition at line 14 of file doAna.cpp.

double calib_factor = 5.6e-3

Definition at line 18 of file doAna.cpp.

TFile* ef =new TFile("/dune/data/users/higuera/data/SCE_DataDriven_180kV_v3.root")

Definition at line 25 of file doAna.cpp.

double normalisation_factor =1.0

Definition at line 17 of file doAna.cpp.

double recombination =0.6417

Definition at line 19 of file doAna.cpp.

double Rho = 1.383

Definition at line 13 of file doAna.cpp.

double Wion = 23.6e-6

Definition at line 16 of file doAna.cpp.

TH3F* xneg =(TH3F*)ef->Get("Reco_ElecField_X_Neg")

Definition at line 26 of file doAna.cpp.

TH3F* xpos =(TH3F*)ef->Get("Reco_ElecField_X_Pos")

Definition at line 29 of file doAna.cpp.

TH3F* yneg =(TH3F*)ef->Get("Reco_ElecField_Y_Neg")

Definition at line 27 of file doAna.cpp.

TH3F* ypos =(TH3F*)ef->Get("Reco_ElecField_Y_Pos")

Definition at line 30 of file doAna.cpp.

TH3F* zneg =(TH3F*)ef->Get("Reco_ElecField_Z_Neg")

Definition at line 28 of file doAna.cpp.

TH3F* zpos =(TH3F*)ef->Get("Reco_ElecField_Z_Pos")

Definition at line 31 of file doAna.cpp.

Generated by   doxygen 1.8.11