yangyong/zSelection/ecalGap.cc



// EB gap positions
double _barrelCGap[169][360][2];
double _barrelSGap[33][180][2];
double _barrelTGap[33][180][2];
double _barrelMGap[7][18][2];

// EE crystal existence and gap positions
bool   _endcapCrystal[100][100];
double _endcapCGap[2][7080][2];
double _endcapSGap[2][264][2];
double _endcapMGap[2][1][2];

double _endcapCGapxy[2][7080][2];
double _endcapSGapxy[2][264][2];
double _endcapMGapxy[2][1][2];


// Actual data for each instantiated object
unsigned _be,_hl;
double _e,_eta,_phi,_r9;
double _aC,_aS,_aM,_bC,_bS,_bM;
double _aT,_bT;
double _xyaC,_xyaS,_xyaM,_xybC,_xybS,_xybM;


void loadecalGapCoordinates(){
  
  TChain *feg = new TChain("ecalGap");

  if(isRealData){
    feg->Add("GapGeometryFile/ecalGap.cmssw425.GeometryDB.FT_R_42_V17A.root"); ///data 
  }else{
    feg->Add("GapGeometryFile/ecalGap.cmssw425.GeometryDB.START42_V12_new.root"); //MC 
  }
  

  feg->SetBranchAddress("_barrelCGap", _barrelCGap);
  feg->SetBranchAddress("_barrelSGap", _barrelSGap);
  feg->SetBranchAddress("_barrelTGap", _barrelTGap);
  feg->SetBranchAddress("_barrelMGap", _barrelMGap);
  feg->SetBranchAddress("_endcapCGap", _endcapCGap);
  feg->SetBranchAddress("_endcapSGap", _endcapSGap);
  feg->SetBranchAddress("_endcapMGap", _endcapMGap);
  feg->SetBranchAddress("_endcapCGapxy", _endcapCGapxy);
  feg->SetBranchAddress("_endcapSGapxy", _endcapSGapxy);
  feg->SetBranchAddress("_endcapMGapxy", _endcapMGapxy);

  
  feg->GetEntry(0);
    
  cout<<"ecalGap Positions loaded " << isRealData <<" "<< endl; 
  
}

const double _onePi = acos(-1);
const double _twoPi = 2.0 * acos(-1);


double dPhi(double f0, double f1) {
  double df(f0-f1);
  if(df> _onePi) df-=_twoPi;
  if(df<-_onePi) df+=_twoPi;
  return df;
}

double aPhi(double f0, double f1) {
  double af(0.5*(f0+f1));
  if(fabs(dPhi(af,f0))>0.5*_onePi) {
    if(af>=0.0) af-=_onePi;
    else        af+=_onePi;
  }
  
  assert(fabs(dPhi(af,f0))<0.5*_onePi);
  assert(fabs(dPhi(af,f1))<0.5*_onePi);
  
  return af;
}


double xZ(double eta, double phi){
  assert(_be==1);
  return asinh(cos(phi)/sinh(eta));
}

double yZ(double eta, double phi) {
  assert(_be==1);
  return asinh(sin(phi)/sinh(eta));
}


