yangyong/combineICEE/common_functions.cc



////152, not xyz of xtal ,read from 167 root file
void get_xyzEBrechits(){
  TChain *ch = new TChain("Analysis");
  //  ch->Add("/afs/cern.ch/user/y/yangyong/w0/cmssw167.xyzEB.root");
  ch->Add("/afs/cern.ch/user/y/yangyong/w0/xyzECAL.GR_R_42_V20.root");
  //ch->Add("/uscms_data/d2/yongy/calib/backup/cmssw217.xyzECAL.root");
  
  
  ch->SetBranchAddress("xEBAll",xEBAll);
  ch->SetBranchAddress("yEBAll",yEBAll);
  ch->SetBranchAddress("zEBAll",zEBAll);
  ch->SetBranchAddress("etaEBAll",etaEBAll);
  ch->SetBranchAddress("phiEBAll",phiEBAll);

  
  ch->SetBranchAddress("dxEBAll",dxEBAll);
  ch->SetBranchAddress("dyEBAll",dyEBAll);
  ch->SetBranchAddress("dzEBAll",dzEBAll);
  
  
  ch->SetBranchAddress("xEEAll",xEEAll);
  ch->SetBranchAddress("yEEAll",yEEAll);
  ch->SetBranchAddress("zEEAll",zEEAll);
  
  ch->SetBranchAddress("dxEEAll",dxEEAll);
  ch->SetBranchAddress("dyEEAll",dyEEAll);
  ch->SetBranchAddress("dzEEAll",dzEEAll);
  
  
  ch->SetBranchAddress("etaEEAll",etaEEAll);
  ch->SetBranchAddress("phiEEAll",phiEEAll);
  
  
  ch->GetEntry(0);

  ch->Delete();
  
  cout<<"got xyzEcal.."<<endl; 
  
  
}


/////////////////////////////////////////////////////////////////////////////////////////////////

void convxtalid(Int_t &nphi,Int_t &neta)
{
  // Changed to what Yong's convention; output will give just two indices
  // phi is unchanged; only eta now runs from
  //
  // 03/01/2008 changed to the new definition in CMSSW. The output is still the same...
  // Barrel only
  // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
  // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.

     if(neta > 0) neta -= 1;
     if(nphi > 359) nphi=nphi-360;

     // final check
   if(nphi >359 || nphi <0 || neta< -85 || neta > 84)
     {
       cout <<" output not in range: "<<  nphi <<  " " << neta <<  " " <<endl;
    exit(1);
     }
} //end of convxtalid


// Calculate the distance in xtals taking into account possibly different sides
// change to coincide with yongs definition
Int_t diff_neta(Int_t neta1, Int_t neta2){
    Int_t mdiff;
    mdiff=abs(neta1-neta2);
    return mdiff;
}
 
// Calculate the absolute distance in xtals taking into account the periodicity of the Barrel
Int_t diff_nphi(Int_t nphi1,Int_t nphi2) {
  Int_t mdiff;
   mdiff=abs(nphi1-nphi2);
   if (mdiff > (360-abs(nphi1-nphi2))) mdiff=(360-abs(nphi1-nphi2));
   return mdiff;
}

// Calculate the distance in xtals taking into account possibly different sides
// Then the distance would be from the 1st to the 2nd argument
// _s means that it gives the sign; the only difference from the above !
// also changed to coincide with Yong's definition
Int_t diff_neta_s(Int_t neta1, Int_t neta2){
  Int_t mdiff;
  mdiff=(neta1-neta2);
  return mdiff;
}
 
// Calculate the distance in xtals taking into account the periodicity of the Barrel
Int_t diff_nphi_s(Int_t nphi1,Int_t nphi2) {
   Int_t mdiff;
   if(abs(nphi1-nphi2) < (360-abs(nphi1-nphi2))) {
     mdiff=nphi1-nphi2;
   }
   else {
   mdiff=360-abs(nphi1-nphi2);
   if(nphi1>nphi2) mdiff=-mdiff;
   }
   return mdiff;
}


void sortPtvector(vector<float> pt, vector<int> &ind){
  
  ind.clear();

  int n = int( pt.size());

  vector<float> pt0; 
  for(int j=0; j< n ; j++){
    pt0.push_back(pt[j]);
  }
  
  sort(pt.begin(),pt.end());
  
  for(int j = n -1; j>=0; j--){
    
    float dptmin = 0.1; 
    int jmin = -1; 
    for(int k=0; k< n; k++){
      
      float dpt = fabs( pt[j] - pt0[k]); 
      vector<int>::iterator it = find( ind.begin(),ind.end(),k);
      if( it != ind.end()) continue; 
      
      if( dpt < dptmin){
        dptmin = dpt; 
        jmin = k; 
      }

    }
    ind.push_back(jmin);
    
  }
  
  if(int(ind.size()) != n){
    cout<<"erorr_sortPtvector.."<<endl; 
    exit(1);
  }
  
}


