Pi0Calibration/Common/PositionCalc.cc



////to make the cmssw official-like cluster position


void calcClusterLocationBarrel(float eTot,float energy[],int ietaHit[],int iphiHit[],int nhit,float posi[]){
  
  float param_T0_barl_ = 7.4; 
  // float param_T0_endc_ = 3.1; 
  //  float param_t0_endcPresh = 1.2; 
  float param_W0_ = 4.2; 
  float param_X0_ = 0.89; 
  
  ///no preshower..
  ///bool preshower = false;
  
  if(eTot<=0 || nhit <1){
    cout<<"fatal error calcClusterLocation ..negative energy "<<eTot<<" "<<nhit<<endl;
    exit(1);
  }
  
  
  int nmaxhit = 0; 
  float eMax = 0; 

  for(int j=0; j<nhit; j++){
    if(eMax < energy[j]){
      eMax = energy[j];
      nmaxhit = j; 
    }
  }
  

  int ieta = ietaHit[nmaxhit] + 85; 
  int iphi = iphiHit[nmaxhit]; 
  iphi = getIndphixyzEBAll(iphi);
  
  TVector3 vmaxHit(xEBAll[ieta][iphi],yEBAll[ieta][iphi],zEBAll[ieta][iphi]);
    
  float T0 = param_T0_barl_;
  
  if(fabs(vmaxHit.Eta())>1.482){
    cout<<"error..barrel..?"<<" "<< vmaxHit.Eta()<<" "<<ieta <<" "<<iphi <<" "<<etaEBAll[ieta][iphi]<<endl; 
    exit(1);
  }
  
  float depth =  param_X0_ * (T0 + log(eTot));
  
  TVector3 dvmaxHit(depth*dxEBAll[ieta][iphi],depth*dyEBAll[ieta][iphi],depth*dzEBAll[ieta][iphi]);
  
  TVector3 center_pos = vmaxHit + dvmaxHit;
  
  
  // Loop over hits and get weights
  double weight = 0;
  double total_weight = 0;
  
  double center_phi = center_pos.Phi();
  double center_theta = center_pos.Theta();
  
  double delta_theta = 0;
  double delta_phi = 0;

  float w_min = 0.;
  
  double dphi = 0;
  
  const double k_PI = acos(-1);
  
  
  for(int j=0; j<nhit; j++){
    
    double e_j = energy[j];
    
    if(e_j<0) continue; 
    
    
    weight = param_W0_ +log ( e_j/eTot); 
    if(weight<w_min) weight = w_min;  
    
    total_weight += weight;

    ieta = ietaHit[j] + 85; 
    iphi = iphiHit[j]; 
    iphi = getIndphixyzEBAll(iphi);
    TVector3 vHit(xEBAll[ieta][iphi],yEBAll[ieta][iphi],zEBAll[ieta][iphi]);
    TVector3 dvHit(depth*dxEBAll[ieta][iphi],depth*dyEBAll[ieta][iphi],depth*dzEBAll[ieta][iphi]);
    TVector3 jth_pos = vHit + dvHit;
    
    
    delta_theta += weight * (jth_pos.Theta() - center_theta);
    dphi = (jth_pos.Phi() - center_phi);
    
    // Check the 2*pi problem for delta_phi
    if (dphi > k_PI)
      dphi -= 2.*k_PI;
    if (dphi < -k_PI)
      dphi += 2.*k_PI;
    
    delta_phi += dphi*weight;    
  }
  
  delta_theta /= total_weight;
  delta_phi /= total_weight;
  
  double cluster_theta = center_theta + delta_theta;
  double cluster_phi = center_phi + delta_phi;

  // Check the 2*pi problem for cluster_phi
  if (cluster_phi > k_PI)
    cluster_phi -= 2.*k_PI;
  if (cluster_phi < -k_PI)
    cluster_phi += 2.*k_PI;
  
  double radius = sqrt(center_pos.x()*center_pos.x()
                       + center_pos.y()*center_pos.y()
                       + center_pos.z()*center_pos.z());
  
  float xpos = radius*cos(cluster_phi)*sin(cluster_theta);
  float ypos = radius*sin(cluster_phi)*sin(cluster_theta);
  float zpos = radius*cos(cluster_theta);
  
  
  posi[0] = xpos; 
  posi[1] = ypos; 
  posi[2] = zpos; 
  
  
}