void getGapCoordinates(double eta, double phi){
  // Check constants have been set up
  //assert(_initialised);
  
  // Determine if EB or EE
  _be=(fabs(eta)<1.482?0:1);
  
  //  // Determine if high or low R9
  //   if(_be==0) _hl=(_r9>=0.94?0:1);
  //   else       _hl=(_r9>=0.95?0:1);
  
  // Coordinates relative to cracks
  double r2Min;
  if(_be==0) {
    
    r2Min=1.0e6;
    for(unsigned i(0);i<169;i++) {
      for(unsigned j(0);j<360;j++) {
        double de(eta-_barrelCGap[i][j][0]);
        double df(dPhi(phi,_barrelCGap[i][j][1]));
        double r2(de*de+df*df);
        
        if(r2<r2Min) {
          r2Min=r2;
          if(i>=84) {
            _aC= de;
            _bC=-df;
          } else {
            _aC=-de;
            _bC= df;
          }
        }
      }
    }
    
    r2Min=1.0e6;
    for(unsigned i(0);i<33;i++) {
      for(unsigned j(0);j<180;j++) {
        double de(eta-_barrelSGap[i][j][0]);
        double df(dPhi(phi,_barrelSGap[i][j][1]));
        double r2(de*de+df*df);
        
        if(r2<r2Min) {
          r2Min=r2;
          if(i>=16) {
            _aS= de;
            _bS=-df;
            
          } else {
            _aS=-de;
            _bS= df;
          }
        }
      }
    }
    
    r2Min=1.0e6;
    for(unsigned i(0);i<7;i++) {
      for(unsigned j(0);j<18;j++) {
        double de(eta-_barrelMGap[i][j][0]);
        double df(dPhi(phi,_barrelMGap[i][j][1]));
        double r2(de*de+df*df);
        
        if(r2<r2Min) {
          r2Min=r2;
          if(i>=3) {
            _aM= de;
            _bM=-df;
          } else {
            _aM=-de;
            _bM= df;
          }
        }
      }
    }
    
    
    r2Min=1.0e6;
    for(unsigned i(0);i<33;i++) {
      for(unsigned j(0);j<72;j++) {
        double de(eta-_barrelTGap[i][j][0]);
        double df(dPhi(phi,_barrelTGap[i][j][1]));
        double r2(de*de+df*df);
        
        if(r2<r2Min) {
          r2Min=r2;
          if(i>=16) {
            _aT= de;
            _bT=-df;
          } else {
            _aT=-de;
            _bT= df;
          }
        }
      }
    }
    
  } else {
    unsigned iz(eta>=0.0?0:1);
    double r[2]={xZ(eta,phi),yZ(eta,phi)};
    
    r2Min=1.0e6;
    for(unsigned i(0);i<7080;i++) {
      double dx(r[0]-_endcapCGap[iz][i][0]);
      double dy(r[1]-_endcapCGap[iz][i][1]);
      double r2(dx*dx+dy*dy);

      if(r2<r2Min) {
        r2Min=r2;
        if(r[0]>0.0) _aC= dx;
        else         _aC=-dx;
        if(r[1]>0.0) _bC= dy;
        else         _bC=-dy;
      }
    }
    
    r2Min=1.0e6;
    for(unsigned i(0);i<264;i++) {
      double dx(r[0]-_endcapSGap[iz][i][0]);
      double dy(r[1]-_endcapSGap[iz][i][1]);
      double r2(dx*dx+dy*dy);

      if(r2<r2Min) {
        r2Min=r2;
        if(r[0]>0.0) _aS= dx;
        else         _aS=-dx;
        if(r[1]>0.0) _bS= dy;
        else         _bS=-dy;
      }
    }
    
    r2Min=1.0e6;
    for(unsigned i(0);i<1;i++) {
      double dx(r[0]-_endcapMGap[iz][i][0]);
      double dy(r[1]-_endcapMGap[iz][i][1]);
      double r2(dx*dx+dy*dy);

      if(r2<r2Min) {
        r2Min=r2;
        if(iz==0) {_aM= dx;_bM= dy;}
        else      {_aM=-dx;_bM=-dy;}
      }
    }
  }
}


