yangyong/zSelection/physUtils.cc



///to access xEBAll[ieta][iphi]
/// input ieta -85,0,84
int getIndetaxyzEBAll(int ieta){
  return ieta+85; 
}

////something not consistent with 167,152?


///input 0, 359 after convxtalid
int getIndphixyzEBAll(int iphi){
  
  iphi = iphi-1; 
  if(iphi<0) return 359; 
  else return iphi;
  
}

double DeltaPhi(double phi1, double phi2){
  
  //double diff = fabs(phi2 - phi1);
  double diff = phi2 - phi1;
  
  while (diff >acos(-1)) diff -= 2*acos(-1);
  while (diff <= -acos(-1)) diff += 2*acos(-1);
  
  return diff; 
  
}


double GetDeltaR(double eta1, double eta2, double phi1, double phi2){
  
  return sqrt( (eta1-eta2)*(eta1-eta2) 
               + DeltaPhi(phi1, phi2)*DeltaPhi(phi1, phi2) );
  
}


Float_t getcosd(Float_t eta1, Float_t phi1, Float_t eta2, Float_t phi2) {
  Float_t theta1 = 2*atan(exp(-eta1));
  Float_t theta2 = 2*atan(exp(-eta2));
  Float_t cosd;
  Float_t dphi = DeltaPhi(phi1,phi2);
  cosd = cos(theta1)*cos(theta2)+sin(theta1)*sin(theta2)*cos(dphi);  //opening angle 
  return cosd;
}
void separation(Float_t sceta1, Float_t scphi1, Float_t sceta2, Float_t scphi2, Float_t &dr)
{
  float dphi=fabs(scphi1-scphi2);
  if(dphi > (2*acos(-1)-fabs(scphi1-scphi2))) dphi = ( 2*acos(-1)-fabs(scphi1-scphi2));
  dr=sqrt((sceta1- sceta2)*(sceta1- sceta2)+dphi*dphi);
}


void calcPairObjects(int pid1, int pid2,float en[],float eta[],float phi[],float res[]){
  
  TLorentzVector vpht[2];
  
  TLorentzVector vpair; 
  
  float mass[2]={0,0};
  
  if(pid1==11) mass[0] = 0.511*0.001; 
  else if(pid1==13) mass[0] = 0.105658;
  else if(pid1==22) mass[0] = 0;
  else{
    cout<<"calcPairMuPht, mass.pid1 "<<endl;
    exit(1);
  }
  
  if(pid2==11) mass[1] = 0.511*0.001; 
  else if(pid2==13) mass[1] = 0.105658;
  else if(pid2 ==22) mass[1] = 0; 
  else {
    cout<<"calcPairMuPht, mass. pid2"<<endl;
    exit(1);
  }
  

  if( en[0] < mass[0] || en[1] < mass[1]) {
    cout<<"warning calcPairObjects E<m? "<< en[0]<<" "<<en[1] <<" "<< pid1<<" "<<pid2<<endl; 
    exit(1); 
  }
  
  for( int j= 0; j<2; j++){
    
    float e = sqrt( en[j] * en[j] - mass[j]*mass[j]);
    float px = e * sin ( 2*atan(exp(-eta[j]))) * cos(phi[j]);
    float py = e * sin ( 2*atan(exp(-eta[j]))) * sin(phi[j]);
    float pz = e * cos ( 2*atan(exp(-eta[j]))) ;
    
    vpht[j].SetXYZT(px,py,pz,e);
    
  }
  
  vpair = vpht[0] + vpht[1];

  res[0] = vpair.M();
  res[1] = vpair.Eta();
  res[2] = vpair.Phi();
  res[3] = vpair.Pt();
  
  
}


void calcPairPhoton(float en[],float eta[],float phi[],float res[]){
  
  TLorentzVector vpht[2];
  
  TLorentzVector vpair; 
  
  
  for( int j= 0; j<2; j++){
    
    float e = en[j];
    float px = e * sin ( 2*atan(exp(-eta[j]))) * cos(phi[j]);
    float py = e * sin ( 2*atan(exp(-eta[j]))) * sin(phi[j]);
    float pz = e * cos ( 2*atan(exp(-eta[j]))) ;
    
    vpht[j].SetXYZT(px,py,pz,e);
    
  }
  
  vpair = vpht[0] + vpht[1];

  res[0] = vpair.M();
  res[1] = vpair.Eta();
  res[2] = vpair.Phi();
  res[3] = vpair.Pt();
  
  
}

void phinorm(Float_t & PHI)
{
  while (PHI<0)  PHI= PHI + 2*acos(-1);
  if(PHI>2*acos(-1))  PHI= PHI - 2*acos(-1);
  
}


////change to [-pi,pi];
float phinorm2(float phi){
  while( phi > acos(-1) ) phi -= 2*acos(-1);
  while(phi< -acos(-1)) phi += 2*acos(-1);

  return phi;


}