void calcClusterLocationEndcap(float eTot,float energy[],int ixHit[],int iyHit[],int izHit[],int nhit,float posi[]){
  
  //float param_T0_barl_ = 7.4; 
  float param_T0_endc_ = 3.1; 
  float param_t0_endcPresh = 1.2; 
  float param_W0_ = 4.2; 
  float param_X0_ = 0.89; 
  
  ///no preshower..
  bool preshower = false;
  
  if(eTot<=0 || nhit <1){
    cout<<"fatal error calcClusterLocationEndcap ..negative energy "<<eTot<<" "<<nhit<<endl;
    exit(1);
  }
  
  
  int nmaxhit = 0; 
  float eMax = 0; 

  for(int j=0; j<nhit; j++){
    if(eMax < energy[j]){
      eMax = energy[j];
      nmaxhit = j; 
    }
  }
  

  int ix = ixHit[nmaxhit];
  int iy = iyHit[nmaxhit]; 
  int iz = izHit[nmaxhit]; 
  int indz = iz > 0; 
  
  TVector3 vmaxHit(xEEAll[indz][ix][iy],yEEAll[indz][ix][iy],zEEAll[indz][ix][iy]);
  
  
  float T0 = param_T0_endc_;
  if(preshower){
    T0 = param_t0_endcPresh;
  }  
    
  if(fabs(vmaxHit.Eta())<1.482){
    cout<<"error..encap..?"<<vmaxHit.Eta()<<endl;
    exit(1);
  }
  
  float depth =  param_X0_ * (T0 + log(eTot));
  
  TVector3 dvmaxHit(depth*dxEEAll[indz][ix][iy],depth*dyEEAll[indz][ix][iy],depth*dzEEAll[indz][ix][iy]);
  
  TVector3 center_pos = vmaxHit + dvmaxHit;
  
  
  // Loop over hits and get weights
  double weight = 0;
  double total_weight = 0;

  double center_phi = center_pos.Phi();
  double center_theta = center_pos.Theta();
  
  double delta_theta = 0;
  double delta_phi = 0;

  float w_min = 0.;
  
  double dphi = 0;
  
  const double k_PI = acos(-1);
    

  for(int j=0; j<nhit; j++){
    
    double e_j = energy[j];
    
    if(e_j<0) continue; 
    
    
    weight = param_W0_ +log ( e_j/eTot); 
    if(weight<w_min) weight = w_min;  
    
    total_weight += weight;

    //    ieta = ietaHit[j] + 85; 
    /// iphi = iphiHit[j]; 

    ix = ixHit[j];
    iy = iyHit[j];
    iz = izHit[nmaxhit]; 
    indz = iz > 0;    
    
    TVector3 vHit(xEEAll[indz][ix][iy],yEEAll[indz][ix][iy],zEEAll[indz][ix][iy]);
    TVector3 dvHit(depth*dxEEAll[indz][ix][iy],depth*dyEEAll[indz][ix][iy],depth*dzEEAll[indz][ix][iy]);
    

    TVector3 jth_pos = vHit + dvHit;
    
    
    delta_theta += weight * (jth_pos.Theta() - center_theta);
    dphi = (jth_pos.Phi() - center_phi);
    
    // Check the 2*pi problem for delta_phi
    if (dphi > k_PI)
      dphi -= 2.*k_PI;
    if (dphi < -k_PI)
      dphi += 2.*k_PI;
    
    delta_phi += dphi*weight;    
  }
  
  delta_theta /= total_weight;
  delta_phi /= total_weight;
  
  double cluster_theta = center_theta + delta_theta;
  double cluster_phi = center_phi + delta_phi;

  // Check the 2*pi problem for cluster_phi
  if (cluster_phi > k_PI)
    cluster_phi -= 2.*k_PI;
  if (cluster_phi < -k_PI)
    cluster_phi += 2.*k_PI;
  
  double radius = sqrt(center_pos.x()*center_pos.x()
                       + center_pos.y()*center_pos.y()
                       + center_pos.z()*center_pos.z());
  
  float xpos = radius*cos(cluster_phi)*sin(cluster_theta);
  float ypos = radius*sin(cluster_phi)*sin(cluster_theta);
  float zpos = radius*cos(cluster_theta);
  
  
  posi[0] = xpos; 
  posi[1] = ypos; 
  posi[2] = zpos; 
  
    
}
Revision:	1.1
Committed:	Thu Jul 26 08:48:21 2012 UTC (12 years, 9 months ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-00-01g, V00-00-01f, V00-00-01e, V00-00-01d, V00-00-01c, V00-00-01b, V00-00-01a, pi0calibrator_v01a, HEAD
Error occurred while calculating annotation data.
Log Message:	first commit
#	Content
1
2
3	////to make the cmssw official-like cluster position
4
5
6	void calcClusterLocationBarrel(float eTot,float energy[],int ietaHit[],int iphiHit[],int nhit,float posi[]){
7
8	float param_T0_barl_ = 7.4;
9	// float param_T0_endc_ = 3.1;
10	// float param_t0_endcPresh = 1.2;
11	float param_W0_ = 4.2;
12	float param_X0_ = 0.89;
13
14	///no preshower..
15	///bool preshower = false;
16
17	if(eTot<=0 \|\| nhit <1){
18	cout<<"fatal error calcClusterLocation ..negative energy "<<eTot<<" "<<nhit<<endl;
19	exit(1);
20	}
21
22
23	int nmaxhit = 0;
24	float eMax = 0;
25
26	for(int j=0; j<nhit; j++){
27	if(eMax < energy[j]){
28	eMax = energy[j];
29	nmaxhit = j;
30	}
31	}
32
33
34	int ieta = ietaHit[nmaxhit] + 85;
35	int iphi = iphiHit[nmaxhit];
36	iphi = getIndphixyzEBAll(iphi);
37
38	TVector3 vmaxHit(xEBAll[ieta][iphi],yEBAll[ieta][iphi],zEBAll[ieta][iphi]);
39
40	float T0 = param_T0_barl_;
41
42	if(fabs(vmaxHit.Eta())>1.482){
43	cout<<"error..barrel..?"<<" "<< vmaxHit.Eta()<<" "<<ieta <<" "<<iphi <<" "<<etaEBAll[ieta][iphi]<<endl;
44	exit(1);
45	}
46
47	float depth = param_X0_ * (T0 + log(eTot));
48
49	TVector3 dvmaxHit(depthdxEBAll[ieta][iphi],depthdyEBAll[ieta][iphi],depth*dzEBAll[ieta][iphi]);
50
51	TVector3 center_pos = vmaxHit + dvmaxHit;
52
53
54
55	// Loop over hits and get weights
56	double weight = 0;
57	double total_weight = 0;
58
59	double center_phi = center_pos.Phi();
60	double center_theta = center_pos.Theta();
61
62	double delta_theta = 0;
63	double delta_phi = 0;
64
65	float w_min = 0.;
66
67	double dphi = 0;
68
69	const double k_PI = acos(-1);
70
71
72	for(int j=0; j<nhit; j++){
73
74	double e_j = energy[j];
75
76	if(e_j<0) continue;
77
78
79	weight = param_W0_ +log ( e_j/eTot);
80	if(weight<w_min) weight = w_min;
81
82	total_weight += weight;
83
84	ieta = ietaHit[j] + 85;
85	iphi = iphiHit[j];
86	iphi = getIndphixyzEBAll(iphi);
87	TVector3 vHit(xEBAll[ieta][iphi],yEBAll[ieta][iphi],zEBAll[ieta][iphi]);
88	TVector3 dvHit(depthdxEBAll[ieta][iphi],depthdyEBAll[ieta][iphi],depth*dzEBAll[ieta][iphi]);
89	TVector3 jth_pos = vHit + dvHit;
90
91
92	delta_theta += weight * (jth_pos.Theta() - center_theta);
93	dphi = (jth_pos.Phi() - center_phi);
94
95	// Check the 2*pi problem for delta_phi
96	if (dphi > k_PI)
97	dphi -= 2.*k_PI;
98	if (dphi < -k_PI)
99	dphi += 2.*k_PI;
100
101	delta_phi += dphi*weight;
102	}
103
104	delta_theta /= total_weight;
105	delta_phi /= total_weight;
106
107	double cluster_theta = center_theta + delta_theta;
108	double cluster_phi = center_phi + delta_phi;
109
110	// Check the 2*pi problem for cluster_phi
111	if (cluster_phi > k_PI)
112	cluster_phi -= 2.*k_PI;
113	if (cluster_phi < -k_PI)
114	cluster_phi += 2.*k_PI;
115
116	double radius = sqrt(center_pos.x()*center_pos.x()
117	+ center_pos.y()*center_pos.y()
118	+ center_pos.z()*center_pos.z());
119
120	float xpos = radiuscos(cluster_phi)sin(cluster_theta);
121	float ypos = radiussin(cluster_phi)sin(cluster_theta);
122	float zpos = radius*cos(cluster_theta);
123
124
125	posi[0] = xpos;
126	posi[1] = ypos;
127	posi[2] = zpos;
128
129
130
131
132	}
133
134
135
136
137	void calcClusterLocationEndcap(float eTot,float energy[],int ixHit[],int iyHit[],int izHit[],int nhit,float posi[]){
138
139	//float param_T0_barl_ = 7.4;
140	float param_T0_endc_ = 3.1;
141	float param_t0_endcPresh = 1.2;
142	float param_W0_ = 4.2;
143	float param_X0_ = 0.89;
144
145	///no preshower..
146	bool preshower = false;
147
148	if(eTot<=0 \|\| nhit <1){
149	cout<<"fatal error calcClusterLocationEndcap ..negative energy "<<eTot<<" "<<nhit<<endl;
150	exit(1);
151	}
152
153
154	int nmaxhit = 0;
155	float eMax = 0;
156
157	for(int j=0; j<nhit; j++){
158	if(eMax < energy[j]){
159	eMax = energy[j];
160	nmaxhit = j;
161	}
162	}
163
164
165	int ix = ixHit[nmaxhit];
166	int iy = iyHit[nmaxhit];
167	int iz = izHit[nmaxhit];
168	int indz = iz > 0;
169
170	TVector3 vmaxHit(xEEAll[indz][ix][iy],yEEAll[indz][ix][iy],zEEAll[indz][ix][iy]);
171
172
173	float T0 = param_T0_endc_;
174	if(preshower){
175	T0 = param_t0_endcPresh;
176	}
177
178	if(fabs(vmaxHit.Eta())<1.482){
179	cout<<"error..encap..?"<<vmaxHit.Eta()<<endl;
180	exit(1);
181	}
182
183	float depth = param_X0_ * (T0 + log(eTot));
184
185	TVector3 dvmaxHit(depthdxEEAll[indz][ix][iy],depthdyEEAll[indz][ix][iy],depth*dzEEAll[indz][ix][iy]);
186
187	TVector3 center_pos = vmaxHit + dvmaxHit;
188
189
190
191
192	// Loop over hits and get weights
193	double weight = 0;
194	double total_weight = 0;
195
196	double center_phi = center_pos.Phi();
197	double center_theta = center_pos.Theta();
198
199	double delta_theta = 0;
200	double delta_phi = 0;
201
202	float w_min = 0.;
203
204	double dphi = 0;
205
206	const double k_PI = acos(-1);
207
208
209	for(int j=0; j<nhit; j++){
210
211	double e_j = energy[j];
212
213	if(e_j<0) continue;
214
215
216	weight = param_W0_ +log ( e_j/eTot);
217	if(weight<w_min) weight = w_min;
218
219	total_weight += weight;
220
221	// ieta = ietaHit[j] + 85;
222	/// iphi = iphiHit[j];
223
224	ix = ixHit[j];
225	iy = iyHit[j];
226	iz = izHit[nmaxhit];
227	indz = iz > 0;
228
229	TVector3 vHit(xEEAll[indz][ix][iy],yEEAll[indz][ix][iy],zEEAll[indz][ix][iy]);
230	TVector3 dvHit(depthdxEEAll[indz][ix][iy],depthdyEEAll[indz][ix][iy],depth*dzEEAll[indz][ix][iy]);
231
232
233	TVector3 jth_pos = vHit + dvHit;
234
235
236	delta_theta += weight * (jth_pos.Theta() - center_theta);
237	dphi = (jth_pos.Phi() - center_phi);
238
239	// Check the 2*pi problem for delta_phi
240	if (dphi > k_PI)
241	dphi -= 2.*k_PI;
242	if (dphi < -k_PI)
243	dphi += 2.*k_PI;
244
245	delta_phi += dphi*weight;
246	}
247
248	delta_theta /= total_weight;
249	delta_phi /= total_weight;
250
251	double cluster_theta = center_theta + delta_theta;
252	double cluster_phi = center_phi + delta_phi;
253
254	// Check the 2*pi problem for cluster_phi
255	if (cluster_phi > k_PI)
256	cluster_phi -= 2.*k_PI;
257	if (cluster_phi < -k_PI)
258	cluster_phi += 2.*k_PI;
259
260	double radius = sqrt(center_pos.x()*center_pos.x()
261	+ center_pos.y()*center_pos.y()
262	+ center_pos.z()*center_pos.z());
263
264	float xpos = radiuscos(cluster_phi)sin(cluster_theta);
265	float ypos = radiussin(cluster_phi)sin(cluster_theta);
266	float zpos = radius*cos(cluster_theta);
267
268
269	posi[0] = xpos;
270	posi[1] = ypos;
271	posi[2] = zpos;
272
273
274	}