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
#	User	Rev	Content
1	yangyong	1.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			}