Skip to content

Commit f4d623e

Browse files
kmtsuikmtsui
authored
Merge pull request #10 from hyperk/angular_response_gain
PMT QE correction by angular response function
2 parents f5d8576 + cacdb70 commit f4d623e

10 files changed

Lines changed: 234 additions & 138 deletions

‎app/application/Makefile‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,10 @@ CXXFLAGS += -Wall -std=c++11 -g $(shell root-config --cflags)\
66
-I${MDTROOT}/cpp/include\
77
-I${WCRDROOT}/include
88

9-
LDFLAGS += $(shell root-config --ldflags) $(shell root-config --libs) -lTreePlayer\
9+
LDFLAGS += -L${MDTROOT}/cpp -lMDT\
10+
-L${WCRDROOT} -lWCRData\
1011
-L$(WCSIM_BUILD_DIR)/lib -lWCSimRoot\
11-
-L${MDTROOT}/cpp -lMDT\
12-
-L${WCRDROOT} -lWCRData
12+
$(shell root-config --ldflags) $(shell root-config --libs) -lTreePlayer
1313

1414
.PHONY: clean Execs
1515

‎app/application/PMTResponse3inchR12199_02.cc‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@ void PMTResponse3inchR12199_02::Initialize(int seed, const string &pmtname)
3232
s["SPECDFFile"] = "SPECDFFile";
3333
s["PMTDE"] = "PMTDE";
3434
s["PMTTime"] = "PMTTime";
35+
s["AngularResponse"] = "AngularResponse";
3536
if( fPMTType!="" )
3637
{
3738
map<string, string>::iterator i;
@@ -47,9 +48,11 @@ void PMTResponse3inchR12199_02::Initialize(int seed, const string &pmtname)
4748
Conf->GetValue<string>(s["SPECDFFile"], fTxtFileSPECDF);
4849
Conf->GetValue<string>(s["PMTDE"], fPMTDEFile);
4950
Conf->GetValue<string>(s["PMTTime"], fPMTTFile);
51+
Conf->GetValue<string>(s["AngularResponse"], fARFile);
5052
this->LoadCDFOfSPE(fTxtFileSPECDF);
5153
this->LoadPMTDE(fPMTDEFile);
5254
this->LoadPMTTime(fPMTTFile);
55+
this->LoadAngularResponse(fARFile);
5356
}
5457

5558
float PMTResponse3inchR12199_02::HitTimeSmearing(float Q)

‎cpp/include/HitDigitizer.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,8 @@ class HitDigitizer_mPMT : public HitDigitizer
4848
private:
4949
TH1F* hWF;
5050
float fAmplitudeSigma;
51+
float fADCMax;
52+
int fADCOverflow;
5153
float fDt;
5254
float fDv;
5355
float fWaveformOffset;

‎cpp/include/PMTResponse.h‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
#include <string>
44
#include <TGraph.h>
5+
#include <TFile.h>
56
#include "MTRandom.h"
67
#include "HitTube.h"
78

@@ -48,12 +49,16 @@ class GenericPMTResponse : public PMTResponse
4849
string fTxtFileSPECDF;
4950
void LoadPMTDE(const string &s);
5051
int fLoadDE;
51-
std::vector<double> fDE;
52+
std::vector<TGraph*> fDE;
5253
string fPMTDEFile;
5354
void LoadPMTTime(const string &s);
5455
int fLoadT;
5556
std::vector<double> fT;
5657
string fPMTTFile;
58+
void LoadAngularResponse(const string &s);
59+
int fLoadAR;
60+
std::vector<double> fAR;
61+
string fARFile;
5762

5863
//private:
5964
protected:

‎cpp/src/HitDigitizer.cc‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -162,6 +162,8 @@ HitDigitizer_mPMT::HitDigitizer_mPMT(int seed) : HitDigitizer(seed)
162162
fHitInsensitivityPeriod = 8;
163163
fAmplitudeThreshold = 20.;
164164
fAmplitudeSigma = 0.37;
165+
fADCMax = 999999;
166+
fADCOverflow = 0;
165167

166168
Configuration *Conf = Configuration::GetInstance();
167169
Conf->GetValue<float>("SamplingInterval", fDt);
@@ -178,6 +180,8 @@ HitDigitizer_mPMT::HitDigitizer_mPMT(int seed) : HitDigitizer(seed)
178180

179181
this->LoadWaveform(Conf->GetValue<string>("WaveformFile"));
180182
Conf->GetValue<float>("AmplitudeSigma", fAmplitudeSigma);
183+
Conf->GetValue<float>("ADCMax", fADCMax);
184+
Conf->GetValue<int>("ADCOverflow", fADCOverflow);
181185
}
182186