///input phi : [-pi,pi];
////this is the not exactly the same as; 
////RecoEcal/EgammaCoreTools/src/LogPositionCalc.cc
///need x,y,z information of each crystal to do that. 

/////center iphi = 190 or 191 phi change signs. 3.133, -3.133

void simpleLogWeightedEtaPhi(int nxt, float esum, float energy[],float eta[],float phi[],int phimax, float res[]){
  
  
  float etasum = 0; 
  float phisum = 0; 
  float wtsum = 0; 
  
  if(phimax==190 || phimax==191){
    for( int j=0; j<nxt; j++){
      phinorm(phi[j]);
    }
  }
    
  
  for( int j=0; j<nxt; j++){
    float mw=4.2+log(fabs(energy[j])/esum);
    if(mw < 0.) mw=0.;
    wtsum += mw; 
    etasum += mw * eta[j];
    phisum += mw * phi[j];
    
  }
  
  
  etasum /= wtsum; 
  phisum /= wtsum; 
  
  ///change to [-pi,pi]
  phinorm2(phisum);
  
  
  res[0] = etasum; 
  res[1] = phisum; 
  
  ////cluster shape

  
}

// // ///cluster shape SigmaEtaEta, SigmaEtaPhi,SigmaPhiPhi
// // ///


void Calculate_ClusterCovariance(int nxt, float esum,float ceta,float cphi,float en[],float eta[],float phi[],float res[]){
  
  double numeratorEtaEta = 0;
  double numeratorEtaPhi = 0;
  double numeratorPhiPhi = 0;
  double denominator     = 0;
  
  for( int j=0; j<nxt; j++){
    
    float dPhi = phi[j] - cphi; 

    if(dPhi > acos(-1)) dPhi = 2*acos(-1) - dPhi; 
    if(dPhi <-acos(-1)) dPhi = 2*acos(-1) + dPhi; 
    
    float dEta = eta[j] - ceta; 
    
    
    float w=4.2+log(fabs(en[j])/esum);
    if(w < 0.) w=0.;
    
    
    denominator += w;
    numeratorEtaEta += w * dEta * dEta;
    numeratorEtaPhi += w * dEta * dPhi;
    numeratorPhiPhi += w * dPhi * dPhi;
    
  }
  
     
  res[0] = numeratorEtaEta / denominator;
  res[1] = numeratorEtaPhi / denominator;
  res[2] = numeratorPhiPhi / denominator;
  
    
}
   

float Calculate_LAT(int nxt, float xclus, float yclus, float zclus,float en[],float x[],float y[],float z[]){
  
  if( nxt <3) return 10; 
  
  TVector3 clVect(xclus,yclus,zclus);
  
  TVector3 clDir = clVect; 
  clDir *= 1.0/clDir.Mag();
  
  float redmoment = 0; 
  
  float e12 = en[0] + en[1]; 
  TVector3 gblPos; 
  for( int j=2; j< nxt; j++){
    gblPos.SetXYZ(x[j],y[j],z[j]);
    TVector3 diff = gblPos - clVect;

    TVector3 DigiVect =  diff - diff.Dot(clDir)*clDir;
    float r = DigiVect.Mag();
    redmoment += r*r*fabs(en[j]);
  }
  
  float lat = redmoment/(redmoment + 2.19*2.19*e12);

  return lat; 
  

}


///change [-85,84] to bin 1, to bin 170
///change [0,359] to bin 1, 260; 

void getBinEtaPhi(int eta, int phi, float res[]){
  
  if( eta <0) {
    res[0] = eta + 85 +1+0.001; 
    
    res[2] = eta; 
    
  }
  else {
    res[0] = eta + 85 +2+0.001; 
    res[2] = eta +1; 
  }
  
  if( phi==0) {
    res[1] = 360; 
    res[3] = 360; 
  }
  else {
    res[1] = phi+0.001; 
    res[3] = phi; 
  }
  
  
}


////check crystal at boader of em obj boader

int IsatBoaderEMObjPhi(int iphi){
  
  // if( iphi ==0 || iphi ==1) return 1; 

  if( iphi%20 ==0 || (iphi-1)%20 ==0) return 1; 
  
  return 0; 

}

int IsatBoaderEMObjEta(int ieta){
  
  if( ieta ==0 || ieta == -1) return 1; 
  if( ieta ==19 || ieta==20) return 1; 
  if( ieta ==39 || ieta==40) return 1; 
  if( ieta == 59 || ieta ==60) return 1; 
  
  
  if( ieta==-20 || ieta==-21) return 1; 
  if( ieta==-40 || ieta==-41) return 1;
  if( ieta==-60 || ieta==-61) return 1;
    
  
  return 0; 

}


//convert ietaTT , iphiTT to ieta,iphi of central crystal

