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
Log Message:	first commit
#	User	Rev	Content
1	yangyong	1.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			}