Skip to content

Commit feab9cb

Browse files
Sergey MartynenkoSergey Martynenko
Sergey Martynenko
authored and
Sergey Martynenko
committed
Change code to work with new MC with crosstalk
add new attinuation, new MC file, new wway of reading true hits
1 parent 2e2b8dd commit feab9cb

File tree

3 files changed

+298
-13
lines changed

3 files changed

+298
-13
lines changed

‎CInputRead_CT.hxx

+264
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,264 @@
1+
#ifndef CInputRead_hxx_seen
2+
#define CInputRead_hxx_seen
3+
4+
#include "TROOT.h"
5+
#include "TObject.h"
6+
#include "TMath.h"
7+
#include "CHit2D.hxx"
8+
#include <iostream>
9+
#include <vector>
10+
11+
#include "TH2F.h"
12+
#include "TPad.h"
13+
#include "TH1.h"
14+
#include "TVector3.h"
15+
16+
bool SortHits3DTrue(const TVector3& lhs, const TVector3& rhs){
17+
18+
if(lhs.X()!=rhs.X()){
19+
return lhs.X()<rhs.X();
20+
}else if(lhs.X()==rhs.X()){
21+
if(lhs.Z()!=rhs.Z()){
22+
return lhs.Z()<rhs.Z();
23+
}else if(lhs.Z()==rhs.Z()){
24+
return lhs.Y()<rhs.Y();
25+
}
26+
}
27+
return false;
28+
29+
};
30+
31+
32+
33+
void CInputRead(TTree *tree,int& nevts,int& skip){
34+
35+
36+
std::vector<CHit2D> hitsXY;
37+
std::vector<CHit2D> hitsXZ;
38+
std::vector<CHit2D> hitsYZ;
39+
40+
int NHITS=4928;
41+
int XMIN=-110;
42+
int XMAX=110;
43+
int YMIN=-110;
44+
int YMAX=110;
45+
int ZMIN=-100;
46+
int ZMAX=100;
47+
double dimention=1;
48+
49+
50+
int event;
51+
float hitLocation[NHITS][3];
52+
int hitPE[NHITS][3],hitT[NHITS][3];
53+
54+
// tree->SetBranchAddress("event",&event);
55+
tree->SetBranchAddress("hitLocation",&hitLocation);
56+
tree->SetBranchAddress("hitPE",&hitPE);
57+
tree->SetBranchAddress("hitTime",&hitT);
58+
59+
TFile* hfile2 = new TFile("FileWith2DHits.root","RECREATE");
60+
TTree *newTree = new TTree("treeWith2DHits","tree with 2D hits");
61+
62+
newTree->Branch("2DHitsXY","vector<CHit2D>",&hitsXY,64000,1);
63+
newTree->Branch("2DHitsXZ","vector<CHit2D>",&hitsXZ,64000,1);
64+
newTree->Branch("2DHitsYZ","vector<CHit2D>",&hitsYZ,64000,1);
65+
66+
int nevents = tree->GetEntries();
67+
if(nevts>nevents){
68+
std::cout<<"process max possible number of events : "<<nevents<<std::endl;
69+
nevts=nevents;
70+
}
71+
72+
73+
for(int evt=skip;evt<nevts;++evt){
74+
std::cout<<"Process event "<<evt<<"/"<<nevts<<std::endl;
75+
tree->GetEntry(evt);
76+
77+
78+
TH2F* chargeXYplot = new TH2F("chargeXY","chargeXY",fabs(XMIN)+fabs(XMAX),XMIN,XMAX,fabs(YMIN)+fabs(YMAX),YMIN,YMAX);
79+
chargeXYplot->GetXaxis()->SetTitle("X,cm");
80+
chargeXYplot->GetYaxis()->SetTitle("Y,cm");
81+
TH2F* chargeXZplot = new TH2F("chargeXZ","chargeXZ",fabs(XMIN)+fabs(XMAX),XMIN,XMAX,fabs(ZMIN)+fabs(ZMAX),ZMIN,ZMAX);
82+
chargeXZplot->GetXaxis()->SetTitle("X,cm");
83+
chargeXZplot->GetYaxis()->SetTitle("Z,cm");
84+
TH2F* chargeYZplot = new TH2F("chargeYZ","chargeYZ",fabs(YMIN)+fabs(YMAX),YMIN,YMAX,fabs(ZMIN)+fabs(ZMAX),ZMIN,ZMAX);
85+
chargeYZplot->GetXaxis()->SetTitle("Y,cm");
86+
chargeYZplot->GetYaxis()->SetTitle("Z,cm");
87+
88+
TH2F* timeXYplot = new TH2F("timeXY","timeXY",fabs(XMIN)+fabs(XMAX),XMIN,XMAX,fabs(YMIN)+fabs(YMAX),YMIN,YMAX);
89+
timeXYplot->GetXaxis()->SetTitle("X,cm");
90+
timeXYplot->GetYaxis()->SetTitle("Y,cm");
91+
TH2F* timeXZplot = new TH2F("timeXZ","timeXZ",fabs(XMIN)+fabs(XMAX),XMIN,XMAX,fabs(ZMIN)+fabs(ZMAX),ZMIN,ZMAX);
92+
timeXZplot->GetXaxis()->SetTitle("X,cm");
93+
timeXZplot->GetYaxis()->SetTitle("Z,cm");
94+
TH2F* timeYZplot = new TH2F("timeYZ","timeYZ",fabs(YMIN)+fabs(YMAX),YMIN,YMAX,fabs(ZMIN)+fabs(ZMAX),ZMIN,ZMAX);
95+
timeYZplot->GetXaxis()->SetTitle("Y,cm");
96+
timeYZplot->GetYaxis()->SetTitle("Z,cm");
97+
98+
std::vector<std::pair<int,std::shared_ptr<std::vector<int>>>> trueInfoXY;
99+
std::vector<std::pair<int,std::shared_ptr<std::vector<int>>>> trueInfoXZ;
100+
std::vector<std::pair<int,std::shared_ptr<std::vector<int>>>> trueInfoYZ;
101+
102+
for(std::size_t i=0;i<(fabs(XMIN)+fabs(XMAX))*(fabs(YMIN)+fabs(YMAX))+10;++i){
103+
std::shared_ptr<std::vector<int>> vect1(new std::vector<int>);
104+
trueInfoXY.push_back(std::make_pair(i,vect1));
105+
}
106+
for(std::size_t i=0;i<(fabs(XMIN)+fabs(XMAX))*(fabs(ZMIN)+fabs(ZMAX))+10;++i){
107+
108+
std::shared_ptr<std::vector<int>> vect2(new std::vector<int>);
109+
trueInfoXZ.push_back(std::make_pair(i,vect2));
110+
}
111+
for(std::size_t i=0;i<(fabs(YMIN)+fabs(YMAX))*(fabs(ZMIN)+fabs(ZMAX))+10;++i){
112+
113+
std::shared_ptr<std::vector<int>> vect3(new std::vector<int>);
114+
trueInfoYZ.push_back(std::make_pair(i,vect3));
115+
}
116+
117+
for(int i=0;i<NHITS;i++){
118+
if(hitPE[i][0]>-1){
119+
if(hitPE[i][0]>9999 || hitPE[i][1]>9999 || hitPE[i][2]>9999)continue;
120+
if(hitLocation[i][0]>9999 || hitLocation[i][1]>9999 || hitLocation[i][2]>9999 )continue;
121+
if(isnan(hitLocation[i][0]) || isnan(hitLocation[i][1]) || isnan(hitLocation[i][2]) )continue;
122+
TVector3 position(hitLocation[i][0],hitLocation[i][1],hitLocation[i][2]);
123+
124+
125+
double chargeXY = (double)hitPE[i][2];
126+
double chargeXZ = (double)hitPE[i][1];
127+
double chargeYZ = (double)hitPE[i][0];
128+
129+
130+
double timeXY = hitT[i][2];
131+
double timeXZ = hitT[i][1];
132+
double timeYZ = hitT[i][0];
133+
int ibinXY=chargeXYplot->FindBin(position.X()/dimention,position.Y()/dimention);
134+
chargeXYplot->AddBinContent(chargeXYplot->FindBin(position.X()/dimention,position.Y()/dimention),chargeXY);
135+
double pT1 = timeXYplot->GetBinContent(timeXYplot->FindBin(position.X()/dimention,position.Y()/dimention));
136+
if(pT1!=0 && pT1>timeXY)timeXYplot->SetBinContent(timeXYplot->FindBin(position.X()/dimention,position.Y()/dimention),timeXY);
137+
if(pT1==0)timeXYplot->SetBinContent(timeXYplot->FindBin(position.X()/dimention,position.Y()/dimention),timeXY);
138+
139+
for(std::size_t j=0;j<trueInfoXY.size();++j){
140+
if(ibinXY==trueInfoXY[j].first){
141+
trueInfoXY[j].second->push_back(i);
142+
break;
143+
}
144+
}
145+
146+
147+
148+
int ibinXZ=chargeXZplot->FindBin(position.X()/dimention,position.Z()/dimention);
149+
chargeXZplot->AddBinContent(chargeXZplot->FindBin(position.X()/dimention,position.Z()/dimention),chargeXZ);
150+
double pT2 = timeXZplot->GetBinContent(timeXZplot->FindBin(position.X()/dimention,position.Z()/dimention));
151+
if(pT2!=0 && pT2>timeXZ)timeXZplot->SetBinContent(timeXZplot->FindBin(position.X()/dimention,position.Z()/dimention),timeXZ);
152+
if(pT2==0)timeXZplot->SetBinContent(timeXZplot->FindBin(position.X()/dimention,position.Z()/dimention),timeXZ);
153+
154+
for(std::size_t j=0;j<trueInfoXZ.size();++j){
155+
if(ibinXZ==trueInfoXZ[j].first){
156+
trueInfoXZ[j].second->push_back(i);
157+
break;
158+
}
159+
}
160+
161+
int ibinYZ=chargeYZplot->FindBin(position.Y()/dimention,position.Z()/dimention);
162+
chargeYZplot->AddBinContent(chargeYZplot->FindBin(position.Y()/dimention,position.Z()/dimention),chargeYZ);
163+
double pT3 = timeYZplot->GetBinContent(timeYZplot->FindBin(position.Y()/dimention,position.Z()/dimention));
164+
if(pT3!=0 && pT3>timeYZ)timeYZplot->SetBinContent(timeYZplot->FindBin(position.Y()/dimention,position.Z()/dimention),timeYZ);
165+
if(pT3==0)timeYZplot->SetBinContent(timeYZplot->FindBin(position.Y()/dimention,position.Z()/dimention),timeYZ);
166+
167+
for(std::size_t j=0;j<trueInfoYZ.size();++j){
168+
if(ibinYZ==trueInfoYZ[j].first){
169+
trueInfoYZ[j].second->push_back(i);
170+
break;
171+
}
172+
}
173+
}
174+
}
175+
int nXY=0;
176+
int nXZ=0;
177+
int nYZ=0;
178+
for(int i=1; i<(fabs(XMIN)+fabs(XMAX))*(fabs(YMIN)+fabs(YMAX))+2;++i){
179+
if(chargeXYplot->GetBinContent(trueInfoXY[i-1].first)>0){
180+
CHit2D hit;
181+
int nx = chargeXYplot->GetXaxis()->GetNbins()+2;
182+
int ny = chargeXYplot->GetYaxis()->GetNbins()+2;
183+
int binx = trueInfoXY[i-1].first%nx;
184+
int biny = ((trueInfoXY[i-1].first-binx)/nx)%ny;
185+
hit.SetRow(chargeXYplot->GetXaxis()->GetBinLowEdge(binx)+0.5);
186+
hit.SetColumn(chargeXYplot->GetYaxis()->GetBinLowEdge(biny)+0.5);
187+
hit.SetCharge(chargeXYplot->GetBinContent(trueInfoXY[i-1].first));
188+
hit.SetTime(timeXYplot->GetBinContent(trueInfoXY[i-1].first));
189+
hit.SetConstituents(*trueInfoXY[i-1].second);
190+
hit.SetPlane(2);
191+
hit.SetId(nXY);
192+
hitsXY.push_back(hit);
193+
nXY++;
194+
}
195+
}
196+
for(int i=1; i<(fabs(XMIN)+fabs(XMAX))*(fabs(ZMIN)+fabs(ZMAX))+2;++i){
197+
198+
if(chargeXZplot->GetBinContent(trueInfoXZ[i-1].first)>0){
199+
CHit2D hit;
200+
int nx = chargeXZplot->GetXaxis()->GetNbins()+2;
201+
int ny = chargeXZplot->GetYaxis()->GetNbins()+2;
202+
int binx = trueInfoXZ[i-1].first%nx;
203+
int biny = ((trueInfoXZ[i-1].first-binx)/nx)%ny;
204+
hit.SetRow(chargeXZplot->GetXaxis()->GetBinLowEdge(binx)+0.5);
205+
hit.SetColumn(chargeXZplot->GetYaxis()->GetBinLowEdge(biny)+0.5);
206+
hit.SetCharge(chargeXZplot->GetBinContent(trueInfoXZ[i-1].first));
207+
hit.SetTime(timeXZplot->GetBinContent(trueInfoXZ[i-1].first));
208+
hit.SetConstituents(*trueInfoXZ[i-1].second);
209+
hit.SetPlane(1);
210+
hit.SetId(nXZ);
211+
hitsXZ.push_back(hit);
212+
nXZ++;
213+
}
214+
}
215+
for(int i=1; i<(fabs(YMIN)+fabs(YMAX))*(fabs(ZMIN)+fabs(ZMAX))+2;++i){
216+
if(chargeYZplot->GetBinContent(trueInfoYZ[i-1].first)>0){
217+
CHit2D hit;
218+
int nx = chargeYZplot->GetXaxis()->GetNbins()+2;
219+
int ny = chargeYZplot->GetYaxis()->GetNbins()+2;
220+
int binx = trueInfoYZ[i-1].first%nx;
221+
int biny = ((trueInfoYZ[i-1].first-binx)/nx)%ny;
222+
hit.SetRow(chargeYZplot->GetXaxis()->GetBinLowEdge(binx)+0.5);
223+
hit.SetColumn(chargeYZplot->GetYaxis()->GetBinLowEdge(biny)+0.5);
224+
hit.SetCharge(chargeYZplot->GetBinContent(trueInfoYZ[i-1].first));
225+
hit.SetTime(timeYZplot->GetBinContent(trueInfoYZ[i-1].first));
226+
hit.SetConstituents(*trueInfoYZ[i-1].second);
227+
hit.SetPlane(0);
228+
hit.SetId(nYZ);
229+
hitsYZ.push_back(hit);
230+
nYZ++;
231+
}
232+
}
233+
234+
235+
delete chargeXZplot;
236+
delete chargeXYplot;
237+
delete chargeYZplot;
238+
delete timeXYplot;
239+
delete timeXZplot;
240+
delete timeYZplot;
241+
242+
newTree->Fill();
243+
hfile2->Write();
244+
245+
hitsXY.clear();
246+
hitsXZ.clear();
247+
hitsYZ.clear();
248+
249+
trueInfoXY.clear();
250+
trueInfoXZ.clear();
251+
trueInfoYZ.clear();
252+
253+
254+
}
255+
256+
257+
hfile2->Close();
258+
delete hfile2;
259+
260+
261+
};
262+
263+
264+
#endif

