Project

General

Profile

Documentation of libsbsdig » GMnBeambkgdSample.C

Eric Fuchey, 09/18/2020 01:03 PM

 
1
#include "TMatrixD.h"
2
#include "TVectorD.h"
3
#include "TFile.h"
4
#include "TChain.h"
5
#include <iostream>
6
#include <fstream>
7
#include <vector>
8
#include "TString.h"
9
#include "TDecompSVD.h"
10
#include "TCut.h"
11
#include "TEventList.h"
12
#include "G4SBSRunData.hh"
13
//#include "g4sbs_gmn_tree.C"
14
#include "gmn_deftree.C"
15
//#include "gmn_nope_tree.C"
16
#include "TH1D.h"
17
#include "TH2D.h"
18
#include "TProfile.h"
19
#include "TVector3.h"
20
#include "TMath.h"
21
#include "TF2.h"
22
#include "HistLoader.h"
23
#include "TChainElement.h"
24
#include "TRandom3.h"
25
#include "TGraph.h"
26
#include "TSystem.h"
27

    
28
/*
29
//Return an simplified PID using the G4 PID for the most common particles: e+-, mu+-, pi+-0, K+-, n, p/
30
int findReducedPID(int G4pid, double *M = 0)
31
{
32
  int sPID;
33
  switch(G4pid){
34
  case(22):
35
    sPID = 0;
36
    if(M)*M = 0.0;
37
    break;
38
  case(11):
39
    sPID = -1;
40
    if(M)*M = 0.000511;
41
    break;
42
  case(-11):
43
    sPID = 1;
44
    if(M)*M = 0.000511;
45
    break;
46
  case(13):
47
    sPID = -2;
48
    if(M)*M = 0.105658;
49
    break;
50
  case(-13):
51
    sPID = 2;
52
    if(M)*M = 0.105658;
53
    break;
54
  case(211):
55
    sPID = 3;
56
    if(M)*M = 0.139570;
57
    break;
58
  case(-211):
59
    sPID = -3;
60
    if(M)*M = 0.139570;
61
    break;
62
  case(321):
63
    sPID = 4;
64
    if(M)*M = 0.493677;
65
    break;
66
  case(-321):
67
    sPID = -4;
68
    if(M)*M = 0.493677;
69
    break;
70
  case(2112):
71
    sPID = 5;
72
    if(M)*M = 0.939565;
73
    break;
74
  case(111):
75
    sPID = -5;
76
    if(M)*M = 0.134977;
77
    break;
78
  case(2212):
79
    sPID = 6;
80
    if(M)*M = 0.938272;
81
    break;
82
  default:
83
    sPID = -10;
84
    if(M)*M = 10000.0;
85
    break;
86
  }
87
  return(sPID);
88
}
89

    
90
int retrieveG4PID(int sPID){
91
  int G4PID[12] = {111, -321, -211, 13, 11, 22, -11, -13, 211, 321, 2112, 2212};
92
  if(sPID>-10){
93
    return(G4PID[sPID+5]);
94
  }else{
95
    return(420000001);
96
  }
97
}
98
*/
99
void BeamBackground(const char *inputfilename, 
100
		    const char *root_outfile)
101
		    
