BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
DedxCorrecSvc.cxx
Go to the documentation of this file.
1#include "DedxCorrecSvc/DedxCorrecSvc.h"
2#include "DedxCorrecSvc/DedxCorrParameters.h"
3#include "GaudiKernel/Kernel.h"
4#include "GaudiKernel/IInterface.h"
5#include "GaudiKernel/StatusCode.h"
6//#include "GaudiKernel/ISvcFactory.h"
7#include "GaudiKernel/SvcFactory.h"
8#include "GaudiKernel/MsgStream.h"
9#include "GaudiKernel/SmartDataPtr.h"
10#include "GaudiKernel/DataSvc.h"
11
12#include "GaudiKernel/IIncidentSvc.h"
13#include "GaudiKernel/Incident.h"
14#include "GaudiKernel/IIncidentListener.h"
15
16#include "Identifier/Identifier.h"
17#include "Identifier/MdcID.h"
18#include "CLHEP/Vector/ThreeVector.h"
19#include "CalibData/Dedx/DedxCalibData.h"
20#include "CalibData/CalibModel.h"
21#include <math.h>
22#include <vector>
23#include <iostream>
24#include <fstream>
25#include <string>
26
27#include "TFile.h"
28#include "TTree.h"
29#include "TLeafD.h"
30
31using namespace std;
32using CLHEP::Hep3Vector;
33
34//static SvcFactory<DedxCorrecSvc> s_factory;
35//const ISvcFactory& DedxCorrecSvcFactory = s_factory;
36//double Iner_Stepdoca = (Iner_DocaMax-Iner_DocaMin)/NumSlicesDoca;
37//double Out_Stepdoca = (Out_DocaMax-Out_DocaMin)/NumSlicesDoca;
38
39DECLARE_COMPONENT(DedxCorrecSvc)
40DedxCorrecSvc::DedxCorrecSvc( const string& name, ISvcLocator* svcloc) :
41 geosvc(0),base_class(name, svcloc) {
42 declareProperty("Run",m_run=1);
43 declareProperty("init",m_init=1);
44 declareProperty("par_flag",m_par_flag=0);
45 declareProperty("Debug",m_debug=false);
46 declareProperty("DebugI",m_debug_i=1);
47 declareProperty("DebugJ",m_debug_j=1);
48 m_initConstFlg = false;
49 // cout<<"DedxCorrecSvc::DedxCorrecSvc"<<endl;
50 }
51
54
55
56StatusCode DedxCorrecSvc::initialize(){
57 // cout<<"DedxCorrecSvc::initialize"<<endl;
58 MsgStream log(messageService(), name());
59 log << MSG::INFO << name() << "DedxCorrecSvc::initialize()" << endreq;
60
61 StatusCode sc = Service::initialize();
62 if( sc.isFailure() ) return sc;
63
64 IIncidentSvc* incsvc;
65 sc = service("IncidentSvc", incsvc);
66 int priority = 100;
67 if( sc.isSuccess() ){
68 //incsvc -> addListener(this, "BeginEvent", priority);
69 incsvc -> addListener(this, "NewRun", priority);
70 }
71
72 StatusCode scgeo = service("MdcGeomSvc", geosvc);
73 if(scgeo.isFailure() ) {
74 log << MSG::ERROR << "GeoSvc failing!!!!!!!" << endreq;
75 return scgeo;
76 }
77
78 StatusCode scmgn = service ("MagneticFieldSvc",m_pIMF);
79 if(scmgn!=StatusCode::SUCCESS) {
80 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
81 }
82
83 return StatusCode::SUCCESS;
84}
85
86StatusCode DedxCorrecSvc::finalize(){
87 // cout<<"DedxCorrecSvc::finalize()"<<endl;
88 MsgStream log(messageService(), name());
89 log << MSG::INFO << name() << "DedxCorrecSvc::finalize()" << endreq;
90 return StatusCode::SUCCESS;
91}
92
93void DedxCorrecSvc::handle(const Incident& inc){
94 // cout<<"DedxCorrecSvc::handle"<<endl;
95 MsgStream log( messageService(), name() );
96 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
97
98 if ( inc.type() == "NewRun" ){
99 log << MSG::DEBUG << "New Run" << endreq;
100
101 if( ! m_initConstFlg )
102 {
103 if( m_init == 0 ) { init_param(); }
104 else if( m_init ==1 ) { init_param_svc(); }
105 /* if( ! init_param() ){
106 log << MSG::ERROR
107 << "can not initilize Mdc Calib Constants" << endreq;
108 }*/
109 }
110 }
111}
112
113double DedxCorrecSvc::f_larcos(double x, int nbin) {
114 if(nbin!=80) return x; // temporally for only nbin = 80
115 double m_absx(fabs(x));
116 double m_charge(x/m_absx);
117 if(m_absx>0.925 && m_absx<=0.950) return 0.9283*m_charge; // these numbers are from the mean values
118 if(m_absx>0.900 && m_absx<=0.925) return 0.9115*m_charge; // of each bin in cos(theta) distribution
119 if(m_absx>0.875 && m_absx<=0.900) return 0.8877*m_charge;
120 if(m_absx>0.850 && m_absx<=0.875) return 0.8634*m_charge;
121 if(m_absx>0.825 && m_absx<=0.850) return 0.8460*m_charge;
122 if(m_absx>0.800 && m_absx<=0.825) return 0.8089*m_charge;
123}
124
125void
126DedxCorrecSvc::init_param_svc() {
127 // cout<<"DedxCorrecSvc::init_param_svc()"<<endl;
128 MsgStream log(messageService(), name());
129 IDataProviderSvc* pCalibDataSvc;
130 StatusCode sc = service("CalibDataSvc", pCalibDataSvc, true);
131 if ( !sc.isSuccess() ) {
132 log << MSG::ERROR
133 << "Could not get IDataProviderSvc interface of CalibXmlCnvSvc"
134 << endreq;
135 return ;
136 } else {
137 log << MSG::DEBUG
138 << "Retrieved IDataProviderSvc interface of CalibXmlCnvSvc"
139 << endreq;
140 }
141
142 Iner_Stepdoca = (Iner_DocaMax-Iner_DocaMin)/NumSlicesDoca;
143 Out_Stepdoca = (Out_DocaMax-Out_DocaMin)/NumSlicesDoca;
144 //cout<<"Out_DocaMax: "<<Out_DocaMax<<" Out_DocaMin: "<<Out_DocaMin<<" Out_Stepdoca : "<<Out_Stepdoca<<" NumSlicesDoca : "<<NumSlicesDoca<<endl;
145 std::string fullPath = "/Calib/DedxCal";
146
147 SmartDataPtr<CalibData::DedxCalibData> test(pCalibDataSvc, fullPath);
148
149 //rung par----------------------------------------------------------
150 N_run = test->getrunNO(); // not the real runNO, a few events blocks may in a run
151 int run_temp = 0;
152 cout<<"N_run = "<<N_run<<endl;
153 for(int j=0;j<100000;j++)
154 {
155 if(j<N_run)
156 {
157 for(int i=0;i<6;i++)
158 {
159 m_rung[i][j] = test->getrung(i,j);
160 if(i==2 && m_rung[2][j]>run_temp) run_temp = m_rung[2][j];
161 }
162 }
163 else
164 {
165 run_temp++;
166 m_rung[0][j] = 1;
167 m_rung[1][j] = 550;
168 m_rung[2][j] = run_temp;
169 m_rung[3][j] = 26.8;
170 m_rung[4][j] = 0;
171 m_rung[5][j] = 1000000000;
172 }
173 }
174
175 //for(int i=0;i<4;i++)
176 //{
177 // for(int j=0;j<10000;j++)
178 // {
179 // std::cout<<"rung["<<i<<"]["<<j<<"]= "<<m_rung[i][j]<<std::endl;
180 // }
181 //}
182
183 //TFile* fruncalib=new TFile("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/rung_calib/first_calib/DedxConst.root");
184 /*TFile* fruncalib=new TFile("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/rung_calib/fourth_calib/temp/DedxConst.root");
185
186 double rungain, runmean, runresol;
187 int runno;
188 TTree *rungtree= (TTree*) fruncalib-> Get("runcalib");
189 rungtree -> SetBranchAddress("rungain", &rungain);
190 rungtree -> SetBranchAddress("runmean", &runmean);
191 rungtree -> SetBranchAddress("runno", &runno);
192 rungtree -> SetBranchAddress("runresol", &runresol);
193 rungtree -> GetEntry(0);
194 int NRUN = rungtree -> GetEntries();
195 N_run = rungtree -> GetEntries();
196 //cout<<"N_run = "<<N_run<<endl;
197 //cout<<"NRUN = "<<NRUN<<endl;
198 for(int runid =0; runid<NRUN; runid++) {
199 rungtree->GetEntry(runid);
200 m_rung[0][runid] = rungain;
201 m_rung[1][runid] = runmean;
202 m_rung[2][runid] = runno;
203 m_rung[3][runid] = runresol;
204 //cout<<"runid :"<<runid<<" run_no = "<<m_rung[2][runid]<<" run_mean = "<<m_rung[1][runid]<<" run_gain = "<<m_rung[0][runid]<<" runresol = "<<m_rung[3][runid]<<endl;
205 }
206
207 //cout<<"run by run const ------------------------------------------------- "<<endl;
208 for(int i=0;i<4;i++) {
209 for(int j=0;j<NRUN;j++) {
210 //std::cout<<"rung["<<i<<"]["<<j<<"]="<<m_rung[i][j]<<endl;
211 //std::cout<<"\n";
212 }
213 }*/
214
215 for(int i=0;i<4;i++) {
216 for(int j=0;j<43;j++) {
217 m_ddg[i][j] = test->getddg(i,j);
218 m_ggs[i][j] = test->getggs(i,j);
219 m_zdep[i][j] = test->getzdep(i,j);
220 // m_enta[i][j] = test->getenta(i,j);
221 //std::cout<<"ddg["<<i<<"]["<<j<<"]="<<m_ddg[i][j];
222 //std::cout<<" ggs["<<i<<"]["<<j<<"]="<<m_ggs[i][j];
223 //std::cout<<" zdep["<<i<<"]["<<j<<"]="<<m_zdep[i][j];
224 //std::cout<<"\n";
225 }
226 }
227
228 m_dedx_gain = test->getgain();
229 std::cout<<"gain="<<m_dedx_gain<<"\n";
230
231 m_dedx_resol = test->getresol();
232 std::cout<<"resol="<<m_dedx_resol<<"\n";
233
234 for(int i=0;i<43;i++) {
235 m_layer_g[i] = test -> getlayerg(i);
236 if(m_layer_g[i]>2.0 || m_layer_g[i]<0.5) m_layer_g[i] =1;
237 //std::cout<<"layerg["<<i<<"]="<<m_layer_g[i]<<endl;
238 }
239
240 /*const int n = 6796;
241 double wireg[n];
242 TFile* fwiregcalib=new TFile("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/wireg_calib/wiregcalib.root");
243 TTree *wiregtree= (TTree*)fwiregcalib -> Get("wiregcalib");
244 wiregtree -> SetBranchAddress("wireg", &wireg);
245 wiregtree -> GetEntry(0);
246 //cout<<"wire gain ------------------------------------------------- "<<endl;
247 for(int i=0; i<n; i++){
248 m_wire_g[i] = wireg[i];
249 //std::cout<<"wireg["<<i<<"] = "<<m_wire_g[i]<<"\n";
250 }*/
251
252 for(int i=0;i<6796;i++) {
253 m_wire_g[i] = test -> getwireg(i);
254 //std::cout<<"wireg["<<i<<"] = "<<m_wire_g[i]<<"\n";
255 }
256
257 //int kk=0;
258 //if(N_run<11397) kk=0;
259 //else if(N_run<12739) kk=1;
260 //else if(N_run<14395) kk=2;
261 //else kk=3;
262 int kk=3;
263 for(int i=0;i<6796;i++) {
264 m_valid[i] = 1;
265 /*if(i<=483) {m_valid[i] =0;}*/
266 for(int j=0; j<25; j++){
267 if( i == dead_wire[kk][j] )
268 {m_valid[i] = 0;
269 //cout<<"N_run: "<<N_run<<" kk: "<<kk<<endl;
270 //cout<<"m_valid["<<i<<"]: "<<m_valid[i]<<endl;
271 }
272 }
273 //std::cout<<"valid["<<i<<"]="<<m_valid[i]<<"\n";
274 }
275
276 //tempoary way to get costheta constants
277 cos_k.clear();
278 cos_b.clear();
279 //double cos1[80];
280 if(true){
281 const int nbin=80;
282 vector<double> cos;
283 /*TFile* fcosthecalib=new TFile("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/costhe_calib/costhetacalib2.root");
284 TTree *costhetree= (TTree*)fcosthecalib -> Get("costhetacalib");
285 TBranch* b_costheta;
286 double costheta;
287 costhetree->SetBranchAddress("costheta",&costheta,&b_costheta);
288 //const int nbin=tc->GetEntries();
289 //cout<<"costheta gain ------------------------------------------------- "<<endl;
290 for(int i=0; i<nbin; i++){
291 costhetree->GetEntry(i);
292 cos.push_back(costheta);
293 //cout<<"i : "<<i<<" costheta gain : "<<costheta<<endl;
294 }*/
295
296 for(int i=0; i<nbin; i++){
297 cos.push_back(test -> get_costheta(i));
298 }
299 for(int i=0;i<nbin-1;i++){
300 double k=0;
301 double b=0;
302 if(cos[i]>0.0001){ // not empty bin corresponding to large angle
303 if(cos[i+1]>0.0001){
304 double x1=-1.00+(0.5+i)*(2.00/nbin);
305 double y1=cos[i];
306 double x2=-1.00+(0.5+i+1)*(2.00/nbin);
307 double y2=cos[i+1];
308 // correct around absolute large cos(theta) angle
309 if(fabs(x1)>0.8) x1 = f_larcos(x1, nbin);
310 if(fabs(x2)>0.8) x2 = f_larcos(x2, nbin);
311 k=(y2-y1)/(x2-x1); // k is the slope
312 b=y1-k*x1;
313 //cout << "x1 " << x1 << " x2 " << x2 << " y1 " << y1 << " y2 " << y2 << endl;
314 }
315 else if(i>0){
316 if(cos[i-1]>0.0001){
317 double x1=-1.00+(0.5+i-1)*(2.00/nbin);
318 double y1=cos[i-1];
319 double x2=-1.00+(0.5+i)*(2.00/nbin);
320 double y2=cos[i];
321 // correct around absolute large cos(theta) angle
322 if(fabs(x1)>0.8) x1 = f_larcos(x1, nbin);
323 if(fabs(x2)>0.8) x2 = f_larcos(x2, nbin);
324 k=(y2-y1)/(x2-x1);
325 b=y1-k*x1;
326 //cout << "x1 " << x1 << " x2 " << x2 << " y1 " << y1 << " y2 " << y2 << endl;
327 }
328 }
329 }
330 else {
331 if(i>0){
332 if(cos[i-1]>0.0001 && cos[i+1]>0.0001){
333 double x1=-1.00+(0.5+i-1)*(2.00/nbin);
334 double y1=cos[i-1];
335 double x2=-1.00+(0.5+i+1)*(2.00/nbin);
336 double y2=cos[i+1];
337 // correct around absolute large cos(theta) angle
338 if(fabs(x1)>0.8) x1 = f_larcos(x1, nbin);
339 if(fabs(x2)>0.8) x2 = f_larcos(x2, nbin);
340 k=(y2-y1)/(x2-x1);
341 b=y1-k*x1;
342 //cout << "x1 " << x1 << " x2 " << x2 << " y1 " << y1 << " y2 " << y2 << endl;
343 }
344 }
345 }
346 //cout<<"i : "<<i<<" costheta gain : "<<cos[i] << " k=" << k << " b=" << b <<endl;
347 cos_k.push_back(k);
348 cos_b.push_back(b);
349 }
350 }
351
352 t0_k.clear();
353 t0_b.clear();
354 if(true){
355 //const int nbin=35;
356 /*vector<double> xbin;
357 vector<double> ybin;
358 TFile* ft0calib=new TFile("/ihepbatch/besd13/chunlei/calib/t0calib.root");
359 TTree *tree_t0=(TTree*)ft0calib->Get("t0calib");
360 TBranch* b_t0;
361 TBranch* b_dedx;
362 double t0, dedx;
363 tree_t0->SetBranchAddress("t0",&t0,&b_t0);
364 tree_t0->SetBranchAddress("dedx",&dedx,&b_dedx);
365 const int nbin=tree_t0->GetEntries();
366 //cout<<"costheta t0 ------------------------------------------------- "<<endl;
367 for(int i=0; i<tree_t0->GetEntries(); i++){
368 tree_t0->GetEntry(i);
369 xbin.push_back(t0);
370 ybin.push_back(dedx);
371 //cout<<"xbin = "<<t0<<" ybin = "<<dedx<<endl;
372 }*/
373
374 const int nbin=35;
375 vector<double> xbin;
376 vector<double> ybin;
377 for(int i=0; i<35; i++){
378 xbin.push_back(test ->get_t0(i));
379 ybin.push_back(test ->get_dedx(i));
380 //cout<<"xbin = "<<test ->get_t0(i)<<" ybin = "<<test ->get_dedx(i)<<endl;
381 }
382 for(int i=0;i<nbin-1;i++){
383 double k=0;
384 double b=0;
385 if(ybin[i]>0){
386 if(ybin[i+1]>0){
387 double x1=xbin[i];
388 double y1=ybin[i];
389 double x2=xbin[i+1];
390 double y2=ybin[i+1];
391 k=(y2-y1)/(x2-x1);
392 b=(y1*x2-y2*x1)/(x2-x1);
393 }
394 else if(i>0){
395 if(ybin[i-1]>0){
396 double x1=xbin[i-1];
397 double y1=ybin[i-1];
398 double x2=xbin[i];
399 double y2=ybin[i];
400 k=(y2-y1)/(x2-x1);
401 b=(y1*x2-y2*x1)/(x2-x1);
402 }
403 }
404 }
405 else {
406 if(i>0){
407 if(ybin[i-1]>0 && ybin[i+1]>0){
408 double x1=xbin[i-1];
409 double y1=ybin[i-1];
410 double x2=xbin[i+1];
411 double y2=ybin[i+1];
412 k=(y2-y1)/(x2-x1);
413 b=(y1*x2-y2*x1)/(x2-x1);
414 }
415 }
416 }
417 t0_k.push_back(k);
418 t0_b.push_back(b);
419 }
420 }
421
422 if(true){
423 /*TFile fconst9("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/doca_eangle_calib/check/doca_eangle_kal_doca110.root");
424 const int n = 1600;
425 double Iner_gain[n], Iner_chi[n], Iner_hits[n];
426 double Out_gain[n], Out_chi[n], Out_hits[n];
427 double Id_doca[n], Ip_eangle[n];
428
429 TTree *tree_docasin= (TTree*)fconst9.Get("docasincalib");
430 tree_docasin -> SetBranchAddress("Iner_gain", &Iner_gain);
431 tree_docasin -> SetBranchAddress("Iner_chi", &Iner_chi);
432 tree_docasin -> SetBranchAddress("Iner_hits", &Iner_hits);
433 tree_docasin -> SetBranchAddress("Out_gain", &Out_gain);
434 tree_docasin -> SetBranchAddress("Out_chi", &Out_chi);
435 tree_docasin -> SetBranchAddress("Out_hits", &Out_hits);
436 tree_docasin -> SetBranchAddress("Id_doca", &Id_doca);
437 tree_docasin -> SetBranchAddress("Ip_eangle", &Ip_eangle);
438 tree_docasin -> GetEntry(0);
439 double docaeangle_gain[n];
440 double docaeangle_chisq[n];
441 double docaeangle_entries[n];
442 double cell[n];
443 //cout<<"doca eangle gain ------------------------------------------------- "<<endl;
444 for(int i=0; i<n; i++){
445 tree_docasin->GetEntry(i);
446 cell[i] = i;
447 docaeangle_gain[i] = Out_gain[i];
448 docaeangle_chisq[i] = Out_chi[i];
449 docaeangle_entries[i] = Out_hits[i];
450 int Id_Doca_temp = Id_doca[i];
451 int Ip_Eangle_temp = Ip_eangle[i];
452 m_docaeangle[Id_Doca_temp][Ip_Eangle_temp] = Out_gain[i];
453 //cout<<i<<" "<<Id_Doca_temp<<" "<<Ip_Eangle_temp<<" "<<docaeangle_gain[i]<<endl;
454 //m_docaeangle[Id_Doca_temp][Ip_Eangle_temp] = docaeangle_gain[i];
455 //cout<<i<<" "<<Id_doca[i]<<" "<<Ip_eangle[i]<<" Entries : "<<docaeangle_entries[i]<<" Gain value : "<<docaeangle_gain[i]<<" chisq : "<<docaeangle_chisq[i] <<endl;
456
457 }
458 m_docaeangle[38][28]=0.72805;*/
459
460 const int n = 1600;
461 double docaeangle_gain[n];
462 double docaeangle_chisq[n];
463 double docaeangle_entries[n];
464 double cell[n];
465 for(int i=0; i<n; i++){
466 cell[i] = i;
467 docaeangle_gain[i] = test -> get_out_gain(i);
468 docaeangle_chisq[i] = test -> get_out_chi(i);
469 docaeangle_entries[i] = test -> get_out_hits(i);
470 int Id_Doca_temp = test -> get_id_doca(i);
471 int Ip_Eangle_temp = test -> get_ip_eangle(i);
472 m_docaeangle[Id_Doca_temp][Ip_Eangle_temp] = docaeangle_gain[i];
473 if(m_debug && (Id_Doca_temp==m_debug_i) && (Ip_Eangle_temp==m_debug_j)) std::cout<<"docaeangle_gain["<<Id_Doca_temp<<"]["<<Ip_Eangle_temp<<"]="<<m_docaeangle[m_debug_i][m_debug_j] << std::endl;
474 //cout<<i<<" "<<Id_Doca_temp<<" "<<Ip_Eangle_temp<<" "<<docaeangle_gain[i]<<endl;
475 }
476 }
477
478
479 //get 1d entrance angle constants
480 int onedsize=test->get_enanglesize();
481 m_venangle.clear();
482 for(int i=0; i< onedsize; i++){
483 m_venangle.push_back(test->get_enangle(i));
484 }
485
486
487 //temporary way to get hadron saturation constants
488 //for Psai hadron saturation
489
490 //only valide for jpsi data now !!!!!
491 //have to find ways pass constans for different data
492 /* m_alpha= 1.35630e-02;
493 m_gamma= -6.78907e-04;
494 m_delta= 1.18037e-02;
495 m_power= 1.66268e+00;
496 m_ratio= 9.94775e-01;*/
497
498 const int hadronNo=test -> get_hadronNo();
499 // cout<<"@@@hadronNO===="<<hadronNo<<endl;
500 m_alpha=test -> get_hadron(0);
501 // cout<<"@@@@m_alpha===="<<m_alpha<<endl;
502 m_gamma=test -> get_hadron(1);
503 // cout<<"@@@@m_gamma===="<<m_gamma<<endl;
504 m_delta=test -> get_hadron(2);
505 // cout<<"@@@@m_delta====="<<m_delta<<endl;
506 m_power=test -> get_hadron(3);
507 // cout<<"@@@@m_power====="<<m_power<<endl;
508 m_ratio=test -> get_hadron(4);
509 // cout<<"@@@m_ratio======"<<m_ratio<<endl;
510
511
512
513 /*
514 m_delta=0.1;
515 m_alpha=0.0282;
516 m_gamma=-0.030;
517 m_power=1.0;
518 m_ratio=1.0;
519 //m_delta=0.1;
520 //m_alpha=0.0382;
521 //m_gamma=-0.0438;
522 */
523
524 //m_initConstFlg =true;
525}
526
527double DedxCorrecSvc::WireGainCorrec(int wireid, double ex) const{
528 if( m_wire_g[wireid] > 0 ){
529 double ch_dedx = (ex/m_wire_g[wireid])*m_valid[wireid];
530 return ch_dedx;
531 }
532 else if( m_wire_g[wireid] == 0 ){
533 return ex;
534 }
535 else return 0;
536}
537
538double
539DedxCorrecSvc::SaturCorrec( int layer, double costheta, double dedx ) const {
540 //cout<<"DedxCorrecSvc::SaturCorrec"<<endl;
541 double dedx_ggs;
542 //cout<<"costheta vaule is = "<<costheta<<endl;
543 costheta = fabs(costheta);
544 if(m_par_flag == 1) {
545 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer]*costheta +
546 m_ggs[2][layer]*pow(costheta,2) + m_ggs[3][layer]*pow(costheta,3);
547 } else if(m_par_flag == 0) {
548 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer]*T1(costheta) +
549 m_ggs[2][layer]*T2(costheta) + m_ggs[3][layer]*T3(costheta);
550 }
551 //cout<<"dedx_ggs is = "<<dedx_ggs<<endl;
552 dedx_ggs = dedx/dedx_ggs;
553 if(dedx_ggs>0.0) return dedx_ggs;
554 else return dedx;
555}
556
557
558double
559DedxCorrecSvc::CosthetaCorrec( double costheta, double dedx ) const {
560 //cout<<"DedxCorrecSvc::CosthetaCorrec"<<endl;
561 double dedx_cos;
562 //cout<<"costheta vaule is = "<<costheta<< " dedx = " << dedx << endl;
563 if(fabs(costheta)>1.0) return dedx;
564
565 const int nbin=cos_k.size()+1;
566 const double step=2.00/nbin;
567 int n= (int)((costheta+1.00-0.5*step)/step);
568 if(costheta>(1.00-0.5*step))
569 n=nbin-2;
570
571 if(costheta>-0.5*step && costheta<=0)
572 dedx_cos=cos_k[n-1]*costheta+cos_b[n-1];
573 else if(costheta>0 && costheta<0.5*step)
574 dedx_cos=cos_k[n+1]*costheta+cos_b[n+1];
575 else
576 dedx_cos=cos_k[n]*costheta+cos_b[n];
577
578 //cout << "cos_k[n]="<<cos_k[n] << " cos_b[n]=" << cos_b[n] <<
579 // " dedx_cos=" << dedx_cos << " final dedx=" << dedx/dedx_cos << endl;
580
581 //conside physical edge around 0.93
582 if(nbin==80){ // temporally for only nbin = 80
583 if(costheta>0.80 && costheta<0.95){
584 n = 77;
585 if(costheta<0.9282) n--;
586 if(costheta<0.9115) n--;
587 if(costheta<0.8877) n--;
588 if(costheta<0.8634) n--;
589 if(costheta<0.8460) n--;
590 if(costheta<0.8089) n--;
591 dedx_cos=cos_k[n]*costheta+cos_b[n];
592 }
593 else if(costheta<-0.80 && costheta>-0.95){
594 n = 2;
595 if(costheta>-0.9115) n++;
596 if(costheta>-0.8877) n++;
597 if(costheta>-0.8634) n++;
598 if(costheta>-0.8460) n++;
599 if(costheta>-0.8089) n++;
600 dedx_cos=cos_k[n]*costheta+cos_b[n];
601 }
602 }
603
604 if(dedx_cos>0){
605 dedx_cos = dedx/dedx_cos;
606 return dedx_cos;
607 }
608 else return dedx;
609}
610
611
612double
613DedxCorrecSvc::T0Correc( double t0, double dedx ) const {
614 // cout<<"DedxCorrecSvc::T0Correc"<<endl;
615 double dedx_t0;
616 const int nbin=t0_k.size()+1;
617 if(nbin!=35)
618 cout<<"there should be 35 bins for t0 correction, check it!"<<endl;
619
620 int n=0;
621 if(t0<575)
622 n=0;
623 else if(t0<610)
624 n=1;
625 else if(t0>=610 && t0<1190){
626 n=(int)( (t0-610.0)/20.0 ) + 2;
627 }
628 else if(t0>=1190 && t0<1225)
629 n=31;
630 else if(t0>=1225 && t0<1275)
631 n=32;
632 else if(t0>=1275)
633 n=33;
634
635 dedx_t0=t0_k[n]*t0+t0_b[n];
636
637 if(dedx_t0>0){
638 dedx_t0 = dedx/dedx_t0;
639 return dedx_t0;
640 }
641 else return dedx;
642}
643
644double
645DedxCorrecSvc::HadronCorrec( double costheta, double dedx ) const {
646 // cout<<"DedxCorrecSvc::HadronCorrec"<<endl;
647 //constant 1.00 in the following function is the mean dedx of normalized electrons.
648 dedx=dedx/550.00;
649 return D2I(costheta, I2D(costheta,1.00)/1.00*dedx)*550;
650}
651double
652DedxCorrecSvc::LayerGainCorrec( int layid, double dedx ) const {
653 // cout<<"DedxCorrecSvc::LayerGainCorrec"<<endl;
654 if( m_layer_g[layid] > 0 ){
655 double ch_dedx = dedx/m_layer_g[layid];
656 return ch_dedx;
657 }
658 else if( m_layer_g[layid] == 0 ){
659 return dedx;
660 }
661 else return 0;
662}
663
664
665double
666DedxCorrecSvc::RungCorrec( int runNO, int evtNO, double dedx ) const{
667 //cout<<"DedxCorrecSvc::RungCorrec"<<endl;
668 double dedx_rung;
669 int run_ture =0;
670 // cout << " runNO : "<<runNO << " evtNO " << evtNO << endl;
671
672 for(int j=0;j<100000;j++) {
673 // for(int j=0;j<10;j++) {
674 if((runNO == m_rung[2][j]) && (evtNO >= m_rung[4][j]) && (evtNO <= m_rung[5][j])) {
675 dedx_rung = dedx/m_rung[0][j];
676 // cout <<"j " << j << "runno " << m_rung[2][j] << " evtNO " << evtNO << " evtfrom " << m_rung[4][j] << " evtto " << m_rung[5][j] << " rung " << m_rung[0][j] <<endl;
677 run_ture =1;
678 return dedx_rung;
679 }
680 }
681 if(run_ture ==0)
682 {
683 cout<<"Warning!!! in DedxCorrecSvc::RungCorrec(): no rungain to "<<runNO<<endl;
684 exit(1);
685 }
686}
687
688
689double
690DedxCorrecSvc::DriftDistCorrec( int layer, double dd, double dedx ) const {
691 // cout<<"DedxCorrecSvc::DriftDistCorrec"<<endl;
692 double dedx_ddg;
693 //cout<<"m_par_flag = "<<m_par_flag<<endl;
694 //cout<<"dd vaule is = "<<dd<<endl;
695 dd = fabs(dd);
696 if(m_par_flag == 1) {
697 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer]*dd +
698 m_ddg[2][layer]*dd*dd + m_ddg[3][layer]*pow(dd,3);
699 } else if(m_par_flag == 0) {
700 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer]*T1(dd) +
701 m_ddg[2][layer]*T2(dd) + m_ddg[3][layer]*T3(dd);
702 }
703 //cout<<"dedx_ddg is = "<<dedx_ddg<<endl;
704 dedx_ddg = dedx/dedx_ddg;
705 if (dedx_ddg>0.0) return dedx_ddg;
706 else return dedx;
707}
708
709
710double
711DedxCorrecSvc::EntaCorrec( int layer, double eangle, double dedx ) const {
712 // cout<<"DedxCorrecSvc::EntaCorrec"<<endl;
713// double dedx_enta;
714// if(m_par_flag == 1) {
715// dedx_enta = m_enta[0][layer] + m_enta[1][layer]*sinenta +
716// m_enta[2][layer]*sinenta*sinenta + m_enta[3][layer]*pow(sinenta,3);
717// } else if(m_par_flag == 0) {
718// dedx_enta = m_enta[0][layer] + m_enta[1][layer]*T1(sinenta) +
719// m_enta[2][layer]*T2(sinenta) + m_enta[3][layer]*T3(sinenta);
720// }
721// dedx_enta = dedx/dedx_enta;
722// if (dedx_enta>0.0) return dedx_enta;
723// else return dedx;
724
725 if(eangle>-0.25 && eangle<0.25){
726 double stepsize= 0.5/m_venangle.size();
727 int nsize= m_venangle.size();
728 double slope=0;
729 double offset=1;
730 double y1=0,y2=0,x1=0,x2=0;
731
732 if(eangle>=-0.25+0.5*stepsize && eangle<=0.25-0.5*stepsize){
733 int bin = (int)( (eangle-(-0.25+0.5*stepsize))/stepsize);
734 y1=m_venangle[bin];
735 x1=-0.25+0.5*stepsize+bin*stepsize;
736 y2=m_venangle[bin+1];
737 x2=-0.25+1.5*stepsize+bin*stepsize;
738 }
739 else if(eangle<=-0.25+0.5*stepsize){
740 y1=m_venangle[0];
741 x1=-0.25+0.5*stepsize;
742 y2=m_venangle[1];
743 x2=-0.25+1.5*stepsize;
744 }
745 else if( eangle>=0.25-0.5*stepsize){
746 y1=m_venangle[nsize-2];
747 x1=0.25-1.5*stepsize;
748 y2=m_venangle[nsize-1];
749 x2=0.25-0.5*stepsize;
750 }
751 double kcof= (y2-y1)/(x2-x1);
752 double bcof= y1-kcof*x1;
753 double ratio = kcof*eangle+bcof;
754 dedx=dedx/ratio;
755 }
756
757 return dedx;
758}
759
760
761double
762DedxCorrecSvc::ZdepCorrec( int layer, double z, double dedx ) const {
763 // cout<<"DedxCorrecSvc::ZdepCorrec"<<endl;
764 double dedx_zdep;
765 if(m_par_flag == 1) {
766 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer]*z +
767 m_zdep[2][layer]*z*z + m_zdep[3][layer]*pow(z,3);
768 } else if(m_par_flag == 0) {
769 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer]*T1(z) +
770 m_zdep[2][layer]*T2(z) + m_zdep[3][layer]*T3(z);
771 }
772 dedx_zdep = dedx/dedx_zdep;
773 if (dedx_zdep>0.0) return dedx_zdep;
774 else return dedx;
775
776 //return dedx;
777}
778
779
780double
781DedxCorrecSvc::DocaSinCorrec( int layer,double doca, double eangle, double dedx ) const {
782 if(m_debug) cout<<"DedxCorrecSvc::DocaSinCorrec"<<endl;
783 // cout<<"doca = "<<doca<<" eangle = "<<eangle<<endl;
784
785 if(eangle>0.25) eangle = eangle -0.5;
786 else if(eangle < -0.25) eangle = eangle +0.5;
787 int iDoca;
788 int iEAng = 0;
789 doca = doca/doca_norm[layer];
790 iDoca =(Int_t)floor((doca-Out_DocaMin)/Out_Stepdoca);
791 if(doca<Out_DocaMin) iDoca=0;
792 if(doca>Out_DocaMax) iDoca=NumSlicesDoca-1;
793 if(iDoca >=(Int_t)floor((Out_DocaMax-Out_DocaMin)/Out_Stepdoca) )
794 iDoca = (Int_t)floor((Out_DocaMax-Out_DocaMin)/Out_Stepdoca)-1;
795 if(iDoca<0) iDoca = 0;
796 if(m_debug) cout<<"iDoca : "<<iDoca<<" doca : "<<doca<<" Out_DocaMin : "<<Out_DocaMin<<" Out_Stepdoca : "<<Out_Stepdoca<<endl;
797
798 //iEAng = (Int_t)floor((eangle-Phi_Min)/IEangle_step);
799 for(int i =0; i<14; i++){
800 //iEAng =0;
801 if(eangle <Eangle_value_cut[i] || eangle > Eangle_value_cut[i+1]) continue;
802 else if(eangle>= Eangle_value_cut[i] && eangle < Eangle_value_cut[i+1]) {
803 for(int k =0; k<i; k++){
804 iEAng+= Eangle_cut_bin[k];
805 }
806 double eangle_bin_cut_step = (Eangle_value_cut[i+1] - Eangle_value_cut[i])/Eangle_cut_bin[i];
807 int temp_bin = int((eangle-Eangle_value_cut[i])/eangle_bin_cut_step);
808 iEAng+= temp_bin;
809 }
810 }
811 //cout<<iDoca <<" "<<iEAng<<endl;
812 if(m_docaeangle[iDoca][iEAng]>0) {
813 if(m_debug && (iDoca==m_debug_i) && (iEAng==m_debug_j))
814 cout << "doca: " << doca << " eang" << eangle << " dedx before: " << dedx << endl;
815 dedx = dedx/m_docaeangle[iDoca][iEAng];
816 if(m_debug && (iDoca==m_debug_i) && (iEAng==m_debug_j))
817 cout << "gain: " << m_docaeangle[iDoca][iEAng] << " dedx after: " << dedx << endl;
818 }
819 return dedx;
820}
821
822
823double
824DedxCorrecSvc::DipAngCorrec(int layer, double costheta, double dedx) const {
825}
826
827double
828DedxCorrecSvc::GlobalCorrec( double dedx) const{
829 if( m_mdcg_flag == 0 ) return dedx;
830 double calib_ex = dedx;
831 if( m_dedx_gain > 0 ) calib_ex /= m_dedx_gain;
832 return calib_ex;
833}
834
835
836
837double
838DedxCorrecSvc::CellCorrec( int ser, double adc, double dd, double sinenta,
839 double z, double costheta ) const {
840
841 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0
842 && m_layerg_flag == 0 && m_zdep_flag == 0 &&
843 m_ggs_flag == 0) return adc;
844 adc = m_valid[ser]*m_wire_g[ser]*adc;
845 // int lyr = Wire2Lyr(ser);
846 int lyr = ser;
847 double ex = adc;
848 double correct1_ex , correct2_ex, correct3_ex, correct4_ex, correct5_ex;
849
850 if( m_ggs_flag) {
851 correct1_ex = SaturCorrec( lyr, costheta, adc );
852 ex = correct1_ex;
853 } else {
854 correct1_ex = adc;
855 }
856
857 if( m_ddg_flag) {
858 correct2_ex = DriftDistCorrec( lyr, dd, correct1_ex );
859 ex = correct2_ex;
860 } else {
861 correct2_ex =correct1_ex;
862 }
863 if( m_enta_flag) {
864 correct3_ex = DriftDistCorrec( lyr, sinenta, correct2_ex );
865 ex = correct3_ex;
866 } else {
867 correct3_ex =correct2_ex;
868 }
869
870 if( m_zdep_flag) {
871 correct4_ex = ZdepCorrec( lyr, z, correct3_ex );
872 ex = correct4_ex;
873 } else {
874 correct4_ex = correct3_ex;
875 }
876
877 if( m_layerg_flag) {
878 correct5_ex = LayerGainCorrec( lyr, correct4_ex );
879 ex = correct5_ex;
880 } else {
881 correct5_ex = correct4_ex;
882 }
883 return ex;
884
885}
886
887double
888DedxCorrecSvc::LayerCorrec( int layer, double z, double costheta, double ex ) const{
889 //cout<<"DedxCorrecSvc::LayerCorrec"<<endl;
890
891 if( m_zdep_flag == 0 && m_ggs_flag == 0 ) return ex;
892
893 double calib_ex = ZdepCorrec( layer, z, ex );
894 if( m_ggs_flag != 0 ) calib_ex = SaturCorrec( layer, costheta, calib_ex );
895 return calib_ex;
896
897}
898
899double
900DedxCorrecSvc::TrkCorrec( double costheta, double dedx ) const{
901 // cout<<"DedxCorrecSvc::TrkCorrec"<<endl;
902 if( m_mdcg_flag == 0 ) return dedx;
903 ///dEdx calib. for each track
904 double calib_ex = dedx;
905
906 double run_const = m_dedx_gain;
907 if( run_const > 0 && m_mdcg_flag != 0 ) calib_ex /= run_const;
908
909 return calib_ex;
910
911}
912
913
914double
915DedxCorrecSvc::StandardCorrec( int runFlag, int ntpFlag, int runNO, int evtNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double z, double costheta ) const {
916 // cout<<"DedxCorrecSvc::StandardCorrec"<<endl;
917 //int layid = MdcID::layer(mdcid);
918 //int localwid = MdcID::wire(mdcid);
919 //int w0id = geosvc->Layer(layid)->Wirst();
920 //int wid = w0id + localwid;
921 //double pathl = PathL(ntpFlag, hel, layid, localwid, z);
922 ////pathl= PathLCosmic(hel, layid, localwid, z, sigmaz );
923 double ex = adc;
924 if( runNO<0 ) return ex;
925 ex = ex*1.5/pathl;
926 //double ex = adc/pathl;
927 //if( runNO<0 ) return ex;
928 if( ntpFlag ==0) return ex;
929 //double ex = adc/pathl;
930 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_dcasin_flag==0 && m_ddg_flag == 0
931 && m_layerg_flag == 0 && m_zdep_flag == 0 &&
932 m_ggs_flag == 0) return ex;
933
934 if(m_rung_flag) {
935 ex = RungCorrec( runNO, evtNO, ex );
936 }
937
938 if( m_wireg_flag) {
939 ex = WireGainCorrec(wid, ex);
940 }
941
942 if( m_dcasin_flag) {
943 if(runFlag<3) {
944 ex = DriftDistCorrec( layid, dd, ex );
945 }
946 else{ ex = DocaSinCorrec(layid, dd, eangle, ex);}
947 }
948
949 if(m_enta_flag && runFlag>=3){
950 ex = EntaCorrec(layid, eangle, ex);
951 }
952
953 // ddg is not use anymore, it's replaced by DocaSin
954 if(m_ddg_flag) {
955 ex = DriftDistCorrec( layid, dd, ex );
956 }
957
958 if(m_ggs_flag) {
959 if(runFlag <3 ){
960 ex = SaturCorrec( layid, costheta, ex );
961 }
962 else{ ex = CosthetaCorrec( costheta, ex );}
963 // Staur is OLD, Costheta is NEW.
964 }
965
966 if( m_sat_flag) {
967 ex = HadronCorrec( costheta, ex );
968 }
969
970 if( m_zdep_flag) {
971 ex = ZdepCorrec( layid, z, ex );
972 }
973
974 if( m_layerg_flag) {
975 ex = LayerGainCorrec( layid, ex );
976 }
977
978 if( m_dip_flag){
979 ex = DipAngCorrec(layid, costheta, ex);
980 }
981
982 if( m_mdcg_flag) {
983 ex = GlobalCorrec( ex );
984 }
985 return ex;
986}
987
988double
989DedxCorrecSvc::StandardHitCorrec(int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, int evtNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta ) const {
990 if(m_debug) cout<<"DedxCorrecSvc::StandardHitCorrec"<<endl;
991 //int layid = MdcID::layer(mdcid);
992 //int localwid = MdcID::wire(mdcid);
993 //int w0id = geosvc->Layer(layid)->Wirst();
994 //int wid = w0id + localwid;
995 //double pathl = PathL(ntpFlag, hel, layid, localwid, z);
996 //cout<<"DedxCorrecSvc::StandardHitCorrec pathl= "<<pathl<<endl;
997 //pathl= PathLCosmic(hel, layid, localwid, z, sigmaz );
998 double ex = adc;
999 if( runNO<0 ) return ex;
1000 ex = ex*1.5/pathl;
1001 //if( runNO<0 ) return ex;
1002 //double ex = adc/pathl;
1003 if( ntpFlag ==0) return ex;
1004 //if(ntpFlag>0) cout<<"dE/dx value after path correction: "<<ex<<endl;
1005 //double ex = adc/pathl;
1006 //cout<<m_rung_flag << " "<<m_wireg_flag<<" "<<m_ddg_flag<<" "<<m_ggs_flag<<endl;
1007 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0
1008 && m_layerg_flag == 0 && m_zdep_flag == 0 && m_dcasin_flag==0
1009 && m_ggs_flag == 0 && m_enta_flag==0) return ex;
1010
1011 if(m_rung_flag) {
1012 ex = RungCorrec( runNO, evtNO, ex );
1013 }
1014 //if(ntpFlag>0) cout<<"dE/dx value after run by run correction: "<<ex<<endl;
1015
1016 if( m_wireg_flag) {
1017 ex = WireGainCorrec(wid, ex);
1018 }
1019 //if(ntpFlag>0) cout<<"dE/dx value after wire gain correction: "<<ex<<endl;
1020
1021 if( m_dcasin_flag) {
1022 if(runFlag<3){
1023 ex = DriftDistCorrec( layid, dd, ex );
1024 }
1025 else{
1026 //cout<<layid<<" "<<dd<<" "<<eangle<<" "<<ex<<endl;
1027 ex = DocaSinCorrec(layid, dd, eangle, ex);
1028 }
1029 }
1030
1031 // 1D entrance angle correction
1032 if(m_enta_flag && runFlag>=3){
1033 ex = EntaCorrec(layid, eangle, ex);
1034 }
1035
1036 // ddg is not used anymore
1037 if( m_ddg_flag) {
1038 ex = DriftDistCorrec( layid, dd, ex );
1039 }
1040 //if(ntpFlag>0) cout<<"dE/dx value after ddg correction: "<<ex<<endl;
1041
1042 if(m_ggs_flag) {
1043 if(runFlag <3 ){
1044 ex = SaturCorrec( layid, costheta, ex );
1045 }
1046 else{ ex = ex;} // do not do the cos(theta) correction at hit level
1047 }
1048 //if(ntpFlag>0) cout<<"dE/dx value after costheta correction: "<<ex<<endl;
1049 return ex;
1050}
1051
1052
1053double
1054DedxCorrecSvc::StandardTrackCorrec(int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, int evtNO, double ex, double costheta, double t0 ) const {
1055 if(m_debug) cout<<"DedxCorrecSvc::StandardTrackCorrec"<<endl;
1056 if( runNO<0 ) return ex;
1057 if( ntpFlag ==0) return ex;
1058 if( runFlag <3) return ex;
1059 if( calib_rec_Flag ==1){
1060 ex = CosthetaCorrec( costheta, ex );
1061 if(m_t0_flag) ex= T0Correc(t0, ex);
1062 return ex;
1063 }
1064
1065 //if(ntpFlag>0) cout<<"trcak value before costheta correction: "<<ex<<endl;
1066 if( m_ggs_flag) {
1067 ex = CosthetaCorrec( costheta, ex );
1068 }
1069 //if(ntpFlag>0) cout<<"trcak value after costheta correction: "<<ex<<endl;
1070 if(m_t0_flag){
1071 // if(runFlag==3) {ex= T0Correc(t0, ex);}
1072 // else if(runFlag>3) {ex=ex;}
1073 // do t0 correction for all data samples
1074 ex= T0Correc(t0, ex);
1075 }
1076 //if(ntpFlag>0) cout<<"trcak value after all correction: "<<ex<<endl;
1077 if( m_sat_flag) {
1078 ex = HadronCorrec( costheta, ex );
1079 }
1080
1081 if(m_rung_flag && calib_rec_Flag ==2){
1082 double rungain_temp =RungCorrec( runNO, evtNO, ex )/ex;
1083 ex = ex/rungain_temp;
1084 }
1085
1086 //if(ntpFlag>0) cout<<"trcak value no run gain correction: "<<ex<<endl;
1087 return ex;
1088
1089}
1090
1091
1092
1093void
1094DedxCorrecSvc::init_param() {
1095
1096 //cout<<"DedxCorrecSvc::init_param()"<<endl;
1097
1098 for( int i = 0; i<6796 ; i++) {
1099 m_valid[i] = 1.0;
1100 m_wire_g[i] = 1.0;
1101 }
1102 m_dedx_gain = 1.0;
1103 m_dedx_resol = 1.0;
1104 for(int j = 0; j<100; j++){
1105 m_rung[0][j] =1;
1106 }
1107 for(int j = 0; j<43; j++) {
1108 m_layer_g[j] = 1.0;
1109 m_ggs[0][j] = 1.0;
1110 m_ddg[0][j] = 1.0;
1111 m_enta[0][j] = 1.0;
1112 m_zdep[0][j] = 1.0;
1113 for(int k = 1; k<4; k++ ) {
1114 m_ggs[k][j] = 0.0;
1115 m_ddg[k][j] = 0.0;
1116 m_enta[k][j] = 0.0;
1117 m_zdep[k][j] = 0.0;
1118 }
1119 }
1120
1121 std::cout<<"DedxCorrecSvc::init_param()!!!"<<std::endl;
1122 std::cout<<"hello,init_param"<<endl;
1123 m_initConstFlg =true;
1124
1125}
1126
1127
1128void
1129DedxCorrecSvc::set_flag( int calib_F ) {
1130 // cout<<"DedxCorrecSvc::set_flag"<<endl;
1131 cout<<"calib_F is = "<<calib_F<<endl;
1132 if(calib_F<0){ m_enta_flag = 0; calib_F = abs(calib_F); }
1133 else m_enta_flag = 1;
1134 m_rung_flag = ( calib_F & 1 )? 1 : 0;
1135 m_wireg_flag = ( calib_F & 2 )? 1 : 0;
1136 m_dcasin_flag = ( calib_F & 4 )? 1 : 0;
1137 m_ggs_flag = ( calib_F & 8 )? 1 : 0;
1138 m_ddg_flag = 0;
1139 //m_ddg_flag = ( calib_F & 8 )? 1 : 0;
1140 m_t0_flag = ( calib_F & 16 )? 1 : 0;
1141 m_sat_flag = ( calib_F & 32 )? 1 : 0;
1142 m_layerg_flag = ( calib_F & 64 )? 1 : 0;
1143 //m_dcasin_flag = ( calib_F & 128 )? 1 : 0;
1144 m_dip_flag = ( calib_F & 256 )? 1 : 0;
1145 m_mdcg_flag = ( calib_F & 512 )? 1 : 0;
1146}
1147
1148
1149
1150//funcation to calculate the path pength of Cosmic data in the layer
1151/*double
1152 DedxCorrecSvc::PathLCosmic(const Dedx_Helix& hel, int layer, int cellid, double z,double sigmaz ) const {
1153////// calculate the cooridate of P0
1154HepPoint3D piv(hel.pivot());
1155double dr = hel.a()[0];
1156double phi0 = hel.a()[1];
1157double tanl = hel.a()[4];
1158double dz = hel.a()[3];
1159
1160double dr = -19.1901;
1161double phi0 = 3.07309;
1162double tanl = -0.466654;
1163double dz = 64.8542;
1164
1165double csf0 = cos(phi0);
1166double snf0 = sin(phi0);
1167double x_c = dr*csf0;
1168double y_c = dr*snf0;
1169double z_c = hel.a()[3];
1170
1171////// calculate the cooridate of hits
1172////calculate the track length in r_phi plane
1173
1174double m_crio[2];
1175double m_zb, m_zf, Calpha;
1176// double sintheta = sin(M_PI_2-atan(hel.a()[4]));
1177double sintheta = sin(M_PI_2-atan(tanl));
1178// retrieve Mdc geometry information
1179double rcsiz1= geosvc->Layer(layer)->RCSiz1();
1180double rcsiz2= geosvc->Layer(layer)->RCSiz2();
1181double rlay= geosvc->Layer(layer)->Radius();
1182int ncell= geosvc->Layer(layer)->NCell();
1183double phioffset= geosvc->Layer(layer)->Offset();
1184float shift= geosvc->Layer(layer)->Shift();
1185double slant= geosvc->Layer(layer)->Slant();
1186double length = geosvc->Layer(layer)->Length();
1187int type = geosvc->Layer(layer)->Sup()->Type();
1188// shift = shift*type;
1189
1190//conversion of the units(mm=>cm)
1191length = 0.1*length;
1192rcsiz1 = 0.1*rcsiz1;
1193rcsiz2 = 0.1*rcsiz2;
1194rlay = 0.1*rlay;
1195m_zb = 0.5*length;
1196m_zf = -0.5*length;
1197m_crio[0] = rlay - rcsiz1;
1198m_crio[1] = rlay + rcsiz2;
1199
1200if(z == 0)
1201{ if(rlay<= fabs(dr)) rlay = rlay + rcsiz2;
1202if(rlay<fabs(dr)) return -1;
1203double t_digi = sqrt(rlay*rlay - dr*dr);
1204double x_t_digi = x_c - t_digi*snf0;
1205double y_t_digi = y_c + t_digi*csf0;
1206double z_t_digi = z_c + t_digi*tanl;
1207double x__t_digi = x_c + t_digi*snf0;
1208double y__t_digi = y_c - t_digi*csf0;
1209double z__t_digi = z_c - t_digi*tanl;
1210double phi_t_digi = atan2(y_t_digi,x_t_digi);
1211double phi__t_digi = atan2(y__t_digi,x__t_digi);
1212phi_t_digi = fmod( phi_t_digi+4*M_PI,2*M_PI );
1213phi__t_digi = fmod( phi__t_digi+4*M_PI,2*M_PI );
1214double phibin_digi = 2.0*M_PI/ncell;
1215double phi_cellid_digi = phioffset + shift*phibin_digi*0.5 + cellid *phibin_digi;
1216if(fabs(phi_cellid_digi - phi_t_digi)<fabs(phi_cellid_digi - phi__t_digi))
1217z = z_t_digi;
1218else if (fabs(phi_cellid_digi - phi_t_digi)>fabs(phi_cellid_digi - phi__t_digi))
1219z = z__t_digi;
1220else z = z_t_digi;
1221}
1222
1223int sign = 1; ///assumed the first half circle
1224int epflag[2];
1225Hep3Vector iocand;
1226Hep3Vector cell_IO[2];
1227double dphi, downin;
1228Hep3Vector zvector;
1229if( type ) {
1230 downin = (z*z-m_zb*m_zb)*pow(tan(slant),2);
1231 m_crio[0] = sqrt(m_crio[0]*m_crio[0]+downin);
1232 m_crio[1] = sqrt(m_crio[1]*m_crio[1]+downin);
1233}
1234
1235// calculate t value
1236double t[2];
1237if(m_crio[1]<fabs(dr)) return -1;
1238else if(m_crio[0]<fabs(dr)) {
1239 t[0] = sqrt(m_crio[1]*m_crio[1] - dr*dr);
1240 t[1] = -t[0];
1241}
1242else{
1243 for( int i = 0; i < 2; i++ ) {
1244 if(m_crio[i]<fabs(dr)) return -1;
1245 t[i] = sqrt(m_crio[i]*m_crio[i] - dr*dr);
1246 }
1247}
1248
1249// calculate the cooridate of hits
1250double x_t[2],y_t[2],z_t[2];
1251double x__t[2],y__t[2],z__t[2];
1252double x_p[2],y_p[2],z_p[2];
1253for( int i = 0; i < 2; i++){
1254 x_t[i] = x_c - t[i]*snf0;
1255 y_t[i] = y_c + t[i]*csf0;
1256 z_t[i] = z_c + t[i]*tanl;
1257 x__t[i] = x_c + t[i]*snf0;
1258 y__t[i] = y_c - t[i]*csf0;
1259 z__t[i] = z_c - t[i]*tanl;
1260}
1261
1262double phi_in_t,phi_in__t, phi_out_t,phi_out__t,phi_t,phi__t;
1263double phibin = 2.0*M_PI/ncell;
1264double phi_cellid = phioffset + shift*phibin*(0.5-z/length) + cellid *phibin;
1265phi_in_t = atan2(y_t[0],x_t[0]);
1266phi_out_t = atan2(y_t[1],x_t[1]);
1267phi_in__t = atan2(y__t[0],x__t[0]);
1268phi_out__t = atan2(y__t[1],x__t[1]);
1269phi_t = atan2(((y_t[0]+y_t[1])/2),((x_t[0]+x_t[1])/2));
1270phi__t = atan2(((y__t[0]+y__t[1])/2),((x__t[0]+x__t[1])/2));
1271
1272phi_in_t = fmod( phi_in_t+4*M_PI,2*M_PI );
1273phi_out_t = fmod( phi_out_t+4*M_PI,2*M_PI );
1274phi_in__t = fmod( phi_in__t+4*M_PI,2*M_PI );
1275phi_out__t = fmod( phi_out__t+4*M_PI,2*M_PI );
1276phi_t = fmod( phi_t+4*M_PI,2*M_PI );
1277phi__t = fmod( phi__t+4*M_PI,2*M_PI );
1278phi_cellid = fmod( phi_cellid+4*M_PI,2*M_PI );
1279
1280if(fabs(phi_cellid - phi_t)<fabs(phi_cellid - phi__t))
1281{
1282 for(int i =0; i<2; i++ )
1283 { x_p[i] = x_t[i];
1284 y_p[i] = y_t[i];
1285 z_p[i] = z_t[i];}
1286}
1287else if (fabs(phi_cellid - phi_t)>fabs(phi_cellid - phi__t))
1288{
1289 for(int i =0; i<2; i++ )
1290 { x_p[i] = x__t[i];
1291 y_p[i] = y__t[i];
1292 z_p[i] = z__t[i];}
1293}
1294else
1295{
1296 for(int i =0; i<2; i++ )
1297 { x_p[i] = x_t[i];
1298 y_p[i] = y_t[i];
1299 z_p[i] = z_t[i];}
1300}
1301
1302//calculate path length in r_phi plane and all length in this layer
1303double ch_ltrk_rp, ch_ltrk;
1304ch_ltrk_rp = sqrt((x_p[0]-x_p[1])*(x_p[0]-x_p[1])+(y_p[0]-y_p[1])*(y_p[0]-y_p[1]));
1305ch_ltrk = ch_ltrk_rp/sintheta;
1306//give cellid of in and out of this layer
1307double phi_in,phi_out;
1308phi_in = atan2(y_p[0],x_p[0]);
1309phi_out = atan2(y_p[1],x_p[1]);
1310phi_in = fmod( phi_in+4*M_PI,2*M_PI );
1311phi_out = fmod( phi_out+4*M_PI,2*M_PI );
1312
1313//determine the in/out cellid
1314double inlow, inup, outlow, outup, gap, phi1, phi2, phi_mid, phi_midin, phi_midout;
1315int cid_in, cid_out;
1316// int cid_in_t,cid_in__t,cid_out_t,cid_out__t;
1317//cache sampl in each cell for this layer
1318std::vector<double> phi_entrance;
1319std::vector<double> sampl;
1320sampl.resize(ncell);
1321phi_entrance.resize(ncell);
1322for(int k=0; k<ncell; k++) {
1323 sampl[k] = -1.0;
1324 phi_entrance[k] = 0;
1325}
1326
1327cid_in = cid_out = -1;
1328phibin = 2.0*M_PI/ncell;
1329//determine the in/out cell id
1330double stphi[2], phioff[2];
1331stphi[0] = shift*phibin*(0.5-z_p[0]/length);
1332stphi[1] = shift*phibin*(0.5-z_p[1]/length);
1333phioff[0] = phioffset+stphi[0];
1334phioff[1] = phioffset+stphi[1];
1335for(int i=0; i<ncell; i++) {
1336 // for stereo layer
1337 // phi[i] = phioff[0]+phibin*i;
1338 inlow = phioff[0]+phibin*i-phibin*0.5;
1339 inup = phioff[0]+phibin*i+phibin*0.5;
1340 outlow = phioff[1]+phibin*i-phibin*0.5;
1341 outup = phioff[1]+phibin*i+phibin*0.5;
1342 inlow = fmod( inlow+4*M_PI,2*M_PI );
1343 inup = fmod( inup+4*M_PI,2*M_PI );
1344 outlow = fmod( outlow+4*M_PI,2*M_PI );
1345 outup = fmod( outup+4*M_PI,2*M_PI );
1346#ifdef YDEBUG
1347 cout << "shift " << shift<< " cellid " << i
1348 <<" phi_in " << phi_in << " phi_out " << phi_out
1349 << " inup "<< inup << " inlow " << inlow
1350 << " outup "<< outup << " outlow " << outlow << endl;
1351#endif
1352 if(phi_in>=inlow&&phi_in<inup) cid_in = i;
1353 if(phi_out>=outlow&&phi_out<outup) cid_out = i;
1354 if ( inlow>inup) {
1355 if((phi_in>=inlow&&phi_in<2.0*M_PI)||(phi_in>=0.0&&phi_in<inup)) cid_in = i;
1356 }
1357 if ( outlow>outup) {
1358 if((phi_out>=outlow&&phi_out<2.0*M_PI)||(phi_out>=0.0&&phi_out<outup)) cid_out = i;
1359 }
1360}
1361
1362// judge how many cells are traversed and deal with different situations
1363phi_midin = phi_midout = phi1 = phi2 = -999.0;
1364gap = -999.0;
1365//only one cell traversed
1366if( cid_in == cid_out) {
1367 sampl[cid_in] = ch_ltrk;
1368 // return ch_ltrk;
1369
1370} else if(cid_in < cid_out) {
1371 //in cell id less than out cell id
1372 //deal with the special case crossing x axis
1373 if( cid_out-cid_in>ncell/2 ) {
1374 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
1375 phi_midin = phioff[0]+phibin*cid_in-phibin*0.5;
1376 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5;
1377 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1378 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1379 phi1 = phi_midout-phi_out;
1380 phi2 = phi_in-phi_midin;
1381 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1382 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1383 gap = phi1+phi2+(ncell-1-cid_out+cid_in)*phibin;
1384 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1385 sampl[cid_in]=phi2/gap*ch_ltrk;
1386 sampl[cid_out]=phi1/gap*ch_ltrk;
1387
1388 for( int jj = cid_out+1; jj<ncell; jj++) {
1389 sampl[jj]=phibin/gap*ch_ltrk;
1390 }
1391 for( int jj = 0; jj<cid_in; jj++) {
1392 sampl[jj]=phibin/gap*ch_ltrk;
1393 }
1394
1395 } else {
1396 //normal case
1397 phi_midin = phioff[0]+phibin*cid_in+phibin*0.5;
1398 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5;
1399 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1400 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1401 phi1 = phi_midin-phi_in;
1402 phi2 = phi_out-phi_midout;
1403 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1404 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1405 gap = phi1+phi2+(cid_out-cid_in-1)*phibin;
1406 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1407 sampl[cid_in]=phi1/gap*ch_ltrk;
1408 sampl[cid_out]=phi2/gap*ch_ltrk;
1409 phi_entrance[cid_in] = phi_in;
1410 phi_entrance[cid_out] = phi_midout;
1411 for( int jj = cid_in+1; jj<cid_out; jj++) {
1412 sampl[jj]=phibin/gap*ch_ltrk;
1413 }
1414 }
1415} else if(cid_in > cid_out) {
1416 //in cell id greater than out cell id
1417 //deal with the special case crossing x axis
1418 if( cid_in-cid_out>ncell/2 ) {
1419 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
1420 phi_midin = phioff[0]+phibin*cid_in+phibin*0.5;
1421 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5;
1422 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1423 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1424 phi1 = phi_midin-phi_in;
1425 phi2 = phi_out-phi_midout;
1426 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1427 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1428 gap = phi1+phi2+(ncell-1-cid_in+cid_out)*phibin;
1429 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1430 sampl[cid_out]=phi2/gap*ch_ltrk;
1431 sampl[cid_in]=phi1/gap*ch_ltrk;
1432
1433 for( int jj = cid_in+1; jj<ncell; jj++) {
1434 sampl[jj]=phibin/gap*ch_ltrk;
1435 }
1436 for( int jj = 0; jj<cid_out; jj++) {
1437 sampl[jj]=phibin/gap*ch_ltrk;
1438 }
1439 } else {
1440 //normal case
1441 phi_midin = phioff[0]+phibin*cid_in-phibin*0.5;
1442 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5;
1443 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1444 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1445 phi1 = phi_midout-phi_out;
1446 phi2 = phi_in-phi_midin;
1447 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1448 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1449 gap = phi1+phi2+(cid_in-cid_out-1)*phibin;
1450 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1451 sampl[cid_in]=phi2/gap*ch_ltrk;
1452 sampl[cid_out]=phi1/gap*ch_ltrk;
1453 for( int jj = cid_out+1; jj<cid_in; jj++) {
1454 sampl[jj]=phibin/gap*ch_ltrk;
1455 }
1456 }
1457}
1458
1459#ifdef YDEBUG
1460if(sampl[cellid]<0.0) {
1461 if(cid_in!=cid_out) cout<<"?????????"<<endl;
1462 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
1463 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
1464
1465 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
1466 << phi_out << " phi_midout " << phi_midout <<endl;
1467 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
1468 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
1469
1470
1471 for(int l=0; l<ncell; l++) {
1472 if(sampl[l]>=0.0)
1473 cout<<"picked cellid " << l << " sampling length "<< sampl[l]<< endl;
1474 }
1475}
1476#endif
1477return sampl[cellid];
1478
1479}
1480
1481*/
1482
1483
1484// function to calculate the path length in the layer
1485
1486double
1487DedxCorrecSvc::PathL( int ntpFlag, const Dedx_Helix& hel, int layer, int cellid, double z ) const {
1488
1489 double length = geosvc->Layer(layer)->Length();
1490 int genlay = geosvc->Layer(layer)->Gen();
1491 double East_lay_X = geosvc->GeneralLayer(genlay)->SxEast();
1492 double East_lay_Y = geosvc->GeneralLayer(genlay)->SyEast();
1493 double East_lay_Z = geosvc->GeneralLayer(genlay)->SzEast();
1494 double West_lay_X = geosvc->GeneralLayer(genlay)->SxWest();
1495 double West_lay_Y = geosvc->GeneralLayer(genlay)->SyWest();
1496 double West_lay_Z = geosvc->GeneralLayer(genlay)->SzWest();
1497
1498 HepPoint3D East_origin(East_lay_X, East_lay_Y, East_lay_Z);
1499 HepPoint3D West_origin(West_lay_X, West_lay_Y, West_lay_Z);
1500 Hep3Vector wire = (CLHEP::Hep3Vector)East_origin - (CLHEP::Hep3Vector)West_origin;
1501 HepPoint3D piovt_z =(z*10+length/2 )/length * wire + West_origin;
1502 piovt_z = piovt_z*0.1; // conversion of the units(mm=>cm)
1503
1504
1505 //-------------------------------temp -------------------------------//
1506 HepPoint3D piv(hel.pivot());
1507 //-------------------------------temp -------------------------------//
1508
1509 double dr0 = hel.a()[0];
1510 double phi0 = hel.a()[1];
1511 double kappa = hel.a()[2];
1512 double dz0 = hel.a()[3];
1513 double tanl = hel.a()[4];
1514
1515 // Choose the local field !!
1516 double Bz = 1000*(m_pIMF->getReferField());
1517 double ALPHA_loc = 1000/(2.99792458*Bz);
1518
1519 // Choose the local field !!
1520 //double Bz = 1.0; // will be obtained from magnetic field service
1521 //double ALPHA_loc = 1000/(2.99792458*Bz);
1522 //double ALPHA_loc = 333.564095;
1523 int charge = ( kappa >= 0 )? 1 : -1;
1524 double rho = ALPHA_loc/kappa;
1525 double pt = fabs( 1.0/kappa );
1526 double lambda = atan( tanl );
1527 double theta = M_PI_2- lambda;
1528 double sintheta = sin(M_PI_2-atan(tanl));
1529 // theta = 2.0*M_PI*theta/360.;
1530 double phi = fmod(phi0 + M_PI*4, M_PI*2);
1531 double csf0 = cos(phi);
1532 double snf0 = (1. - csf0) * (1. + csf0);
1533 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
1534 if(phi > M_PI) snf0 = - snf0;
1535 //if(ntpFlag>0)
1536 //cout<<"rho = "<<rho<<" hel.dr() + rho = "<<hel.dr() + rho<<endl;
1537
1538 //-------------------------------temp -------------------------------//
1539 double x_c = piv.x() + ( hel.dr() + rho )*csf0;
1540 double y_c = piv.y() + ( hel.dr() + rho )*snf0;
1541 double z_c = piv.z() + hel.dz();
1542 HepPoint3D ccenter(x_c, y_c, 0.0);
1543 double m_c_perp(ccenter.perp());
1544 Hep3Vector m_c_unit((HepPoint3D)ccenter.unit());
1545 //-------------------------------temp -------------------------------//
1546
1547 ////change to boost coordinate system
1548 double x_c_boost = x_c - piovt_z.x();
1549 double y_c_boost = y_c - piovt_z.y();
1550 double z_c_boost = z_c - piovt_z.z();
1551 HepPoint3D ccenter_boost(x_c_boost, y_c_boost, 0.0);
1552 double m_c_perp_boost(ccenter_boost.perp());
1553 //if(ntpFlag>0)
1554 //cout<<" ccenter = "<<ccenter<<" m_c_perp ="<<m_c_perp<<endl;
1555 Hep3Vector m_c_unit_boost((HepPoint3D)ccenter_boost.unit());
1556
1557 //phi region estimation
1558 double phi_io[2];
1559 HepPoint3D IO = (-1)*charge*m_c_unit;
1560 double dphi0 = fmod( IO.phi()+4*M_PI,2*M_PI ) - phi;
1561 double IO_phi = fmod( IO.phi()+4*M_PI,2*M_PI );
1562 //double dphi0_bak = IO_phi - phi;
1563 //if(dphi0 != 0)
1564 //cout<<"dphi value is : "<<dphi0<<" phi:"<<phi<<" IO.phi:"<<IO_phi<<endl;
1565 if( dphi0 > M_PI ) dphi0 -= 2*M_PI;
1566 else if( dphi0 < -M_PI ) dphi0 += 2*M_PI;
1567 //cout<<"dphi value is : "<<dphi0<<" phi:"<<phi<<" IO.phi:"<<IO_phi<<endl;
1568 //cout<<"charge is = "<<charge<<endl;
1569 phi_io[0] = -(1+charge)*M_PI_2 - charge*dphi0;
1570 //phi_io[0] = -(1+charge)*M_PI_2 + charge*dphi0;
1571 phi_io[1] = phi_io[0]+1.5*M_PI;
1572 //cout<<"phi_io[0] is : "<<phi_io[0]<<" phi_io[1]:"<<phi_io[1]<<endl;
1573 double m_crio[2];
1574 double m_zb, m_zf, Calpha;
1575
1576 // retrieve Mdc geometry information
1577 double rcsiz1= geosvc->Layer(layer)->RCSiz1();
1578 double rcsiz2= geosvc->Layer(layer)->RCSiz2();
1579 double rlay= geosvc->Layer(layer)->Radius();
1580 int ncell= geosvc->Layer(layer)->NCell();
1581 double phioffset= geosvc->Layer(layer)->Offset();
1582 float shift= geosvc->Layer(layer)->Shift();
1583 double slant= geosvc->Layer(layer)->Slant();
1584 //double length = geosvc->Layer(layer)->Length();
1585 int type = geosvc->Layer(layer)->Sup()->Type();
1586
1587 ///***********-------------------***************------------------****************//
1588 //check the value if same //
1589 ///***********-------------------***************------------------****************//
1590 int w0id = geosvc->Layer(layer)->Wirst();
1591 int wid = w0id + cellid;
1592 HepPoint3D backkward = geosvc->Wire(wid)->BWirePos();
1593 HepPoint3D forward = geosvc->Wire(wid)->FWirePos();
1594 double x_lay_backward = geosvc->Wire(layer, cellid)->Backward().x();
1595 double y_lay_backward = geosvc->Wire(layer, cellid)->Backward().y();
1596 double x_lay_forward = geosvc->Wire(layer, cellid)->Forward().x();
1597 double y_lay_forward = geosvc->Wire(layer, cellid)->Forward().y();
1598 double r_lay_backward = sqrt(x_lay_backward*x_lay_backward+y_lay_backward*y_lay_backward);
1599 double r_lay_forward = sqrt(x_lay_forward*x_lay_forward+y_lay_forward*y_lay_forward);
1600 double r_lay_use = ((z*10+length/2)/length)*(r_lay_backward-r_lay_forward) + r_lay_forward;
1601 /*for(int i=0; i<43; i++){
1602 double r_up= geosvc->Layer(i)->RCSiz1();
1603 double r_down= geosvc->Layer(i)->RCSiz2();
1604 cout<<i<<" "<<r_up<<" "<<r_down<<endl;
1605 }*/
1606 // shift = shift*type;
1607 // cout<< "type "<< type <<endl;
1608 r_lay_use = 0.1*r_lay_use;
1609 rcsiz1 = 0.1*rcsiz1;
1610 rcsiz2 = 0.1*rcsiz2;
1611 rlay = 0.1*rlay;
1612 length = 0.1*length;
1613 m_zb = 0.5*length;
1614 m_zf = -0.5*length;
1615 m_crio[0] = rlay - rcsiz1;
1616 m_crio[1] = rlay + rcsiz2;
1617 //conversion of the units(mm=>cm)
1618 int sign = -1; ///assumed the first half circle
1619 int epflag[2];
1620 Hep3Vector iocand[2];
1621 Hep3Vector cell_IO[2];
1622 double dphi, downin;
1623 Hep3Vector zvector;
1624 //if (ntpFlag>0) cout<<"z = "<<z<<" , m_zb = "<<m_zb<<" , m_zf = "<<m_zf<<endl;
1625 if( type ) {
1626 downin = (z*z-m_zb*m_zb)*pow(tan(slant),2);
1627 m_crio[0] = sqrt(m_crio[0]*m_crio[0]+downin);
1628 m_crio[1] = sqrt(m_crio[1]*m_crio[1]+downin);
1629 }
1630
1631 //int stced[2];
1632
1633 for( int i = 0; i < 2; i++ ) {
1634 double cos_alpha = m_c_perp_boost*m_c_perp_boost + m_crio[i]*m_crio[i] - rho*rho;
1635 cos_alpha = 0.5*cos_alpha/( m_c_perp_boost*m_crio[i] );
1636 if(fabs(cos_alpha)>1&&i==0) return(-1.0);
1637 if(fabs(cos_alpha)>1&&i==1) {
1638 cos_alpha = m_c_perp_boost*m_c_perp_boost + m_crio[0]*m_crio[0] - rho*rho;
1639 cos_alpha = 0.5*cos_alpha/( m_c_perp_boost*m_crio[0] );
1640 Calpha = 2.0*M_PI-acos( cos_alpha );
1641 } else {
1642 Calpha = acos( cos_alpha );
1643 }
1644 epflag[i] = 0;
1645 iocand[i] = m_c_unit_boost;
1646 iocand[i].rotateZ( charge*sign*Calpha );
1647 iocand[i]*= m_crio[i];
1648 //if (ntpFlag>0) cout<<"iocand corridate is : "<<iocand[i]<<endl;
1649
1650 ///***********-------------------***************------------------****************//
1651 //compare with standard coordinate system results //
1652 ///***********-------------------***************------------------****************//
1653 //-------------------------------temp-------------------------//
1654 iocand[i] = iocand[i]+ piovt_z; //change to standard coordinate system
1655 //iocand[i].y() = iocand[i].y() + piovt_z.y();
1656 //if (ntpFlag>0) cout<<i<<setw(10)<<iocand[i].x()<<setw(10)<<iocand[i].y()<<endl;
1657 //------------------------------temp -------------------------//
1658
1659 double xx = iocand[i].x() - x_c;
1660 double yy = iocand[i].y() - y_c;
1661 //dphi = atan2(yy, xx) - phi0 - M_PI_2*(1+charge);
1662 dphi = atan2(yy, xx) - phi0 - M_PI_2*(1-charge);
1663 dphi = fmod( dphi + 8.0*M_PI, 2*M_PI );
1664
1665 if( dphi < phi_io[0] ) {
1666 dphi += 2*M_PI;
1667 }
1668 else if( phi_io[1] < dphi ) {
1669 dphi -= 2*M_PI;
1670 }
1671 ///in the Local Helix case, phi must be small
1672
1673 //Hep3Vector zvector( 0., 0., z_c-rho*dphi*tanl);
1674 Hep3Vector zvector( 0., 0., z_c-rho*dphi*tanl-piovt_z.z());
1675 //if (ntpFlag>0) cout<<"z_c-rho*dphi*tanl : "<<z_c-rho*dphi*tanl<<endl;
1676 cell_IO[i] = iocand[i];
1677 cell_IO[i] += zvector;
1678 //---------------------------------temp ---------------------------------//
1679 //cell_IO[i] = hel.x(dphi);//compare with above results
1680 //---------------------------------temp ---------------------------------//
1681
1682 double xcio, ycio, phip;
1683 ///// z region check active volume is between zb+2. and zf-2. [cm]
1684 /*
1685 float delrin = 2.0;
1686 if( m_zf-delrin > zvector.z() ){
1687 phip = z_c - m_zb + delrin;
1688 phip = phip/( rho*tanl );
1689 phip = phip + phi0;
1690 xcio = x_c - rho*cos( phip );
1691 ycio = y_c - rho*sin( phip );
1692 cell_IO[i].setX( xcio );
1693 cell_IO[i].setY( ycio );
1694 cell_IO[i].setZ( m_zb+delrin );
1695 epflag[i] = 1;
1696 }
1697 else if( m_zb+delrin < zvector.z() ){
1698 phip = z_c - m_zf-delrin;
1699 phip = phip/( rho*tanl );
1700 phip = phip + phi0;
1701 xcio = x_c - rho*cos( phip );
1702 ycio = y_c - rho*sin( phip );
1703 cell_IO[i].setX( xcio );
1704 cell_IO[i].setY( ycio );
1705 cell_IO[i].setZ( m_zf-delrin );
1706 epflag[i] = 1;
1707 }
1708 else{
1709 */
1710 // cell_IO[i] = iocand;
1711 // cell_IO[i] += zvector;
1712 // }
1713 //dphi = dphi -M_PI;
1714 cell_IO[i] = hel.x(dphi);
1715 //if (ntpFlag>0) { cout<<"cell_IO corridate : "<<cell_IO[i]<<endl;}
1716 // if(i==0) cout<<"zhit value is = "<<z<<endl;
1717 }
1718
1719 // path length estimation
1720 // phi calculation
1721 Hep3Vector cl = cell_IO[1] - cell_IO[0];
1722
1723 //float ch_tphi, ch_tdphi;
1724 double ch_theta;
1725 double ch_dphi;
1726 double ch_ltrk = 0;
1727 double ch_ltrk_rp = 0;
1728 ch_dphi = cl.perp()*0.5/(ALPHA_loc*pt);
1729 ch_dphi = 2.0 * asin( ch_dphi );
1730 ch_ltrk_rp = ch_dphi*ALPHA_loc*pt;
1731 double rpi_path = sqrt(cl.x()*cl.x()+cl.y()*cl.y());
1732 ch_ltrk = sqrt( ch_ltrk_rp*ch_ltrk_rp + cl.z()*cl.z() );
1733 double path = ch_ltrk_rp/ sintheta;
1734 ch_theta = cl.theta();
1735 /* if (ntpFlag>0)
1736 {
1737 // if((ch_ltrk_rp-rpi_path)>0.001 || (ch_ltrk-path)>0.001)
1738 cout<<"ch_ltrk_rp : "<<ch_ltrk_rp<<" rpi_path: "<<rpi_path<<endl;
1739 cout<<"ch_ltrk : "<<ch_ltrk<<" path:"<<path<<endl;
1740 }
1741 */
1742 /*
1743 if( ch_theta < theta*0.85 || 1.15*theta < ch_theta)
1744 ch_ltrk *= -1; /// miss calculation
1745
1746 if( epflag[0] == 1 || epflag [1] == 1 )
1747 ch_ltrk *= -1; /// invalid region for dE/dx or outside of Mdc
1748 */
1749 // judge how many cells are traversed and deal with different situations
1750 double phibin, phi_in, phi_out, inlow, inup, outlow, outup, gap, phi1, phi2, phi_mid, phi_midin, phi_midout;
1751 int cid_in, cid_out;
1752 double inlow_z, inup_z, outlow_z, outup_z, gap_z, phi1_z, phi2_z, phi_mid_z, phi_midin_z, phi_midout_z;
1753 //cache sampl in each cell for this layer
1754
1755 std::vector<double> sampl;
1756 sampl.resize(ncell);
1757 for(int k=0; k<ncell; k++) {
1758 sampl[k] = -1.0;
1759 }
1760
1761 cid_in = cid_out = -1;
1762 phi_in = cell_IO[0].phi();
1763 phi_out = cell_IO[1].phi();
1764 //phi_in = iocand[0].phi();
1765 //phi_out = iocand[1].phi();
1766 phi_in = fmod( phi_in+4*M_PI,2*M_PI );
1767 phi_out = fmod( phi_out+4*M_PI,2*M_PI );
1768 phibin = 2.0*M_PI/ncell;
1769 // gap = fabs(phi_out-phi_in);
1770
1771 //determine the in/out cell id
1772 Hep3Vector cell_mid=0.5*(cell_IO[0]+cell_IO[1]);
1773 //Hep3Vector cell_mid=0.5*(iocand[0]+iocand[0]);
1774 //cout<<cell_mid<<endl;
1775 //double stereophi = shift*phibin*(0.5-cell_mid.z()/length);
1776 //phioffset = phioffset+stereophi;
1777 double stphi[2], phioff[2];
1778 stphi[0] = shift*phibin*(0.5-cell_IO[0].z()/length);
1779 stphi[1] = shift*phibin*(0.5-cell_IO[1].z()/length);
1780 //stphi[0] = shift*phibin*(0.5-z/length);
1781 //stphi[1] = shift*phibin*(0.5-z/length);
1782 phioff[0] = phioffset+stphi[0];
1783 phioff[1] = phioffset+stphi[1];
1784
1785 for(int i=0; i<ncell; i++) {
1786
1787 double x_lay_backward_cell = geosvc->Wire(layer, i)->Backward().x()*0.1;
1788 double y_lay_backward_cell = geosvc->Wire(layer, i)->Backward().y()*0.1;
1789 double x_lay_forward_cell = geosvc->Wire(layer, i)->Forward().x()*0.1;
1790 double y_lay_forward_cell = geosvc->Wire(layer, i)->Forward().y()*0.1;
1791 //if(ntpFlag>0) cout<<layer<<setw(10)<<i<<setw(10)<<x_lay_forward<<setw(10)<<y_lay_forward<<setw(10)<<x_lay_backward<<setw(10)<<y_lay_backward<<setw(10)<<r_lay_forward<<setw(10)<<endl;
1792 //Hep3Vector lay_backward, lay_forward;
1793 Hep3Vector lay_backward(x_lay_backward_cell, y_lay_backward_cell, 0);
1794 Hep3Vector lay_forward(x_lay_forward_cell, y_lay_forward_cell, 0);
1795 Hep3Vector Cell_z[2];
1796 Cell_z[0] = ((cell_IO[0].z()+length/2)/length)*(lay_backward - lay_forward) + lay_forward;
1797 Cell_z[1] = ((cell_IO[1].z()+length/2)/length)*(lay_backward - lay_forward) + lay_forward;
1798 double z_phi[2];
1799 z_phi[0] = Cell_z[0].phi();
1800 z_phi[1] = Cell_z[1].phi();
1801 z_phi[0] = fmod( z_phi[0]+4*M_PI,2*M_PI );
1802 z_phi[1] = fmod( z_phi[1]+4*M_PI,2*M_PI );
1803 //double backward_phi = lay_backward.phi();
1804 //double forward_phi = lay_forward.phi();
1805 //backward_phi = fmod( backward_phi+4*M_PI,2*M_PI );
1806 //forward_phi = fmod( forward_phi+4*M_PI,2*M_PI );
1807 //if(ntpFlag>0) cout<<"backward_phi cal : "<<atan2(y_lay_backward,x_lay_backward)<<" forward_phi : "<<atan2(y_lay_forward,x_lay_forward)<<endl;
1808 //if(ntpFlag>0) cout<<layer<<" backward_phi : "<<backward_phi<<" forward_phi : "<<forward_phi<<endl;
1809
1810 //double z_phi[2];
1811 //z_phi[0] = ((cell_IO[0].z()+length/2)/length)*(backward_phi - forward_phi) + forward_phi;
1812 //z_phi[1] = ((cell_IO[1].z()+length/2)/length)*(backward_phi - forward_phi) + forward_phi;
1813 //if(ntpFlag>0) cout<<"phi z : "<<z_phi[0]<<" "<<z_phi[1]<<endl;
1814 inlow_z = z_phi[0] - phibin*0.5;
1815 inup_z = z_phi[0] + phibin*0.5;
1816 outlow_z = z_phi[1] - phibin*0.5;
1817 outup_z = z_phi[1] + phibin*0.5;
1818 inlow_z = fmod( inlow_z+4*M_PI,2*M_PI );
1819 inup_z = fmod( inup_z+4*M_PI,2*M_PI );
1820 outlow_z = fmod( outlow_z+4*M_PI,2*M_PI );
1821 outup_z = fmod( outup_z+4*M_PI,2*M_PI );
1822
1823
1824 // for stereo layer
1825 inlow = phioff[0]+phibin*i-phibin*0.5;
1826 inup = phioff[0]+phibin*i+phibin*0.5;
1827 outlow = phioff[1]+phibin*i-phibin*0.5;
1828 outup = phioff[1]+phibin*i+phibin*0.5;
1829 inlow = fmod( inlow+4*M_PI,2*M_PI );
1830 inup = fmod( inup+4*M_PI,2*M_PI );
1831 outlow = fmod( outlow+4*M_PI,2*M_PI );
1832 outup = fmod( outup+4*M_PI,2*M_PI );
1833
1834#ifdef YDEBUG
1835 if(ntpFlag > 0) cout << "shift " << shift
1836 <<" phi_in " << phi_in << " phi_out " << phi_out
1837 << " inup "<< inup << " inlow " << inlow
1838 << " outup "<< outup << " outlow " << outlow << endl;
1839#endif
1840
1841 /*if(phi_in>=inlow&&phi_in<inup) cid_in = i;
1842 if(phi_out>=outlow&&phi_out<outup) cid_out = i;
1843 if ( inlow>inup) {
1844 if((phi_in>=inlow&&phi_in<2.0*M_PI)||(phi_in>=0.0&&phi_in<inup)) cid_in = i;
1845 }
1846 if ( outlow>outup) {
1847 if((phi_out>=outlow&&phi_out<2.0*M_PI)||(phi_out>=0.0&&phi_out<outup)) cid_out = i;
1848 }*/
1849
1850 if(phi_in>=inlow_z&&phi_in<inup_z) cid_in = i;
1851 if(phi_out>=outlow_z&&phi_out<outup_z) cid_out = i;
1852 if ( inlow_z>inup_z) {
1853 if((phi_in>=inlow_z&&phi_in<2.0*M_PI)||(phi_in>=0.0&&phi_in<inup_z)) cid_in = i;
1854 }
1855 if ( outlow_z>outup_z) {
1856 if((phi_out>=outlow_z&&phi_out<2.0*M_PI)||(phi_out>=0.0&&phi_out<outup_z)) cid_out = i;
1857 }
1858 }
1859
1860 phi_midin = phi_midout = phi1 = phi2 = -999.0;
1861 gap = -999.0;
1862 //only one cell traversed
1863 //if(ntpFlag > 0) cout<<"cid_in = "<<cid_in<<" cid_out = "<<cid_out<<endl;
1864 if(cid_in == -1 || cid_out == -1) return -1;
1865
1866 if( cid_in == cid_out) {
1867 sampl[cid_in]= ch_ltrk;
1868 //return ch_ltrk;
1869 return sampl[cellid];
1870 } else if(cid_in < cid_out) {
1871 //in cell id less than out cell id
1872 //deal with the special case crossing x axis
1873 if( cid_out-cid_in>ncell/2 ) {
1874
1875 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
1876 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1;
1877 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1;
1878 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1;
1879 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1;
1880 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
1881 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
1882 Hep3Vector Cell_z[2];
1883 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
1884 double phi_cin_z = Cell_z[0].phi();
1885 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI );
1886 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1;
1887 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1;
1888 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1;
1889 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1;
1890 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
1891 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
1892 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
1893 double phi_cout_z = Cell_z[1].phi();
1894 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI );
1895
1896 phi_midin_z = phi_cin_z-phibin*0.5;
1897 phi_midout_z = phi_cout_z+phibin*0.5;
1898 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI );
1899 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI );
1900 phi1_z = phi_midout_z-phi_out;
1901 phi2_z = phi_in-phi_midin_z;
1902 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI);
1903 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI);
1904 gap_z = phi1_z+phi2_z+(ncell-1-cid_out+cid_in)*phibin;
1905 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI);
1906 sampl[cid_in]=phi2_z/gap_z*ch_ltrk;
1907 sampl[cid_out]=phi1_z/gap_z*ch_ltrk;
1908 for( int jj = cid_out+1; jj<ncell; jj++) {
1909 sampl[jj]=phibin/gap_z*ch_ltrk;
1910 }
1911 for( int jj = 0; jj<cid_in; jj++) {
1912 sampl[jj]=phibin/gap_z*ch_ltrk;
1913 }
1914
1915 /*
1916 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
1917 phi_midin = phioff[0]+phibin*cid_in-phibin*0.5;
1918 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5;
1919 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1920 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1921 phi1 = phi_midout-phi_out;
1922 phi2 = phi_in-phi_midin;
1923 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1924 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1925 gap = phi1+phi2+(ncell-1-cid_out+cid_in)*phibin;
1926 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1927 sampl[cid_in]=phi2/gap*ch_ltrk;
1928 sampl[cid_out]=phi1/gap*ch_ltrk;
1929 for( int jj = cid_out+1; jj<ncell; jj++) {
1930 sampl[jj]=phibin/gap*ch_ltrk;
1931 }
1932 for( int jj = 0; jj<cid_in; jj++) {
1933 sampl[jj]=phibin/gap*ch_ltrk;
1934 }*/
1935 /*
1936 cout<<"LLLLLLLLLLLLL"<<endl;
1937 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
1938 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
1939 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
1940 << phi_out << " phi_midout " << phi_midout <<endl;
1941 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
1942 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
1943 */
1944 } else {
1945 //normal case
1946 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1;
1947 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1;
1948 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1;
1949 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1;
1950 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
1951 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
1952 Hep3Vector Cell_z[2];
1953 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
1954 double phi_cin_z = Cell_z[0].phi();
1955 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI );
1956 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1;
1957 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1;
1958 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1;
1959 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1;
1960 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
1961 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
1962 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
1963 double phi_cout_z = Cell_z[1].phi();
1964 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI );
1965
1966 phi_midin_z = phi_cin_z+phibin*0.5;
1967 phi_midout_z = phi_cout_z-phibin*0.5;
1968 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI );
1969 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI );
1970 phi1_z = phi_midin_z-phi_in;
1971 phi2_z = phi_out-phi_midout_z;
1972 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI);
1973 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI);
1974 gap_z = phi1_z+phi2_z+(cid_out-cid_in-1)*phibin;
1975 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI);
1976 sampl[cid_in]=phi1_z/gap_z*ch_ltrk;
1977 sampl[cid_out]=phi2_z/gap_z*ch_ltrk;
1978 for( int jj = cid_in+1; jj<cid_out; jj++) {
1979 sampl[jj]=phibin/gap_z*ch_ltrk;
1980 }
1981 //normal case
1982 /*phi_midin = phioff[0]+phibin*cid_in+phibin*0.5;
1983 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5;
1984 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1985 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1986 phi1 = phi_midin-phi_in;
1987 phi2 = phi_out-phi_midout;
1988 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1989 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1990 gap = phi1+phi2+(cid_out-cid_in-1)*phibin;
1991 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1992 sampl[cid_in]=phi1/gap*ch_ltrk;
1993 sampl[cid_out]=phi2/gap*ch_ltrk;
1994 for( int jj = cid_in+1; jj<cid_out; jj++) {
1995 sampl[jj]=phibin/gap*ch_ltrk;
1996 }*/
1997 }
1998
1999 } else if(cid_in > cid_out) {
2000 //in cell id greater than out cell id
2001 //deal with the special case crossing x axis
2002 if( cid_in-cid_out>ncell/2 ) {
2003 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1;
2004 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1;
2005 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1;
2006 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1;
2007 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
2008 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
2009 Hep3Vector Cell_z[2];
2010 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
2011 double phi_cin_z = Cell_z[0].phi();
2012 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI );
2013 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1;
2014 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1;
2015 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1;
2016 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1;
2017 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
2018 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
2019 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
2020 double phi_cout_z = Cell_z[1].phi();
2021 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI );
2022
2023 phi_midin_z = phi_cin_z+phibin*0.5;
2024 phi_midout_z = phi_cout_z-phibin*0.5;
2025 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI );
2026 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI );
2027 phi1_z = phi_midin_z-phi_in;
2028 phi2_z = phi_out-phi_midout_z;
2029 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI);
2030 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI);
2031 gap_z = phi1_z+phi2_z+(ncell-1-cid_in+cid_out)*phibin;
2032 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI);
2033 sampl[cid_out]=phi2_z/gap_z*ch_ltrk;
2034 sampl[cid_in]=phi1_z/gap_z*ch_ltrk;
2035 for( int jj = cid_in+1; jj<ncell; jj++) {
2036 sampl[jj]=phibin/gap_z*ch_ltrk;
2037 }
2038 for( int jj = 0; jj<cid_out; jj++) {
2039 sampl[jj]=phibin/gap_z*ch_ltrk;
2040 }
2041
2042 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
2043 /*phi_midin = phioff[0]+phibin*cid_in+phibin*0.5;
2044 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5;
2045 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
2046 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
2047 phi1 = phi_midin-phi_in;
2048 phi2 = phi_out-phi_midout;
2049 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
2050 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
2051 gap = phi1+phi2+(ncell-1-cid_in+cid_out)*phibin;
2052 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
2053 sampl[cid_out]=phi2/gap*ch_ltrk;
2054 sampl[cid_in]=phi1/gap*ch_ltrk;
2055 for( int jj = cid_in+1; jj<ncell; jj++) {
2056 sampl[jj]=phibin/gap*ch_ltrk;
2057 }
2058 for( int jj = 0; jj<cid_out; jj++) {
2059 sampl[jj]=phibin/gap*ch_ltrk;
2060 }*/
2061 /*
2062 cout<<"LLLLLLLLLLLLL"<<endl;
2063 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
2064 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
2065 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
2066 << phi_out << " phi_midout " << phi_midout <<endl;
2067 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
2068 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
2069 */
2070 } else {
2071 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1;
2072 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1;
2073 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1;
2074 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1;
2075 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
2076 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
2077 Hep3Vector Cell_z[2];
2078 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
2079 double phi_cin_z = Cell_z[0].phi();
2080 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI );
2081 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1;
2082 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1;
2083 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1;
2084 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1;
2085 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
2086 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
2087 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
2088 double phi_cout_z = Cell_z[1].phi();
2089 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI );
2090
2091 phi_midin_z = phi_cin_z-phibin*0.5;
2092 phi_midout_z = phi_cout_z+phibin*0.5;
2093 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI );
2094 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI );
2095 phi1_z = phi_midout_z-phi_out;
2096 phi2_z = phi_in-phi_midin_z;
2097 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI);
2098 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI);
2099 gap_z = phi1_z+phi2_z+(cid_in-cid_out-1)*phibin;
2100 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI);
2101 sampl[cid_in]=phi2_z/gap_z*ch_ltrk;
2102 sampl[cid_out]=phi1_z/gap_z*ch_ltrk;
2103 for( int jj = cid_out+1; jj<cid_in; jj++) {
2104 sampl[jj]=phibin/gap_z*ch_ltrk;
2105 }
2106
2107 //normal case
2108 /*phi_midin = phioff[0]+phibin*cid_in-phibin*0.5;
2109 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5;
2110 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
2111 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
2112 phi1 = phi_midout-phi_out;
2113 phi2 = phi_in-phi_midin;
2114 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
2115 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
2116 gap = phi1+phi2+(cid_in-cid_out-1)*phibin;
2117 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
2118 sampl[cid_in]=phi2/gap*ch_ltrk;
2119 sampl[cid_out]=phi1/gap*ch_ltrk;
2120 for( int jj = cid_out+1; jj<cid_in; jj++) {
2121 sampl[jj]=phibin/gap*ch_ltrk;
2122 }*/
2123 }
2124 }
2125
2126#ifdef YDEBUG
2127 if(sampl[cellid]<0.0) {
2128 if(cid_in!=cid_out) cout<<"?????????"<<endl;
2129 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
2130 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
2131
2132 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
2133 << phi_out << " phi_midout " << phi_midout <<endl;
2134 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
2135 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
2136
2137
2138 for(int l=0; l<ncell; l++) {
2139 if(sampl[l]>=0.0)
2140 cout<<"picked cellid " << l << " sampling length "<< sampl[l]<< endl;
2141 }
2142 }
2143#endif
2144 return sampl[cellid];
2145}
2146
2147
2148
2149
2150// function to convert measured ionization (D) to actural ionization (I)
2151double DedxCorrecSvc::D2I(const double& cosTheta, const double& D) const
2152{
2153 //cout<<"DedxCorrecSvc::D2I"<<endl;
2154 double absCosTheta = fabs(cosTheta);
2155 double projection = pow(absCosTheta,m_power) + m_delta; // Projected length on wire
2156 double chargeDensity = D/projection;
2157 double numerator = 1 + m_alpha*chargeDensity;
2158 double denominator = 1 + m_gamma*chargeDensity;
2159 double I = D;
2160
2161 //if(denominator > 0.1)
2162
2163 I = D * m_ratio* numerator/denominator;
2164// cout << "m_alpha " << m_alpha << endl;
2165// cout << "m_gamma " << m_gamma << endl;
2166// cout << "m_delta " << m_delta << endl;
2167// cout << "m_power " << m_power << endl;
2168// cout << "m_ratio " << m_ratio << endl;
2169 return I;
2170}
2171
2172// Convert actural ionization (I) to measured ionization (D)
2173double DedxCorrecSvc::I2D(const double& cosTheta, const double& I) const
2174{
2175 // cout<<" DedxCorrecSvc::I2D"<<endl;
2176 double absCosTheta = fabs(cosTheta);
2177 double projection = pow(absCosTheta,m_power) + m_delta; // Projected length on wire
2178
2179 // Do the quadratic to compute d of the electron
2180 double a = m_alpha / projection;
2181 double b = 1 - m_gamma / projection*(I/m_ratio);
2182 double c = -I/m_ratio;
2183
2184 if(b==0 && a==0){
2185 cout<<"both a and b coefficiants for hadron correction are 0"<<endl;
2186 return I;
2187 }
2188
2189 double D = a != 0? (-b + sqrt(b*b - 4.0*a*c))/(2.0*a) : -c/b;
2190 if(D<0){
2191 cout<<"D is less 0! will try another solution"<<endl;
2192 D=a != 0? (-b - sqrt(b*b + 4.0*a*c))/(2.0*a) : -c/b;
2193 if(D<0){
2194 cout<<"D is still less 0! just return uncorrectecd value"<<endl;
2195 return I;
2196 }
2197 }
2198
2199 return D;
2200}
2201
2202
tan
double tan(const BesAngle a)
Definition BesAngle.h:216
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
test
std::string test
Definition CalibModel.cxx:43
n
const Int_t n
Definition DataBase/tau_mode.c:65
phi2
Double_t phi2
Definition DataBase/tau_mode.c:7
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
phi1
Double_t phi1
Definition DataBase/tau_mode.c:7
NumSlicesDoca
#define NumSlicesDoca
Definition DedxCalibParameters.h:62
Iner_DocaMax
#define Iner_DocaMax
Definition DedxCorrParameters.h:52
Out_DocaMin
#define Out_DocaMin
Definition DedxCorrParameters.h:53
Iner_DocaMin
#define Iner_DocaMin
Definition DedxCorrParameters.h:51
Out_DocaMax
#define Out_DocaMax
Definition DedxCorrParameters.h:54
I
const DifComplex I
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition FoamA.h:85
M_PI
#define M_PI
Definition TConstant.h:4
DedxCorrecSvc::CosthetaCorrec
double CosthetaCorrec(double costheta, double ex) const
Definition DedxCorrecSvc.cxx:559
DedxCorrecSvc::D2I
double D2I(const double &cosTheta, const double &D) const
Definition DedxCorrecSvc.cxx:2151
DedxCorrecSvc::f_larcos
double f_larcos(double x, int nbin)
Definition DedxCorrecSvc.cxx:113
DedxCorrecSvc::finalize
virtual StatusCode finalize()
Definition DedxCorrecSvc.cxx:86
DedxCorrecSvc::LayerCorrec
double LayerCorrec(int layer, double z, double costheta, double ex) const
Definition DedxCorrecSvc.cxx:888
DedxCorrecSvc::PathL
double PathL(int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const
Definition DedxCorrecSvc.cxx:1487
DedxCorrecSvc::StandardTrackCorrec
double StandardTrackCorrec(int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, int evtNO, double ex, double costheta, double t0) const
Definition DedxCorrecSvc.cxx:1054
DedxCorrecSvc::handle
void handle(const Incident &)
Definition DedxCorrecSvc.cxx:93
DedxCorrecSvc::T0Correc
double T0Correc(double t0, double ex) const
Definition DedxCorrecSvc.cxx:613
DedxCorrecSvc::LayerGainCorrec
double LayerGainCorrec(int layid, double ex) const
Definition DedxCorrecSvc.cxx:652
DedxCorrecSvc::EntaCorrec
double EntaCorrec(int layid, double enta, double ex) const
Definition DedxCorrecSvc.cxx:711
DedxCorrecSvc::SaturCorrec
double SaturCorrec(int layid, double costheta, double ex) const
Definition DedxCorrecSvc.cxx:539
DedxCorrecSvc::StandardHitCorrec
double StandardHitCorrec(int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, int evtNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const
Definition DedxCorrecSvc.cxx:989
DedxCorrecSvc::CellCorrec
double CellCorrec(int ser, double adc, double dd, double enta, double z, double costheta) const
Definition DedxCorrecSvc.cxx:838
DedxCorrecSvc::initialize
virtual StatusCode initialize()
Definition DedxCorrecSvc.cxx:56
DedxCorrecSvc::WireGainCorrec
double WireGainCorrec(int wireid, double ex) const
Definition DedxCorrecSvc.cxx:527
DedxCorrecSvc::HadronCorrec
double HadronCorrec(double costheta, double dedx) const
Definition DedxCorrecSvc.cxx:645
DedxCorrecSvc::TrkCorrec
double TrkCorrec(double costheta, double dedx) const
Definition DedxCorrecSvc.cxx:900
DedxCorrecSvc::RungCorrec
double RungCorrec(int runNO, int evtNO, double ex) const
Definition DedxCorrecSvc.cxx:666
DedxCorrecSvc::DriftDistCorrec
double DriftDistCorrec(int layid, double ddrift, double ex) const
Definition DedxCorrecSvc.cxx:690
DedxCorrecSvc::DocaSinCorrec
double DocaSinCorrec(int layid, double doca, double enta, double ex) const
Definition DedxCorrecSvc.cxx:781
DedxCorrecSvc::DipAngCorrec
double DipAngCorrec(int layid, double costheta, double ex) const
Definition DedxCorrecSvc.cxx:824
DedxCorrecSvc::ZdepCorrec
double ZdepCorrec(int layid, double zhit, double ex) const
Definition DedxCorrecSvc.cxx:762
DedxCorrecSvc::I2D
double I2D(const double &cosTheta, const double &I) const
Definition DedxCorrecSvc.cxx:2173
DedxCorrecSvc::StandardCorrec
double StandardCorrec(int runFlag, int ntpFlag, int runNO, int evtNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double z, double costheta) const
Definition DedxCorrecSvc.cxx:915
DedxCorrecSvc::GlobalCorrec
double GlobalCorrec(double dedx) const
Definition DedxCorrecSvc.cxx:828
DedxCorrecSvc::set_flag
void set_flag(int calib_F)
Definition DedxCorrecSvc.cxx:1129
Dedx_Helix
Helix parameter class.
Definition Dedx_Helix.h:33
Dedx_Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition Dedx_Helix.h:238
Dedx_Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Dedx_Helix.cxx:209
Dedx_Helix::dz
double dz(void) const
Definition Dedx_Helix.h:220
Dedx_Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Dedx_Helix.h:184
Dedx_Helix::dr
double dr(void) const
returns an element of parameters.
Definition Dedx_Helix.h:202
IMagneticFieldSvc::getReferField
virtual double getReferField()=0
IMdcGeomSvc::Layer
virtual const MdcGeoLayer *const Layer(unsigned id)=0
IMdcGeomSvc::Wire
virtual const MdcGeoWire *const Wire(unsigned id)=0
IMdcGeomSvc::GeneralLayer
virtual const MdcGeoGeneral *const GeneralLayer(unsigned id)=0
MdcGeoGeneral::SzWest
double SzWest(void) const
Definition MdcGeoGeneral.h:220
MdcGeoGeneral::SxEast
double SxEast(void) const
Definition MdcGeoGeneral.h:211
MdcGeoGeneral::SyWest
double SyWest(void) const
Definition MdcGeoGeneral.h:219
MdcGeoGeneral::SxWest
double SxWest(void) const
Definition MdcGeoGeneral.h:218
MdcGeoGeneral::SzEast
double SzEast(void) const
Definition MdcGeoGeneral.h:213
MdcGeoGeneral::SyEast
double SyEast(void) const
Definition MdcGeoGeneral.h:212
MdcGeoLayer::Radius
double Radius(void) const
Definition MdcGeoLayer.h:160
MdcGeoLayer::Slant
double Slant(void) const
Definition MdcGeoLayer.h:158
MdcGeoLayer::Gen
int Gen(void) const
Definition MdcGeoLayer.h:175
MdcGeoLayer::RCSiz2
double RCSiz2(void) const
Definition MdcGeoLayer.h:163
MdcGeoLayer::Length
double Length(void) const
Definition MdcGeoLayer.h:161
MdcGeoLayer::Sup
MdcGeoSuper * Sup(void) const
Definition MdcGeoLayer.h:174
MdcGeoLayer::Shift
double Shift(void) const
Definition MdcGeoLayer.h:167
MdcGeoLayer::RCSiz1
double RCSiz1(void) const
Definition MdcGeoLayer.h:162
MdcGeoLayer::NCell
int NCell(void) const
Definition MdcGeoLayer.h:165
MdcGeoLayer::Wirst
int Wirst(void) const
Definition MdcGeoLayer.h:157
MdcGeoLayer::Offset
double Offset(void) const
Definition MdcGeoLayer.h:166
MdcGeoSuper::Type
int Type(void) const
Definition MdcGeoSuper.h:35
MdcGeoWire::FWirePos
HepPoint3D FWirePos(void) const
Definition MdcGeoWire.h:131
MdcGeoWire::BWirePos
HepPoint3D BWirePos(void) const
Definition MdcGeoWire.h:130
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition MdcGeoWire.h:129
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition MdcGeoWire.h:128
m_debug
bool m_debug
Definition draw_charge_space_coarse_II.cxx:3
costheta
float costheta
Definition prepare_data.cxx:28
charge
float charge
Definition prepare_data.cxx:28
slope
void slope()
Definition slope.cxx:12
b
const double b
Definition slope.cxx:9