doAna.cpp
Go to the documentation of this file.
1 //README
2 //This macro is an example of how to loop over the
3 //PandoraBeam TTree that was create by the LArsoft model ProtoDUNEelectronAnaTree_module.cc
4 //It would make a selection sample based on the tech note
5 //to run do $make -f Makefile, then execute the binary ./doAna mydatasample.root myoutputfile.root
6 
7 #include <iostream>
8 #include "include.h"
9 #include "AnalysisTree.h"
10 using namespace std;
11 
12 /***************modified box model parameters and function*****************/
13 double Rho = 1.383;//g/cm^3 (liquid argon density at a pressure 18.0 psia)
14 double betap = 0.212;//(kV/cm)(g/cm^2)/MeV
15 double alpha = 0.93;//parameter from ArgoNeuT experiment at 0.481kV/cm
16 double Wion = 23.6e-6;//parameter from ArgoNeuT experiment at 0.481kV/cm
18 double calib_factor = 5.6e-3;
19 double recombination =0.6417;
20 
21 double dEdx(float dqdx, float Efield) {
22  return (exp(dqdx*(betap/(Rho*Efield)*Wion))-alpha)/(betap/(Rho*Efield));
23 }
24 
25  TFile *ef=new TFile("/dune/data/users/higuera/data/SCE_DataDriven_180kV_v3.root");
26  TH3F *xneg=(TH3F*)ef->Get("Reco_ElecField_X_Neg");
27  TH3F *yneg=(TH3F*)ef->Get("Reco_ElecField_Y_Neg");
28  TH3F *zneg=(TH3F*)ef->Get("Reco_ElecField_Z_Neg");
29  TH3F *xpos=(TH3F*)ef->Get("Reco_ElecField_X_Pos");
30  TH3F *ypos=(TH3F*)ef->Get("Reco_ElecField_Y_Pos");
31  TH3F *zpos=(TH3F*)ef->Get("Reco_ElecField_Z_Pos");
32 
33 double tot_Ef(double xval,double yval,double zval){
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 }
47 
48 
49 int EventSelection( const char *filename, const char *outfile ){
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);
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 }
216 
217 int main( int argc, char *argv[] ){
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 }
TH3F * yneg
Definition: doAna.cpp:27
Double_t primaryShower_hit_Z[4847]
Definition: AnalysisTree.h:82
Double_t primaryShower_hit_X[4847]
Definition: AnalysisTree.h:80
TH3F * zneg
Definition: doAna.cpp:28
double betap
Definition: doAna.cpp:14
int EventSelection(const char *filename, const char *outfile)
Definition: doAna.cpp:49
TH3F * xneg
Definition: doAna.cpp:26
TH3F * xpos
Definition: doAna.cpp:29
Int_t primaryIsshower
Definition: AnalysisTree.h:54
Int_t primaryShower_nHits
Definition: AnalysisTree.h:76
constexpr T pow(T x)
Definition: pow.h:72
Int_t cerenkovStatus[2]
Definition: AnalysisTree.h:35
STL namespace.
Double_t primary_calZ[4847]
Definition: AnalysisTree.h:93
Double_t primarydQdx[4847]
Definition: AnalysisTree.h:90
Double_t primary_cal_pitch[4847]
Definition: AnalysisTree.h:94
string filename
Definition: train.py:213
TH3F * ypos
Definition: doAna.cpp:30
double dEdx(float dqdx, float Efield)
Definition: doAna.cpp:21
TH3F * zpos
Definition: doAna.cpp:31
Double_t primary_calY[4847]
Definition: AnalysisTree.h:92
double Rho
Definition: doAna.cpp:13
double alpha
Definition: doAna.cpp:15
double recombination
Definition: doAna.cpp:19
Double_t primaryShower_hit_Y[4847]
Definition: AnalysisTree.h:81
Double_t primaryStartDirection[3]
Definition: AnalysisTree.h:65
double calib_factor
Definition: doAna.cpp:18
Double_t beamtrackMomentum
Definition: AnalysisTree.h:38
int main(int argc, char *argv[])
Definition: doAna.cpp:217
Double_t beamtrackDir[3]
Definition: AnalysisTree.h:42
double normalisation_factor
Definition: doAna.cpp:17
E
Definition: 018_def.c:13
virtual Int_t GetEntry(Long64_t entry)
list x
Definition: train.py:276
double tot_Ef(double xval, double yval, double zval)
Definition: doAna.cpp:33
Vector cross(Vector const &a, Vector const &b)
Return cross product of two vectors.
double Wion
Definition: doAna.cpp:16
Int_t primaryNHits
Definition: AnalysisTree.h:56
Int_t primarynCal
Definition: AnalysisTree.h:89
Double_t primaryLength
Definition: AnalysisTree.h:59
Double_t primaryShower_hit_q[4847]
Definition: AnalysisTree.h:77
QTextStream & endl(QTextStream &s)
TFile * ef
Definition: doAna.cpp:25
Double_t primary_calX[4847]
Definition: AnalysisTree.h:91