183187
HitDigitizer_mPMT::~HitDigitizer_mPMT()
@@ -369,7 +373,13 @@ TH1F HitDigitizer_mPMT::BuildWavetrain(const vector<TrueHit*> PEs, double wavefo
369373

370374
for (int i=1;i<=hWT.GetNbinsX();i++)
371375
{
372-
double val = (int)(hWT.GetBinContent(i)/fDv);
376+
double val = hWT.GetBinContent(i);
377+
if (val>fADCMax)
378+
{
379+
val = fADCMax;
380+
if (fADCOverflow>0) val = -fADCMax;
381+
}
382+
val = (int)(val/fDv);
373383
val *= fDv;
374384
hWT.SetBinContent(i,val);
375385
hWT.SetBinError(i,fDv);

‎cpp/src/PMTResponse.cc‎

Lines changed: 93 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ GenericPMTResponse::GenericPMTResponse(int seed, const string &pmtname)
2020
fSclFacTTS = 1.0;
2121
fLoadDE = 0;
2222
fLoadT = 0;
23+
fLoadAR = 0;
2324
this->Initialize(seed, pmtname);
2425
}
2526

@@ -28,6 +29,7 @@ GenericPMTResponse::GenericPMTResponse()
2829
fSclFacTTS = 1.0;
2930
fLoadDE = 0;
3031
fLoadT = 0;
32+
fLoadAR = 0;
3133
}
3234

3335
GenericPMTResponse::~GenericPMTResponse()
@@ -50,6 +52,7 @@ void GenericPMTResponse::Initialize(int seed, const string &pmtname)
5052
s["SPECDFFile"] = "SPECDFFile";
5153
s["PMTDE"] = "PMTDE";
5254
s["PMTTime"] = "PMTTime";
55+
s["AngularResponse"] = "AngularResponse";
5356

5457
if( fPMTType!="" )
5558
{
@@ -67,36 +70,59 @@ void GenericPMTResponse::Initialize(int seed, const string &pmtname)
6770
Conf->GetValue<string>(s["SPECDFFile"], fTxtFileSPECDF);
6871
Conf->GetValue<string>(s["PMTDE"], fPMTDEFile);
6972
Conf->GetValue<string>(s["PMTTime"], fPMTTFile);
73+
Conf->GetValue<string>(s["AngularResponse"], fARFile);
7074

7175
this->LoadCDFOfSPE(fTxtFileSPECDF);
7276
this->LoadPMTDE(fPMTDEFile);
7377
this->LoadPMTTime(fPMTTFile);
78+
this->LoadAngularResponse(fARFile);
7479
}
7580

7681

7782
double GenericPMTResponse::GetRawSPE(const TrueHit* th, const HitTube* ht)
7883
{
84+
double ar = 1;
85+
if (fLoadAR>0 && ht)
86+
{
87+
int tubeID = ht->GetTubeID();
88+
if (tubeID>=fLoadAR)
89+
{
90+
cout<<" GenericPMTResponse::GetRawSPE" <<endl;
91+
cout<<" - tubeID = " << tubeID << " >= fLoadAR = " << fLoadAR << endl;
92+
cout<<" -> EXIT" <<endl;
93+
exit(-1);
94+
}
95+
if (th->GetParentId()>=0)
96+
{
97+
double costh = -ht->GetOrientation(0)*th->GetDirection(0)-ht->GetOrientation(1)*th->GetDirection(1)-ht->GetOrientation(2)*th->GetDirection(2);
98+
ar = 1 + (costh-0.5)*2*fAR[tubeID];
99+
}
100+
}
79101
int i;
80102
double random1=fRand->Rndm();
81103
for(i = 0; i < 501; i++){
82104
if( random1<=*(fqpe0+i) ){ break; }
83105
}
84-
return (double(i-50) + fRand->Rndm())/22.83;
106+
return (double(i-50) + fRand->Rndm())/22.83*ar;
85107
}
86108