‎CRecon.C

+8-6
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@
2121
#include "CHit3D.hxx"
2222
#include "CBond3D.hxx"
2323
#include "CCluster3D.hxx"
24-
#include "CInputRead.hxx"
24+
#include "CInputRead_CT.hxx"
2525
#include "CCreate3DHits.hxx"
2626
#include "CSharedCharge.hxx"
2727
#include "CCluster3DHits.hxx"
@@ -46,11 +46,13 @@ int CRecon(){
4646
}
4747

4848

49-
TFile* hfile = new TFile("testEvent_3DST+emptyECAL_event_222_sampleT_v2.root","READ");
50-
TTree* tree = (TTree*)hfile->Get("EDepSimTree");
49+
/*TFile* hfile = new TFile("full3DST.neutrino.eleSim.file0_qe.root","READ");
50+
TTree* tree = (TTree*)hfile->Get("EDepSimTree");*/
51+
TFile* hfile = new TFile("/Users/sergey/Desktop/DUNE/work/Reconstruction/Cesar/rec/MC_output.root","READ");
52+
TTree* tree = (TTree*)hfile->Get("AllEvents");
5153

52-
int nprocessEvents = 1;
53-
int skip=0;
54+
int nprocessEvents = 300;
55+
int skip=1;
5456
nprocessEvents+=skip;
5557
CInputRead(tree,nprocessEvents,skip);
5658

