| 1 |
yangyong |
1.1 |
|
| 2 |
|
|
|
| 3 |
|
|
int convertIetaIphiToSMNumber(int ieta, int iphi){
|
| 4 |
|
|
|
| 5 |
|
|
if( ! (abs(ieta)>=1 && abs(ieta)<=85 && iphi>=1 && iphi<=360)) {
|
| 6 |
|
|
cout<<" convertSMNumber " << ieta <<" "<<iphi<<endl;
|
| 7 |
|
|
exit(1);
|
| 8 |
|
|
}
|
| 9 |
|
|
int ism = (iphi-1)/20+1;
|
| 10 |
|
|
if( ieta<0) ism += 18;
|
| 11 |
|
|
return ism;
|
| 12 |
|
|
|
| 13 |
|
|
}
|
| 14 |
|
|
|
| 15 |
|
|
// Declaration of leaf types
|
| 16 |
|
|
double xEBAll[170][360];
|
| 17 |
|
|
double yEBAll[170][360];
|
| 18 |
|
|
double zEBAll[170][360];
|
| 19 |
|
|
|
| 20 |
|
|
double xctEBAll[170][360];
|
| 21 |
|
|
double yctEBAll[170][360];
|
| 22 |
|
|
double zctEBAll[170][360];
|
| 23 |
|
|
|
| 24 |
|
|
|
| 25 |
|
|
|
| 26 |
|
|
double etaEBAll[170][360];
|
| 27 |
|
|
double phiEBAll[170][360];
|
| 28 |
|
|
double coxEBAll[170][360][8];
|
| 29 |
|
|
double coyEBAll[170][360][8];
|
| 30 |
|
|
double cozEBAll[170][360][8];
|
| 31 |
|
|
double xEEAll[2][101][101];
|
| 32 |
|
|
double yEEAll[2][101][101];
|
| 33 |
|
|
double zEEAll[2][101][101];
|
| 34 |
|
|
|
| 35 |
|
|
//in the 2d plane
|
| 36 |
|
|
double coxnewEBAll[170][360][8];
|
| 37 |
|
|
double coynewEBAll[170][360][8];
|
| 38 |
|
|
|
| 39 |
|
|
|
| 40 |
|
|
double xctEEAll[2][101][101];
|
| 41 |
|
|
double yctEEAll[2][101][101];
|
| 42 |
|
|
double zctEEAll[2][101][101];
|
| 43 |
|
|
|
| 44 |
|
|
|
| 45 |
|
|
double etaEEAll[2][101][101];
|
| 46 |
|
|
double phiEEAll[2][101][101];
|
| 47 |
|
|
|
| 48 |
|
|
double coxEEAll[2][101][101][8];
|
| 49 |
|
|
double coyEEAll[2][101][101][8];
|
| 50 |
|
|
double cozEEAll[2][101][101][8];
|
| 51 |
|
|
|
| 52 |
|
|
double coxnewEEAll[2][101][101][8];
|
| 53 |
|
|
double coynewEEAll[2][101][101][8];
|
| 54 |
|
|
|
| 55 |
|
|
|
| 56 |
|
|
//information on the crystal's x/y/z needed to compute cluster position X/Y/Z after
|
| 57 |
|
|
//each inter-calibration step.
|
| 58 |
|
|
void get_xyzEBrechits(){
|
| 59 |
|
|
TChain *ch = new TChain("Analysis");
|
| 60 |
|
|
ch->Add("xyzECAL.root");
|
| 61 |
|
|
|
| 62 |
|
|
ch->SetBranchAddress("xEBAll",xEBAll);
|
| 63 |
|
|
ch->SetBranchAddress("yEBAll",yEBAll);
|
| 64 |
|
|
ch->SetBranchAddress("zEBAll",zEBAll);
|
| 65 |
|
|
|
| 66 |
|
|
ch->SetBranchAddress("coxEBAll",coxEBAll);
|
| 67 |
|
|
ch->SetBranchAddress("coyEBAll",coyEBAll);
|
| 68 |
|
|
ch->SetBranchAddress("cozEBAll",cozEBAll);
|
| 69 |
|
|
|
| 70 |
|
|
|
| 71 |
|
|
ch->SetBranchAddress("etaEBAll",etaEBAll);
|
| 72 |
|
|
ch->SetBranchAddress("phiEBAll",phiEBAll);
|
| 73 |
|
|
|
| 74 |
|
|
// ch->SetBranchAddress("dxEBAll",dxEBAll);
|
| 75 |
|
|
// ch->SetBranchAddress("dyEBAll",dyEBAll);
|
| 76 |
|
|
// ch->SetBranchAddress("dzEBAll",dzEBAll);
|
| 77 |
|
|
|
| 78 |
|
|
ch->SetBranchAddress("xEEAll",xEEAll);
|
| 79 |
|
|
ch->SetBranchAddress("yEEAll",yEEAll);
|
| 80 |
|
|
ch->SetBranchAddress("zEEAll",zEEAll);
|
| 81 |
|
|
|
| 82 |
|
|
|
| 83 |
|
|
ch->SetBranchAddress("coxEEAll",coxEEAll);
|
| 84 |
|
|
ch->SetBranchAddress("coyEEAll",coyEEAll);
|
| 85 |
|
|
ch->SetBranchAddress("cozEEAll",cozEEAll);
|
| 86 |
|
|
|
| 87 |
|
|
|
| 88 |
|
|
// ch->SetBranchAddress("dxEEAll",dxEEAll);
|
| 89 |
|
|
// ch->SetBranchAddress("dyEEAll",dyEEAll);
|
| 90 |
|
|
// ch->SetBranchAddress("dzEEAll",dzEEAll);
|
| 91 |
|
|
|
| 92 |
|
|
|
| 93 |
|
|
|
| 94 |
|
|
ch->SetBranchAddress("etaEEAll",etaEEAll);
|
| 95 |
|
|
ch->SetBranchAddress("phiEEAll",phiEEAll);
|
| 96 |
|
|
|
| 97 |
|
|
|
| 98 |
|
|
|
| 99 |
|
|
|
| 100 |
|
|
|
| 101 |
|
|
ch->GetEntry(0);
|
| 102 |
|
|
|
| 103 |
|
|
ch->Delete();
|
| 104 |
|
|
|
| 105 |
|
|
cout<<"got xyzEcal.."<<endl;
|
| 106 |
|
|
|
| 107 |
|
|
|
| 108 |
|
|
|
| 109 |
|
|
}
|
| 110 |
|
|
|
| 111 |
|
|
/////////////////////////////////////////////////////////////////////////////////////////////////
|
| 112 |
|
|
|
| 113 |
|
|
void convxtalid(Int_t &nphi,Int_t &neta)
|
| 114 |
|
|
{
|
| 115 |
|
|
// Changed to what Yong's convention; output will give just two indices
|
| 116 |
|
|
// phi is unchanged; only eta now runs from
|
| 117 |
|
|
//
|
| 118 |
|
|
// 03/01/2008 changed to the new definition in CMSSW. The output is still the same...
|
| 119 |
|
|
// Barrel only
|
| 120 |
|
|
// Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
|
| 121 |
|
|
// neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
|
| 122 |
|
|
|
| 123 |
|
|
if(neta > 0) neta -= 1;
|
| 124 |
|
|
if(nphi > 359) nphi=nphi-360;
|
| 125 |
|
|
|
| 126 |
|
|
// final check
|
| 127 |
|
|
if(nphi >359 || nphi <0 || neta< -85 || neta > 84)
|
| 128 |
|
|
{
|
| 129 |
|
|
cout <<" output not in range: "<< nphi << " " << neta << " " <<endl;
|
| 130 |
|
|
exit(1);
|
| 131 |
|
|
}
|
| 132 |
|
|
} //end of convxtalid
|
| 133 |
|
|
|
| 134 |
|
|
|
| 135 |
|
|
// Calculate the distance in xtals taking into account possibly different sides
|
| 136 |
|
|
// change to coincide with yongs definition
|
| 137 |
|
|
Int_t diff_neta(Int_t neta1, Int_t neta2){
|
| 138 |
|
|
Int_t mdiff;
|
| 139 |
|
|
mdiff=abs(neta1-neta2);
|
| 140 |
|
|
return mdiff;
|
| 141 |
|
|
}
|
| 142 |
|
|
|
| 143 |
|
|
// Calculate the absolute distance in xtals taking into account the periodicity of the Barrel
|
| 144 |
|
|
Int_t diff_nphi(Int_t nphi1,Int_t nphi2) {
|
| 145 |
|
|
Int_t mdiff;
|
| 146 |
|
|
mdiff=abs(nphi1-nphi2);
|
| 147 |
|
|
if (mdiff > (360-abs(nphi1-nphi2))) mdiff=(360-abs(nphi1-nphi2));
|
| 148 |
|
|
return mdiff;
|
| 149 |
|
|
}
|
| 150 |
|
|
|
| 151 |
|
|
// Calculate the distance in xtals taking into account possibly different sides
|
| 152 |
|
|
// Then the distance would be from the 1st to the 2nd argument
|
| 153 |
|
|
// _s means that it gives the sign; the only difference from the above !
|
| 154 |
|
|
// also changed to coincide with Yong's definition
|
| 155 |
|
|
Int_t diff_neta_s(Int_t neta1, Int_t neta2){
|
| 156 |
|
|
Int_t mdiff;
|
| 157 |
|
|
mdiff=(neta1-neta2);
|
| 158 |
|
|
return mdiff;
|
| 159 |
|
|
}
|
| 160 |
|
|
|
| 161 |
|
|
// Calculate the distance in xtals taking into account the periodicity of the Barrel
|
| 162 |
|
|
Int_t diff_nphi_s(Int_t nphi1,Int_t nphi2) {
|
| 163 |
|
|
Int_t mdiff;
|
| 164 |
|
|
if(abs(nphi1-nphi2) < (360-abs(nphi1-nphi2))) {
|
| 165 |
|
|
mdiff=nphi1-nphi2;
|
| 166 |
|
|
}
|
| 167 |
|
|
else {
|
| 168 |
|
|
mdiff=360-abs(nphi1-nphi2);
|
| 169 |
|
|
if(nphi1>nphi2) mdiff=-mdiff;
|
| 170 |
|
|
}
|
| 171 |
|
|
return mdiff;
|
| 172 |
|
|
}
|
| 173 |
|
|
|
| 174 |
|
|
|
| 175 |
|
|
///to access xEBAll[ieta][iphi]
|
| 176 |
|
|
/// input ieta -85,0,84
|
| 177 |
|
|
int getIndetaxyzEBAll(int ieta){
|
| 178 |
|
|
return ieta+85;
|
| 179 |
|
|
}
|
| 180 |
|
|
|
| 181 |
|
|
////something not consistent with 167,152?
|
| 182 |
|
|
|
| 183 |
|
|
|
| 184 |
|
|
///input 0, 359 after convxtalid
|
| 185 |
|
|
int getIndphixyzEBAll(int iphi){
|
| 186 |
|
|
|
| 187 |
|
|
iphi = iphi-1;
|
| 188 |
|
|
if(iphi<0) return 359;
|
| 189 |
|
|
else return iphi;
|
| 190 |
|
|
|
| 191 |
|
|
}
|
| 192 |
|
|
|
| 193 |
|
|
double DeltaPhi(double phi1, double phi2){
|
| 194 |
|
|
|
| 195 |
|
|
double diff = fabs(phi2 - phi1);
|
| 196 |
|
|
|
| 197 |
|
|
while (diff >acos(-1)) diff -= 2*acos(-1);
|
| 198 |
|
|
while (diff <= -acos(-1)) diff += 2*acos(-1);
|
| 199 |
|
|
|
| 200 |
|
|
return diff;
|
| 201 |
|
|
|
| 202 |
|
|
}
|
| 203 |
|
|
|
| 204 |
|
|
|
| 205 |
|
|
double GetDeltaR(double eta1, double phi1,double eta2, double phi2){
|
| 206 |
|
|
|
| 207 |
|
|
return sqrt( (eta1-eta2)*(eta1-eta2)
|
| 208 |
|
|
+ DeltaPhi(phi1, phi2)*DeltaPhi(phi1, phi2) );
|
| 209 |
|
|
|
| 210 |
|
|
}
|
| 211 |
|
|
|
| 212 |
|
|
|
| 213 |
|
|
|
| 214 |
|
|
|
| 215 |
|
|
Float_t getcosd(Float_t eta1, Float_t phi1, Float_t eta2, Float_t phi2) {
|
| 216 |
|
|
Float_t theta1 = 2*atan(exp(-eta1));
|
| 217 |
|
|
Float_t theta2 = 2*atan(exp(-eta2));
|
| 218 |
|
|
Float_t cosd;
|
| 219 |
|
|
Float_t dphi = DeltaPhi(phi1,phi2);
|
| 220 |
|
|
cosd = cos(theta1)*cos(theta2)+sin(theta1)*sin(theta2)*cos(dphi); //opening angle
|
| 221 |
|
|
return cosd;
|
| 222 |
|
|
}
|
| 223 |
|
|
void separation(Float_t sceta1, Float_t scphi1, Float_t sceta2, Float_t scphi2, Float_t &dr)
|
| 224 |
|
|
{
|
| 225 |
|
|
float dphi=fabs(scphi1-scphi2);
|
| 226 |
|
|
if(dphi > (2*acos(-1)-fabs(scphi1-scphi2))) dphi = ( 2*acos(-1)-fabs(scphi1-scphi2));
|
| 227 |
|
|
dr=sqrt((sceta1- sceta2)*(sceta1- sceta2)+dphi*dphi);
|
| 228 |
|
|
}
|
| 229 |
|
|
|
| 230 |
|
|
|
| 231 |
|
|
|
| 232 |
|
|
void phinorm(Float_t & PHI)
|
| 233 |
|
|
{
|
| 234 |
|
|
while (PHI<0) PHI= PHI + 2*acos(-1);
|
| 235 |
|
|
if(PHI>2*acos(-1)) PHI= PHI - 2*acos(-1);
|
| 236 |
|
|
|
| 237 |
|
|
}
|
| 238 |
|
|
|
| 239 |
|
|
|
| 240 |
|
|
////change to [-pi,pi];
|
| 241 |
|
|
float phinorm2(float phi){
|
| 242 |
|
|
while( phi > acos(-1) ) phi -= 2*acos(-1);
|
| 243 |
|
|
while(phi< -acos(-1)) phi += 2*acos(-1);
|
| 244 |
|
|
|
| 245 |
|
|
return phi;
|
| 246 |
|
|
|
| 247 |
|
|
|
| 248 |
|
|
}
|
| 249 |
|
|
|
| 250 |
|
|
|
| 251 |
|
|
|
| 252 |
|
|
void convertindTTindCrystal(int ietaT,int iphiT, int &ieta, int &iphi){
|
| 253 |
|
|
|
| 254 |
|
|
if( iphiT >= 71) iphi = (iphiT-71)*5+3;
|
| 255 |
|
|
else if( iphiT >=1 && iphiT <=70){
|
| 256 |
|
|
iphi = (iphiT+1)*5+3;
|
| 257 |
|
|
}else{
|
| 258 |
|
|
cout<<"fatal error ..iphiT "<<iphiT<<endl;
|
| 259 |
|
|
exit(1);
|
| 260 |
|
|
}
|
| 261 |
|
|
|
| 262 |
|
|
ieta = abs(ietaT)*5-2;
|
| 263 |
|
|
|
| 264 |
|
|
if(ietaT<0) ieta *= -1;
|
| 265 |
|
|
|
| 266 |
|
|
|
| 267 |
|
|
if( ieta >85 || ieta <-85 || ieta ==0){
|
| 268 |
|
|
cout<<"convertindTTindCrystal... "<<ieta<<" "<<ietaT<<endl;
|
| 269 |
|
|
exit(1);
|
| 270 |
|
|
}
|
| 271 |
|
|
|
| 272 |
|
|
if( iphi <0 || iphi >360){
|
| 273 |
|
|
cout<<"convertindTTindCrystal... "<<iphi<<" "<<iphiT<<endl;
|
| 274 |
|
|
exit(1);
|
| 275 |
|
|
}
|
| 276 |
|
|
|
| 277 |
|
|
}
|
| 278 |
|
|
|
| 279 |
|
|
|
| 280 |
|
|
double ecalEta(double EtaParticle ,double Zvertex, double RhoVertex){
|
| 281 |
|
|
|
| 282 |
|
|
|
| 283 |
|
|
// const Double_t PI = 3.1415927;
|
| 284 |
|
|
double PI = acos(-1);
|
| 285 |
|
|
|
| 286 |
|
|
//---Definitions for ECAL
|
| 287 |
|
|
double R_ECAL = 136.5;
|
| 288 |
|
|
double Z_Endcap = 328.0;
|
| 289 |
|
|
double etaBarrelEndcap = 1.479;
|
| 290 |
|
|
|
| 291 |
|
|
if (EtaParticle!= 0.)
|
| 292 |
|
|
{
|
| 293 |
|
|
double Theta = 0.0 ;
|
| 294 |
|
|
double ZEcal = (R_ECAL-RhoVertex)*sinh(EtaParticle)+Zvertex;
|
| 295 |
|
|
|
| 296 |
|
|
if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
|
| 297 |
|
|
if(Theta<0.0) Theta = Theta+PI;
|
| 298 |
|
|
|
| 299 |
|
|
double ETA = - log(tan(0.5*Theta));
|
| 300 |
|
|
|
| 301 |
|
|
if( fabs(ETA) > etaBarrelEndcap )
|
| 302 |
|
|
{
|
| 303 |
|
|
double Zend = Z_Endcap ;
|
| 304 |
|
|
if(EtaParticle<0.0 ) Zend = -Zend ;
|
| 305 |
|
|
double Zlen = Zend - Zvertex ;
|
| 306 |
|
|
double RR = Zlen/sinh(EtaParticle);
|
| 307 |
|
|
Theta = atan((RR+RhoVertex)/Zend);
|
| 308 |
|
|
if(Theta<0.0) Theta = Theta+PI;
|
| 309 |
|
|
ETA = - log(tan(0.5*Theta));
|
| 310 |
|
|
}
|
| 311 |
|
|
return ETA;
|
| 312 |
|
|
}
|
| 313 |
|
|
else
|
| 314 |
|
|
{
|
| 315 |
|
|
return EtaParticle;
|
| 316 |
|
|
}
|
| 317 |
|
|
}
|
| 318 |
|
|
|
| 319 |
|
|
|
| 320 |
|
|
|
| 321 |
|
|
double ecalPhi(double phi,double x0,double y0){
|
| 322 |
|
|
|
| 323 |
|
|
//double R_ECAL = 136.5; ///cm
|
| 324 |
|
|
double r = 136.5;
|
| 325 |
|
|
|
| 326 |
|
|
double r0 = sqrt(x0*x0 + y0*y0);
|
| 327 |
|
|
|
| 328 |
|
|
if(r0<1E-5) return phi;
|
| 329 |
|
|
|
| 330 |
|
|
if( r0 >= r){
|
| 331 |
|
|
cout<<"warning. ecalPhi vtx outside ecal return input" << r0 <<" "<< r <<endl;
|
| 332 |
|
|
return phi;
|
| 333 |
|
|
}
|
| 334 |
|
|
|
| 335 |
|
|
double theta0 ;
|
| 336 |
|
|
if(fabs(y0)>0) theta0= y0/fabs(y0) * acos(x0/r0);
|
| 337 |
|
|
else theta0 = acos(x0/r0);
|
| 338 |
|
|
|
| 339 |
|
|
/// cout<<theta0<<" "<<phi<<endl;
|
| 340 |
|
|
|
| 341 |
|
|
double theta = phi + asin( r0/r *sin(theta0-phi));
|
| 342 |
|
|
|
| 343 |
|
|
//phinorm2(theta);
|
| 344 |
|
|
double PI = acos(-1);
|
| 345 |
|
|
while ( theta < -PI) theta += PI;
|
| 346 |
|
|
while ( theta > PI) theta -= PI;
|
| 347 |
|
|
|
| 348 |
|
|
return theta;
|
| 349 |
|
|
|
| 350 |
|
|
|
| 351 |
|
|
}
|
| 352 |
|
|
|
| 353 |
|
|
|
| 354 |
|
|
|
| 355 |
|
|
int simpleEcalID(int ieta, int iphi,int iz=0){
|
| 356 |
|
|
if(iz==0) return ieta*360+iphi;
|
| 357 |
|
|
else{
|
| 358 |
|
|
return iz*(100000 + ieta * 100 + iphi);
|
| 359 |
|
|
}
|
| 360 |
|
|
|
| 361 |
|
|
}
|
| 362 |
|
|
|
| 363 |
|
|
|
| 364 |
|
|
//return 5x5 crystals around maximum
|
| 365 |
|
|
vector<int> get5x5CrystalSim(float geta,float gphi,float vx,float vy,float vz){
|
| 366 |
|
|
////EB;-85,85 , 360 ; ieta*360+iphi
|
| 367 |
|
|
/// EE: 100*100; iz*(100000 + ix * iy)
|
| 368 |
|
|
|
| 369 |
yangyong |
1.2 |
|
| 370 |
|
|
float etag = ecalEta(geta,vz,sqrt(vx*vx+vy*vy));
|
| 371 |
yangyong |
1.1 |
float phig = ecalPhi(gphi,vx,vy);
|
| 372 |
|
|
|
| 373 |
yangyong |
1.2 |
if(debug) cout<<" eta/phi " << geta<<" "<<gphi<<" ecal "<<etag <<" "<<phig<<" v "<<vx<<" "<<vy<<" "<<vz<<" e "<< ptGenPht[0]/sin(2*atan(exp(-geta)))<<endl;
|
| 374 |
yangyong |
1.1 |
|
| 375 |
|
|
vector<int> vj;
|
| 376 |
|
|
if( fabs(etag)<1.5){
|
| 377 |
|
|
|
| 378 |
|
|
float drmin = 0.1;
|
| 379 |
|
|
float emax = 0;
|
| 380 |
|
|
int jmax = -1;
|
| 381 |
|
|
for(int j=0; j< nsimEB;j++){
|
| 382 |
|
|
|
| 383 |
|
|
int ieta1 = ietasimEB[j];
|
| 384 |
|
|
int iphi1 = iphisimEB[j];
|
| 385 |
|
|
convxtalid(iphi1,ieta1);
|
| 386 |
|
|
int ieta = ieta1+85;
|
| 387 |
|
|
int iphi = getIndphixyzEBAll(iphi1);
|
| 388 |
|
|
TVector3 vv(xEBAll[ieta][iphi],yEBAll[ieta][iphi],zEBAll[ieta][iphi]);
|
| 389 |
|
|
float dr = GetDeltaR(etag,phig,vv.Eta(),vv.Phi());
|
| 390 |
|
|
if(dr>drmin) continue;
|
| 391 |
|
|
float e = esumsimEB[j];
|
| 392 |
|
|
if(emax <e) {
|
| 393 |
|
|
emax = e;
|
| 394 |
|
|
jmax = j;
|
| 395 |
|
|
}
|
| 396 |
|
|
}
|
| 397 |
|
|
|
| 398 |
|
|
if(jmax<0) return vj;
|
| 399 |
|
|
|
| 400 |
|
|
vj.push_back(jmax);
|
| 401 |
|
|
int ietam = ietasimEB[jmax];
|
| 402 |
|
|
int iphim = iphisimEB[jmax];
|
| 403 |
|
|
convxtalid(iphim,ietam);
|
| 404 |
|
|
|
| 405 |
|
|
float esum = esumsimEB[jmax];
|
| 406 |
|
|
|
| 407 |
|
|
for(int j=0; j< nsimEB;j++){
|
| 408 |
|
|
if(j==jmax) continue;
|
| 409 |
|
|
int ieta1 = ietasimEB[j];
|
| 410 |
|
|
int iphi1 = iphisimEB[j];
|
| 411 |
|
|
convxtalid(iphi1,ieta1);
|
| 412 |
|
|
|
| 413 |
|
|
if(diff_neta(ietam,ieta1)<=2 && diff_nphi(iphim,iphi1)<=2){
|
| 414 |
|
|
vj.push_back(j);
|
| 415 |
|
|
esum += esumsimEB[j];
|
| 416 |
|
|
}
|
| 417 |
|
|
}
|
| 418 |
|
|
if(debug) cout<<"esum " << esum<<endl;
|
| 419 |
|
|
|
| 420 |
|
|
}else{
|
| 421 |
|
|
|
| 422 |
|
|
float drmin = 0.1;
|
| 423 |
|
|
float emax = 0;
|
| 424 |
|
|
int jmax = -1;
|
| 425 |
|
|
for(int j=0; j< nsimEE;j++){
|
| 426 |
|
|
|
| 427 |
|
|
int ieta = ixsimEE[j];
|
| 428 |
|
|
int iphi = iysimEE[j];
|
| 429 |
|
|
int indz = izsimEE[j]>0;
|
| 430 |
|
|
TVector3 vv(xEEAll[indz][ieta][iphi],yEEAll[indz][ieta][iphi],zEEAll[indz][ieta][iphi]);
|
| 431 |
|
|
float dr = GetDeltaR(etag,phig,vv.Eta(),vv.Phi());
|
| 432 |
|
|
if(dr>drmin) continue;
|
| 433 |
|
|
float e = esumsimEE[j];
|
| 434 |
|
|
if(emax <e) {
|
| 435 |
|
|
emax = e;
|
| 436 |
|
|
jmax = j;
|
| 437 |
|
|
}
|
| 438 |
|
|
}
|
| 439 |
|
|
|
| 440 |
|
|
if(jmax<0) return vj;
|
| 441 |
|
|
|
| 442 |
|
|
vj.push_back(jmax);
|
| 443 |
|
|
int ietam = ixsimEE[jmax];
|
| 444 |
|
|
int iphim = iysimEE[jmax];
|
| 445 |
|
|
|
| 446 |
|
|
float esum = esumsimEE[jmax];
|
| 447 |
|
|
|
| 448 |
|
|
|
| 449 |
|
|
for(int j=0; j< nsimEE;j++){
|
| 450 |
|
|
if(j==jmax) continue;
|
| 451 |
|
|
int ieta1 = ixsimEE[j];
|
| 452 |
|
|
int iphi1 = iysimEE[j];
|
| 453 |
|
|
if(diff_neta(ietam,ieta1)<=2 && diff_nphi(iphim,iphi1)<=2){
|
| 454 |
|
|
vj.push_back(j);
|
| 455 |
|
|
esum += esumsimEE[j];
|
| 456 |
|
|
}
|
| 457 |
|
|
}
|
| 458 |
|
|
if(debug) cout<<"esum " << esum <<endl;
|
| 459 |
|
|
|
| 460 |
|
|
}
|
| 461 |
|
|
|
| 462 |
|
|
return vj;
|
| 463 |
|
|
|
| 464 |
|
|
}
|
| 465 |
|
|
|
| 466 |
|
|
//return 5x5 crystals around maximum
|
| 467 |
|
|
vector<int> get5x5CrystalStep(float geta,float gphi,float vx,float vy,float vz){
|
| 468 |
|
|
////EB;-85,85 , 360 ; ieta*360+iphi
|
| 469 |
|
|
/// EE: 100*100; iz*(100000 + ix * iy)
|
| 470 |
|
|
|
| 471 |
yangyong |
1.2 |
float etag = ecalEta(geta,vz,sqrt(vx*vx+vy*vy));
|
| 472 |
yangyong |
1.1 |
float phig = ecalPhi(gphi,vx,vy);
|
| 473 |
|
|
|
| 474 |
yangyong |
1.2 |
if(debug) cout<<" eta/phi " << geta<<" "<<gphi<<" ecal "<<etag <<" "<<phig<<" v "<<vx<<" "<<vy<<" "<<vz<<" e "<< ptGenPht[0]/sin(2*atan(exp(-geta)))<<endl;
|
| 475 |
yangyong |
1.1 |
|
| 476 |
|
|
vector<int> vj;
|
| 477 |
|
|
if( fabs(etag)<1.5){
|
| 478 |
|
|
|
| 479 |
|
|
float drmin = 0.1;
|
| 480 |
|
|
float emax = 0;
|
| 481 |
|
|
int jmax = -1;
|
| 482 |
|
|
if(debug) cout<<"ng4EB " << ng4EB<<endl;
|
| 483 |
|
|
|
| 484 |
|
|
for(int j=0; j< ng4EB;j++){
|
| 485 |
|
|
|
| 486 |
|
|
int ieta1 = ietag4EB[j];
|
| 487 |
|
|
int iphi1 = iphig4EB[j];
|
| 488 |
|
|
convxtalid(iphi1,ieta1);
|
| 489 |
|
|
int ieta = ieta1+85;
|
| 490 |
|
|
int iphi = getIndphixyzEBAll(iphi1);
|
| 491 |
|
|
TVector3 vv(xEBAll[ieta][iphi],yEBAll[ieta][iphi],zEBAll[ieta][iphi]);
|
| 492 |
|
|
float dr = GetDeltaR(etag,phig,vv.Eta(),vv.Phi());
|
| 493 |
|
|
if(dr>drmin) continue;
|
| 494 |
|
|
float e = esumg4EB[j];
|
| 495 |
|
|
if(emax <e) {
|
| 496 |
|
|
emax = e;
|
| 497 |
|
|
jmax = j;
|
| 498 |
|
|
}
|
| 499 |
|
|
}
|
| 500 |
|
|
|
| 501 |
|
|
if(jmax<0) return vj;
|
| 502 |
|
|
|
| 503 |
|
|
if(debug) cout<<"jmax "<< jmax <<endl;
|
| 504 |
|
|
|
| 505 |
|
|
vj.push_back(jmax);
|
| 506 |
|
|
int ietam = ietag4EB[jmax];
|
| 507 |
|
|
int iphim = iphig4EB[jmax];
|
| 508 |
|
|
convxtalid(iphim,ietam);
|
| 509 |
|
|
float esum = esumg4EB[jmax];
|
| 510 |
|
|
|
| 511 |
|
|
// vector<float> ev = eg4EB->at(jmax);
|
| 512 |
|
|
// float testesum = 0;
|
| 513 |
|
|
// for(int j=0; j<int(ev.size()); j++){
|
| 514 |
|
|
// testesum += ev[j];
|
| 515 |
|
|
// }
|
| 516 |
|
|
|
| 517 |
|
|
if(debug) cout<<"emaxeb; " << esum <<endl;
|
| 518 |
|
|
|
| 519 |
|
|
for(int j=0; j< ng4EB;j++){
|
| 520 |
|
|
if(j==jmax) continue;
|
| 521 |
|
|
int ieta1 = ietag4EB[j];
|
| 522 |
|
|
int iphi1 = iphig4EB[j];
|
| 523 |
|
|
convxtalid(iphi1,ieta1);
|
| 524 |
|
|
|
| 525 |
|
|
if(diff_neta(ietam,ieta1)<=2 && diff_nphi(iphim,iphi1)<=2){
|
| 526 |
|
|
vj.push_back(j);
|
| 527 |
|
|
esum += esumg4EB[j];
|
| 528 |
|
|
}
|
| 529 |
|
|
}
|
| 530 |
|
|
if(debug) cout<<"esum " << esum<<endl;
|
| 531 |
|
|
|
| 532 |
|
|
}else{
|
| 533 |
|
|
|
| 534 |
|
|
float drmin = 0.1;
|
| 535 |
|
|
float emax = 0;
|
| 536 |
|
|
int jmax = -1;
|
| 537 |
|
|
for(int j=0; j< ng4EE;j++){
|
| 538 |
|
|
|
| 539 |
|
|
int ieta = ixg4EE[j];
|
| 540 |
|
|
int iphi = iyg4EE[j];
|
| 541 |
|
|
int indz = izg4EE[j]>0;
|
| 542 |
|
|
TVector3 vv(xEEAll[indz][ieta][iphi],yEEAll[indz][ieta][iphi],zEEAll[indz][ieta][iphi]);
|
| 543 |
|
|
float dr = GetDeltaR(etag,phig,vv.Eta(),vv.Phi());
|
| 544 |
|
|
if(dr>drmin) continue;
|
| 545 |
|
|
float e = esumg4EE[j];
|
| 546 |
|
|
if(emax <e) {
|
| 547 |
|
|
emax = e;
|
| 548 |
|
|
jmax = j;
|
| 549 |
|
|
}
|
| 550 |
|
|
}
|
| 551 |
|
|
|
| 552 |
|
|
if(jmax<0) return vj;
|
| 553 |
|
|
|
| 554 |
|
|
vj.push_back(jmax);
|
| 555 |
|
|
int ietam = ixg4EE[jmax];
|
| 556 |
|
|
int iphim = iyg4EE[jmax];
|
| 557 |
|
|
|
| 558 |
|
|
float esum = esumg4EE[jmax];
|
| 559 |
|
|
if(debug) cout<<"emaxee; " << esum <<endl;
|
| 560 |
|
|
|
| 561 |
|
|
for(int j=0; j< ng4EE;j++){
|
| 562 |
|
|
if(j==jmax) continue;
|
| 563 |
|
|
int ieta1 = ixg4EE[j];
|
| 564 |
|
|
int iphi1 = iyg4EE[j];
|
| 565 |
|
|
if(diff_neta(ietam,ieta1)<=2 && diff_nphi(iphim,iphi1)<=2){
|
| 566 |
|
|
vj.push_back(j);
|
| 567 |
|
|
esum += esumg4EE[j];
|
| 568 |
|
|
}
|
| 569 |
|
|
}
|
| 570 |
|
|
if(debug)cout<<"esum " << esum <<endl;
|
| 571 |
|
|
|
| 572 |
|
|
}
|
| 573 |
|
|
|
| 574 |
|
|
return vj;
|
| 575 |
|
|
|
| 576 |
|
|
}
|