102
{
103
  cout << "reading input files" << endl;
104
  
105
  ifstream inputfile(inputfilename);
106
  //TFile *fout = new TFile( outputfilename, "RECREATE" );
107
  TChain *C = new TChain("T");
108

    
109
  set<TString> files;
110
  
111
  TString currentline;
112
  while( currentline.ReadLine(inputfile) && !currentline.BeginsWith("endlist") ){
113
    if( !currentline.BeginsWith("#") ){
114
      C->Add(currentline.Data());
115

    
116
      files.insert( currentline );
117
    }
118
  }
119
  
120
  //const double Thr_Hodo = 8.0e-3;// threshold to ensure the slat with max energy deposit is recorded
121
  
122
  //string part_str[12] = {"#pi^{0}", "K^{-}", "#pi^{-}", "#mu^{-}", "e^{-}", "#gamma",
123
  //"e^{+}", "#mu^{+}", "#pi^{+}", "K^{+}", "n", "p"};
124
  
125
  //log binning for momentum: 
126
  double pbins_log10[102];
127
  for(int k = 0; k<102; k++){
128
    pbins_log10[k] = pow(10.0, (double)k*0.05-4.0);
129
  }
130
  
131
  double edepbins_log10[121];
132
  for(int k = 0; k<121; k++){
133
    edepbins_log10[k] = pow(10.0, (double)k*0.05-6.0);
134
  }
135

    
136
  double pebins_log10[81];
137
  for(int k = 0; k<81; k++){
138
    pebins_log10[k] = pow(10.0, (double)k*0.05);
139
  }
140
  
141
  TFile *fout = new TFile( root_outfile, "RECREATE" );
142
  
143
  TH1D *h1_Ntries = new TH1D("h1_Ntries","number of total generated events",1, 0, 1);
144
  
145
  //GEMs
146
  TH1D* h1_BBGEM_nhits_[5];
147
  
148
  TH2D* h1_BBGEM_yVsx_[5];
149
  TH2D* h1_BBGEM_dyVsdx_[5];
150
  TH1D* h1_BBGEM_Edep_[5];
151
  TH1D* h1_BBGEM_Edep_log_[5];
152
  
153
  for(int m = 0; m<5; m++){
154
    h1_BBGEM_nhits_[m] = new TH1D(Form("h1_BBGEM_nhits_%d",m), 
155
				  Form("BB GEM plane %d energy deposit", m+1),
156
				  500, 0.0, 500);
157
    h1_BBGEM_yVsx_[m] = new TH2D(Form("h1_BBGEM_yVsx_%d",m), 
158
  				 Form("BB GEM plane %d y vs x coordinate", m+1),
159
  				 202, -1.01, 1.01, 62, -0.31, 0.31);
160
    h1_BBGEM_dyVsdx_[m] = new TH2D(Form("h1_BBGEM_dyVsdx_%d",m), 
161
				   Form("BB GEM plane %d dy vs dx coordinate", m+1),
162
				   100, -0.1, 0.1, 100, -0.1, 0.1);
163
    h1_BBGEM_Edep_[m] = new TH1D(Form("h1_BBGEM_Edep_%d",m), 
164
  				 Form("BB GEM plane %d energy deposit;MeV", m+1),
165
  				 100, 0.0, 1.e-1);
166
    h1_BBGEM_Edep_log_[m] = new TH1D(Form("h1_BBGEM_Edep_log_%d",m), 
167
				     Form("BB GEM plane %d energy deposit;MeV", m+1),
168
				     110, edepbins_log10);
169
  }
170

    
171
  //HCal
172
  TH2D *h1_HCal_nhitsVsChan = new TH2D("h1_HCal_nhitsVsChan", "", 288, 0, 288, 100, 0.0, 100);
173
  
174
  TH2D *h1_HCal_EdepHitVsChan = new TH2D("h1_HCal_EdepHitVsChan", ";GeV", 288, 0, 288, 1000, 0.0, 2.0);
175
  TH2D *h1_HCal_EdepHitVsChan_log = new TH2D("h1_HCal_EdepHitVsChan_log", ";GeV", 288, 0, 288, 110, edepbins_log10);
176
  TH2D *h1_HCal_zHitVsChan = new TH2D("h1_HCal_zHitVsChan", ";m", 288, 0, 288, 1000, 2.2, 3.2);
177
  
178
  TH2D *h1_HCal_EdepTotVsChan = new TH2D("h1_HCal_EdepTotVsChan;GeV", ";GeV", 288, 0, 288, 1000, 0.0, 2.0);
179
  TH2D *h1_HCal_EdepTotVsChan_log = new TH2D("h1_HCal_EdepTotVsChan_log", ";GeV", 288, 0, 288, 110, edepbins_log10);
180
  
181
  //Hodoscope
182
  TH2D *h1_BBHodo_nhitsVsSlat = new TH2D("h1_BBHodo_nhitsVsSlat", "", 90, 0, 90, 100, 0, 100);
183
  
184
  TH2D *h1_BBHodo_xhitVsSlat = new TH2D("h1_BBHodo_xhitVsSlat", "", 90, 0, 90, 60, -0.3, 0.3);
185
  
186
  TH2D *h1_BBHodo_EdepHitVsSlat = new TH2D("h1_BBHodo_EdepHitVsSlat", ";GeV", 90, 0, 90, 250, 0.0, 0.5);
187
  TH2D *h1_BBHodo_EdepHitVsSlat_log = new TH2D("h1_BBHodo_EdepHitVsSlat_log", ";GeV", 90, 0, 90, 110, edepbins_log10);
188
  
189
  TH2D *h1_BBHodo_EdepTotVsSlat = new TH2D("h1_BBHodo_EdepTotVsSlat", ";GeV", 90, 0, 90, 250, 0.0, 0.5);
190
  TH2D *h1_BBHodo_EdepTotVsSlat_log = new TH2D("h1_BBHodo_EdepTotVsSlat_log", ";GeV", 90, 0, 90, 110, edepbins_log10);
191

    
192
  //PS
193
  TH2D *h1_BBPS_nhitsVsChan = new TH2D("h1_BBPS_nhitsVsChan", "", 52, 0, 52, 150, 0, 150);
194
  
195
  TH2D *h1_BBPS_EdepHitVsChan = new TH2D("h1_BBPS_EdepHitVsChan", ";GeV", 52, 0, 52, 250, 0.0, 0.5);
196
  TH2D *h1_BBPS_EdepHitVsChan_log = new TH2D("h1_BBPS_EdepHitVsChan_log", ";GeV", 52, 0, 52, 110, edepbins_log10);
197
  
198
  //SH
199
  TH2D *h1_BBSH_nhitsVsChan = new TH2D("h1_BBSH_nhitsVsChan", "", 189, 0, 189, 100, 0, 100);
200
  
201
  TH2D *h1_BBSH_EdepHitVsChan = new TH2D("h1_BBSH_EdepHitVsChan", ";GeV", 189, 0, 189, 250, 0.0, 0.5);
202
  TH2D *h1_BBSH_EdepHitVsChan_log = new TH2D("h1_BBSH_EdepHitVsChan_log", ";GeV", 189, 0, 189, 110, edepbins_log10);
203
  
204
  //GRINCH
205
  TH2D *h1_GRINCH_nhitsVsChan = new TH2D("h1_GRINCH_nhitsVsChan", "", 510, 0, 510, 20, 0, 20);
206
   
207
  TH2D *h1_GRINCH_NpeVsChan = new TH2D("h1_GRINCH_NpeVsChan", "", 510, 0, 510, 100, 0, 100);
208
  
209
  int FileCounter = 0;
210

    
211
  int nhits_GEM[5];
212
  int nhits_HCal[288];
213
  int nhits_Hodo[90];
214
  int nhits_PS[52];
215
  int nhits_SH[189];
216
  int nhits_GRINCH[510];
217
  
218
  double HCal_Edep_blocks_100ns[288];
219
  double BBHodo_Edep_slats_100ns[90];
220
  
221
  memset(nhits_GEM, 0, 5*sizeof(int));
222
  memset(nhits_HCal, 0, 288*sizeof(int));
223
  memset(nhits_Hodo, 0, 90*sizeof(int));
224
  memset(nhits_PS, 0, 52*sizeof(int));
225
  memset(nhits_SH, 0, 189*sizeof(int));
226
  memset(nhits_GRINCH, 0, 510*sizeof(int));
227

    
228
  memset(HCal_Edep_blocks_100ns, 0, 288*sizeof(double));
229
  memset(BBHodo_Edep_slats_100ns, 0, 90*sizeof(double));
230

    
231
  Long64_t N1_tot = 0;
232
  int nplane;
233
  int chan;
234

    
235
  double theta_sbs, d_hcal;
236
  double x_ref, z_ref;
237
  double z_hit;
238
  
239
  TObjArray *fileElements=C->GetListOfFiles();
240
  TIter next(fileElements);
241
  TChainElement *chEl=0;
242
  while (( chEl=(TChainElement*)next() )) {
243
    //TFile *f = new TFile(chEl->GetTitle());
244
    TFile f(chEl->GetTitle());
245
     
246
    G4SBSRunData *RD = (G4SBSRunData*)f.Get("run_data");
247
    //N1_tot+= (double)RD->fNtries;
248
    theta_sbs = RD->fSBStheta;
249
    d_hcal = RD->fHCALdist;
250
    x_ref = d_hcal*sin(theta_sbs);
251
    z_ref = d_hcal*cos(theta_sbs);
252

    
253
    TChain *C1 = (TChain*)f.Get("T");
254
    
255
    gmn_deftree *T1 = new gmn_deftree(C1);
256
    
257
    FileCounter++;
258
  
259
    Long64_t MaxEvt = C1->GetEntries();
260
    cout << chEl->GetTitle() << ": " << MaxEvt << " entries" << endl;
261

    
262
    for(Long64_t nevent = 0; nevent<MaxEvt; nevent++){
263
      //while( T1->GetEntry(nevent++) && nevent < 10000 ){
264
      if( nevent%1000 == 0){// && nevent!=0){
265
	cout << nevent << endl;
266
      }
267
      
268
      T1->GetEntry(nevent);
269
      
270
      if(T1->Harm_HCalScint_hit_nhits){
271
	for(int i = 0; i<T1->Harm_HCalScint_hit_nhits; i++){
272
	  nhits_HCal[T1->Harm_HCalScint_hit_cell->at(i)]++;
273
	  HCal_Edep_blocks_100ns[T1->Harm_HCalScint_hit_cell->at(i)]+= T1->Harm_HCalScint_hit_sumedep->at(i);
274
	  
275
	  h1_HCal_EdepHitVsChan->Fill(T1->Harm_HCalScint_hit_cell->at(i), T1->Harm_HCalScint_hit_sumedep->at(i));
276
	  h1_HCal_EdepHitVsChan_log->Fill(T1->Harm_HCalScint_hit_cell->at(i), T1->Harm_HCalScint_hit_sumedep->at(i));
277
	  
278
	  z_hit = -(T1->Harm_HCalScint_hit_xhitg->at(i)-x_ref)*sin(theta_sbs)+(T1->Harm_HCalScint_hit_zhitg->at(i)-z_ref)*cos(theta_sbs);
279
	  
280
	  h1_HCal_zHitVsChan->Fill(T1->Harm_HCalScint_hit_cell->at(i), z_hit);
281
	  // Edep_HCal+= T1->Harm_HCalScint_hit_sumedep->at(i);
282
	  // h1_HCal_blocks_Edep_rates->Fill(T1->Harm_HCalScint_hit_cell->at(i), T1->Harm_HCalScint_hit_sumedep->at(i));  
283
	  // //in MeV
284
	  // Edep_tot_50ns+= T1->Harm_HCalScint_hit_sumedep->at(i)*1.e3;
285
	  // Edep_tot_blocks_50ns[T1->Harm_HCalScint_hit_cell->at(i)]+= T1->Harm_HCalScint_hit_sumedep->at(i)*1.e3;
286
	  
287
	  // h1_HCal_Edep_tot_XC->Fill(T1->Harm_HCalScint_hit_cell->at(i), T1->Harm_HCalScint_hit_sumedep->at(i));
288
	}
289
      }
290
      
291
      //cout << "frourourtre" << endl;	
292
      if(T1->Earm_BBHodoScint_hit_nhits){
293
	for(int i = 0; i<T1->Earm_BBHodoScint_hit_nhits; i++){
294
	  nhits_Hodo[T1->Earm_BBHodoScint_hit_cell->at(i)]++;
295
	  BBHodo_Edep_slats_100ns[T1->Earm_BBHodoScint_hit_cell->at(i)]+= T1->Earm_BBHodoScint_hit_sumedep->at(i);
296
	  
297
	  h1_BBHodo_xhitVsSlat->Fill(T1->Earm_BBHodoScint_hit_cell->at(i), T1->Earm_BBHodoScint_hit_xhit->at(i));
298
	  h1_BBHodo_EdepHitVsSlat->Fill(T1->Earm_BBHodoScint_hit_cell->at(i), T1->Earm_BBHodoScint_hit_sumedep->at(i));
299
	  h1_BBHodo_EdepHitVsSlat_log->Fill(T1->Earm_BBHodoScint_hit_cell->at(i), T1->Earm_BBHodoScint_hit_sumedep->at(i));
300
	  // h1_BBHodo_slat_rates_nothr->Fill(T1->Earm_BBHodoScint_hit_cell->at(i));
301
	  // if(T1->Earm_BBHodoScint_hit_sumedep->at(i)>=Thr_Hodo){
302
	  //   h1_BBHodo_slat_rates->Fill(T1->Earm_BBHodoScint_hit_cell->at(i));
303
	  // }
304
	  // h1_BBHodo_slat_Edep_rates->Fill(T1->Earm_BBHodoScint_hit_cell->at(i), T1->Earm_BBHodoScint_hit_sumedep->at(i));
305
	  
306
	  // Nph_hodo_1 = Npe_E_GeV_bbhs*T1->Earm_BBHodoScint_hit_sumedep->at(i)*LCE_0*exp(-(0.3+(T1->Earm_BBHodoScint_hit_xhit->at(i)-T1->Earm_BBHodoScint_hit_xcellg->at(i))/cos(theta_BB))/LCE_lambda);
307
	  // Nph_hodo_2 = Npe_E_GeV_bbhs*T1->Earm_BBHodoScint_hit_sumedep->at(i)*LCE_0*exp(-(0.3-(T1->Earm_BBHodoScint_hit_xhit->at(i)-T1->Earm_BBHodoScint_hit_xcellg->at(i))/cos(theta_BB))/LCE_lambda);
308
	  
309
	  // Npe_hodo_1 = R.Poisson(Nph_hodo_1*0.24);
310
	  // //0.28 is the QE, 0.1 is a facto to take into account of the variable photon path due to the solid angle (very conservative)
311
	  // Npe_hodo_2 = R.Poisson(Nph_hodo_2*0.24);
312
	  
313
	  // if(Npe_hodo_1)h1_BBHodo_slat_pe_rates->Fill(2*T1->Earm_BBHodoScint_hit_cell->at(i), Npe_hodo_1);
314
	  // if(Npe_hodo_2)h1_BBHodo_slat_pe_rates->Fill(2*T1->Earm_BBHodoScint_hit_cell->at(i)+1, Npe_hodo_2);
315
	}
316
      }
317
      
318
      //cout << "BBGEM: " << T1->Earm_BBGEM_hit_nhits << endl;
319
      if(T1->Earm_BBGEM_hit_nhits){
320
	for(int i = 0; i<T1->Earm_BBGEM_hit_nhits; i++){
321
	  nplane = T1->Earm_BBGEM_hit_plane->at(i)-1;
322
	  
323
	  nhits_GEM[nplane]++;
324
	  // 
325
	  
326
	  // //cout << T1->Earm_BBGEM_hit_tx->at(i) << " " << T1->Earm_BBGEM_hit_ty->at(i) << endl;
327
	  
328
	  h1_BBGEM_yVsx_[nplane]->Fill(T1->Earm_BBGEM_hit_xin->at(i), 
329
	   			       T1->Earm_BBGEM_hit_yin->at(i));
330
	  h1_BBGEM_dyVsdx_[nplane]->Fill(T1->Earm_BBGEM_hit_xout->at(i)-T1->Earm_BBGEM_hit_xin->at(i), 
331
					 T1->Earm_BBGEM_hit_yout->at(i)-T1->Earm_BBGEM_hit_yin->at(i));
332
	  
333
	  h1_BBGEM_Edep_[nplane]->Fill(T1->Earm_BBGEM_hit_edep->at(i)*1.e3);
334
	  h1_BBGEM_Edep_log_[nplane]->Fill(T1->Earm_BBGEM_hit_edep->at(i)*1.e3);
335
	  
336
	}
337
      }
338
      
339
      if(T1->Earm_BBPSTF1_hit_nhits){
340
	for(int i = 0; i<T1->Earm_BBPSTF1_hit_nhits; i++){
341
	  nhits_PS[T1->Earm_BBPSTF1_hit_cell->at(i)]++;
342
	  
343
	  h1_BBPS_EdepHitVsChan->Fill(T1->Earm_BBPSTF1_hit_cell->at(i), T1->Earm_BBPSTF1_hit_sumedep->at(i));
344
	  h1_BBPS_EdepHitVsChan_log->Fill(T1->Earm_BBPSTF1_hit_cell->at(i), T1->Earm_BBPSTF1_hit_sumedep->at(i));
345
	  // h1_BBPS_Edep->Fill(T1->Earm_BBPSTF1_hit_row->at(i), 
346
	  // 		     T1->Earm_BBPSTF1_hit_col->at(i), 
347
	  // 		     T1->Earm_BBPSTF1_hit_sumedep->at(i));
348
	  // h1_BBPS_EdepSpectrum->Fill(T1->Earm_BBPSTF1_hit_cell->at(i), T1->Earm_BBPSTF1_hit_sumedep->at(i));
349
	}
350
      }
351
      //cout << "BBSHTF1: " << T1->Earm_BBSHTF1_hit_nhits << endl;
352
      if(T1->Earm_BBSHTF1_hit_nhits){
353
	for(int i = 0; i<T1->Earm_BBSHTF1_hit_nhits; i++){
354
	  nhits_SH[T1->Earm_BBSHTF1_hit_cell->at(i)]++;
355
	  
356
	  h1_BBSH_EdepHitVsChan->Fill(T1->Earm_BBSHTF1_hit_cell->at(i), T1->Earm_BBSHTF1_hit_sumedep->at(i));
357
	  h1_BBSH_EdepHitVsChan_log->Fill(T1->Earm_BBSHTF1_hit_cell->at(i), T1->Earm_BBSHTF1_hit_sumedep->at(i));
358
	  // h1_BBSH_Edep->Fill(T1->Earm_BBSHTF1_hit_row->at(i), 
359
	  // 		     T1->Earm_BBSHTF1_hit_col->at(i), 
360
	  // 		     T1->Earm_BBSHTF1_hit_sumedep->at(i));
361
	  // h1_BBSH_EdepSpectrum->Fill(T1->Earm_BBSHTF1_hit_cell->at(i), T1->Earm_BBSHTF1_hit_sumedep->at(i));
362
	}
363
      }
364
      
365
       if(T1->Earm_GRINCH_hit_nhits){
366
	for(int i = 0; i<T1->Earm_GRINCH_hit_nhits; i++){
367
	  chan = int(T1->Earm_GRINCH_hit_PMT->at(i)/5);
368
	  nhits_GRINCH[chan]++;
369

    
370
	  h1_GRINCH_NpeVsChan->Fill(chan, T1->Earm_GRINCH_hit_NumPhotoelectrons->at(i));
371
	  // h1_BBSH_Edep->Fill(T1->Earm_BBSHTF1_hit_row->at(i), 
372
	  // 		     T1->Earm_BBSHTF1_hit_col->at(i), 
373
	  // 		     T1->Earm_BBSHTF1_hit_sumedep->at(i));
374
	  // h1_BBSH_EdepSpectrum->Fill(T1->Earm_BBSHTF1_hit_cell->at(i), T1->Earm_BBSHTF1_hit_sumedep->at(i));
375
	}
376
      }
377
     
378
      
379
    }//end for(Long_64t...)
380
    
381
    //if(FileCounter==5){
382
    if(FileCounter%4==0){
383
      cout << "File Counter = 4 => 20 M evts beam-on-target" << endl;
384
      for(int j = 0; j<288; j++){
385
    	h1_HCal_nhitsVsChan->Fill(j, nhits_HCal[j]);
386
    	h1_HCal_EdepTotVsChan->Fill(j, HCal_Edep_blocks_100ns[j]);
387
    	h1_HCal_EdepTotVsChan_log->Fill(j, HCal_Edep_blocks_100ns[j]);
388
      }
389
      
390
      for(int j = 0; j<90; j++){
391
    	h1_BBHodo_nhitsVsSlat->Fill(j, nhits_Hodo[j]);
392
    	h1_BBHodo_EdepTotVsSlat->Fill(j, BBHodo_Edep_slats_100ns[j]);
393
    	h1_BBHodo_EdepTotVsSlat_log->Fill(j, BBHodo_Edep_slats_100ns[j]);
394
      }
395

    
396
      for(int j = 0; j<52; j++){
397
    	h1_BBPS_nhitsVsChan->Fill(j, nhits_PS[j]);
398
      }
399
      
400
      for(int j = 0; j<189; j++){
401
    	h1_BBSH_nhitsVsChan->Fill(j, nhits_SH[j]);
402
      }
403
      
404
      for(int j = 0; j<510; j++){
405
    	h1_GRINCH_nhitsVsChan->Fill(j, nhits_GRINCH[j]);
406
      }
407
      memset(nhits_HCal, 0, 288*sizeof(int));
408
      memset(nhits_Hodo, 0, 90*sizeof(int));
409
      memset(nhits_PS, 0, 52*sizeof(int));
410
      memset(nhits_SH, 0, 189*sizeof(int));
411
      memset(nhits_GRINCH, 0, 510*sizeof(int));
412
      
413
      memset(HCal_Edep_blocks_100ns, 0, 288*sizeof(double));
414
      memset(BBHodo_Edep_slats_100ns, 0, 90*sizeof(double));
415
    }
416
    if(FileCounter%15==0){
417
       cout << "File Counter = 15 => 75 M evts beam-on-target" << endl;
418
     for(int j = 0; j<5; j++){
419
	h1_BBGEM_nhits_[j]->Fill(nhits_GEM[j]);
420
      }
421
      memset(nhits_GEM, 0, 5*sizeof(int));
422
    }
423
    
424
    C1->Delete();
425
    T1->~gmn_deftree();
426
    //T1->~gmn_nope_tree();
427
    
428
  }//end while (( chEl= )) 
429
  //h1_Ntries->Fill(0.5, N1_tot);
430
  
431
  
432
  //h1_HCal_Edep_rates->Scale(I_exp/1.602e-19/N1_tot);
433

    
434
  C->Delete();
435
  fout->Write();
436
  //return;
437
}
(4-4/7)