87109
bool GenericPMTResponse::ApplyDE(const TrueHit* th, const HitTube *ht)
88110
{
89111
if (fLoadDE>0 && ht)
90112
{
91-
int tubeID = ht->GetTubeID();
92-
if (tubeID>=fLoadDE)
113+
int mPMTID = ht->GetmPMTID();
114+
if (mPMTID>=fLoadDE)
93115
{
94116
cout<<" GenericPMTResponse::ApplyDE" <<endl;
95-
cout<<" - tubeID = " << tubeID << " >= fLoadDE = " << fLoadDE << endl;
117+
cout<<" - mPMTID = " << mPMTID << " >= fLoadDE = " << fLoadDE << endl;
96118
cout<<" -> EXIT" <<endl;
97119
exit(-1);
98120
}
99-
return fRand->Rndm() < fDE[tubeID];
121+
122+
double costh = 0;
123+
for (int i=0;i<3;i++) costh -= th->GetDirection(i)*ht->GetOrientation(i);
124+
125+
return fRand->Rndm() < fDE[mPMTID]->Eval(costh);
100126
}
101127

102128
return true;
@@ -165,29 +191,40 @@ void GenericPMTResponse::LoadPMTDE(const string &filename)
165191
{
166192
fLoadDE = 0;
167193
fDE.clear();
168-
ifstream ifs(filename.c_str());
169-
if (!ifs)
194+
195+
if (filename.size()==0)
170196
{
171197
cout<<" GenericPMTResponse::LoadPMTDE" <<endl;
172198
cout<<" - No PMT QE file: " << filename <<endl;
173199
cout<<" - Do not apply individual PMT DE " << endl;
200+
return;
174201
}
175-
string aLine;
176-
while( std::getline(ifs, aLine) )
202+
TFile f(filename.c_str());
203+
if (!f.IsOpen())
177204
{
178-
if( aLine[0] == '#' ){ continue; }
179-
stringstream ssline(aLine);
180-
string item;
181-
while (getline(ssline, item, ssline.widen(' ')))
182-
{
183-
fDE.push_back( atof(item.c_str()) );
184-
}
205+
cout<<" GenericPMTResponse::LoadPMTDE" <<endl;
206+
cout<<" - Cannot open PMT QE file: " << filename <<endl;
207+
cout<<" -> EXIT" <<endl;
208+
exit(-1);
185209
}
186-
ifs.close();
187-
188-
fLoadDE = fDE.size();
210+
TGraph* graph;
211+
while ( (graph=(TGraph*)f.Get(Form("mPMT%i",fLoadDE))) )
212+
{
213+
graph->SetBit(TGraph::kIsSortedX);
214+
fDE.push_back(graph);
215+
fLoadDE++;
216+
}
217+
f.Close();
189218

190-
if (fLoadDE>0)
219+
if (fLoadDE==0)
220+
{
221+
cout<<" GenericPMTResponse::LoadPMTDE" <<endl;
222+
cout<<" - Load PMT QE file: " << filename <<endl;
223+
cout<<" - No graph loaded!!! " << endl;
224+
cout<<" -> EXIT" <<endl;
225+
exit(-1);
226+
}
227+
else
191228
{
192229
cout<<" GenericPMTResponse::LoadPMTDE" <<endl;
193230
cout<<" - Load PMT QE file: " << filename <<endl;
@@ -228,6 +265,39 @@ void GenericPMTResponse::LoadPMTTime(const string &filename)
228265
}
229266
}
230267

268+
void GenericPMTResponse::LoadAngularResponse(const string &filename)
269+
{
270+
fLoadAR = 0;
271+
fAR.clear();
272+
ifstream ifs(filename.c_str());
273+
if (!ifs)
274+
{
275+
cout<<" GenericPMTResponse::LoadAngularResponse" <<endl;
276+
cout<<" - No Angular Response file: " << filename <<endl;
277+
cout<<" - Do not apply Angular Response weights " << endl;
278+
}
279+
string aLine;
280+
while( std::getline(ifs, aLine) )
281+
{
282+
if( aLine[0] == '#' ){ continue; }
283+
stringstream ssline(aLine);
284+
string item;
285+
while (getline(ssline, item, ssline.widen(' ')))
286+
{
287+
fAR.push_back( atof(item.c_str()) );
288+
}
289+
}
290+
ifs.close();
291+
292+
fLoadAR = fAR.size();
293+
if (fLoadAR>0)
294+
{
295+
cout<<" GenericPMTResponse::LoadAngularResponse" <<endl;
296+
cout<<" - Load Anuglar Response file: " << filename <<endl;
297+
cout<<" - # Entries = " << fLoadAR << endl;
298+
}
299+
}
300+
231301
///////////////////////////////////////////////////////////////////
232302
///////////////////////////////////////////////////////////////////
233303
const double ResponseBoxandLine20inchHQE::ksig_param[4] = {0.6314, 0.06260, 0.5711,23.96};
@@ -390,6 +460,7 @@ void Response3inchR14374_WCTE::Initialize(int seed, const string &pmtname)
390460
s["SPECDFFile"] = "SPECDFFile";
391461
s["PMTDE"] = "PMTDE";
392462
s["PMTTime"] = "PMTTime";
463+
s["AngularResponse"] = "AngularResponse";
393464
if( fPMTType!="" )
394465
{
395466
map<string, string>::iterator i;
@@ -403,9 +474,11 @@ void Response3inchR14374_WCTE::Initialize(int seed, const string &pmtname)
403474
Conf->GetValue<string>(s["SPECDFFile"], fTxtFileSPECDF);
404475
Conf->GetValue<string>(s["PMTDE"], fPMTDEFile);
405476
Conf->GetValue<string>(s["PMTTime"], fPMTTFile);
477+
Conf->GetValue<string>(s["AngularResponse"], fARFile);
406478
this->LoadCDFOfSPE(fTxtFileSPECDF);
407479
this->LoadPMTDE(fPMTDEFile);
408480
this->LoadPMTTime(fPMTTFile);
481+
this->LoadAngularResponse(fARFile);
409482
}
410483

411484
float Response3inchR14374_WCTE::HitTimeSmearing(float Q)

0 commit comments

Comments
 (0)