void convertindTTindCrystal(int ietaT,int iphiT, int &ieta, int &iphi){

  if( iphiT >= 71) iphi = (iphiT-71)*5+3;
  else if( iphiT >=1 && iphiT <=70){
    iphi = (iphiT+1)*5+3;
  }else{
    cout<<"fatal error ..iphiT "<<iphiT<<endl;
    exit(1);
  }

  ieta = abs(ietaT)*5-2;

  if(ietaT<0) ieta *= -1;
  

  if( ieta >85 || ieta <-85 || ieta ==0){
    cout<<"convertindTTindCrystal...  "<<ieta<<" "<<ietaT<<endl;
    exit(1);
  }
  
  if( iphi <0 || iphi >360){
    cout<<"convertindTTindCrystal...  "<<iphi<<" "<<iphiT<<endl;
    exit(1);
  }
  
}


////invariant mass of two obejcts given, pid each. 
void calcPairPtEtaPhi(float pt[2],float eta[2],float phi[2],int pid1,int pid2,float res[]){
  //void calcPairMuPht(double pt[2],double eta[2],double phi[2],int pid1,int pid2,double res[]){
  
  TLorentzVector vpht[2];
  
  TLorentzVector vpair; 
  
  double mass[2]={0,0};
  
  if(pid1==11) mass[0] = 0.511*0.001; 
  else if(pid1==13) mass[0] = 0.105658;
  else if(pid1==22) mass[0] = 0;
  else{
    cout<<"calcPairMuPht, mass.pid1 "<<endl;
    exit(1);
  }
  
  if(pid2==11) mass[1] = 0.511*0.001; 
  else if(pid2==13) mass[1] = 0.105658;
  else if(pid2 ==22) mass[1] = 0; 
  else {
    cout<<"calcPairMuPht, mass. pid2"<<endl;
    exit(1);
  }
  
  
  for( int j= 0; j<2; j++){
    
    double p = pt[j]/sin(2*atan(exp(-eta[j])));
    double e = sqrt(p*p+mass[j]*mass[j]);
    double px = p * sin ( 2*atan(exp(-eta[j]))) * cos(phi[j]);
    double py = p * sin ( 2*atan(exp(-eta[j]))) * sin(phi[j]);
    double pz = p * cos ( 2*atan(exp(-eta[j]))) ;
    
    vpht[j].SetXYZT(px,py,pz,e);
    
  }
  
  vpair = vpht[0] + vpht[1];
  
  res[0] = vpair.M();
  res[1] = vpair.Eta();
  res[2] = vpair.Phi();
  res[3] = vpair.Pt();
  res[4] = vpht[0].DeltaR(vpht[1]);
  
  
  ///cout<<"dot: "<< vpht[0].E() <<" "<< vpht[1].E()<< " "<< vpht[0].X() <<" "<< vpht[1].X()<<" "<< vpht[0].Y() <<" "<< vpht[1].Y() <<" "<< vpht[0].Z() <<" "<< vpht[1].Z()<<endl; 
  
  
  // res[5] = vpht[0].Dot(vpht[1])/vpht[0].P()/vpht[1].P();
  
  
}


////invariant mass of two obejcts given, pid each. 
//void calcMass3Objects(float pt[3],float eta[3],float phi[3],int pid1,int pid2,int pid3,float res[]){
void calcMass3Objects(double pt[3],double eta[3],double phi[3],int pid1,int pid2,int pid3,double res[]){
  
  TLorentzVector vpht[3];
  
  TLorentzVector vpair; 
  
  float mass[3]={0,0,0};
  
  if(pid1==11) mass[0] = 0.511*0.001; 
  else if(pid1==13) mass[0] = 0.105658;
  else if(pid1==22) mass[0] = 0;
  else{
    cout<<"calcPairMuPht, mass.pid1 "<<endl;
    exit(1);
  }
  
  if(pid2==11) mass[1] = 0.511*0.001; 
  else if(pid2==13) mass[1] = 0.105658;
  else if(pid2 ==22) mass[1] = 0; 
  else {
    cout<<"calcPairMuPht, mass. pid2"<<endl;
    exit(1);
  }
  
  if(pid3==11) mass[2] = 0.511*0.001; 
  else if(pid3==13) mass[2] = 0.105658;
  else if(pid3 ==22) mass[2] = 0; 
  else {
    cout<<"calcPairMuPht, mass. pid3"<<endl;
    exit(1);
  }
  
  
  for( int j= 0; j<3; j++){
    
    float p = pt[j]/sin(2*atan(exp(-eta[j])));
    float e = sqrt(p*p+mass[j]*mass[j]);
    float px = p * sin ( 2*atan(exp(-eta[j]))) * cos(phi[j]);
    float py = p * sin ( 2*atan(exp(-eta[j]))) * sin(phi[j]);
    float pz = p * cos ( 2*atan(exp(-eta[j]))) ;
    
    vpht[j].SetXYZT(px,py,pz,e);
    
  }
  
  vpair = vpht[0] + vpht[1] + vpht[2];
  
  res[0] = vpair.M();
  res[1] = vpair.Eta();
  res[2] = vpair.Phi();
  res[3] = vpair.Pt();
  res[4] = vpht[0].DeltaR(vpht[1]);
  
  
}


