BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcTunningSvc.cc
Go to the documentation of this file.
1#include <math.h>
2#include <iostream>
3#include <iomanip>
4#include <fstream>
5#include <iostream>
6
7#include "MdcTunningSvc.h"
8#include "GaudiKernel/Kernel.h"
9#include "GaudiKernel/IInterface.h"
10#include "GaudiKernel/StatusCode.h"
11#include "GaudiKernel/SvcFactory.h"
12#include "GaudiKernel/MsgStream.h"
13
14#include "GaudiKernel/IIncidentSvc.h"
15#include "GaudiKernel/Incident.h"
16#include "GaudiKernel/IIncidentListener.h"
17
18#include "GaudiKernel/ISvcLocator.h"
19#include "GaudiKernel/Bootstrap.h"
20
21
22#include "EventModel/EventModel.h"
23#include "EventModel/Event.h"
24#include "EventModel/EventHeader.h"
25#include "GaudiKernel/SmartDataPtr.h"
26
27#include "DatabaseSvc/DatabaseSvc.h"
28
29using namespace std;
30
31MdcTunningSvc::MdcTunningSvc( const string& name, ISvcLocator* svcloc) :
32 Service (name, svcloc){
33 m_BesMdcRes=0;
34 // declare properties
35 declareProperty("UseDatabase",m_dbFlag = false);
36 declareProperty("UseEndcapTuning",m_EndcapTuning = 1); // 0:no endcap Tuning , 1: using endcap Tuning
37 declareProperty("EffFile", m_effFile = std::string("no path"));
38 declareProperty("ResFile", m_resFile = std::string("no path"));
39 declareProperty("EffFile_endcap", m_effFile_endcap = std::string("no path"));
40 declareProperty("ResFile_endcap", m_resFile_endcap = std::string("no path"));
41 declareProperty("path_mdc", m_path = std::string("no path"));
42 declareProperty("Host" , host = std::string("bes3db2.ihep.ac.cn"));
43 declareProperty("DbName" , dbName = std::string("offlinedb"));
44 declareProperty("UserName" , userName = std::string("guest"));
45 declareProperty("Password" , password = std::string("guestpass"));
46 declareProperty("SerialNo" , serialNo = 0);
47 declareProperty("fromDB", m_fromDB = true);
48 declareProperty("ParBossVer", m_ParBossVer = std::string("unknown"));
49
50 int no[43]={40,44,48,56,64,72,80,80,76,76,88,88,100,100,112,112,128,128,140,140,160,160,160,160,176,176,176,176,208,208,208,208,240,240,240,240,256,256,256,256,288,288,288};
51
52 for(int i=0;i<43;i++){
53 cellNo[i]=no[i];
54 }
55
56}
57
58MdcTunningSvc::~MdcTunningSvc(){
59 if(m_BesMdcRes) delete m_BesMdcRes;
60}
61
62StatusCode MdcTunningSvc::queryInterface(const InterfaceID& riid, void** ppvInterface){
63 if( IID_IMdcTunningSvc.versionMatch(riid) ){
64 *ppvInterface = static_cast<IMdcTunningSvc*> (this);
65 } else{
66 return Service::queryInterface(riid, ppvInterface);
67 }
68 return StatusCode::SUCCESS;
69}
70
71
72StatusCode MdcTunningSvc::initialize(){
73 MsgStream log(messageService(), name());
74 log << MSG::INFO << "========== MdcTunningSvc::initialize() ==========" << endreq;
75
76 m_ParBossVer=getenv("BES_RELEASE");
77
78 StatusCode sc = Service::initialize();
79 if( sc.isFailure() ) return sc;
80
81 if(m_fromDB)
82 cout << " MdcTunningSvc read from database. " << endl;
83 else if(!m_fromDB)
84 cout << " MdcTunningSvc read from localfile. " << endl;
85
86 if(m_fromDB==true){
87
88 sc = serviceLocator()->service("DatabaseSvc", m_dbsvc, true);
89 if (sc .isFailure() ) {
90 log << MSG::ERROR << " ERROR: unable to find DatabaseSvc " << endreq;
91 return sc;
92 }
93
94 // MYSQL *conn;
95 // char *opt_host_name = "202.122.33.53";
96 // char *opt_user_name = "maqm";
97 // char *opt_password = "maqm_offline";
98 //// unsigned int opt_port_num = 3306;
99 //// char *opt_socket_name = NULL;
100 // char *opt_db_name = "offlinedb";
101 ////conn = mysql_init(NULL);
102 ////unsigned int opt_flags = 0;
103 ////mysql_real_connect(conn, host.c_str(), userName.c_str(), password.c_str(),
104 //// dbName.c_str(), opt_port_num, opt_socket_name, opt_flags);
105 // mysql_real_connect(conn, opt_host_name, opt_user_name, opt_password,
106 // opt_db_name, opt_port_num, opt_socket_name, opt_flags);
107
108
109 IIncidentSvc* incsvc;
110 sc = service("IncidentSvc", incsvc);
111 int priority = 100;
112 if( sc.isSuccess() ){
113 incsvc -> addListener(this, "NewRun", priority);
114 }
115 sc = serviceLocator()->service("EventDataSvc", m_eventSvc, true);
116 if (sc .isFailure() ) {
117 log << MSG::ERROR << " ERROR: unable to find EventDataSvc " << endreq;
118 return sc;
119 }
120 }
121
122 if(m_fromDB!=true){
123 bool initStat = initTuningConst();
124 // if(m_path!=std::string("no path")) setMdcRes(m_path);
125 if(!initStat){
126 cout << "========== MdcTunningSvc::initialize() failure ! ==========" << endl;
127 return StatusCode::FAILURE;
128 }
129 }
130
131 return StatusCode::SUCCESS;
132}
133
134StatusCode MdcTunningSvc::finalize(){
135 MsgStream log(messageService(), name());
136 log << MSG::INFO << "========== MdcTunningSvc::finalize() ==========" << endreq;
137 //// mysql_close(conn);
138
139 return StatusCode::SUCCESS;
140}
141
142void MdcTunningSvc::handle(const Incident& inc){
143 cout << "========== MdcTunningSvc::handle() ==========" << endl;
144
145 MsgStream log( messageService(), name() );
146 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
147
148 if ( inc.type() == "NewRun" ){
149 log << MSG::DEBUG << "NewRun" << endreq;
150 if(m_fromDB==true) {
151 cout << " Start getMdcTuningTableInfo. " << endl;
152 StatusCode sc = getMdcTuningTableInfo();
153 if(sc.isFailure()) {
154 cout << " ERROR: can not get MdcTuning data from the database. " << endl;
155 log << MSG::ERROR << " ERROR: can not get MdcTuning data from the database. " << endreq;
156 exit(1);
157 }
158 }
159 }
160}
161
162bool MdcTunningSvc::initTuningConst(){
163 std::string FilePath = getenv("MDCTUNNINGSVCROOT");
164
165 if(m_effFile==std::string("no path")){
166 if(!m_fromDB) cout << " ERROR: no mdc tuning eff file, please check the input! " << endl;
167 //m_effFile=FilePath+"/dat/mc_eff.dat";
168 //m_effFile = FilePath+"/dat/mc_eff_psipp.dat";
169 return false;
170 }
171 bool setMcEffStat = setMcEff(m_effFile);
172
173 if(m_resFile==std::string("no path")){
174 if(!m_fromDB) cout << " ERROR: no mdc tuning res file, please check the input! " << endl;
175 //m_resFile=FilePath+"/dat/mc_res.dat";
176 //m_resFile = FilePath+"/dat/mc_res_psipp.dat";
177 return false;
178 }
179 //setMcRes2(m_resFile);
180 bool setMcRes3Stat = setMcRes3(m_resFile);
181
182 if(!(setMcRes3Stat&&setMcEffStat))
183 return false;
184
185 return true;
186}
187
188bool MdcTunningSvc::setMcEff(std::string eff_con){
189 int i,j;
190 string line;
191 double lay,bin,expect,hit;
192 //ifstream readMCEff(eff_con.c_str());
193 std::istringstream readMCEff;
194 // (eff_con);
195 if(m_fromDB){
196 readMCEff.str(eff_con);
197 }
198 if(!m_fromDB){
199 ifstream in(eff_con.c_str());
200 //char ch;
201 //string hhh;
202 //while(ift.get(ch)) hhh.putback(ch);
203 //std::cout<<"hhh:"<<hhh<<std::endl;
204 // stringstream strBuf ;
205 istreambuf_iterator<char> iter(in) ;
206 string strCache = string( iter, (istreambuf_iterator<char>()) );
207 readMCEff.str(strCache);
208 //std::cout<<"strCache:"<<strCache<<std::endl;
209 }
210
211
212 ifstream fin(eff_con.c_str());
213 if(!m_fromDB)
214 if(!fin){
215 cout << " ERROR: the mdc tunning eff file " << m_effFile << " not exist, please check the input! " << endl;
216 return false;
217 }
218 // ifstream readMCEff(eff_con);
219 if(!readMCEff.good()){
220 cout << " ERROR: mdc tuning eff file " << m_effFile << " not exist. " << endl;
221 return false;
222 }else{
223 if(!m_fromDB)
224 if(fin)
225 cout << " Open mdc tuning eff file: " << m_effFile << endl;
226 for(i=0;i<43;i++){
227 readMCEff>>lay;
228 getline(readMCEff,line);
229 for(j=0;j<docaNo;j++){
230 readMCEff>>bin>>docaEff[i][j]>>expect>>hit;
231 }
232 readMCEff>>lay;
233 getline(readMCEff,line);
234 for(j=0;j<thetaNo;j++){
235 readMCEff>>bin>>thetaEff[i][j]>>expect>>hit;
236 }
237 readMCEff>>lay;
238 getline(readMCEff,line);
239 for(j=0;j<cellNo[i];j++){
240 readMCEff>>bin>>cellEff[i][j]>>expect>>hit;
241
242 }
243 }
244 for(i=0;i<43;i++){
245 readMCEff>>lay;
246 getline(readMCEff,line);
247 for(j=0;j<docaNo;j++){
248 readMCEff>>bin>>docaEff_2[i][j]>>expect>>hit;
249 }
250 readMCEff>>lay;
251 getline(readMCEff,line);
252 for(j=0;j<thetaNo;j++){
253 readMCEff>>bin>>thetaEff_2[i][j]>>expect>>hit;
254
255 }
256 readMCEff>>lay;
257 getline(readMCEff,line);
258 for(j=0;j<cellNo[i];j++){
259 readMCEff>>bin>>cellEff_2[i][j]>>expect>>hit;
260
261 }
262 }
263 }
264 return true;
265}
266
267bool MdcTunningSvc::setMcRes(){
268
269
270 int i,j;
271 string line;
272 double lay,bin;
273 ifstream readMCRes(m_resFile.c_str());
274 if(!readMCRes.good()){
275 cout << " ERROR: mdc tuning file: " << m_resFile << " not exist. " << endl;
276 return false;
277 }else{
278 cout << " Open mdc tuning file: " << m_resFile << endl;
279 for(i=0;i<43;i++){
280 readMCRes>>lay;
281 getline(readMCRes,line);
282 getline(readMCRes,line);
283 for(j=0;j<docaNo;j++){
284 readMCRes>>bin>>docaRes[i][j][0][0]>>docaRes[i][j][0][1]; //entranceAngle<0
285 }
286 readMCRes>>lay;
287 getline(readMCRes,line);
288 getline(readMCRes,line);
289 for(j=0;j<docaNo;j++){
290 readMCRes>>bin>>docaRes[i][j][1][0]>>docaRes[i][j][1][1]; //entranceAngle>0
291 /*
292 if(i==0||i==42){
293 cout<<"lay "<<i<<" docaNo "<<j<<" <0 mean "<<docaRes[i][j][0][0]<<" sigma "<<docaRes[i][j][0][1]<<" >0 mean "<<docaRes[i][j][1][0]<<" sigma "<<docaRes[i][j][1][1]<<endl;
294 }
295 */
296 }
297 }
298 }
299 return true;
300}
301
302bool MdcTunningSvc::setMcRes2(std::string res_con){
303
304 int i,j;
305 string line;
306 double lay,bin;
307 //ifstream readMCRes(m_resFile.c_str());
308 std::istringstream readMCRes;
309 if(m_fromDB)
310 {
311 readMCRes.str(res_con);
312 }
313 if(!m_fromDB)
314 {
315 ifstream in(res_con.c_str());
316 istreambuf_iterator<char> iter(in) ;
317 string strCache = string( iter, (istreambuf_iterator<char>()) );
318 readMCRes.str(strCache);
319 }
320 if(!readMCRes.good()){
321 cout << " ERROR: mdc tuning file: " << m_resFile << " not exist. " << endl;
322 return false;
323 }else{
324 cout << " MdcTunningSvc::setMcRes2() Open mdc tuning resfile " << m_resFile << endl;
325 for(i=0;i<43;i++){
326 readMCRes>>lay;
327 getline(readMCRes,line);
328 getline(readMCRes,line);
329 for(j=0;j<docaNo;j++){
330 readMCRes>>bin>>docaF[i][j][0]>>docaMean1[i][j][0]>>docaSigma1[i][j][0]>>docaMean2[i][j][0]>>docaSigma2[i][j][0]; //entranceAngle<0
331 }
332 readMCRes>>lay;
333 getline(readMCRes,line);
334 getline(readMCRes,line);
335 for(j=0;j<docaNo;j++){
336 readMCRes>>bin>>docaF[i][j][1]>>docaMean1[i][j][1]>>docaSigma1[i][j][1]>>docaMean2[i][j][1]>>docaSigma2[i][j][1]; //entranceAngle>0
337 }
338 }
339 for(i=0;i<43;i++){
340 readMCRes>>lay;
341 getline(readMCRes,line);
342 getline(readMCRes,line);
343 for(j=0;j<docaNo;j++){
344 readMCRes>>bin>>docaF_2[i][j][0]>>docaMean1_2[i][j][0]>>docaSigma1_2[i][j][0]>>docaMean2_2[i][j][0]>>docaSigma2_2[i][j][0]; //entranceAngle<0
345 }
346 readMCRes>>lay;
347 getline(readMCRes,line);
348 getline(readMCRes,line);
349 for(j=0;j<docaNo;j++){
350 readMCRes>>bin>>docaF_2[i][j][1]>>docaMean1_2[i][j][1]>>docaSigma1_2[i][j][1]>>docaMean2_2[i][j][1]>>docaSigma2_2[i][j][1]; //entranceAngle>0
351 }
352 }
353 }
354 return true;
355}
356
357bool MdcTunningSvc::setMcRes3(std::string res_con){
358
359 int i,j;
360 string line;
361 double lay,bin;
362 //ifstream readMCRes(m_resFile.c_str());
363 std::istringstream readMCRes;
364 if(m_fromDB){
365 readMCRes.str(res_con);
366 }
367 if(!m_fromDB){
368 ifstream in(res_con.c_str());
369 istreambuf_iterator<char> iter(in) ;
370 string strCache = string( iter, (istreambuf_iterator<char>()) );
371 readMCRes.str(strCache);
372 }
373
374 ifstream fin(res_con.c_str());
375
376 if(!m_fromDB)
377 if(!fin){
378 cout << " ERROR: the mdc tunning res file " << m_resFile << " not exist, please check the input!" << endl;
379 return false;
380 }
381 if(!readMCRes.good()){
382 cout << " ERROR: the mdc tuning res file: " << m_resFile << " not exist, please check the input! " << endl;
383 return false;
384 }else{
385 if(!m_fromDB)
386 if(fin)
387 cout << " MdcTunningSvc::setMcRes3() Open mdc tuning resfile: " << m_resFile << endl;
388 for(i=0;i<43;i++){
389 readMCRes>>lay;
390 getline(readMCRes,line);
391 getline(readMCRes,line);
392 for(j=0;j<docaNo;j++){
393 readMCRes>>bin>>docaF[i][j][0]>>docaMean1[i][j][0]>>docaSigma1[i][j][0]>>docaMean2[i][j][0]>>docaSigma2[i][j][0]>>resLargest[i][j][0]>>resSmallest[i][j][0]>>resRatio[i][j][0]; //entranceAngle<0
394 }
395 readMCRes>>lay;
396 getline(readMCRes,line);
397 getline(readMCRes,line);
398 for(j=0;j<docaNo;j++){
399 readMCRes>>bin>>docaF[i][j][1]>>docaMean1[i][j][1]>>docaSigma1[i][j][1]>>docaMean2[i][j][1]>>docaSigma2[i][j][1]>>resLargest[i][j][1]>>resSmallest[i][j][1]>>resRatio[i][j][1]; //entranceAngle>0
400
401 }
402 }
403 for(i=0;i<43;i++){
404 readMCRes>>lay;
405 getline(readMCRes,line);
406 getline(readMCRes,line);
407 for(j=0;j<docaNo;j++){
408 readMCRes>>bin>>docaF_2[i][j][0]>>docaMean1_2[i][j][0]>>docaSigma1_2[i][j][0]>>docaMean2_2[i][j][0]>>docaSigma2_2[i][j][0]>>resLargest_2[i][j][0]>>resSmallest_2[i][j][0]>>resRatio_2[i][j][0]; //entranceAngle<0
409 }
410 readMCRes>>lay;
411 getline(readMCRes,line);
412 getline(readMCRes,line);
413 for(j=0;j<docaNo;j++){
414
415 readMCRes>>bin>>docaF_2[i][j][1]>>docaMean1_2[i][j][1]>>docaSigma1_2[i][j][1]>>docaMean2_2[i][j][1]>>docaSigma2_2[i][j][1]>>resLargest_2[i][j][1]>>resSmallest_2[i][j][1]>>resRatio_2[i][j][1]; //entranceAngle>0
416
417
418 }
419 }
420 }
421
422 return true;
423}
424
425
426
427BesMdcRes* MdcTunningSvc::getMdcRes(){
428 return m_BesMdcRes;
429}
430
431void MdcTunningSvc::setMdcRes(std::string path){
432 if(m_BesMdcRes) delete m_BesMdcRes;
433 m_BesMdcRes = new BesMdcRes(path);
434}
435
436double MdcTunningSvc::NewSig(int layerId,double driftD){
437 int bindD = 0;
438 double mindD = 0.0 ;
439 double maxdD = 9.0 ;
440 int maxbin =8.0 ;
441
442 if((driftD<mindD)||(driftD>maxdD)){
443 bindD = maxbin ;
444 }else {
445 for(int kk = 0; kk < 9; kk++){
446 if((driftD>=(double)kk)&&(driftD<(double)(kk+1))){
447 bindD = kk ; }
448 }
449 }
450
451 double sigma1 = 0 ;
452
453 sigma1 = (m_BesMdcRes -> getD_dD(layerId ,bindD)) ;
454
455 return sigma1;
456}
457
458
459double MdcTunningSvc::DeldriftD(int layerId,double driftD){
460 int bindD = 0;
461 int maxbin =8 ;
462 double mindD = 0.0 ;
463 double maxdD = 9.0 ;
464
465 for(int jj =0;jj<9;jj++){
466 if((driftD<mindD)||(driftD>maxdD)){
467 bindD = maxbin;
468 }else if((driftD>=jj)&&(driftD<(jj+1))){
469 bindD = jj ;
470 }
471 }
472 double y0D = (m_BesMdcRes -> getD_dD(layerId ,bindD)) ;
473 double y1D = (m_BesMdcRes -> getD_dD(layerId, bindD+1)) ;
474 double yD = y0D + (y1D-y0D)*(driftD - bindD); // calculate the data
475 double y0M = (m_BesMdcRes -> getM_dD(layerId ,bindD)) ;
476 double y1M = (m_BesMdcRes -> getM_dD(layerId ,bindD+1));
477 double yM = y0M + (y1M-y0M)*(driftD - bindD); // calculate the mc data
478 double dely = yD - yM ;
479
480
481 /*if((bindD>=1)&&(bindD<=4)){
482 // cout<<"layerId :"<<layerId<<" dritfD :"<<driftD<<" dely*0.618:"<<dely*0.618<<endl;
483 return dely*0.618 ;
484 }else if((bindD>=5)&&(bindD<7)||(bindD ==8)){
485 // cout<<"layerId :"<<layerId<<" dritfD :"<<driftD<<" dely*0.418:"<<dely*0.418<<endl;
486 return dely*0.418 ;
487 }else {
488 // cout<<"layerId :"<<layerId<<" dritfD :"<<driftD<<" dely:"<<dely<<endl;
489 return dely;
490 }*/
491 return dely;
492
493}
494
495double MdcTunningSvc::Delcostta(int layerId,double costta){
496 int binTa = 0;
497 int maxTa = 15;
498 double minCos = -0.8 ;
499 double minCos2 = -0.7 ;
500 double maxCos = 0.8 ;
501
502 for(int ii = 0; ii <16; ii++){
503 if((costta<minCos)||(costta>maxCos)){
504 binTa = maxTa;
505 }else if((costta>=(minCos + ii*0.1))&&(costta<(minCos2 + ii*0.1))){
506 binTa = ii;
507 }
508 }
509
510 double y0D = (m_BesMdcRes -> getD_theta(layerId ,binTa));
511 double y1D = (m_BesMdcRes -> getD_theta(layerId ,binTa+1));
512 double y0M = (m_BesMdcRes -> getM_theta(layerId,binTa));
513 double y1M = (m_BesMdcRes -> getM_theta(layerId,binTa+1));
514
515 double yD[16],yM[16],Del[16];
516 for(int ll =0;ll<16;ll++){
517 yD[ll] = y0D + (y1D - y0D)/0.1*(costta - (minCos + ll*0.1));
518 yM[ll] = y0M + (y1M - y0D)/0.1*(costta - (minCos + ll*0.1));
519 Del[ll] = yD[ll] - yM[ll] ;
520 }
521
522 double delTha = 0;
523
524 if((binTa>=0)&&(binTa<=5)){
525 delTha = Del[binTa] * 0.118 ;
526 }else if((binTa>5)&&(binTa<=10)){
527 delTha = Del[binTa] * 0.518 ;
528 }else if((binTa>10)&&(binTa<=15)){
529 delTha = Del[binTa] * 0.218 ;
530 }
531
532 return delTha ;
533
534}
535
536
537double MdcTunningSvc::GetEff(int layerId,int cellId,double driftD,double cosTheta,int posFlag){
538 driftD=fabs(driftD);
539 if(driftD>12){
540 //std::cout<<"MdcTuningSvc: driftD overflow "<<driftD<<std::endl;
541 driftD=12;
542 }
543 if(posFlag==0)driftD *= -1;
544
545 if(layerId<0 || layerId>42)std::cout<<" MdcTuningSvc:wrong LayerId "<<layerId<<std::endl;
546 if(cellId<0 || cellId>=cellNo[layerId])std::cout<<"MdcTuningSvc:wrong cellId "<<cellId<<std::endl;
547
548 if(fabs(cosTheta)>1){
549 std::cout<<"MdcTuningSvc:wrong coseTheta "<<cosTheta<<std::endl;
550 cosTheta=1;
551 }
552 double eff;
553 int thetaBin=(int)floor((cosTheta+1)*thetaNo/2.);
554 //debug
555 // std::cout<<"cosTheta "<<cosTheta<<" caled "<<(cosTheta+1)*thetaNo/2.<<" floor "<<thetaBin<<std::endl;
556 int docaBin=(int)floor((driftD+12)*docaNo/24.);
557 if(m_EndcapTuning==0)
558 eff = docaEff[layerId][docaBin] * thetaEff[layerId][thetaBin] * cellEff[layerId][cellId];
559 else {
560 if(fabs(cosTheta)<=0.83)
561 eff = docaEff[layerId][docaBin] * thetaEff[layerId][thetaBin] * cellEff[layerId][cellId];
562 else
563 eff = docaEff_2[layerId][docaBin] * thetaEff_2[layerId][thetaBin] * cellEff_2[layerId][cellId];
564
565 }
566 //debug
567 //std::cout<<"tuning svc layer "<<layerId<<"doca"<<docaBin<<" theta "<<thetaBin<<"cell"<<cellId<<" eff "<<eff<<std::endl;
568
569 return eff;
570}
571
572
573
574double MdcTunningSvc::GetRes(int layerId,int cellId,double driftD,double cosTheta,int posFlag,double entranceAngle,double& mean,double& sigma){
575
576 driftD=fabs(driftD);
577 if(driftD>12){
578 //std::cout<<"MdcTuningSvc: driftD overflow "<<driftD<<std::endl;
579 driftD=12;
580 }
581 if(posFlag==0)driftD *= -1;
582
583 if(layerId<0 || layerId>42)std::cout<<" MdcTuningSvc:wrong LayerId "<<layerId<<std::endl;
584 if(cellId<0 || cellId>=cellNo[layerId])std::cout<<"MdcTuningSvc:wrong cellId "<<cellId<<std::endl;
585
586 if(fabs(cosTheta)>1){
587 std::cout<<"MdcTuningSvc:wrong cosTheta "<<cosTheta<<std::endl;
588 cosTheta=1;
589 }
590
591 //// int thetaBin=(int)floor((cosTheta+1)*thetaNo/2.);
592 //debug
593 // std::cout<<"cosTheta "<<cosTheta<<" caled "<<(cosTheta+1)*thetaNo/2.<<" floor "<<thetaBin<<std::endl;
594 int docaBin=(int)floor((driftD+12)*docaNo/24.);
595 if(entranceAngle<0){
596 mean=docaRes[layerId][docaBin][0][0];
597 sigma=docaRes[layerId][docaBin][0][1];
598 }else{
599 mean=docaRes[layerId][docaBin][1][0];
600 sigma=docaRes[layerId][docaBin][1][1];
601 }
602
603 //debug
604 //std::cout<<"tuning svc layer "<<layerId<<" driftD "<<driftD<<" posFlag "<<posFlag<<" doca "<<docaBin<<" theta "<<thetaBin<<" angle "<<entranceAngle<<" mean "<<mean<<" sigma "<<sigma<<std::endl;
605
606 return 1;
607}
608
609double MdcTunningSvc::GetRes2(int layerId,int cellId,double driftD,double cosTheta,int posFlag,double entranceAngle,double& f,double& mean1,double& sigma1,double& mean2,double& sigma2){
610
611 driftD=fabs(driftD);
612 if(driftD>12){
613 //std::cout<<"MdcTuningSvc: driftD overflow "<<driftD<<std::endl;
614 driftD=12;
615 }
616 if(posFlag==0)driftD *= -1;
617
618 if(layerId<0 || layerId>42)std::cout<<" MdcTuningSvc:wrong LayerId "<<layerId<<std::endl;
619 if(cellId<0 || cellId>=cellNo[layerId])std::cout<<"MdcTuningSvc:wrong cellId "<<cellId<<std::endl;
620
621 if(fabs(cosTheta)>1){
622 std::cout<<"MdcTuningSvc:wrong cosTheta "<<cosTheta<<std::endl;
623 cosTheta=1;
624 }
625
626 //// int thetaBin=(int)floor((cosTheta+1)*thetaNo/2.);
627 //debug
628 // std::cout<<"cosTheta "<<cosTheta<<" caled "<<(cosTheta+1)*thetaNo/2.<<" floor "<<thetaBin<<std::endl;
629 int docaBin=(int)floor((driftD+12)*docaNo/24.);
630 if(m_EndcapTuning==0) {
631 if(entranceAngle<0){
632 f=docaF[layerId][docaBin][0];
633 mean1=docaMean1[layerId][docaBin][0];
634 sigma1=docaSigma1[layerId][docaBin][0];
635 mean2=docaMean2[layerId][docaBin][0];
636 sigma2=docaSigma2[layerId][docaBin][0];
637 }else{
638 f=docaF[layerId][docaBin][1];
639 mean1=docaMean1[layerId][docaBin][1];
640 sigma1=docaSigma1[layerId][docaBin][1];
641 mean2=docaMean2[layerId][docaBin][1];
642 sigma2=docaSigma2[layerId][docaBin][1];
643 }
644 }else {
645 if(fabs(cosTheta)<=0.83) {
646 if(entranceAngle<0){
647 f=docaF[layerId][docaBin][0];
648 mean1=docaMean1[layerId][docaBin][0];
649 sigma1=docaSigma1[layerId][docaBin][0];
650 mean2=docaMean2[layerId][docaBin][0];
651 sigma2=docaSigma2[layerId][docaBin][0];
652 }else{
653 f=docaF[layerId][docaBin][1];
654 mean1=docaMean1[layerId][docaBin][1];
655 sigma1=docaSigma1[layerId][docaBin][1];
656 mean2=docaMean2[layerId][docaBin][1];
657 sigma2=docaSigma2[layerId][docaBin][1];
658 }
659 } else {
660 if(entranceAngle<0){
661 f=docaF_2[layerId][docaBin][0];
662 mean1=docaMean1_2[layerId][docaBin][0];
663 sigma1=docaSigma1_2[layerId][docaBin][0];
664 mean2=docaMean2_2[layerId][docaBin][0];
665 sigma2=docaSigma2_2[layerId][docaBin][0];
666 }else{
667 f=docaF_2[layerId][docaBin][1];
668 mean1=docaMean1_2[layerId][docaBin][1];
669 sigma1=docaSigma1_2[layerId][docaBin][1];
670 mean2=docaMean2_2[layerId][docaBin][1];
671 sigma2=docaSigma2_2[layerId][docaBin][1];
672 }
673 }
674 }
675
676 //debug
677 //std::cout<<"tuning svc layer "<<layerId<<" driftD "<<driftD<<" posFlag "<<posFlag<<" doca "<<docaBin<<" theta "<<thetaBin<<" angle "<<entranceAngle<<" f "<<f<<" mean1 "<<mean1<<" sigma1 "<<sigma1<<" mean2 "<<mean2<<" sigma2 "<<sigma2<<std::endl;
678
679 return 1;
680}
681
682double MdcTunningSvc::GetRes3(int layerId,int cellId,double driftD,double cosTheta,int posFlag,double entranceAngle,double& f,double& mean1,double& sigma1,double& mean2,double& sigma2,double& ResLargest,double& ResSmallest,double& ResRatio){
683
684 driftD=fabs(driftD);
685 if(driftD>12){
686 //std::cout<<"MdcTuningSvc: driftD overflow "<<driftD<<std::endl;
687 driftD=12;
688 }
689 if(posFlag==0)driftD *= -1;
690
691 if(layerId<0 || layerId>42)std::cout<<" MdcTuningSvc:wrong LayerId "<<layerId<<std::endl;
692 if(cellId<0 || cellId>=cellNo[layerId])std::cout<<"MdcTuningSvc:wrong cellId "<<cellId<<std::endl;
693
694 if(fabs(cosTheta)>1){
695 std::cout<<"MdcTuningSvc:wrong cosTheta "<<cosTheta<<std::endl;
696 cosTheta=1;
697 }
698
699 ////int thetaBin=(int)floor((cosTheta+1)*thetaNo/2.);
700 //debug
701 // std::cout<<"cosTheta "<<cosTheta<<" caled "<<(cosTheta+1)*thetaNo/2.<<" floor "<<thetaBin<<std::endl;
702 int docaBin=(int)floor((driftD+12)*docaNo/24.);
703 if(m_EndcapTuning==0) {
704 if(entranceAngle<0){
705 f=docaF[layerId][docaBin][0];
706 mean1=docaMean1[layerId][docaBin][0];
707 sigma1=docaSigma1[layerId][docaBin][0];
708 mean2=docaMean2[layerId][docaBin][0];
709 sigma2=docaSigma2[layerId][docaBin][0];
710 ResLargest=resLargest[layerId][docaBin][0];
711 ResSmallest=resSmallest[layerId][docaBin][0];
712 ResRatio=resRatio[layerId][docaBin][0];
713
714 }else{
715 f=docaF[layerId][docaBin][1];
716 mean1=docaMean1[layerId][docaBin][1];
717 sigma1=docaSigma1[layerId][docaBin][1];
718 mean2=docaMean2[layerId][docaBin][1];
719 sigma2=docaSigma2[layerId][docaBin][1];
720 ResLargest=resLargest[layerId][docaBin][1];
721 ResSmallest=resSmallest[layerId][docaBin][1];
722 ResRatio=resRatio[layerId][docaBin][1];
723 }
724 }else {
725 if(fabs(cosTheta)<=0.83) {
726 if(entranceAngle<0){
727 f=docaF[layerId][docaBin][0];
728 mean1=docaMean1[layerId][docaBin][0];
729 sigma1=docaSigma1[layerId][docaBin][0];
730 mean2=docaMean2[layerId][docaBin][0];
731 sigma2=docaSigma2[layerId][docaBin][0];
732 ResLargest=resLargest[layerId][docaBin][0];
733 ResSmallest=resSmallest[layerId][docaBin][0];
734 ResRatio=resRatio[layerId][docaBin][0];
735 }else{
736 f=docaF[layerId][docaBin][1];
737 mean1=docaMean1[layerId][docaBin][1];
738 sigma1=docaSigma1[layerId][docaBin][1];
739 mean2=docaMean2[layerId][docaBin][1];
740 sigma2=docaSigma2[layerId][docaBin][1];
741 ResLargest=resLargest[layerId][docaBin][1];
742 ResSmallest=resSmallest[layerId][docaBin][1];
743 ResRatio=resRatio[layerId][docaBin][1];
744 }
745 } else {
746 if(entranceAngle<0){
747 f=docaF_2[layerId][docaBin][0];
748 mean1=docaMean1_2[layerId][docaBin][0];
749 sigma1=docaSigma1_2[layerId][docaBin][0];
750 mean2=docaMean2_2[layerId][docaBin][0];
751 sigma2=docaSigma2_2[layerId][docaBin][0];
752 ResLargest=resLargest_2[layerId][docaBin][0];
753 ResSmallest=resSmallest_2[layerId][docaBin][0];
754 ResRatio=resRatio_2[layerId][docaBin][0];
755 }else{
756 f=docaF_2[layerId][docaBin][1];
757 mean1=docaMean1_2[layerId][docaBin][1];
758 sigma1=docaSigma1_2[layerId][docaBin][1];
759 mean2=docaMean2_2[layerId][docaBin][1];
760 sigma2=docaSigma2_2[layerId][docaBin][1];
761 ResLargest=resLargest_2[layerId][docaBin][1];
762 ResSmallest=resSmallest_2[layerId][docaBin][1];
763 ResRatio=resRatio_2[layerId][docaBin][1];
764
765 }
766 }
767 }
768
769 //debug
770 //std::cout<<"tuning svc layer "<<layerId<<" driftD "<<driftD<<" posFlag "<<posFlag<<" doca "<<docaBin<<" theta "<<thetaBin<<" angle "<<entranceAngle<<" f "<<f<<" mean1 "<<mean1<<" sigma1 "<<sigma1<<" mean2 "<<mean2<<" sigma2 "<<sigma2<<std::endl;
771 //debug information
772 //std::cout<<"MdcTunningSvc::GetRes3() debug Info."<<endl
773 // <<"layer docaBin thetaBin entranceAngle f mean1 sigma1 mean2 sigma2 largest smallest ratio"<<endl
774 // <<setw(5)<<layerId<<setw(5)<<docaBin<<setw(5)<<thetaBin<<setw(18)<<entranceAngle<<setw(15)<<f<<setw(15)<<mean1<<setw(15)<<sigma1<<setw(15)<<mean2<<setw(15)<<sigma2<<setw(15)<<ResLargest<<setw(15)<<ResSmallest<<ResRatio<<endl;
775
776 return 1;
777}
778
779double MdcTunningSvc::ResvEntr(int layerId,double enterA,int ilr ,double driftD){
780 int bindD =0;
781 int maxbin = 17;
782 int iEntr = 0;
783 int ll = -1;
784 double mindD = -9.;
785 double maxdD = 9.;
786 double sigmaE = 0.0;
787
788 if(enterA < 0){
789 iEntr = 0 ;
790 }else{ iEntr = 1 ;}
791
792 if(ilr == 0 ){
793 driftD = -1.*driftD;
794 }else{ driftD = driftD ;}
795
796 if( (driftD < mindD) || (driftD > maxdD) ){
797 bindD = maxbin ;
798 }else{
799 for(double dd=-9.;dd<9.;dd++){
800 ll++;
801 if( (driftD>= dd ) && (driftD < (dd+1.)) ){
802 bindD = ll ;
803 }
804 }
805 }
806
807 sigmaE = (m_BesMdcRes -> getD_iEntr(layerId,iEntr,bindD) );
808 //cout<<"Svc lay : "<<layerId<<" iEntr : "<<iEntr<<" bindD : "<<bindD<<" sigmaE :"<<sigmaE<<endl;
809 return sigmaE;
810}
811
812double MdcTunningSvc::DelEtr_Sig(int lay,double enterA,int ilr,double driftD){
813 int bindD = 0;
814 int maxbin =17;
815 int iEntr =0;
816 int ll = -1;
817 double mindD = -9.;
818 double maxdD =9.0;
819 //// double delsigE =0;
820
821 if(enterA < 0){
822 iEntr = 0;
823 }else {iEntr = 1;}
824
825 if(ilr == 0 ){
826 driftD = -1.*driftD;
827 }else{ driftD = driftD ;}
828
829 if( (driftD < mindD) || (driftD > maxdD) ){
830 bindD = maxbin;
831 }else {
832 for(double dd =-9.; dd<9.;dd++){
833 ll++;
834 if( (driftD>=dd ) && (driftD < (dd+1.)) ){
835 bindD = ll;
836 dD[bindD] = dd;
837 }
838 }
839 }
840
841 double y0D = (m_BesMdcRes -> getD_iEntr(lay,iEntr,bindD) );
842 double y1D = (m_BesMdcRes -> getD_iEntr(lay,iEntr,bindD+1) );
843 double yD = y0D + (y1D-y0D)*(driftD - dD[bindD]); // calculate the data
844 double y0M = (m_BesMdcRes -> getM_iEntr(lay,iEntr,bindD) );
845 double y1M = (m_BesMdcRes -> getM_iEntr(lay,iEntr,bindD+1));
846 double yM = y0M + (y1M-y0M)*(driftD - dD[bindD]); // calculate the mc data
847 double dely = yD - yM ;
848
849 //cout<<"Svc lay:"<<lay<<" iEntr :"<<iEntr<<" bindD :"<<bindD
850 // <<" dD["<<bindD<<"] : "<<dD[bindD]
851 // <<" y0D : "<<y0D<<" y1D : "<<y1D<<" yD :"<<yD
852 // <<" y0M : "<<y0M<<" y1M : "<<y1M<<" yM :"<<yM<<" dely : "<<dely<<endl;
853
854 return dely;
855
856}
857
858StatusCode MdcTunningSvc::getMdcTuningTableInfo(){
859 MsgStream log(messageService(), name());
860 SmartDataPtr<Event::EventHeader> eventHeader(m_eventSvc,"/Event/EventHeader");
861 int run=eventHeader->runNumber();
862
863 log << MSG::INFO << "MdcTuningSvc::getMdcTuningTableInfo() run =" << run << endreq;
864 if(m_ParBossVer==std::string("unknown"))
865 cout << "MdcTuningSvc::getMdcTuningTableInfo() : ERROR: there is no ParBossVer " << endl;
866 else log << MSG::INFO << "MdcTuningSvc::getMdcTuningTableInfo() : ParBossVer = " << m_ParBossVer << endl;
867
868 run = fabs(run);
869 ////unsigned long *lengths;
870 //// MYSQL_RES *res_set;
871 //// MYSQL_ROW row;
872 int j=0;
873 int cnt=0;
874 DatabaseRecordVector result;
875 for(int i=0;i<1000;i++){
876 char stmt1[200];
877 cnt = i;
878 int run1 = run+i;
879 cout << " ==================== " << endl;
880 //sprintf(stmt1,"select MdcRes,MdcEff,dEdxTuning from MdcTuning where RunFrom <= %d and RunTo >= %d ",run1,run1);
881 sprintf(stmt1,"select MdcRes,MdcEff from MdcTuning where RunFrom <= %d and RunTo >= %d and SftVer = \"%s\"",run1,run1,m_ParBossVer.c_str());
882 cout << " stmt1: " << stmt1 << endl;
883
884 result.clear();
885 int status = m_dbsvc->query("offlinedb",stmt1,result);
886 if(status<0){
887 cout << " ERROR: can not read MdcRes, MdcEff from the MdcTuning table " << endl;
888 log << MSG::ERROR << " ERROR Read MdcRes, MdcEff from the MdcTuning table " << endreq;
889 return StatusCode::FAILURE;
890 }
891
892 if(result.size()>=1){
893
894 //// int aaa = -1;
895 ////aaa = mysql_real_query(conn, stmt1, strlen(stmt1));
896
897 ////std::cout<<" mysql_real_query: "<<aaa<<std::endl;
898 //// if(aaa){
899 // if(mysql_real_query(conn, stmt1, strlen(stmt1)))
900 ////fprintf(stderr, "query error\n");
901 ////return StatusCode::FAILURE;
902 ////}
903 ////res_set = mysql_store_result (conn);
904 ////row = mysql_fetch_row(res_set);
905 ////int row_no = mysql_num_rows(res_set);
906 ////std::cout<<"row_no="<<row_no<<std::endl;
907 ////if (row_no) {}
908 j=1;
909 break;
910 }
911
912 int run2 = run-i;
913 char stmt2[200];
914 //sprintf(stmt2,"select MdcRes,MdcEff.dEdxTuning from MdcTuning where RunFrom <= %d and RunTo >= %d ",run2,run2);
915 sprintf(stmt2,"select MdcRes,MdcEff from MdcTuning where RunFrom <= %d and RunTo >= %d and SftVer = \"%s\"",run2,run2,m_ParBossVer.c_str());
916
917 result.clear();
918 status = m_dbsvc->query("offlinedb",stmt2,result);
919 if(status<0){
920 log << MSG::ERROR << " ERROR Read MdcRes, MdcEff.dEdxTuning from the MdcTuning table " << endreq;
921 return StatusCode::FAILURE;
922 }
923
924
925 if(result.size()>=1){
926 //// mysql_real_query(conn, stmt2, strlen(stmt2));
927 ////res_set = mysql_store_result (conn);
928 ////row = mysql_fetch_row(res_set);
929 ////row_no = mysql_num_rows(res_set);
930 ////std::cout<<"row_no="<<row_no<<std::endl;
931 ////if (row_no) {}
932 j=-1;
933 break;
934 }
935 }
936
937 if(cnt!=0&&cnt!=1000) {
938 log << MSG::INFO << " get MDC tuning data from run " << run + cnt*j << " instead of run " << run<< endreq;
939 //// std::cout<<"get MDC tuning data form run " <<run+cnt*j<<" instread of run"<< run<<std::endl;
940 }
941 cout << " cnt = " << cnt << endl;
942 if(cnt==1000){
943 log << MSG::ERROR << " can not read Data from DB" << endreq;
944 //// mysql_close(conn);
945 return StatusCode::FAILURE;
946 }
947
948 //// mysql_field_seek (res_set, 0);
949
950 //// lengths = mysql_fetch_lengths(res_set);
951 // string fff = std::string(row[0]);
952 // m_BesMdcRes->setMdcRes(fff);
953 ////string ggg = std::string(row[1]);
954 ////string fff = std::string(row[0]);
955
956 //get last row
957 int row = result.size()-1;
958 cout << " row = " << row << endl;
959
960 if(row<0){
961 cout << " ERROR: can not read Data from DB, please check MdcTunningSvc Version. " << endl;
962 return StatusCode::FAILURE;
963 }
964
965 string ggg = result[row]->GetString("MdcEff");
966 string fff = result[row]->GetString("MdcRes");
967
968 log << MSG::DEBUG << " MdcTunning Data: MdcEff: " << ggg << " MdcRes: " << fff << endreq;
969 //cout << " MdcTunning Data: MdcEff: " << ggg << " MdcRes: " << fff << endl;
970
971 bool stRes = setMcRes3(fff);
972 //bool stRes = setMcRes2(fff);
973 bool stEff = setMcEff(ggg);
974
975 if(!stEff) return StatusCode::FAILURE;
976 if(!stRes) return StatusCode::FAILURE;
977 return StatusCode::SUCCESS;
978}
no
const int no
Definition: CheckRung/check/dEdxVsCos.cxx:17
sigma
TTree * sigma
Definition: CurConst_dEdXCal.cxx:6
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
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
IDatabaseSvc::query
virtual int query(const std::string &dbName, const std::string &sql, DatabaseRecordVector &res)=0
MdcTunningSvc::DelEtr_Sig
double DelEtr_Sig(int lay, double enterA, int ilr, double driftD)
Definition: MdcTunningSvc.cc:812
MdcTunningSvc::setMdcRes
void setMdcRes(std::string path)
Definition: MdcTunningSvc.cc:431
MdcTunningSvc::setMcRes2
bool setMcRes2(std::string res_con)
Definition: MdcTunningSvc.cc:302
MdcTunningSvc::setMcRes3
bool setMcRes3(std::string res_con)
Definition: MdcTunningSvc.cc:357
MdcTunningSvc::getMdcTuningTableInfo
StatusCode getMdcTuningTableInfo()
Definition: MdcTunningSvc.cc:858
MdcTunningSvc::Delcostta
double Delcostta(int layerId, double costta)
Definition: MdcTunningSvc.cc:495
MdcTunningSvc::handle
void handle(const Incident &)
Definition: MdcTunningSvc.cc:142
MdcTunningSvc::getMdcRes
BesMdcRes * getMdcRes()
Definition: MdcTunningSvc.cc:427
MdcTunningSvc::ResvEntr
double ResvEntr(int layerId, double enterA, int ilr, double driftD)
Definition: MdcTunningSvc.cc:779
MdcTunningSvc::setMcEff
bool setMcEff(std::string eff_con)
Definition: MdcTunningSvc.cc:188
MdcTunningSvc::NewSig
double NewSig(int layerId, double driftD)
Definition: MdcTunningSvc.cc:436
MdcTunningSvc::DeldriftD
double DeldriftD(int layerId, double driftD)
Definition: MdcTunningSvc.cc:459
MdcTunningSvc::initTuningConst
bool initTuningConst()
Definition: MdcTunningSvc.cc:162
MdcTunningSvc::GetRes
double GetRes(int layerId, int cellId, double driftD, double cosTheta, int posFlag, double entranceAngle, double &mean, double &sigma)
Definition: MdcTunningSvc.cc:574
MdcTunningSvc::GetEff
double GetEff(int layerId, int cellId, double driftD, double cosTheta, int posFlag)
Definition: MdcTunningSvc.cc:537
MdcTunningSvc::MdcTunningSvc
MdcTunningSvc(const std::string &name, ISvcLocator *svcloc)
MdcTunningSvc::queryInterface
virtual StatusCode queryInterface(const InterfaceID &riid, void **ppvUnknown)
Definition: MdcTunningSvc.cc:62
MdcTunningSvc::finalize
virtual StatusCode finalize()
Definition: MdcTunningSvc.cc:134
MdcTunningSvc::GetRes2
double GetRes2(int layerId, int cellId, double driftD, double cosTheta, int posFlag, double entranceAngle, double &f, double &mean1, double &sigma1, double &mean2, double &sigma2)
Definition: MdcTunningSvc.cc:609
MdcTunningSvc::GetRes3
double GetRes3(int layerId, int cellId, double driftD, double cosTheta, int posFlag, double entranceAngle, double &f, double &mean1, double &sigma1, double &mean2, double &sigma2, double &ResLargest, double &ResSmallest, double &ResRatio)
Definition: MdcTunningSvc.cc:682
MdcTunningSvc::initialize
virtual StatusCode initialize()
Definition: MdcTunningSvc.cc:72
sprintf
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")