void getGapCoordinatesXY(double x, double y,double eta){
  // Check constants have been set up
  //assert(_initialised);
  
  // Determine if EB or EE
  ///_be=(fabs(eta)<1.482?0:1);
  
  if( fabs(eta) < 1.482) return;  //for endcap only
  
  unsigned iz(eta>=0.0?0:1);
  double r[2]={x,y};
  
  double r2Min=1.0e6;
  for(unsigned i(0);i<7080;i++) {
    double dx(r[0]-_endcapCGapxy[iz][i][0]);
    double dy(r[1]-_endcapCGapxy[iz][i][1]);
    double r2(dx*dx+dy*dy);

    if(r2<r2Min) {
      r2Min=r2;
      if(r[0]>0.0) _xyaC= dx;
      else         _xyaC=-dx;
      if(r[1]>0.0) _xybC= dy;
      else         _xybC=-dy;
    }
  }
    
  r2Min=1.0e6;
  for(unsigned i(0);i<264;i++) {
    double dx(r[0]-_endcapSGapxy[iz][i][0]);
    double dy(r[1]-_endcapSGapxy[iz][i][1]);
    double r2(dx*dx+dy*dy);

    if(r2<r2Min) {
      r2Min=r2;
      if(r[0]>0.0) _xyaS= dx;
      else         _xyaS=-dx;
      if(r[1]>0.0) _xybS= dy;
      else         _xybS=-dy;
    }
  }
    
  r2Min=1.0e6;
  for(unsigned i(0);i<1;i++) {
    double dx(r[0]-_endcapMGapxy[iz][i][0]);
    double dy(r[1]-_endcapMGapxy[iz][i][1]);
    double r2(dx*dx+dy*dy);

    if(r2<r2Min) {
      r2Min=r2;
      if(iz==0) {_xyaM= dx;_xybM= dy;}
      else      {_xyaM=-dx;_xybM=-dy;}
    }
  }
  
}

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			// EB gap positions
4			double _barrelCGap[169][360][2];
5			double _barrelSGap[33][180][2];
6			double _barrelTGap[33][180][2];
7			double _barrelMGap[7][18][2];
8
9			// EE crystal existence and gap positions
10			bool _endcapCrystal[100][100];
11			double _endcapCGap[2][7080][2];
12			double _endcapSGap[2][264][2];
13			double _endcapMGap[2][1][2];
14
15			double _endcapCGapxy[2][7080][2];
16			double _endcapSGapxy[2][264][2];
17			double _endcapMGapxy[2][1][2];
18
19
20			// Actual data for each instantiated object
21			unsigned _be,_hl;
22			double _e,_eta,_phi,_r9;
23			double _aC,_aS,_aM,_bC,_bS,_bM;
24			double _aT,_bT;
25			double _xyaC,_xyaS,_xyaM,_xybC,_xybS,_xybM;
26
27
28			void loadecalGapCoordinates(){
29
30			TChain *feg = new TChain("ecalGap");
31
32			if(isRealData){
33			feg->Add("GapGeometryFile/ecalGap.cmssw425.GeometryDB.FT_R_42_V17A.root"); ///data
34			}else{
35			feg->Add("GapGeometryFile/ecalGap.cmssw425.GeometryDB.START42_V12_new.root"); //MC
36			}
37
38
39			feg->SetBranchAddress("_barrelCGap", _barrelCGap);
40			feg->SetBranchAddress("_barrelSGap", _barrelSGap);
41			feg->SetBranchAddress("_barrelTGap", _barrelTGap);
42			feg->SetBranchAddress("_barrelMGap", _barrelMGap);
43			feg->SetBranchAddress("_endcapCGap", _endcapCGap);
44			feg->SetBranchAddress("_endcapSGap", _endcapSGap);
45			feg->SetBranchAddress("_endcapMGap", _endcapMGap);
46			feg->SetBranchAddress("_endcapCGapxy", _endcapCGapxy);
47			feg->SetBranchAddress("_endcapSGapxy", _endcapSGapxy);
48			feg->SetBranchAddress("_endcapMGapxy", _endcapMGapxy);
49
50
51			feg->GetEntry(0);
52
53			cout<<"ecalGap Positions loaded " << isRealData <<" "<< endl;
54
55			}
56
57			const double _onePi = acos(-1);
58			const double _twoPi = 2.0 * acos(-1);
59
60
61			double dPhi(double f0, double f1) {
62			double df(f0-f1);
63			if(df> _onePi) df-=_twoPi;
64			if(df<-_onePi) df+=_twoPi;
65			return df;
66			}
67
68			double aPhi(double f0, double f1) {
69			double af(0.5*(f0+f1));
70			if(fabs(dPhi(af,f0))>0.5*_onePi) {
71			if(af>=0.0) af-=_onePi;
72			else af+=_onePi;
73			}
74
75			assert(fabs(dPhi(af,f0))<0.5*_onePi);
76			assert(fabs(dPhi(af,f1))<0.5*_onePi);
77
78			return af;
79			}
80
81
82			double xZ(double eta, double phi){
83			assert(_be==1);
84			return asinh(cos(phi)/sinh(eta));
85			}
86
87			double yZ(double eta, double phi) {
88			assert(_be==1);
89			return asinh(sin(phi)/sinh(eta));
90			}
91
92
93
94			void getGapCoordinates(double eta, double phi){
95			// Check constants have been set up
96			//assert(_initialised);
97
98			// Determine if EB or EE
99			_be=(fabs(eta)<1.482?0:1);
100
101			// // Determine if high or low R9
102			// if(_be==0) _hl=(_r9>=0.94?0:1);
103			// else _hl=(_r9>=0.95?0:1);
104
105			// Coordinates relative to cracks
106			double r2Min;
107			if(_be==0) {
108
109			r2Min=1.0e6;
110			for(unsigned i(0);i<169;i++) {
111			for(unsigned j(0);j<360;j++) {
112			double de(eta-_barrelCGap[i][j][0]);
113			double df(dPhi(phi,_barrelCGap[i][j][1]));
114			double r2(dede+dfdf);
115
116			if(r2<r2Min) {
117			r2Min=r2;
118			if(i>=84) {
119			_aC= de;
120			_bC=-df;
121			} else {
122			_aC=-de;
123			_bC= df;
124			}
125			}
126			}
127			}
128
129			r2Min=1.0e6;
130			for(unsigned i(0);i<33;i++) {
131			for(unsigned j(0);j<180;j++) {
132			double de(eta-_barrelSGap[i][j][0]);
133			double df(dPhi(phi,_barrelSGap[i][j][1]));
134			double r2(dede+dfdf);
135
136			if(r2<r2Min) {
137			r2Min=r2;
138			if(i>=16) {
139			_aS= de;
140			_bS=-df;
141
142			} else {
143			_aS=-de;
144			_bS= df;
145			}
146			}
147			}
148			}
149
150			r2Min=1.0e6;
151			for(unsigned i(0);i<7;i++) {
152			for(unsigned j(0);j<18;j++) {
153			double de(eta-_barrelMGap[i][j][0]);
154			double df(dPhi(phi,_barrelMGap[i][j][1]));
155			double r2(dede+dfdf);
156
157			if(r2<r2Min) {
158			r2Min=r2;
159			if(i>=3) {
160			_aM= de;
161			_bM=-df;
162			} else {
163			_aM=-de;
164			_bM= df;
165			}
166			}
167			}
168			}
169
170
171			r2Min=1.0e6;
172			for(unsigned i(0);i<33;i++) {
173			for(unsigned j(0);j<72;j++) {
174			double de(eta-_barrelTGap[i][j][0]);
175			double df(dPhi(phi,_barrelTGap[i][j][1]));
176			double r2(dede+dfdf);
177
178			if(r2<r2Min) {
179			r2Min=r2;
180			if(i>=16) {
181			_aT= de;
182			_bT=-df;
183			} else {
184			_aT=-de;
185			_bT= df;
186			}
187			}
188			}
189			}
190
191			} else {
192			unsigned iz(eta>=0.0?0:1);
193			double r[2]={xZ(eta,phi),yZ(eta,phi)};
194
195			r2Min=1.0e6;
196			for(unsigned i(0);i<7080;i++) {
197			double dx(r[0]-_endcapCGap[iz][i][0]);
198			double dy(r[1]-_endcapCGap[iz][i][1]);
199			double r2(dxdx+dydy);
200
201			if(r2<r2Min) {
202			r2Min=r2;
203			if(r[0]>0.0) _aC= dx;
204			else _aC=-dx;
205			if(r[1]>0.0) _bC= dy;
206			else _bC=-dy;
207			}
208			}
209
210			r2Min=1.0e6;
211			for(unsigned i(0);i<264;i++) {
212			double dx(r[0]-_endcapSGap[iz][i][0]);
213			double dy(r[1]-_endcapSGap[iz][i][1]);
214			double r2(dxdx+dydy);
215
216			if(r2<r2Min) {
217			r2Min=r2;
218			if(r[0]>0.0) _aS= dx;
219			else _aS=-dx;
220			if(r[1]>0.0) _bS= dy;
221			else _bS=-dy;
222			}
223			}
224
225			r2Min=1.0e6;
226			for(unsigned i(0);i<1;i++) {
227			double dx(r[0]-_endcapMGap[iz][i][0]);
228			double dy(r[1]-_endcapMGap[iz][i][1]);
229			double r2(dxdx+dydy);
230
231			if(r2<r2Min) {
232			r2Min=r2;
233			if(iz==0) {_aM= dx;_bM= dy;}
234			else {_aM=-dx;_bM=-dy;}
235			}
236			}
237			}
238			}
239
240
241
242			void getGapCoordinatesXY(double x, double y,double eta){
243			// Check constants have been set up
244			//assert(_initialised);
245
246			// Determine if EB or EE
247			///_be=(fabs(eta)<1.482?0:1);
248
249			if( fabs(eta) < 1.482) return; //for endcap only
250
251			unsigned iz(eta>=0.0?0:1);
252			double r[2]={x,y};
253
254			double r2Min=1.0e6;
255			for(unsigned i(0);i<7080;i++) {
256			double dx(r[0]-_endcapCGapxy[iz][i][0]);
257			double dy(r[1]-_endcapCGapxy[iz][i][1]);
258			double r2(dxdx+dydy);
259
260			if(r2<r2Min) {
261			r2Min=r2;
262			if(r[0]>0.0) _xyaC= dx;
263			else _xyaC=-dx;
264			if(r[1]>0.0) _xybC= dy;
265			else _xybC=-dy;
266			}
267			}
268
269			r2Min=1.0e6;
270			for(unsigned i(0);i<264;i++) {
271			double dx(r[0]-_endcapSGapxy[iz][i][0]);
272			double dy(r[1]-_endcapSGapxy[iz][i][1]);
273			double r2(dxdx+dydy);
274
275			if(r2<r2Min) {
276			r2Min=r2;
277			if(r[0]>0.0) _xyaS= dx;
278			else _xyaS=-dx;
279			if(r[1]>0.0) _xybS= dy;
280			else _xybS=-dy;
281			}
282			}
283
284			r2Min=1.0e6;
285			for(unsigned i(0);i<1;i++) {
286			double dx(r[0]-_endcapMGapxy[iz][i][0]);
287			double dy(r[1]-_endcapMGapxy[iz][i][1]);
288			double r2(dxdx+dydy);
289
290			if(r2<r2Min) {
291			r2Min=r2;
292			if(iz==0) {_xyaM= dx;_xybM= dy;}
293			else {_xyaM=-dx;_xybM=-dy;}
294			}
295			}
296
297			}
298