void printRatioTwoHist(TH1F *h1, TH1F *h2){
  
  int nbins = h1->GetNbinsX();
  
  float xmin = h1->GetXaxis()->GetXmin() ;
  float xmax = h1->GetXaxis()->GetXmax() ;
  
  TH1F *hh = new TH1F("hh","eff",nbins,xmin,xmax);
  
    
  for(int j=1; j<= nbins; j++){
    
    float y1 = h1->GetBinContent(j);
    float y2 = h2->GetBinContent(j);

    float eff = 0; 
    float effErr = 0; 
    if(y1>0 && y2>0){
      eff = y2/y1; 
      effErr = sqrt(eff*(1-eff)/y1);
      
      hh->SetBinContent(j,eff);
      hh->SetBinError(j,effErr);
    }
  }
  
  // TCanvas *can0 = new TCanvas("can0","c000",200,10,550,500);
  hh->Draw();
  
  
}

double etaTransformation(  float EtaParticle , float Zvertex)  {
  
  //---Definitions
  const float pi = 3.1415927;

  //---Definitions for ECAL
  const float R_ECAL           = 136.5;
  const float Z_Endcap         = 328.0;
  const float etaBarrelEndcap  = 1.479; 
   
  //---ETA correction

  float Theta = 0.0  ; 
  float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;

  if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
  if(Theta<0.0) Theta = Theta+pi ;
  double ETA = - log(tan(0.5*Theta));
         
  if( fabs(ETA) > etaBarrelEndcap )
    {
      float Zend = Z_Endcap ;
      if(EtaParticle<0.0 )  Zend = -Zend ;
      float Zlen = Zend - Zvertex ;
      float RR = Zlen/sinh(EtaParticle); 
      Theta = atan(RR/Zend);
      if(Theta<0.0) Theta = Theta+pi ;
      ETA = - log(tan(0.5*Theta));                    
    } 
  //---Return the result
  return ETA;
  //---end
}

////transform eta ( z, pho), to eta at ecal ( w.r.t 0,0,0,)
double ecalEta(double EtaParticle ,double Zvertex, double RhoVertex){
  
  
  //  const Double_t PI    = 3.1415927;
  double PI    = acos(-1);
  
  //---Definitions for ECAL
  double R_ECAL           = 136.5;
  double Z_Endcap         = 328.0;
  double etaBarrelEndcap  = 1.479; 

  if (EtaParticle!= 0.)
    {
      double Theta = 0.0  ;
      double ZEcal = (R_ECAL-RhoVertex)*sinh(EtaParticle)+Zvertex;
      
      if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
      if(Theta<0.0) Theta = Theta+PI;

      double ETA = - log(tan(0.5*Theta));
      
      if( fabs(ETA) > etaBarrelEndcap )
        {
          double Zend = Z_Endcap ;
          if(EtaParticle<0.0 )  Zend = -Zend ;
          double Zlen = Zend - Zvertex ;
          double RR = Zlen/sinh(EtaParticle);
          Theta = atan((RR+RhoVertex)/Zend);
          if(Theta<0.0) Theta = Theta+PI;
          ETA = - log(tan(0.5*Theta));
        }
      return ETA;
    }
  else
    {
      return EtaParticle;
    }
}


