LIDTrainingDedx.cxx
Go to the documentation of this file.
1 #include "TChain.h"
2 #include "TFile.h"
3 #include "TH1D.h"
4 
5 #include <iostream>
6 #include <fstream>
7 #include <string>
9 #include <math.h>
10 
11 namespace slidt{
12 
13  /////////////////////////////////////////////////////////////////////////////////////////
14 
16  {
17 
18  // if(fActivated)
19  // return;
20 
21  fChain->SetBranchAddress("showerTrueNuCCNC", &showerTrueNuCCNC);
22  fChain->SetBranchAddress("showerTrueNuMode", &showerTrueNuMode);
23  fChain->SetBranchAddress("showerTrueNuPdg", &showerTrueNuPdg);
24  fChain->SetBranchAddress("showerNPlane", &showerNPlane);
25  fChain->SetBranchAddress("showerStart[3]", showerStart);
26  fChain->SetBranchAddress("showerStop[3]", showerStop);
27  fChain->SetBranchAddress("showerEnergy", &showerEnergy);
28  fChain->SetBranchAddress("showerIsFid", &showerIsFid);
29  fChain->SetBranchAddress("showerPlaneDedx[200]", showerPlaneDedx);
30  fChain->SetBranchAddress("showerTransCellDedx[20]", showerTransCellDedx);
31  fChain->SetBranchAddress("showerTrueDang", &showerTrueDang);
32  fChain->SetBranchAddress("showerTruePdg", &showerTruePdg);
33  fChain->SetBranchAddress("showerTrueP4[4]", showerTrueP4);
34  fChain->SetBranchAddress("showerDir[3]", showerDir);
35 
36  fActivated = true;
37 
38  }// end of ActivateSetBranchAddresses
39 
40 
41  /////////////////////////////////////////////////////////////////////////////////////////
42 
43  void SaveHists(const int inPdg,
44  const std::string outDir,
45  const int* inNuMode)
46  {
47  std::string partname, mode = "";
48 
49  switch(inPdg){
50  case 11:
51  partname = "electron";
52  break;
53  case 22:
54  partname = "gamma";
55  break;
56  case 13:
57  partname = "muon";
58  break;
59  case 111:
60  partname = "pi0";
61  break;
62  case 211:
63  partname = "pion";
64  break;
65  case 2212:
66  partname = "proton";
67  break;
68  case 2112:
69  partname = "neutron";
70  break;
71  default:
72  std::cerr<<"The specified PDG particle is not included in PID.\n";
73  std::abort();
74  }
75 
76  if( inNuMode != NULL ){
77  switch(*inNuMode){
78  case 0:
79  mode = "_qe";
80  break;
81  case 1:
82  mode = "_res";
83  break;
84  case 2:
85  mode = "_dis";
86  break;
87  case 3:
88  mode = "_coh";
89  break;
90  default:
91  std::cerr<<"Invalid neutrino interaction mode. Valid choices are 0,1,2,3.\n";
92  std::abort();
93  }
94  }
95 
96  std::string outFileName = outDir + partname + mode + "_dedxhists.root";
97  std::cout<<"output file name is "<<outFileName<<std::endl;
98  TFile fout(outFileName.c_str(), "RECREATE");
99 
100  for( int idet=0; idet < kNumXYRegion; ++idet){
101  for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
102  for( int iplane=0; iplane < kNumLongitudinalPlane; ++iplane){
103  hlongdedx[idet][iebin][iplane]->Write();
104  }
105  }
106  }
107  for( int idet=0; idet < kNumXYRegion; ++idet){
108  for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
109  for( int iplane=0; iplane < kNumTransversePlane; ++iplane){
110  htransdedx[idet][iebin][iplane]->Write();
111  }
112  }
113  }
114  fout.Close();
115  }// End of SaveHists
116 
117 
118  /////////////////////////////////////////////////////////////////////////////////////////
119 
120  void Reset()
121  {
122  for( int idet=0; idet < kNumXYRegion; ++idet){
123  for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
124  for( int iplane=0; iplane < kNumTransversePlane; ++iplane){
125  htransdedx[idet][iebin][iplane]->Reset();
126  }
127  }
128  }
129  for( int idet=0; idet < kNumXYRegion; ++idet){
130  for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
131  for( int iplane=0; iplane < kNumLongitudinalPlane; ++iplane){
132  hlongdedx[idet][iebin][iplane]->Reset();
133  }
134  }
135  }
136 
137  fChain->Reset();
138  }// Reset all histograms
139 
140  /////////////////////////////////////////////////////////////////////////////////////////
141 
143  {
144  if( fInitialized )
145  return;
146 
147  for( int idet=0; idet < kNumXYRegion; ++idet){
148  for( int iebin=0; iebin < kNumEnergyBin; ++iebin){
149 
150  for( int iplane=0; iplane < kNumLongitudinalPlane; ++iplane){
151  char hname[100], htitle[200];
152  sprintf(hname,"hlongdedx_det%d_ebin%d_plane%d", idet, iebin, iplane);
153  sprintf(htitle,"Longitudinal dE/dx: Det Region %d, Energy Bin %d, Plane %d;dE/dX (GeV/cm);Events", idet, iebin, iplane);
154  hlongdedx[idet][iebin][iplane] = new TH1D(hname, htitle, 50, 0, 0.050);
155  }// end loop over longitudinal planes
156 
157  for( int iplane=0; iplane < kNumTransversePlane; ++iplane){
158  char hname[100], htitle[200];
159  sprintf(hname,"htransdedx_det%d_ebin%d_plane%d", idet, iebin, iplane);
160  sprintf(htitle,"Transverse dE/dx: Det Region %d, Energy Bin %d, Plane %d;dE/dX (GeV/cm);Events", idet, iebin, iplane);
161  float hmax = 0.01;
162  if( iplane==0 )
163  hmax = 0.01;
164  if( iplane==1 )
165  hmax = 0.01;
166  if( iplane==2 )
167  hmax = 0.002;
168  if( iplane>=3 && iplane<=5 )
169  hmax = 0.001;
170  if( iplane>=6 && iplane<=9 )
171  hmax = 0.0005;
172  if( iplane>=10 && iplane<=14 )
173  hmax = 0.00010;
174  if( iplane>=15 )
175  hmax = 0.00005;
176  htransdedx[idet][iebin][iplane] = new TH1D(hname, htitle, 500, 0, hmax);
177 
178  }// end loop over transverse planes
179 
180  }// end loop over energy bins
181  }// end loop over detregions
182  fInitialized = true;
183  return;
184  }// end of InitializeHists
185 
186  /////////////////////////////////////////////////////////////////////////////////////////
187 
188  void FillHists( const char* inFnames,
189  const char* inTreename,
190  const char* outDir,
191  const int inPdg,
192  const float inDangLower,
193  const float inDangUpper,
194  const int* inNuPdg,
195  const int* inNuMode,
196  const int* inNuCC,
197  const double* eTrueFrac)
198  {
199 
200  fChain = new TChain(inTreename);
201  fChain->Add(inFnames);
202 
204  InitializeHists();
205 
206  int nentries = fChain->GetEntries();
207  std::cout<<"Number of entries "<<nentries<<std::endl;
208 
209  for(int e = 0; e < nentries; ++e){
210 
211  fChain->GetEntry(e);
212 
213  // skip to next event if any of the required cuts are not
214  // passed
215 
216  if( abs(showerTruePdg) != inPdg ||
217  showerTrueDang <= inDangLower ||
218  showerTrueDang >= inDangUpper ||
219  showerIsFid != 1 )
220  continue;
221 
222  // the following three cuts don't apply to all particles.
223  // see what's been provided.
224  if( inNuPdg != NULL ){
225  if( showerTrueNuPdg != *inNuPdg )
226  continue;
227  }
228 
229  if( inNuMode != NULL ){
230  if( showerTrueNuMode != *inNuMode )
231  continue;
232  }
233 
234  if( inNuCC != NULL ){
235  if( showerTrueNuCCNC != *inNuCC )
236  continue;
237  }
238 
239  if( eTrueFrac != NULL ){
240  if( showerEnergy/showerTrueP4[3] <= *eTrueFrac )
241  continue;
242  }
243 
244  // check which energy bin the shower falls in
245  static const double* eBins = kEnergyBins;
246  int nEBins = kNumEnergyBin;
247  int showerEBin = -1;
248  for(int ibin = 0; ibin <= nEBins; ++ibin){
249  if(showerEnergy > eBins[ibin] && showerEnergy < eBins[ibin+1]){
250  showerEBin = ibin;
251  if( showerEBin == 2)
252  break;
253  }
254  }// end loop over energy bins
255 
256  if( showerEBin == -1)
257  continue;
258 
259  // figure out which of the four detector regions we are in
260  int detRegion = -1;
261  if( (showerStart[0]+showerStop[0])/2.0>0&&(showerStart[1]+showerStop[1])/2.0>0 )
262  detRegion = 0;
263  if( (showerStart[0]+showerStop[0])/2.0>0&&(showerStart[1]+showerStop[1])/2.0<0 )
264  detRegion = 1;
265  if( (showerStart[0]+showerStop[0])/2.0<0&&(showerStart[1]+showerStop[1])/2.0>0 )
266  detRegion = 2;
267  if( (showerStart[0]+showerStop[0])/2.0<0&&(showerStart[1]+showerStop[1])/2.0<0 )
268  detRegion = 3;
269 
271 
272  for( int iplane =0; iplane < nplanes; ++iplane)
273  hlongdedx[detRegion][showerEBin][iplane]->Fill( showerPlaneDedx[int(iplane*fabs(showerDir[2]))] );
274 
275  for( int iplane =0; iplane < kNumTransversePlane; ++iplane)
276  htransdedx[detRegion][showerEBin][iplane]->Fill( showerTransCellDedx[iplane] );
277 
278  }// end loop over entries
279 
280  std::cout<<"the output directory is "<<outDir<<std::endl;
281  std::string outDirStr(outDir);
282 
283  SaveHists( inPdg, outDirStr, inNuMode);
284  Reset();
285  }// end of FillHists
286 
287 }// end of namespace
double showerPlaneDedx[200]
TH1D * htransdedx[kNumXYRegion][kNumEnergyBin][kNumTransversePlane]
void InitializeHists()
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
static const double kEnergyBins[]
void SaveHists(const int inPdg, const std::string outDir, const int *inNuMode)
int showerTruePdg
int showerTrueNuCCNC
static const int kNumTransversePlane
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
double showerTrueDang
double showerDir[3]
std::string outDir
void FillHists(const char *inFnames, const char *inTreename, const char *outDir, const int inPdg, const float inDangLower, const float inDangUpper, const int *inNuPdg, const int *inNuMode, const int *inNuCC, const double *eTrueFrac)
TH1D * hlongdedx[kNumXYRegion][kNumEnergyBin][kNumLongitudinalPlane]
bool fInitialized
Long64_t nentries
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
int showerTrueNuPdg
void Reset()
static const int kNumEnergyBin
TChain * fChain
int showerNPlane
double showerEnergy
void ActivateBranches()
int showerTrueNuMode
OStream cout
Definition: OStream.cxx:6
int nplanes
Definition: geom.C:145
bool fActivated
double showerTrueP4[4]
double showerStop[3]
double showerStart[3]
double showerIsFid
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
double showerTransCellDedx[20]
static const int kNumLongitudinalPlane
static const int kNumXYRegion
enum BeamMode string