@@ -74,7 +76,7 @@ int CRecon(){
7476

7577
hfile3D->Close();
7678
delete hfile3D;
77-
79+
7880

7981
return 0;
8082
}

‎CSharedCharge.hxx

+26-7
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,8 @@ void AttachOutput(std::string output3D) {
7777
// Find the amount of light attenuated in the fiber. The attenuation values
7878
// are taken from Guang's ElecSim.C macro which is tuned for the super-fgd
7979
// prototype at CERN.
80+
81+
#ifdef GUANG
8082
double FiberAttenuation(double ell) {
8183
const double LongCompFrac_FGD = 0.816;
8284
const double LongAtt_FGD = 11926.; //*CLHEP::mm;
@@ -85,7 +87,21 @@ double FiberAttenuation(double ell) {
8587
(1.0-LongCompFrac_FGD)*ell/ShortAtt_FGD;
8688
return std::exp(-arg);
8789
}
88-
90+
#endif
91+
#define SFGD
92+
#ifdef SFGD
93+
double FiberAttenuation(double ell) {
94+
// const double LongCompFrac_FGD = 0.816;
95+
// const double LongAtt_FGD = 11926.; //*CLHEP::mm;
96+
// const double ShortAtt_FGD = 312.; //*CLHEP::mm;
97+
98+
const double LongCompFrac_FGD = 0.77;
99+
const double LongAtt_FGD = 4634.;//*CLHEP::mm;
100+
const double ShortAtt_FGD = 332.;//*CLHEP::mm;
101+
102+
return ( LongCompFrac_FGD*exp((-ell)/LongAtt_FGD) + (1-LongCompFrac_FGD)*exp((-ell)/ShortAtt_FGD) );
103+
}
104+
#endif
89105
// Forward declare the Augmented objects.
90106
struct AugmentedCube;
91107
struct AugmentedDeposit;
@@ -344,8 +360,9 @@ void FillAugmented(const std::vector<CHit3D>& hit3D,
344360
Deposit.Cube = cube;
345361
newCube.Deposits.push_back(Deposit.Index);
346362
theFiber.Deposits.push_back(Deposit.Index);
347-
#define MPPC_POSITION (0.0)
348-
double dist = theCube.GetPosition().X() - MPPC_POSITION;
363+
//-120.5
364+
#define MPPC_POSITION_X (-102)
365+
double dist = (theCube.GetPosition().X() - MPPC_POSITION_X)*10;
349366
Deposit.Attenuation = FiberAttenuation(dist);
350367
Deposit.SetMeasurement(theFiber.Measurement);
351368
gAugmentedDeposits.push_back(Deposit);
@@ -360,8 +377,9 @@ void FillAugmented(const std::vector<CHit3D>& hit3D,
360377
Deposit.Cube = cube;
361378
newCube.Deposits.push_back(Deposit.Index);
362379
theFiber.Deposits.push_back(Deposit.Index);
363-
#define MPPC_POSITION (0.0)
364-
double dist = theCube.GetPosition().Y() - MPPC_POSITION;
380+
//-120.5
381+
#define MPPC_POSITION_Y (-94)
382+
double dist = (theCube.GetPosition().Y() - MPPC_POSITION_Y)*10;
365383
Deposit.Attenuation = FiberAttenuation(dist);
366384
Deposit.SetMeasurement(theFiber.Measurement);
367385
gAugmentedDeposits.push_back(Deposit);
@@ -376,8 +394,9 @@ void FillAugmented(const std::vector<CHit3D>& hit3D,
376394
Deposit.Cube = cube;
377395
newCube.Deposits.push_back(Deposit.Index);
378396
theFiber.Deposits.push_back(Deposit.Index);
379-
#define MPPC_POSITION (0.0)
380-
double dist = theCube.GetPosition().Z() - MPPC_POSITION;
397+
//-100.5
398+
#define MPPC_POSITION_Z (-28)
399+
double dist = (theCube.GetPosition().Z() - MPPC_POSITION_Z)*10;
381400
Deposit.Attenuation = FiberAttenuation(dist);
382401
Deposit.SetMeasurement(theFiber.Measurement);
383402
gAugmentedDeposits.push_back(Deposit);

0 commit comments

Comments
 (0)