float  calculateS4(int nn, float en[],int ieta[], int iphi[]){
  ///shower-shape
  Float_t S4_0 =0.;      Float_t S4_1 =0.; 
  Float_t S4_2 =0.;      Float_t S4_3 =0.; 
  
  
  for(int n=0; n< nn; n++){
    float energy = en[n]; 
    int deta = diff_neta_s(ieta[0],ieta[n]);
    int dphi = diff_nphi_s(iphi[0],iphi[n]);
    if( abs(deta) <=1 && abs(dphi) <= 1){
      if( dphi <= 0 && deta <= 0) S4_0+=energy;
      if( dphi >= 0 && deta <= 0) S4_1+=energy;
      if( dphi <= 0 && deta >= 0) S4_2+=energy;
      if( dphi >= 0 && deta >= 0) S4_3+=energy;
    }
  }
  Float_t S4max=S4_0;
  if(S4_1 > S4max) S4max=S4_1;
  if(S4_2 > S4max) S4max=S4_2;
  if(S4_3 > S4max) S4max=S4_3;
  
  return S4max; 
  
}

void calculateS4_v1(int nn, float en[],int ieta[], int iphi[],float res[]){
  ///shower-shape
  Float_t S4_0 =0.;      Float_t S4_1 =0.; 
  Float_t S4_2 =0.;      Float_t S4_3 =0.; 

  float S9 = 0; 
    
  for(int n=0; n< nn; n++){
    float energy = en[n]; 
    int deta = diff_neta_s(ieta[0],ieta[n]);
    int dphi = diff_nphi_s(iphi[0],iphi[n]);
    if( abs(deta) <=1 && abs(dphi) <= 1){

      S9 += energy; 
      
      if( dphi <= 0 && deta <= 0) S4_0+=energy;
      if( dphi >= 0 && deta <= 0) S4_1+=energy;
      if( dphi <= 0 && deta >= 0) S4_2+=energy;
      if( dphi >= 0 && deta >= 0) S4_3+=energy;
    }
  }
  Float_t S4max=S4_0;
  if(S4_1 > S4max) S4max=S4_1;
  if(S4_2 > S4max) S4max=S4_2;
  if(S4_3 > S4max) S4max=S4_3;
  
  ////  return S4max; 
  
  res[0] = S4max; 
  res[1] = S9; 
  

}


///input ieta iphi originall 