double ecalPhi(double phi,double x0,double y0){
  
  //double R_ECAL = 136.5; ///cm 
  double r = 136.5; 

  double r0 = sqrt(x0*x0 + y0*y0);
  
  if(r0<1E-5) return phi; 
  
  if( r0 >= r){
    cout<<"warning. ecalPhi vtx outside ecal return input" << r0 <<" "<< r <<endl;
    return phi; 
  }
  
  double theta0 ;
  if(fabs(y0)>0) theta0= y0/fabs(y0) * acos(x0/r0);
  else theta0 = acos(x0/r0);
  
  ///  cout<<theta0<<" "<<phi<<endl;
  
  double theta = phi + asin( r0/r *sin(theta0-phi));

  //phinorm2(theta);
  double PI    = acos(-1);
  while ( theta < -PI) theta += PI; 
  while ( theta > PI) theta -= PI; 
  
  return theta; 
  
  
}
Revision:	1.1
Committed:	Fri Sep 30 15:23:48 2011 UTC (13 years, 7 months ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V01-00-01, V01-00-00, HEAD
Log Message:	z selection code
#	Content
1
2
3
4
5	///to access xEBAll[ieta][iphi]
6	/// input ieta -85,0,84
7	int getIndetaxyzEBAll(int ieta){
8	return ieta+85;
9	}
10
11	////something not consistent with 167,152?
12
13
14	///input 0, 359 after convxtalid
15	int getIndphixyzEBAll(int iphi){
16
17	iphi = iphi-1;
18	if(iphi<0) return 359;
19	else return iphi;
20
21	}
22
23	double DeltaPhi(double phi1, double phi2){
24
25	//double diff = fabs(phi2 - phi1);
26	double diff = phi2 - phi1;
27
28	while (diff >acos(-1)) diff -= 2*acos(-1);
29	while (diff <= -acos(-1)) diff += 2*acos(-1);
30
31	return diff;
32
33	}
34
35
36	double GetDeltaR(double eta1, double eta2, double phi1, double phi2){
37
38	return sqrt( (eta1-eta2)*(eta1-eta2)
39	+ DeltaPhi(phi1, phi2)*DeltaPhi(phi1, phi2) );
40
41	}
42
43
44
45
46	Float_t getcosd(Float_t eta1, Float_t phi1, Float_t eta2, Float_t phi2) {
47	Float_t theta1 = 2*atan(exp(-eta1));
48	Float_t theta2 = 2*atan(exp(-eta2));
49	Float_t cosd;
50	Float_t dphi = DeltaPhi(phi1,phi2);
51	cosd = cos(theta1)cos(theta2)+sin(theta1)sin(theta2)*cos(dphi); //opening angle
52	return cosd;
53	}
54	void separation(Float_t sceta1, Float_t scphi1, Float_t sceta2, Float_t scphi2, Float_t &dr)
55	{
56	float dphi=fabs(scphi1-scphi2);
57	if(dphi > (2acos(-1)-fabs(scphi1-scphi2))) dphi = ( 2acos(-1)-fabs(scphi1-scphi2));
58	dr=sqrt((sceta1- sceta2)(sceta1- sceta2)+dphidphi);
59	}
60
61
62
63	void calcPairObjects(int pid1, int pid2,float en[],float eta[],float phi[],float res[]){
64
65	TLorentzVector vpht[2];
66
67	TLorentzVector vpair;
68
69	float mass[2]={0,0};
70
71	if(pid1==11) mass[0] = 0.511*0.001;
72	else if(pid1==13) mass[0] = 0.105658;
73	else if(pid1==22) mass[0] = 0;
74	else{
75	cout<<"calcPairMuPht, mass.pid1 "<<endl;
76	exit(1);
77	}
78
79	if(pid2==11) mass[1] = 0.511*0.001;
80	else if(pid2==13) mass[1] = 0.105658;
81	else if(pid2 ==22) mass[1] = 0;
82	else {
83	cout<<"calcPairMuPht, mass. pid2"<<endl;
84	exit(1);
85	}
86
87
88	if( en[0] < mass[0] \|\| en[1] < mass[1]) {
89	cout<<"warning calcPairObjects E<m? "<< en[0]<<" "<<en[1] <<" "<< pid1<<" "<<pid2<<endl;
90	exit(1);
91	}
92
93	for( int j= 0; j<2; j++){
94
95	float e = sqrt( en[j] * en[j] - mass[j]*mass[j]);
96	float px = e * sin ( 2atan(exp(-eta[j]))) cos(phi[j]);
97	float py = e * sin ( 2atan(exp(-eta[j]))) sin(phi[j]);
98	float pz = e * cos ( 2*atan(exp(-eta[j]))) ;
99
100	vpht[j].SetXYZT(px,py,pz,e);
101
102	}
103
104	vpair = vpht[0] + vpht[1];
105
106	res[0] = vpair.M();
107	res[1] = vpair.Eta();
108	res[2] = vpair.Phi();
109	res[3] = vpair.Pt();
110
111
112
113
114	}
115
116
117	void calcPairPhoton(float en[],float eta[],float phi[],float res[]){
118
119	TLorentzVector vpht[2];
120
121	TLorentzVector vpair;
122
123
124	for( int j= 0; j<2; j++){
125
126	float e = en[j];
127	float px = e * sin ( 2atan(exp(-eta[j]))) cos(phi[j]);
128	float py = e * sin ( 2atan(exp(-eta[j]))) sin(phi[j]);
129	float pz = e * cos ( 2*atan(exp(-eta[j]))) ;
130
131	vpht[j].SetXYZT(px,py,pz,e);
132
133	}
134
135	vpair = vpht[0] + vpht[1];
136
137	res[0] = vpair.M();
138	res[1] = vpair.Eta();
139	res[2] = vpair.Phi();
140	res[3] = vpair.Pt();
141
142
143
144
145	}
146
147	void phinorm(Float_t & PHI)
148	{
149	while (PHI<0) PHI= PHI + 2*acos(-1);
150	if(PHI>2acos(-1)) PHI= PHI - 2acos(-1);
151
152	}
153
154
155	////change to [-pi,pi];
156	float phinorm2(float phi){
157	while( phi > acos(-1) ) phi -= 2*acos(-1);
158	while(phi< -acos(-1)) phi += 2*acos(-1);
159
160	return phi;
161
162
163	}
164
165
166
167
168	///input phi : [-pi,pi];
169	////this is the not exactly the same as;
170	////RecoEcal/EgammaCoreTools/src/LogPositionCalc.cc
171	///need x,y,z information of each crystal to do that.
172
173	/////center iphi = 190 or 191 phi change signs. 3.133, -3.133
174
175	void simpleLogWeightedEtaPhi(int nxt, float esum, float energy[],float eta[],float phi[],int phimax, float res[]){
176
177
178
179	float etasum = 0;
180	float phisum = 0;
181	float wtsum = 0;
182
183	if(phimax==190 \|\| phimax==191){
184	for( int j=0; j<nxt; j++){
185	phinorm(phi[j]);
186	}
187	}
188
189
190	for( int j=0; j<nxt; j++){
191	float mw=4.2+log(fabs(energy[j])/esum);
192	if(mw < 0.) mw=0.;
193	wtsum += mw;
194	etasum += mw * eta[j];
195	phisum += mw * phi[j];
196
197	}
198
199
200	etasum /= wtsum;
201	phisum /= wtsum;
202
203	///change to [-pi,pi]
204	phinorm2(phisum);
205
206
207	res[0] = etasum;
208	res[1] = phisum;
209
210	////cluster shape
211
212
213
214
215	}
216
217	// // ///cluster shape SigmaEtaEta, SigmaEtaPhi,SigmaPhiPhi
218	// // ///
219
220
221	void Calculate_ClusterCovariance(int nxt, float esum,float ceta,float cphi,float en[],float eta[],float phi[],float res[]){
222
223	double numeratorEtaEta = 0;
224	double numeratorEtaPhi = 0;
225	double numeratorPhiPhi = 0;
226	double denominator = 0;
227
228	for( int j=0; j<nxt; j++){
229
230	float dPhi = phi[j] - cphi;
231
232	if(dPhi > acos(-1)) dPhi = 2*acos(-1) - dPhi;
233	if(dPhi <-acos(-1)) dPhi = 2*acos(-1) + dPhi;
234
235	float dEta = eta[j] - ceta;
236
237
238	float w=4.2+log(fabs(en[j])/esum);
239	if(w < 0.) w=0.;
240
241
242	denominator += w;
243	numeratorEtaEta += w * dEta * dEta;
244	numeratorEtaPhi += w * dEta * dPhi;
245	numeratorPhiPhi += w * dPhi * dPhi;
246
247	}
248
249
250	res[0] = numeratorEtaEta / denominator;
251	res[1] = numeratorEtaPhi / denominator;
252	res[2] = numeratorPhiPhi / denominator;
253
254
255
256
257	}
258
259
260
261
262	float Calculate_LAT(int nxt, float xclus, float yclus, float zclus,float en[],float x[],float y[],float z[]){
263
264	if( nxt <3) return 10;
265
266	TVector3 clVect(xclus,yclus,zclus);
267
268	TVector3 clDir = clVect;
269	clDir *= 1.0/clDir.Mag();
270
271	float redmoment = 0;
272
273	float e12 = en[0] + en[1];
274	TVector3 gblPos;
275	for( int j=2; j< nxt; j++){
276	gblPos.SetXYZ(x[j],y[j],z[j]);
277	TVector3 diff = gblPos - clVect;
278
279	TVector3 DigiVect = diff - diff.Dot(clDir)*clDir;
280	float r = DigiVect.Mag();
281	redmoment += rrfabs(en[j]);
282	}
283
284	float lat = redmoment/(redmoment + 2.192.19e12);
285
286	return lat;
287
288
289	}
290
291
292	///change [-85,84] to bin 1, to bin 170
293	///change [0,359] to bin 1, 260;
294
295	void getBinEtaPhi(int eta, int phi, float res[]){
296
297	if( eta <0) {
298	res[0] = eta + 85 +1+0.001;
299
300	res[2] = eta;
301
302	}
303	else {
304	res[0] = eta + 85 +2+0.001;
305	res[2] = eta +1;
306	}
307
308	if( phi==0) {
309	res[1] = 360;
310	res[3] = 360;
311	}
312	else {
313	res[1] = phi+0.001;
314	res[3] = phi;
315	}
316
317
318
319
320	}
321
322
323
324	////check crystal at boader of em obj boader
325
326	int IsatBoaderEMObjPhi(int iphi){
327
328	// if( iphi ==0 \|\| iphi ==1) return 1;
329
330	if( iphi%20 ==0 \|\| (iphi-1)%20 ==0) return 1;
331
332	return 0;
333
334	}
335
336	int IsatBoaderEMObjEta(int ieta){
337
338	if( ieta ==0 \|\| ieta == -1) return 1;
339	if( ieta ==19 \|\| ieta==20) return 1;
340	if( ieta ==39 \|\| ieta==40) return 1;
341	if( ieta == 59 \|\| ieta ==60) return 1;
342
343
344	if( ieta==-20 \|\| ieta==-21) return 1;
345	if( ieta==-40 \|\| ieta==-41) return 1;
346	if( ieta==-60 \|\| ieta==-61) return 1;
347
348
349	return 0;
350
351	}
352
353
354
355	//convert ietaTT , iphiTT to ieta,iphi of central crystal
356
357	void convertindTTindCrystal(int ietaT,int iphiT, int &ieta, int &iphi){
358
359	if( iphiT >= 71) iphi = (iphiT-71)*5+3;
360	else if( iphiT >=1 && iphiT <=70){
361	iphi = (iphiT+1)*5+3;
362	}else{
363	cout<<"fatal error ..iphiT "<<iphiT<<endl;
364	exit(1);
365	}
366
367	ieta = abs(ietaT)*5-2;
368
369	if(ietaT<0) ieta *= -1;
370
371
372	if( ieta >85 \|\| ieta <-85 \|\| ieta ==0){
373	cout<<"convertindTTindCrystal... "<<ieta<<" "<<ietaT<<endl;
374	exit(1);
375	}
376
377	if( iphi <0 \|\| iphi >360){
378	cout<<"convertindTTindCrystal... "<<iphi<<" "<<iphiT<<endl;
379	exit(1);
380	}
381
382	}
383
384
385	////invariant mass of two obejcts given, pid each.
386	void calcPairPtEtaPhi(float pt[2],float eta[2],float phi[2],int pid1,int pid2,float res[]){
387	//void calcPairMuPht(double pt[2],double eta[2],double phi[2],int pid1,int pid2,double res[]){
388
389	TLorentzVector vpht[2];
390
391	TLorentzVector vpair;
392
393	double mass[2]={0,0};
394
395	if(pid1==11) mass[0] = 0.511*0.001;
396	else if(pid1==13) mass[0] = 0.105658;
397	else if(pid1==22) mass[0] = 0;
398	else{
399	cout<<"calcPairMuPht, mass.pid1 "<<endl;
400	exit(1);
401	}
402
403	if(pid2==11) mass[1] = 0.511*0.001;
404	else if(pid2==13) mass[1] = 0.105658;
405	else if(pid2 ==22) mass[1] = 0;
406	else {
407	cout<<"calcPairMuPht, mass. pid2"<<endl;
408	exit(1);
409	}
410
411
412
413
414	for( int j= 0; j<2; j++){
415
416	double p = pt[j]/sin(2*atan(exp(-eta[j])));
417	double e = sqrt(pp+mass[j]mass[j]);
418	double px = p * sin ( 2atan(exp(-eta[j]))) cos(phi[j]);
419	double py = p * sin ( 2atan(exp(-eta[j]))) sin(phi[j]);
420	double pz = p * cos ( 2*atan(exp(-eta[j]))) ;
421
422	vpht[j].SetXYZT(px,py,pz,e);
423
424	}
425
426	vpair = vpht[0] + vpht[1];
427
428	res[0] = vpair.M();
429	res[1] = vpair.Eta();
430	res[2] = vpair.Phi();
431	res[3] = vpair.Pt();
432	res[4] = vpht[0].DeltaR(vpht[1]);
433
434
435	///cout<<"dot: "<< vpht[0].E() <<" "<< vpht[1].E()<< " "<< vpht[0].X() <<" "<< vpht[1].X()<<" "<< vpht[0].Y() <<" "<< vpht[1].Y() <<" "<< vpht[0].Z() <<" "<< vpht[1].Z()<<endl;
436
437
438	// res[5] = vpht[0].Dot(vpht[1])/vpht[0].P()/vpht[1].P();
439
440
441
442	}
443
444
445
446
447
448
449
450
451
452	////invariant mass of two obejcts given, pid each.
453	//void calcMass3Objects(float pt[3],float eta[3],float phi[3],int pid1,int pid2,int pid3,float res[]){
454	void calcMass3Objects(double pt[3],double eta[3],double phi[3],int pid1,int pid2,int pid3,double res[]){
455
456	TLorentzVector vpht[3];
457
458	TLorentzVector vpair;
459
460	float mass[3]={0,0,0};
461
462	if(pid1==11) mass[0] = 0.511*0.001;
463	else if(pid1==13) mass[0] = 0.105658;
464	else if(pid1==22) mass[0] = 0;
465	else{
466	cout<<"calcPairMuPht, mass.pid1 "<<endl;
467	exit(1);
468	}
469
470	if(pid2==11) mass[1] = 0.511*0.001;
471	else if(pid2==13) mass[1] = 0.105658;
472	else if(pid2 ==22) mass[1] = 0;
473	else {
474	cout<<"calcPairMuPht, mass. pid2"<<endl;
475	exit(1);
476	}
477
478	if(pid3==11) mass[2] = 0.511*0.001;
479	else if(pid3==13) mass[2] = 0.105658;
480	else if(pid3 ==22) mass[2] = 0;
481	else {
482	cout<<"calcPairMuPht, mass. pid3"<<endl;
483	exit(1);
484	}
485
486
487
488
489	for( int j= 0; j<3; j++){
490
491	float p = pt[j]/sin(2*atan(exp(-eta[j])));
492	float e = sqrt(pp+mass[j]mass[j]);
493	float px = p * sin ( 2atan(exp(-eta[j]))) cos(phi[j]);
494	float py = p * sin ( 2atan(exp(-eta[j]))) sin(phi[j]);
495	float pz = p * cos ( 2*atan(exp(-eta[j]))) ;
496
497	vpht[j].SetXYZT(px,py,pz,e);
498
499	}
500
501	vpair = vpht[0] + vpht[1] + vpht[2];
502
503	res[0] = vpair.M();
504	res[1] = vpair.Eta();
505	res[2] = vpair.Phi();
506	res[3] = vpair.Pt();
507	res[4] = vpht[0].DeltaR(vpht[1]);
508
509
510
511
512
513	}
514
515
516
517
518
519
520	void printRatioTwoHist(TH1F h1, TH1F h2){
521
522	int nbins = h1->GetNbinsX();
523
524	float xmin = h1->GetXaxis()->GetXmin() ;
525	float xmax = h1->GetXaxis()->GetXmax() ;
526
527	TH1F *hh = new TH1F("hh","eff",nbins,xmin,xmax);
528
529
530	for(int j=1; j<= nbins; j++){
531
532	float y1 = h1->GetBinContent(j);
533	float y2 = h2->GetBinContent(j);
534
535	float eff = 0;
536	float effErr = 0;
537	if(y1>0 && y2>0){
538	eff = y2/y1;
539	effErr = sqrt(eff*(1-eff)/y1);
540
541	hh->SetBinContent(j,eff);
542	hh->SetBinError(j,effErr);
543	}
544	}
545
546	// TCanvas *can0 = new TCanvas("can0","c000",200,10,550,500);
547	hh->Draw();
548
549
550	}
551
552	double etaTransformation( float EtaParticle , float Zvertex) {
553
554	//---Definitions
555	const float pi = 3.1415927;
556
557	//---Definitions for ECAL
558	const float R_ECAL = 136.5;
559	const float Z_Endcap = 328.0;
560	const float etaBarrelEndcap = 1.479;
561
562	//---ETA correction
563
564	float Theta = 0.0 ;
565	float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;
566
567	if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
568	if(Theta<0.0) Theta = Theta+pi ;
569	double ETA = - log(tan(0.5*Theta));
570
571	if( fabs(ETA) > etaBarrelEndcap )
572	{
573	float Zend = Z_Endcap ;
574	if(EtaParticle<0.0 ) Zend = -Zend ;
575	float Zlen = Zend - Zvertex ;
576	float RR = Zlen/sinh(EtaParticle);
577	Theta = atan(RR/Zend);
578	if(Theta<0.0) Theta = Theta+pi ;
579	ETA = - log(tan(0.5*Theta));
580	}
581	//---Return the result
582	return ETA;
583	//---end
584	}
585
586	////transform eta ( z, pho), to eta at ecal ( w.r.t 0,0,0,)
587	double ecalEta(double EtaParticle ,double Zvertex, double RhoVertex){
588
589
590	// const Double_t PI = 3.1415927;
591	double PI = acos(-1);
592
593	//---Definitions for ECAL
594	double R_ECAL = 136.5;
595	double Z_Endcap = 328.0;
596	double etaBarrelEndcap = 1.479;
597
598	if (EtaParticle!= 0.)
599	{
600	double Theta = 0.0 ;
601	double ZEcal = (R_ECAL-RhoVertex)*sinh(EtaParticle)+Zvertex;
602
603	if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
604	if(Theta<0.0) Theta = Theta+PI;
605
606	double ETA = - log(tan(0.5*Theta));
607
608	if( fabs(ETA) > etaBarrelEndcap )
609	{
610	double Zend = Z_Endcap ;
611	if(EtaParticle<0.0 ) Zend = -Zend ;
612	double Zlen = Zend - Zvertex ;
613	double RR = Zlen/sinh(EtaParticle);
614	Theta = atan((RR+RhoVertex)/Zend);
615	if(Theta<0.0) Theta = Theta+PI;
616	ETA = - log(tan(0.5*Theta));
617	}
618	return ETA;
619	}
620	else
621	{
622	return EtaParticle;
623	}
624	}
625
626
627	double ecalPhi(double phi,double x0,double y0){
628
629	//double R_ECAL = 136.5; ///cm
630	double r = 136.5;
631
632	double r0 = sqrt(x0x0 + y0y0);
633
634	if(r0<1E-5) return phi;
635
636	if( r0 >= r){
637	cout<<"warning. ecalPhi vtx outside ecal return input" << r0 <<" "<< r <<endl;
638	return phi;
639	}
640
641	double theta0 ;
642	if(fabs(y0)>0) theta0= y0/fabs(y0) * acos(x0/r0);
643	else theta0 = acos(x0/r0);
644
645	/// cout<<theta0<<" "<<phi<<endl;
646
647	double theta = phi + asin( r0/r *sin(theta0-phi));
648
649	//phinorm2(theta);
650	double PI = acos(-1);
651	while ( theta < -PI) theta += PI;
652	while ( theta > PI) theta -= PI;
653
654	return theta;
655
656
657	}