1 |
pam-fi |
1.1 |
#include <CaloHough.h> |
2 |
|
|
//////////////////////////////////////////////////////////// |
3 |
|
|
/// INITIALIZATIONS & GLOBAL PARAMETERS |
4 |
|
|
//////////////////////////////////////////////////////////// |
5 |
|
|
|
6 |
|
|
CaloHough_parameters * CaloHough_parameters::_parameters = 0; |
7 |
|
|
|
8 |
|
|
/** |
9 |
|
|
* Set external parameters to default values |
10 |
|
|
*/ |
11 |
|
|
void CaloHough_parameters::SetDefault(){ |
12 |
|
|
|
13 |
|
|
cout << "CaloHough --> SET DEFAULT PARAMETERS" << endl; |
14 |
|
|
|
15 |
|
|
isimu = 0; |
16 |
|
|
|
17 |
|
|
// Q_UP = 4.;// mip |
18 |
|
|
// Q_DN = 0.; |
19 |
|
|
// M_UP = 3; |
20 |
|
|
// M_DN = 0; |
21 |
|
|
Q_UP = 10000.; |
22 |
|
|
Q_DN = 100.; |
23 |
|
|
M_UP = 3000; |
24 |
|
|
M_DN = 0; |
25 |
|
|
|
26 |
|
|
h_qtot[0] = new TH2F("houghqtot_0","Hit pattern",96,-0.5,95.5,22,-0.5,21.5); |
27 |
|
|
h_qtot[1] = new TH2F("houghqtot_1","Hit pattern",96,-0.5,95.5,22,-0.5,21.5); |
28 |
|
|
|
29 |
|
|
int n1 = 150; |
30 |
|
|
int n2 = 150; |
31 |
|
|
|
32 |
|
|
h_par[0] = new TH2F("houghpar_0","par1 vs par2",n1,-40.,40.,n2,-10.,10.); |
33 |
|
|
h_par[1] = new TH2F("houghpar_1","par1 vs par2",n1,-40.,40.,n2,-10.,10.); |
34 |
|
|
|
35 |
|
|
h_par_fine[0] = new TH2F("houghparfine_0","par1 vs par2",10*n1,-40.,40.,10*n2,-10.,10.); |
36 |
|
|
h_par_fine[1] = new TH2F("houghparfine_1","par1 vs par2",10*n1,-40.,40.,10*n2,-10.,10.); |
37 |
|
|
|
38 |
|
|
h_par1[0] = new TH1F("houghpar1_0","par1",n1,-10.,10.); |
39 |
|
|
h_par1[1] = new TH1F("houghpar1_1","par1",n1,-10.,10.); |
40 |
|
|
|
41 |
|
|
h_par2[0] = new TH1F("houghpar2_0","par2",n2,-40.,40.); |
42 |
|
|
h_par2[1] = new TH1F("houghpar2_1","par2",n2,-40.,40.); |
43 |
|
|
}; |
44 |
|
|
//////////////////////////////////////////////////////////// |
45 |
|
|
/// CLASS IMPLEMENTATION |
46 |
|
|
//////////////////////////////////////////////////////////// |
47 |
|
|
/** |
48 |
|
|
* Reset event |
49 |
|
|
* it has to be called for each event |
50 |
|
|
*/ |
51 |
|
|
void CaloHough::Reset(){ |
52 |
|
|
return; |
53 |
|
|
} |
54 |
|
|
|
55 |
|
|
/** |
56 |
|
|
* \brief Set event |
57 |
|
|
*/ |
58 |
|
|
bool CaloHough::Set(CaloLevel1 *cl1){ |
59 |
|
|
|
60 |
|
|
// bool debug = CaloHough_parameters::Get()->debug; |
61 |
|
|
float Q_UP = CaloHough_parameters::Get()->Q_UP; |
62 |
|
|
int M_UP = CaloHough_parameters::Get()->M_UP; |
63 |
|
|
float Q_DN = CaloHough_parameters::Get()->Q_DN; |
64 |
|
|
int M_DN = CaloHough_parameters::Get()->M_DN; |
65 |
|
|
|
66 |
|
|
cout << "bool CaloHough::Set(CaloLevel1*)"<<endl; |
67 |
|
|
|
68 |
|
|
if(!cl1)return false; |
69 |
|
|
|
70 |
|
|
TH2F* h_qtot[2]; |
71 |
|
|
TH2F* h_par[2]; |
72 |
|
|
TH2F* h_par_fine[2]; |
73 |
|
|
TH1F* h_par1[2]; |
74 |
|
|
TH1F* h_par2[2]; |
75 |
|
|
for(int iv=0; iv<2; iv++){ |
76 |
|
|
h_qtot[iv] = CaloHough_parameters::Get()->h_qtot[iv]; |
77 |
|
|
h_par[iv] = CaloHough_parameters::Get()->h_par[iv]; |
78 |
|
|
h_par_fine[iv] = CaloHough_parameters::Get()->h_par_fine[iv]; |
79 |
|
|
h_par1[iv] = CaloHough_parameters::Get()->h_par1[iv]; |
80 |
|
|
h_par2[iv] = CaloHough_parameters::Get()->h_par2[iv]; |
81 |
|
|
if(!h_qtot[iv] || !h_par[iv] || !h_par1[iv] || !h_par2[iv] || !h_par_fine[iv])return false; |
82 |
|
|
h_qtot[iv]->Reset(); |
83 |
|
|
h_par[iv]->Reset(); |
84 |
|
|
h_par_fine[iv]->Reset(); |
85 |
|
|
h_par1[iv]->Reset(); |
86 |
|
|
h_par2[iv]->Reset(); |
87 |
|
|
} |
88 |
|
|
|
89 |
|
|
// ----------------------------------- |
90 |
|
|
// set charge collected in each plane |
91 |
|
|
// ----------------------------------- |
92 |
|
|
// |
93 |
|
|
CaloStrip cstrip; |
94 |
|
|
int isimu = CaloHough_parameters::Get()->isimu; |
95 |
|
|
if(isimu==0)cstrip = CaloStrip(cl1,false); |
96 |
|
|
if(isimu==1)cstrip = CaloStrip(cl1,true); |
97 |
|
|
// |
98 |
|
|
vector<int> ids[22][2]; //id strip 0-95 |
99 |
|
|
vector<float> xy[22][2]; //cluster coordinate |
100 |
|
|
int islast = -1; |
101 |
|
|
int iplast = -1; |
102 |
|
|
int mult = 0; |
103 |
|
|
float qsum = 0.; |
104 |
|
|
float xyqsum = 0.; |
105 |
|
|
float issum = 0.; |
106 |
|
|
int ncluster[] = {0,0}; |
107 |
|
|
for(int ih=0;ih<cl1->istrip;ih++){ |
108 |
|
|
// |
109 |
|
|
int iv=-1; |
110 |
|
|
int ip=-1; |
111 |
|
|
int is=-1; |
112 |
|
|
cl1->DecodeEstrip(ih,iv,ip,is); |
113 |
|
|
float xyz[3]; |
114 |
|
|
cstrip.Set(iv,ip,is); |
115 |
|
|
// cout << " * "<<setw(3)<<iv<<setw(3)<<ip<<setw(3)<<is<<endl; |
116 |
|
|
xyz[0] = cstrip.GetX(); |
117 |
|
|
xyz[1] = cstrip.GetY(); |
118 |
|
|
xyz[2] = cstrip.GetZ(); |
119 |
|
|
// |
120 |
|
|
// ...................................... metodo 2 |
121 |
|
|
if( is != islast+1 || ip != iplast ){ // new cluster |
122 |
|
|
if( |
123 |
|
|
qsum < Q_UP && |
124 |
|
|
qsum > Q_DN && |
125 |
|
|
mult < M_UP && |
126 |
|
|
mult > M_DN && |
127 |
|
|
true){ // store the previous cluster |
128 |
|
|
xy[ip][iv].push_back( xyqsum/qsum ); |
129 |
|
|
ids[ip][iv].push_back( (int)(0.5+issum/qsum) ); |
130 |
|
|
h_qtot[iv]->Fill( issum/qsum, (float)(21-ip), qsum); |
131 |
|
|
ncluster[iv]++; |
132 |
|
|
cout << " >> "<<qsum<<endl; |
133 |
|
|
}// ... and reset counters |
134 |
|
|
mult = 0; |
135 |
|
|
qsum = 0.; |
136 |
|
|
xyqsum = 0.; |
137 |
|
|
issum = 0.; |
138 |
|
|
} |
139 |
|
|
|
140 |
|
|
mult++; |
141 |
|
|
qsum += cstrip.GetE(); |
142 |
|
|
xyqsum += xyz[iv] * cstrip.GetE(); |
143 |
|
|
issum += (float)is * cstrip.GetE(); |
144 |
|
|
islast = is; |
145 |
|
|
iplast = ip; |
146 |
|
|
|
147 |
|
|
// ...................................... metodo 1 |
148 |
|
|
// if( cstrip.GetE() < Q_UP ){ |
149 |
|
|
// h_qtot[iv]->Fill( (float)is, (float)(21-ip), cstrip.GetE()); |
150 |
|
|
// ids[ip][iv].push_back(is); |
151 |
|
|
// } |
152 |
|
|
|
153 |
|
|
} |
154 |
|
|
cout << "n.clusters X "<<ncluster[0]<<endl; |
155 |
|
|
cout << "n.clusters Y "<<ncluster[0]<<endl; |
156 |
|
|
// |
157 |
|
|
// cstrip.Set(0,0,0); |
158 |
|
|
// cstrip.Set(0,21,0); |
159 |
|
|
cstrip.Set(0,11,0); |
160 |
|
|
float zref = cstrip.GetZ();//reference plane |
161 |
|
|
// |
162 |
|
|
for(int iv=0; iv<2; iv++){ |
163 |
|
|
for(int ip1=0; ip1<22-1; ip1++){ |
164 |
|
|
// for(int ip2=ip1+1; ip2<22; ip2++){ |
165 |
|
|
for(int ip2=ip1+1; ip2<min(22,ip1+3); ip2++){ |
166 |
|
|
// cout << "ip1 "<<ip1<<" "<<ids[ip1][iv].size()<<endl; |
167 |
|
|
// cout << "ip2 "<<ip2<<" "<<ids[ip2][iv].size()<<endl; |
168 |
mocchiut |
1.2 |
for(int is1=0; is1<(int)ids[ip1][iv].size(); is1++){ |
169 |
|
|
for(int is2=0; is2<(int)ids[ip2][iv].size(); is2++){ |
170 |
pam-fi |
1.1 |
|
171 |
|
|
// ...................................... metodo 1 |
172 |
|
|
// int ids1 = ids[ip1][iv][is1]; |
173 |
|
|
// int ids2 = ids[ip2][iv][is2]; |
174 |
|
|
// float xyz1[3]; |
175 |
|
|
// float xyz2[3]; |
176 |
|
|
// cstrip.Set(iv,ip1,ids1); |
177 |
|
|
// xyz1[0] = cstrip.GetX(); |
178 |
|
|
// xyz1[1] = cstrip.GetY(); |
179 |
|
|
// xyz1[2] = cstrip.GetZ(); |
180 |
|
|
// cstrip.Set(iv,ip2,ids2); |
181 |
|
|
// xyz2[0] = cstrip.GetX(); |
182 |
|
|
// xyz2[1] = cstrip.GetY(); |
183 |
|
|
// xyz2[2] = cstrip.GetZ(); |
184 |
|
|
|
185 |
|
|
// ...................................... metodo 2 |
186 |
mocchiut |
1.2 |
// int ids1 = ids[ip1][iv][is1]; |
187 |
pam-fi |
1.1 |
int ids2 = ids[ip2][iv][is2]; |
188 |
|
|
float xyz1[3]; |
189 |
|
|
float xyz2[3]; |
190 |
|
|
cstrip.Set(iv,ip1,is1); |
191 |
|
|
xyz1[iv] = xy[ip1][iv][is1]; |
192 |
|
|
xyz1[2] = cstrip.GetZ(); |
193 |
|
|
cstrip.Set(iv,ip2,ids2); |
194 |
|
|
xyz2[iv] = xy[ip2][iv][is2]; |
195 |
|
|
xyz2[2] = cstrip.GetZ(); |
196 |
|
|
// |
197 |
|
|
float par1 = 0.; |
198 |
|
|
if(fabs(xyz1[2]-xyz2[2])>0.)par1 = (xyz1[iv]-xyz2[iv])/(xyz1[2]-xyz2[2]); |
199 |
|
|
float par2 = xyz1[iv] - par1*xyz1[2]; |
200 |
|
|
par2 = par1*zref + par2; |
201 |
|
|
h_par[iv]->Fill(par2,par1); // parameter space |
202 |
|
|
h_par_fine[iv]->Fill(par2,par1); // parameter space |
203 |
|
|
h_par1[iv]->Fill(par1); |
204 |
|
|
h_par2[iv]->Fill(par2); |
205 |
|
|
// cout << setw(15)<<par1<< setw(15)<<par2<<endl; |
206 |
|
|
} |
207 |
|
|
} |
208 |
|
|
} |
209 |
|
|
} |
210 |
|
|
// |
211 |
|
|
for( int ibx=1; ibx <= h_par[iv]->GetXaxis()->GetNbins(); ibx++){ |
212 |
|
|
for( int iby=1; iby <= h_par[iv]->GetYaxis()->GetNbins(); iby++){ |
213 |
|
|
// get bin content |
214 |
|
|
float content = h_par[iv]->GetBinContent(ibx,iby); |
215 |
|
|
if( content > 6){ |
216 |
|
|
double xmi = h_par[iv]->GetXaxis()->GetBinLowEdge(ibx); |
217 |
|
|
double xma = h_par[iv]->GetXaxis()->GetBinUpEdge(ibx); |
218 |
|
|
double ymi = h_par[iv]->GetYaxis()->GetBinLowEdge(iby); |
219 |
|
|
double yma = h_par[iv]->GetYaxis()->GetBinUpEdge(iby); |
220 |
|
|
h_par_fine[iv]->GetXaxis()->SetRangeUser(xmi,xma); |
221 |
|
|
h_par_fine[iv]->GetYaxis()->SetRangeUser(ymi,yma); |
222 |
|
|
double par1 = h_par_fine[iv]->GetMean(1); |
223 |
|
|
double par2 = h_par_fine[iv]->GetMean(2); |
224 |
|
|
cout << iv << " <><><><><><><>< "<<setw(15)<<par1<<setw(15)<<par2<<endl; |
225 |
|
|
} |
226 |
|
|
} |
227 |
|
|
} |
228 |
|
|
|
229 |
|
|
// |
230 |
|
|
for(int ip=0; ip<22; ip++){ |
231 |
|
|
ids[ip][iv].clear(); |
232 |
|
|
xy[ip][iv].clear(); |
233 |
|
|
} |
234 |
|
|
}//end loop on views |
235 |
|
|
|
236 |
|
|
return true; |
237 |
|
|
} |
238 |
|
|
ClassImp(CaloHough); |
239 |
|
|
ClassImp(CaloHough_parameters); |
240 |
|
|
|