void calculateS4_v2(int nn, float en[],int ieta[], int iphi[],float res[]){
  ///shower-shape
  Float_t S4_0 =0.;      Float_t S4_1 =0.; 
  Float_t S4_2 =0.;      Float_t S4_3 =0.; 

  float S9 = 0; 

  int ieta_seed = ieta[0];
  int iphi_seed = iphi[0];
  
  convxtalid(iphi_seed,ieta_seed);
  

  for(int n=0; n< nn; n++){
    float energy = en[n]; 

    int eta = ieta[n];
    int phi = iphi[n];
    convxtalid(phi,eta);
    
    int deta = diff_neta_s(ieta_seed,eta);
    int dphi = diff_nphi_s(iphi_seed,phi);
    
    if( abs(deta) <=1 && abs(dphi) <= 1){

      S9 += energy; 
      
      if( dphi <= 0 && deta <= 0) S4_0+=energy;
      if( dphi >= 0 && deta <= 0) S4_1+=energy;
      if( dphi <= 0 && deta >= 0) S4_2+=energy;
      if( dphi >= 0 && deta >= 0) S4_3+=energy;
    }
  }
  Float_t S4max=S4_0;
  if(S4_1 > S4max) S4max=S4_1;
  if(S4_2 > S4max) S4max=S4_2;
  if(S4_3 > S4max) S4max=S4_3;
  
  ////  return S4max; 
  
  res[0] = S4max; 
  res[1] = S9; 
  

}
Revision:	1.1
Committed:	Tue Apr 10 20:03:11 2012 UTC (13 years, 1 month ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	pi0etaee_laser20111122_eefloatalpha
Log Message:	LaserTag: EcalLaserAPDPNRatios_data_20111122_158851_180363 and alpha Tag ( for Endcap Only) EcalLaserAlphas_lto420-620_progr_data_20111122
#	User	Rev	Content
1	yangyong	1.1
2
3
4
5			////152, not xyz of xtal ,read from 167 root file
6			void get_xyzEBrechits(){
7			TChain *ch = new TChain("Analysis");
8			// ch->Add("/afs/cern.ch/user/y/yangyong/w0/cmssw167.xyzEB.root");
9			ch->Add("/afs/cern.ch/user/y/yangyong/w0/xyzECAL.GR_R_42_V20.root");
10			//ch->Add("/uscms_data/d2/yongy/calib/backup/cmssw217.xyzECAL.root");
11
12
13
14			ch->SetBranchAddress("xEBAll",xEBAll);
15			ch->SetBranchAddress("yEBAll",yEBAll);
16			ch->SetBranchAddress("zEBAll",zEBAll);
17			ch->SetBranchAddress("etaEBAll",etaEBAll);
18			ch->SetBranchAddress("phiEBAll",phiEBAll);
19
20
21			ch->SetBranchAddress("dxEBAll",dxEBAll);
22			ch->SetBranchAddress("dyEBAll",dyEBAll);
23			ch->SetBranchAddress("dzEBAll",dzEBAll);
24
25
26			ch->SetBranchAddress("xEEAll",xEEAll);
27			ch->SetBranchAddress("yEEAll",yEEAll);
28			ch->SetBranchAddress("zEEAll",zEEAll);
29
30			ch->SetBranchAddress("dxEEAll",dxEEAll);
31			ch->SetBranchAddress("dyEEAll",dyEEAll);
32			ch->SetBranchAddress("dzEEAll",dzEEAll);
33
34
35
36			ch->SetBranchAddress("etaEEAll",etaEEAll);
37			ch->SetBranchAddress("phiEEAll",phiEEAll);
38
39
40
41
42
43			ch->GetEntry(0);
44
45			ch->Delete();
46
47			cout<<"got xyzEcal.."<<endl;
48
49
50
51			}
52
53
54			/////////////////////////////////////////////////////////////////////////////////////////////////
55
56			void convxtalid(Int_t &nphi,Int_t &neta)
57			{
58			// Changed to what Yong's convention; output will give just two indices
59			// phi is unchanged; only eta now runs from
60			//
61			// 03/01/2008 changed to the new definition in CMSSW. The output is still the same...
62			// Barrel only
63			// Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
64			// neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
65
66			if(neta > 0) neta -= 1;
67			if(nphi > 359) nphi=nphi-360;
68
69			// final check
70			if(nphi >359 \|\| nphi <0 \|\| neta< -85 \|\| neta > 84)
71			{
72			cout <<" output not in range: "<< nphi << " " << neta << " " <<endl;
73			exit(1);
74			}
75			} //end of convxtalid
76
77
78			// Calculate the distance in xtals taking into account possibly different sides
79			// change to coincide with yongs definition
80			Int_t diff_neta(Int_t neta1, Int_t neta2){
81			Int_t mdiff;
82			mdiff=abs(neta1-neta2);
83			return mdiff;
84			}
85
86			// Calculate the absolute distance in xtals taking into account the periodicity of the Barrel
87			Int_t diff_nphi(Int_t nphi1,Int_t nphi2) {
88			Int_t mdiff;
89			mdiff=abs(nphi1-nphi2);
90			if (mdiff > (360-abs(nphi1-nphi2))) mdiff=(360-abs(nphi1-nphi2));
91			return mdiff;
92			}
93
94			// Calculate the distance in xtals taking into account possibly different sides
95			// Then the distance would be from the 1st to the 2nd argument
96			// _s means that it gives the sign; the only difference from the above !
97			// also changed to coincide with Yong's definition
98			Int_t diff_neta_s(Int_t neta1, Int_t neta2){
99			Int_t mdiff;
100			mdiff=(neta1-neta2);
101			return mdiff;
102			}
103
104			// Calculate the distance in xtals taking into account the periodicity of the Barrel
105			Int_t diff_nphi_s(Int_t nphi1,Int_t nphi2) {
106			Int_t mdiff;
107			if(abs(nphi1-nphi2) < (360-abs(nphi1-nphi2))) {
108			mdiff=nphi1-nphi2;
109			}
110			else {
111			mdiff=360-abs(nphi1-nphi2);
112			if(nphi1>nphi2) mdiff=-mdiff;
113			}
114			return mdiff;
115			}
116
117
118
119			void sortPtvector(vector<float> pt, vector<int> &ind){
120
121			ind.clear();
122
123			int n = int( pt.size());
124
125			vector<float> pt0;
126			for(int j=0; j< n ; j++){
127			pt0.push_back(pt[j]);
128			}
129
130			sort(pt.begin(),pt.end());
131
132			for(int j = n -1; j>=0; j--){
133
134			float dptmin = 0.1;
135			int jmin = -1;
136			for(int k=0; k< n; k++){
137
138			float dpt = fabs( pt[j] - pt0[k]);
139			vector<int>::iterator it = find( ind.begin(),ind.end(),k);
140			if( it != ind.end()) continue;
141
142			if( dpt < dptmin){
143			dptmin = dpt;
144			jmin = k;
145			}
146
147			}
148			ind.push_back(jmin);
149
150			}
151
152			if(int(ind.size()) != n){
153			cout<<"erorr_sortPtvector.."<<endl;
154			exit(1);
155			}
156
157			}
158
159
160			float calculateS4(int nn, float en[],int ieta[], int iphi[]){
161			///shower-shape
162			Float_t S4_0 =0.; Float_t S4_1 =0.;
163			Float_t S4_2 =0.; Float_t S4_3 =0.;
164
165
166			for(int n=0; n< nn; n++){
167			float energy = en[n];
168			int deta = diff_neta_s(ieta[0],ieta[n]);
169			int dphi = diff_nphi_s(iphi[0],iphi[n]);
170			if( abs(deta) <=1 && abs(dphi) <= 1){
171			if( dphi <= 0 && deta <= 0) S4_0+=energy;
172			if( dphi >= 0 && deta <= 0) S4_1+=energy;
173			if( dphi <= 0 && deta >= 0) S4_2+=energy;
174			if( dphi >= 0 && deta >= 0) S4_3+=energy;
175			}
176			}
177			Float_t S4max=S4_0;
178			if(S4_1 > S4max) S4max=S4_1;
179			if(S4_2 > S4max) S4max=S4_2;
180			if(S4_3 > S4max) S4max=S4_3;
181
182			return S4max;
183
184			}
185
186			void calculateS4_v1(int nn, float en[],int ieta[], int iphi[],float res[]){
187			///shower-shape
188			Float_t S4_0 =0.; Float_t S4_1 =0.;
189			Float_t S4_2 =0.; Float_t S4_3 =0.;
190
191			float S9 = 0;
192
193			for(int n=0; n< nn; n++){
194			float energy = en[n];
195			int deta = diff_neta_s(ieta[0],ieta[n]);
196			int dphi = diff_nphi_s(iphi[0],iphi[n]);
197			if( abs(deta) <=1 && abs(dphi) <= 1){
198
199			S9 += energy;
200
201			if( dphi <= 0 && deta <= 0) S4_0+=energy;
202			if( dphi >= 0 && deta <= 0) S4_1+=energy;
203			if( dphi <= 0 && deta >= 0) S4_2+=energy;
204			if( dphi >= 0 && deta >= 0) S4_3+=energy;
205			}
206			}
207			Float_t S4max=S4_0;
208			if(S4_1 > S4max) S4max=S4_1;
209			if(S4_2 > S4max) S4max=S4_2;
210			if(S4_3 > S4max) S4max=S4_3;
211
212			//// return S4max;
213
214			res[0] = S4max;
215			res[1] = S9;
216
217
218			}
219
220
221			///input ieta iphi originall
222
223			void calculateS4_v2(int nn, float en[],int ieta[], int iphi[],float res[]){
224			///shower-shape
225			Float_t S4_0 =0.; Float_t S4_1 =0.;
226			Float_t S4_2 =0.; Float_t S4_3 =0.;
227
228			float S9 = 0;
229
230			int ieta_seed = ieta[0];
231			int iphi_seed = iphi[0];
232
233			convxtalid(iphi_seed,ieta_seed);
234
235
236
237			for(int n=0; n< nn; n++){
238			float energy = en[n];
239
240			int eta = ieta[n];
241			int phi = iphi[n];
242			convxtalid(phi,eta);
243
244			int deta = diff_neta_s(ieta_seed,eta);
245			int dphi = diff_nphi_s(iphi_seed,phi);
246
247			if( abs(deta) <=1 && abs(dphi) <= 1){
248
249			S9 += energy;
250
251			if( dphi <= 0 && deta <= 0) S4_0+=energy;
252			if( dphi >= 0 && deta <= 0) S4_1+=energy;
253			if( dphi <= 0 && deta >= 0) S4_2+=energy;
254			if( dphi >= 0 && deta >= 0) S4_3+=energy;
255			}
256			}
257			Float_t S4max=S4_0;
258			if(S4_1 > S4max) S4max=S4_1;
259			if(S4_2 > S4max) S4max=S4_2;
260			if(S4_3 > S4max) S4max=S4_3;
261
262			//// return S4max;
263
264			res[0] = S4max;
265			res[1] = S9;
266
267
268			}