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