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