BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
AbsCor.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2#include "GaudiKernel/AlgFactory.h"
3#include "GaudiKernel/ISvcLocator.h"
4#include "GaudiKernel/SmartDataPtr.h"
5#include "GaudiKernel/IDataProviderSvc.h"
6#include "GaudiKernel/PropertyMgr.h"
7//#
8#include "EventModel/EventModel.h"
9#include "EventModel/Event.h"
10#include "EventModel/EventHeader.h"
11
12#include "EvtRecEvent/EvtRecEvent.h"
13#include "EvtRecEvent/EvtRecTrack.h"
14#include "EmcRecEventModel/RecEmcEventModel.h"
15#include "DstEvent/DstEmcShower.h"
16
17
18#include "TMath.h"
19#include "GaudiKernel/INTupleSvc.h"
20#include "GaudiKernel/NTuple.h"
21#include "GaudiKernel/Bootstrap.h"
22#include "GaudiKernel/IHistogramSvc.h"
23#include "CLHEP/Vector/ThreeVector.h"
24#include "CLHEP/Vector/LorentzVector.h"
25#include "CLHEP/Vector/TwoVector.h"
26using CLHEP::Hep3Vector;
27using CLHEP::Hep2Vector;
28using CLHEP::HepLorentzVector;
29#include "CLHEP/Geometry/Point3D.h"
30#ifndef ENABLE_BACKWARDS_COMPATIBILITY
31 typedef HepGeom::Point3D<double> HepPoint3D;
32#endif
33
34#include "EmcRecEventModel/RecEmcEventModel.h"
35#include "EmcRecGeoSvc/EmcRecBarrelGeo.h"
36#include "EmcRecGeoSvc/EmcRecGeoSvc.h"
37#include "Identifier/TofID.h"
38#include "EvTimeEvent/RecEsTime.h"
39
40#include "AbsCor/AbsCor.h"
41#include "Identifier/Identifier.h"
42#include "TGraphErrors.h"
43#include "TGraph2D.h"
44#include "TGraph2DErrors.h"
45#include <vector>
46#include <string>
47#include <fstream>
48#include <strstream>
49#include "EmcShEnCalibSvc/EmcShEnCalibSvc.h"
50
51/////////////////////////////////////////////////////////////////////////////
52//
53double cbeam[56];
54double ai[4];
55AbsCor::AbsCor(const std::string& name, ISvcLocator* pSvcLocator) :
56 Algorithm(name, pSvcLocator) {
57
58 //Declare the properties
59 ntOut = true;
60 declareProperty("NTupleOut",ntOut=false);
61 declareProperty("McCor", mccor=0);
62 declareProperty("EdgeCor", edgecor=0);
63 declareProperty("HotCellMask", hotcellmask=0);
64 declareProperty("UseTof", usetof=1);
65 declareProperty("DoDataCor", dodatacor = 1);
66 declareProperty("DoPi0Cor",dopi0Cor=1);
67 declareProperty("MCUseTof", MCuseTof=1);
68 declareProperty("MCCorUseFunction", MCCorUseFunction=1);
69 declareProperty("IYear", IYear=2009);
70 declareProperty("ifReadDB", m_ifReadDB=false);
71 declareProperty("CorFunparaPath", m_CorFunparaPath="/afs/ihep.ac.cn/bes3/offline/CalibConst/emc/ShEnCalib/7.0.9/evsetTofCorFunctionPar2021psip.txt");
72 declareProperty("DataPathc3ptof", m_DataPathc3ptof="/afs/ihep.ac.cn/bes3/offline/CalibConst/emc/ShEnCalib/7.0.9/c3ptof2021psip.txt");
73
74
75 runFrom = 0;
76 runTo = 0;
77 m_inFlag=false;
78
79}
80
81
82//TGraphErrors *dt = new TGraphErrors();
83// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
84StatusCode AbsCor::initialize(){
85 MsgStream log(msgSvc(), name());
86 log << MSG::INFO << "in initialize()" << endmsg;
87
88 StatusCode status;
89 if(ntOut == true){
90 NTuplePtr nt1(ntupleSvc(), "FILE1/ec");
91 if ( nt1 ) m_tuple1 = nt1;
92 else {
93 m_tuple1 = ntupleSvc()->book ("FILE1/ec", CLID_ColumnWiseTuple, "ks N-Tuple example");
94 if ( m_tuple1 ) {
95 status = m_tuple1->addItem ("ef", m_ef);
96 status = m_tuple1->addItem ("e5", m_e5);
97 status = m_tuple1->addItem ("ec", m_ec);
98 status = m_tuple1->addItem ("ct", m_ct);
99 status = m_tuple1->addItem ("cedge", m_cedge);
100 }
101 else {
102 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
103 return StatusCode::FAILURE;
104 }
105 }
106 }
107
108 m_index = new int*[56];
109 for (int j=0;j<56;j++ ) {
110 m_index[j] = new int[120];
111 for ( int k=0; k<120; k++) {
112 m_index[j][k]=-1;
113 }
114 }
115
116 m_par = new double*[22];
117 for (int j=0;j<22;j++ ) {
118 m_par[j] = new double[11];
119 for ( int k=0; k<11; k++) {
120 m_par[j][k]=-99;
121 }
122 }
123
124 m_parphi = new double*[22];
125 for (int j=0;j<22;j++ ) {
126 m_parphi[j] = new double[5];
127 for ( int k=0; k<5; k++) {
128 m_parphi[j][k]=-99;
129 }
130 }
131
132 ifstream EnParIn;
133 //EnParIn.exceptions( ifstream::failbit | ifstream::badbit );
134 string EnParInFile;
135 EnParInFile = getenv("ABSCORROOT");
136 EnParInFile += "/share/para.dat";
137 EnParIn.open(EnParInFile.c_str());
138 for(int i=0;i<22;i++) {
139
140 EnParIn>>m_par[i][0]>>m_par[i][1]>>m_par[i][2]>>m_par[i][3]>>m_par[i][4]
141 >>m_par[i][5]>>m_par[i][6]>>m_par[i][7]>>m_par[i][8]>>m_par[i][9]
142 >>m_par[i][10];
143 }
144 EnParIn.close();
145
146 ifstream EnParInphi;
147 //EnParInphi.exceptions( ifstream::failbit | ifstream::badbit );
148 string EnParInFilephi = getenv("ABSCORROOT");
149 EnParInFilephi += "/share/paraphi.dat";
150 EnParInphi.open(EnParInFilephi.c_str());
151 for(int i=0;i<22;i++) {
152 EnParInphi>>m_parphi[i][0]>>m_parphi[i][1]>>m_parphi[i][2]>>m_parphi[i][3]>>m_parphi[i][4];
153
154 }
155 EnParInphi.close();
156
157
158 /* liucx
159 string DataPathevse;
160 DataPathevse = getenv("ABSCORROOT");
161 DataPathevse += "/share/barmccorevse10bin.txt";
162 ifstream inpi0;
163 inpi0.exceptions( ifstream::failbit | ifstream::badbit );
164 inpi0.open(DataPathevse.c_str(),ios::in);
165
166 double epoint[11],peak[11],peakerr[11];
167 for(Int_t i=0;i<11;i++){
168 inpi0>>epoint[i];
169 inpi0>>peak[i];
170 inpi0>>peakerr[i];
171 }
172 for(Int_t i=0;i<11;i++){
173 dt->SetPoint(i, epoint[i],peak[i]);
174 dt->SetPointError(i,0,peakerr[i]);
175 }
176 */
177
178 string DataPathcbeam;
179 DataPathcbeam = getenv("ABSCORROOT");
180 DataPathcbeam += "/share/cbeam.txt";
181 ifstream incbeam;
182 incbeam.exceptions( ifstream::failbit | ifstream::badbit );
183 incbeam.open(DataPathcbeam.c_str(),ios::in);
184 for(int i=0; i<56; i++){
185 double tid,peak,peake,sig,sige;
186 incbeam>>tid;
187 incbeam>>peak;
188 incbeam>>peake;
189 incbeam>>sig;
190 incbeam>>sige;
191 cbeam[i]=1.0/peak;
192 }
193 /* by liucx
194 /////////////////////////////
195 string DataPathc3p;
196 DataPathc3p = getenv("ABSCORROOT");
197 DataPathc3p += "/share/c3p.txt";
198
199 string DataPathc3ptof;
200 DataPathc3ptof = getenv("ABSCORROOT");
201 DataPathc3ptof += "/share/c3ptof.txt";
202 cout<<"mccor = "<<mccor<<endl;
203 //DataPathc3ptof = m_EmcShEnCalibSvc -> getPi0CalibFile();
204
205
206 ifstream inc3p;
207 inc3p.exceptions( ifstream::failbit | ifstream::badbit );
208 if(usetof==0)inc3p.open(DataPathc3p.c_str(),ios::in);
209 if(usetof==1)inc3p.open(DataPathc3ptof.c_str(),ios::in);
210 for(int i=0; i<4; i++){
211 double am,ame;
212 inc3p>>am;
213 inc3p>>ame;
214 ai[i]=am;
215 }
216 */
217 string paraPath = getenv("ABSCORROOT");
218 if(paraPath==""){} //cout<<"AbsCor environment not set!";
219
220 // Read energy correction parameters from PhotonCor/McCor
221 if(MCuseTof==1){
222 paraPath += "/share/evsetTof.txt";
223 }
224 if(MCuseTof==0){
225 paraPath += "/share/evset.txt";
226 }
227
228 ifstream in2;
229 in2.open(paraPath.c_str());
230 assert(in2);
231 double energy,thetaid,peak1,peakerr1,res,reserr;
232 dt = new TGraph2DErrors();
233 dtErr = new TGraph2DErrors();
234 //for(int i=0;i<560;i++){
235 for(int i=0;i<1484;i++){ //53*28
236 in2>>energy;
237 in2>>thetaid;
238 in2>>peak1;
239 in2>>peakerr1;
240 in2>>res;
241 in2>>reserr;
242 dt->SetPoint(i,energy,thetaid,peak1);
243 dt->SetPointError(i,0,0,peakerr1);
244 dtErr->SetPoint(i,energy,thetaid,res);
245 dtErr->SetPointError(i,0,0,reserr);
246 if(i<28) e25min[int(thetaid)]=energy;
247 if(i>=1484-28) e25max[int(thetaid)]=energy;
248 // if(i>=560-28) e25max[int(thetaid)]=energy;
249 }
250 in2.close();
251
252 //-------------------------
253 // Suppression of hot crystals
254 // Reading the map from hotcry.txt (Hajime, Jul 2013)
255 for(int ih=0;ih<10;ih++){
256 hrunstart[ih]=-1;
257 hrunend[ih]=-1;
258 hcell[ih]=-1;
259 }
260 int numhots=4; // numbers of hot crystals
261 int dumring,dumphi,dummod,dumid;
262 string HotList;
263 HotList = getenv("ABSCORROOT");
264 HotList += "/share/hotcry.txt";
265 ifstream hotcrys;
266 hotcrys.exceptions( ifstream::failbit | ifstream::badbit );
267 hotcrys.open(HotList.c_str(),ios::in);
268 for(int il=0; il<numhots; il++){
269 hotcrys>>hrunstart[il];
270 hotcrys>>hrunend[il];
271 hotcrys>>hcell[il];
272 hotcrys>>dumring;
273 hotcrys>>dumphi;
274 hotcrys>>dummod;
275 hotcrys>>dumid;
276 }
277 hotcrys.close();
278 //-------------------------
279 /* by liucx
280 ////////////////////to read the parameters of energy correction function//////////
281 string CorFunparaPath;
282
283 CorFunparaPath = getenv("ABSCORROOT");
284
285 // Read energy correction Function parameters from PhotonCor/McCor
286 if(MCuseTof==1){
287 if (IYear==2009) CorFunparaPath += "/share/evsetTofCorFunctionPar2009Jpsi.txt";
288 if (IYear==2012) CorFunparaPath += "/share/evsetTofCorFunctionPar2012Jpsi.txt";
289 if (IYear==2018) CorFunparaPath += "/share/evsetTofCorFunctionPar2018Jpsi.txt";
290 if (IYear==2019) CorFunparaPath += "/share/evsetTofCorFunctionPar2019Jpsi.txt";
291
292 }
293 if(MCuseTof==0){
294 CorFunparaPath += "/share/evsetCorFunctionPar.txt";
295 }
296
297 ifstream in2corfun;
298 in2corfun.open(CorFunparaPath.c_str());
299 assert(in2corfun);
300
301 for(int i=0;i<28;i++){
302 in2corfun>>m_corFunPar[i][0];
303 in2corfun>>m_corFunPar[i][1];
304 in2corfun>>m_corFunPar[i][2];
305 in2corfun>>m_corFunPar[i][3];
306 in2corfun>>m_corFunPar[i][4];
307 in2corfun>>m_corFunPar[i][5];
308 }
309 in2corfun.close();
310 ///////////////////////
311
312*/
313
314 //////// read the shower correction parameters from database////// by liucx
315 ISvcLocator* svcLocator = Gaudi::svcLocator();
316 StatusCode sc = svcLocator->service("EmcShEnCalibSvc", m_EmcShEnCalibSvc);
317
318 if( sc != StatusCode::SUCCESS){
319 std::cout << "can not use EmcShEnCalibSvc in AbsCor" << endl;
320 }
321 else {
322 std::cout<<"getPi0CalibFile() and getSingleGammaCalibFile() still are empty in initialize()"<< std::endl;
323 std::cout<<"will be assigned a value in execute()"<< std::endl;
324 std::cout << "Test EmcShEnCalibSvc in AbsCor initialize getPi0CalibFile()= "
325 << m_EmcShEnCalibSvc -> getPi0CalibFile()
326 << "Test EmcShEnCalibSvc in AbsCor getSingleGammaCalibFile()= "
327 << m_EmcShEnCalibSvc -> getSingleGammaCalibFile()
328 << std::endl;
329 }
330
331
332
333 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
334 return StatusCode::SUCCESS;
335
336}
337
338// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
339StatusCode AbsCor::execute() {
340 MsgStream log(msgSvc(), name());
341 log << MSG::INFO << "in execute()" << endreq;
342
343 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
344 int runNo=eventHeader->runNumber();
345
346 //// by liucx
347 int runNum = runNo;
348 if (runNo<0) runNum=-runNo; // for MC
349 if (runNum>=runFrom&&runNum<=runTo) {
350 m_inFlag=true;
351 } else {
352 m_inFlag=false;
353 }
354
355 if (m_inFlag==false){
356
357 runFrom=m_EmcShEnCalibSvc -> getRunFrom();
358 runTo=m_EmcShEnCalibSvc -> getRunTo();
359
360 ////to read the parameters of energy correction parameters of single gamma and pi0 calibration//////////
361 // Read energy correction Function parameters
362 cout << "current run=" << runNo<<endl;
363 cout << "in AbsCor open getPi0CalibFile()= " << m_EmcShEnCalibSvc -> getPi0CalibFile()<<endl;
364 cout << "open getSingleGammaCalibFile()= " << m_EmcShEnCalibSvc -> getSingleGammaCalibFile()<<endl;
365 cout <<"getRunFrom()="<< runFrom<<endl;
366 cout <<"getRunTo()="<< runTo<<endl;
367
368
369 string CorFunparaPath;
370 if(MCuseTof==1){
371 if (m_ifReadDB){
372 CorFunparaPath = m_EmcShEnCalibSvc -> getSingleGammaCalibFile();
373 } else {
374 CorFunparaPath = m_CorFunparaPath;
375 cout<<"no read database, but using the set:"<<CorFunparaPath<<endl;
376 }
377 }
378 if(MCuseTof==0){
379 CorFunparaPath = getenv("ABSCORROOT");
380 CorFunparaPath += "/share/evsetCorFunctionPar.txt";
381 }
382
383 ifstream in2corfun;
384 in2corfun.open(CorFunparaPath.c_str());
385 assert(in2corfun);
386
387 for(int i=0;i<28;i++){
388 in2corfun>>m_corFunPar[i][0];
389 in2corfun>>m_corFunPar[i][1];
390 in2corfun>>m_corFunPar[i][2];
391 in2corfun>>m_corFunPar[i][3];
392 in2corfun>>m_corFunPar[i][4];
393 in2corfun>>m_corFunPar[i][5];
394 }
395 in2corfun.close();
396
397 ////////read the parameter of pi0 calibration for data
398 string DataPathc3p;
399 DataPathc3p = getenv("ABSCORROOT");
400 DataPathc3p += "/share/c3p.txt";
401
402 string DataPathc3ptof;
403 if (m_ifReadDB){
404 DataPathc3ptof = m_EmcShEnCalibSvc -> getPi0CalibFile();
405 } else {
406 DataPathc3ptof = m_DataPathc3ptof;
407 cout<<"no read database, but using the set:"<<DataPathc3ptof<<endl;
408 }
409
410 ifstream inc3p;
411 inc3p.exceptions( ifstream::failbit | ifstream::badbit );
412 if(usetof==0)inc3p.open(DataPathc3p.c_str(),ios::in);
413 if(usetof==1)inc3p.open(DataPathc3ptof.c_str(),ios::in);
414 for(int i=0; i<4; i++){
415 double am,ame;
416 inc3p>>am;
417 inc3p>>ame;
418 ai[i]=am;
419 }
420 }
421 ////////////////////////end of reading the correction parameters from Svc and file
422
423 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
424 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
425 if(evtRecEvent->totalTracks()>evtRecTrkCol->size())return SUCCESS;
426 if(evtRecEvent->totalTracks()>50)return SUCCESS;
427
428 IEmcRecGeoSvc* iGeoSvc;
429 ISvcLocator* svcLocator = Gaudi::svcLocator();
430 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
431 if(sc!=StatusCode::SUCCESS) {
432 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
433 }
434
435 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
436 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
437 if(!(*itTrk)->isEmcShowerValid())continue;
438 RecEmcShower *emcTrk = (*itTrk)->emcShower();
439// if(emcTrk->energy()<0.0005)continue;
440
441 /* by liucx
442 //the correction is based on e5x5 energy from the version AbsCor-00-00-29.
443 //It's impossible to correct twice.
444 //So the code is commented out.
445
446 int st = emcTrk->status();
447 if(st>10000)continue;
448 emcTrk->setStatus(st+20000);
449 */
450
451
452
453// cout<<"REC after calib EMC status = "<<emcTrk->status()<<endl;
454
455 //if(emcTrk->e5x5()<0.015)continue;
456
457// if(emcTrk->getTime()<10||emcTrk->getTime()>30)continue;
458/*
459 if(emcTrk->e5x5()<0.02){
460 emcTrk->setEnergy(emcTrk->e5x5()*1.1);
461 continue;
462 }
463*/
464/*
465 int intid = emcTrk->cellId();
466 Identifier id=Identifier(intid);
467 int getthetaid = EmcID::theta_module(id);
468 int getmodule = EmcID::barrel_ec(id);
469 if(getmodule==1)getthetaid=getthetaid+6;
470 if(getmodule==2)getthetaid=55-getthetaid;
471*/
472
473 RecEmcID id = Identifier(emcTrk->cellId());
474
475 unsigned int module, ntheta, nphi;
476 module = EmcID::barrel_ec(id);
477 ntheta = EmcID::theta_module(id);
478 nphi = EmcID::phi_module(id);
479
480 // id=EmcID::crystal_id(module,ntheta,nphi);
481
482 unsigned int thetaModule = EmcID::theta_module(id);
483 unsigned int phiModule = EmcID::phi_module(id);
484
485
486 double e5x5=emcTrk->e5x5();
487 double etof=0;
488
489 if(usetof==1 && (*itTrk)->isTofTrackValid()){
490 SmartRefVector<RecTofTrack> recTofTrackVec=(*itTrk)->tofTrack();
491 if(!recTofTrackVec.empty())etof=recTofTrackVec[0]->energy();
492 if(etof>100.)etof=0;
493 }
494
495 double energyC;
496
497 int thetaId;
498 if (module==0||module==2) thetaId = thetaModule;
499 if (module==1 && thetaModule<=21) thetaId = thetaModule + 6;
500 if (module==1 && thetaModule>21) thetaId = 43 - thetaModule + 6;
501 double DthetaId=double(thetaId);
502
503 if(MCuseTof==1){
504 if (thetaId<6) etof=0.0; //in EMC endcap
505 if (MCCorUseFunction==1){
506 energyC=ECorrFunctionMC(e5x5+etof,DthetaId);
507 } else if (MCCorUseFunction==0){
508 energyC=ECorrMC(e5x5+etof,DthetaId);
509 }
510 }
511 if (MCuseTof==0) {
512 if (MCCorUseFunction==1){
513 energyC=ECorrFunctionMC(e5x5,DthetaId);
514 } else if (MCCorUseFunction==0){
515 energyC=ECorrMC(e5x5,DthetaId);
516 }
517 }
518
519
520 /*
521 if(usetof==1){
522 energyC=ECorrMC(e5x5+etof,thetaId);
523 }
524 if (usetof==0) {
525 energyC=emcTrk->energy();
526 }
527 */
528 // double energyC=emcTrk->energy()+etof;
529/*
530 if(energyC<=0||energyC>4.99)continue;
531 double cor = 1.0;
532 if(runNo>0)cor = dt->Eval(energyC);
533 if(cor<0.001)continue;
534// cout<<cor<<endl;
535 double energyCC= energyC/cor;
536 emcTrk->setEnergy(energyCC);
537*/
538
539 double lnE = std::log(energyC);
540
541 if (energyC>1.0) lnE=std::log(1.0); //gamma energy bin from 0.05 to 0.9GeV for pi0 calibration
542 if (energyC<0.05) lnE=std::log(0.05); //set the range(0.05-1.0GeV) of application for pi0 calibration function
543
544 double lnEcor=1.0;
545 if(dopi0Cor==1){
546 if(runNo>0 && dodatacor==1) {
547 lnEcor = ai[0]+ai[1]*lnE+ai[2]*lnE*lnE+ai[3]*lnE*lnE*lnE;
548 }
549 }
550 if(lnEcor<0.5)continue;
551// cout<<lnEcor<<" "<<etof<<endl;
552
553
554//////////////////////////////////////now no using in the following
555 HepPoint3D pos(emcTrk->position());
556
557 // If it is "hot", return "9999" (Hajime, Jul 2013)
558 double enecor=1.;
559 if(hotcellmask==1){
560 for(int ih=0;ih<10;ih++){
561 if(hrunstart[ih]==-1||hrunend[ih]==-1||hcell[ih]==-1)continue;
562 if(abs(runNo)<hrunstart[ih]||abs(runNo)>hrunend[ih])continue;
563 if ( emcTrk->cellId()==hcell[ih] ) {emcTrk->setStatus(9999);}
564 }
565 }
566
567 if(edgecor==1){
568
569 if(module==1) { //barrel
570 unsigned int theModule;
571 if(thetaModule>21) {
572 theModule = 43 - thetaModule;
573 id = EmcID::crystal_id(module,theModule,phiModule);
574 pos.setZ(-pos.z());
575 } else {
576 theModule = thetaModule;
577 }
578
579 double IRShift;
580 if(theModule==21) {
581 IRShift = 0;
582
583 } else if(theModule==20) {
584 IRShift = 2.5;
585
586 } else {
587 IRShift = 5.;
588
589 }
590
591
592 HepPoint3D IR(0,0,IRShift);
593 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
594 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
595 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
596
597 double theta01,theta23,length,d;
598 theta01 = (center01-IR).theta();
599 theta23 = (center23-IR).theta();
600 length = (center01-IR).mag();
601 d = length*tan(theta01-theta23); //length in x direction
602
603 double xIn;
604 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
605 if (xIn < (-d/2) && theModule!=21){
606
607 theModule=theModule+1 ;
608
609 id = EmcID::crystal_id(module,theModule,phiModule);
610 if(theModule==21) {
611 IRShift = 0;
612
613 } else if(theModule==20) {
614 IRShift = 2.5;
615
616 } else {
617 IRShift = 5.;
618
619 }
620 HepPoint3D IR1(0,0,IRShift);
621 IR=IR1;
622 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
623 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
624 theta01 = (center01-IR).theta();
625 theta23 = (center23-IR).theta();
626 length = (center01-IR).mag();
627 d = length*tan(theta01-theta23); //length in x direction
628
629 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
630
631 } else if (xIn > (d/2) && theModule!=0 ){
632
633 theModule=theModule-1 ;
634 id = EmcID::crystal_id(module,theModule,phiModule);
635 if(theModule==21) {
636 IRShift = 0;
637
638 } else if(theModule==20) {
639 IRShift = 2.5;
640
641 } else {
642 IRShift = 5.;
643
644 }
645
646 HepPoint3D IR1(0,0,IRShift);
647 IR=IR1;
648 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
649 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
650 theta01 = (center01-IR).theta();
651 theta23 = (center23-IR).theta();
652 length = (center01-IR).mag();
653 d = length*tan(theta01-theta23); //length in x direction
654
655 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
656 }
657
658 double yIn;
659 // yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
660 // : pos.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
661
662 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
663 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
664
665 if(nphi==0&&yIn>100) yIn=yIn-360;
666 if(nphi==119&&yIn<-100) yIn=yIn+360;
667
668 if(yIn<-1.5-0.15){
669
670 if(nphi!=0) phiModule=phiModule-1;
671 if(nphi==0) phiModule=119;
672 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
673 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
674 if(phiModule==0&&yIn>100) yIn=yIn-360;
675 if(phiModule==119&&yIn<-100) yIn=yIn+360;
676
677 }
678
679 if(yIn>1.5-0.15){
680
681
682 if(nphi!=119) phiModule=phiModule+1;
683 if(nphi==119) phiModule=0;
684
685 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
686 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
687 if(phiModule==0&&yIn>100) yIn=yIn-360;
688 if(phiModule==119&&yIn<-100) yIn=yIn+360;
689 }
690
691 enecor=m_par[theModule][0]
692 +m_par[theModule][1]*xIn
693 +m_par[theModule][2]*xIn*xIn
694 +m_par[theModule][3]*xIn*xIn*xIn
695 +m_par[theModule][4]*xIn*xIn*xIn*xIn
696 +m_par[theModule][5]*xIn*xIn*xIn*xIn*xIn
697 +m_par[theModule][6]*xIn*xIn*xIn*xIn*xIn*xIn
698 // +m_par[theModule][7]*xIn*xIn*xIn*xIn*xIn*xIn*xIn
699 // +m_par[theModule][8]*xIn*xIn*xIn*xIn*xIn*xIn*xIn*xIn
700 +m_par[theModule][7]*yIn
701 +m_par[theModule][8]*yIn*yIn
702 +m_par[theModule][9]*yIn*yIn*yIn
703 +m_par[theModule][10]*yIn*yIn*yIn*yIn;
704 //////////////////////
705 } else{
706 enecor=1;
707 }
708 }
709//////////////////////////////////////now no using in the above
710
711 double energyCC= energyC/(lnEcor*enecor);
712 // cout<<"Trk Level Corr. and Orig. "<<energyCC<<" "<<emcTrk->energy()<<endl;
713 emcTrk->setEnergy(energyCC);
714
715 if(ntOut==true){
716 m_ef=energyCC;
717 m_e5=emcTrk->e5x5();
718 m_ec=energyC;
719 m_ct=lnEcor ;
720 m_cedge=enecor;
721 m_tuple1->write();
722 }
723
724
725 }
726//==============Modify Dst===============================================================
727 SmartDataPtr<RecEmcShowerCol> recEmcShowerCol(eventSvc(), EventModel::Recon::RecEmcShowerCol);
728 if(!recEmcShowerCol){
729 return( StatusCode::SUCCESS);
730 }
731 SmartDataPtr<DstEmcShowerCol> dstEmcShowerCol(eventSvc(), EventModel::Dst::DstEmcShowerCol);
732 if(!dstEmcShowerCol){
733 return( StatusCode::SUCCESS);
734 }
735
736// cout<<"Rec and Dst Size "<<recEmcShowerCol->size()<<" "<<dstEmcShowerCol->size()<<endl;
737 if(recEmcShowerCol->size()!=dstEmcShowerCol->size())return SUCCESS;
738 for(int i=0;i<recEmcShowerCol->size();i++){
739 RecEmcShowerCol::iterator iter_rec = recEmcShowerCol->begin()+i;
740 DstEmcShowerCol::iterator iter_dst = dstEmcShowerCol->begin()+i;
741// cout<<"Rec and Dst energy "<<(*iter_rec)->energy()<<" "<<(*iter_dst)->energy()<<endl;
742 (*iter_dst)->setEnergy((*iter_rec)->energy());
743 (*iter_dst)->setStatus((*iter_rec)->status());
744// cout<<"DST after calib EMC status = "<<(*iter_dst)->status()<<endl;
745// cout<<"Rec == Dst? "<<(*iter_rec)->energy()<<" "<<(*iter_dst)->energy()<<endl;
746 }
747 return( StatusCode::SUCCESS);
748}
749
750// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
751StatusCode AbsCor::finalize() {
752
753 MsgStream log(msgSvc(), name());
754 log << MSG::INFO << "in finalize()" << endmsg;
755 return StatusCode::SUCCESS;
756}
757
758
759
760double AbsCor::E25min(int n) const
761{
762 return e25min[n];
763}
764double AbsCor::E25max(int n) const
765{
766 return e25max[n];
767}
768
769
770//The following function is copied from PhotonCor/McCor
771double AbsCor::ECorrMC(double eg, double theid) const
772{
773 double Energy5x5=eg;
774 if(eg<E25min(int(theid))) eg=E25min(int(theid));
775 if(eg>E25max(int(theid))) eg=E25max(int(theid))-0.001;
776
777 if(theid<=0)theid=0.001;
778 if(theid>=27)theid=26.999;
779 Float_t einter = eg + 0.00001;
780 Float_t tinter = theid+0.0001;
781 //cout<<"inter="<< einter<<" "<<tinter<<endl;
782 double ecor=dt->Interpolate(einter,tinter);
783 // cout<<"ecor="<<ecor<<endl;
784 if(!(ecor))return Energy5x5;
785 if(ecor<0.5)return Energy5x5;
786 double EnergyCor=Energy5x5/ecor;
787 return EnergyCor;
788}
789
790// Get energy error
791double AbsCor::ErrMC(double eg, double theid) const
792{
793
794 if(eg<E25min(int(theid))) eg=E25min(int(theid));
795 if(eg>E25max(int(theid))) eg=E25max(int(theid))-0.001;
796 if(theid<=0)theid=0.001;
797 if(theid>=27)theid=26.999;
798 Float_t einter = eg + 0.00001;
799 Float_t tinter = theid+0.0001;
800 double err=dtErr->Interpolate(einter,tinter);
801 return err;
802}
803
804
805
806double AbsCor::ECorrFunctionMC(double eg, double theid) const
807{
808 if(theid<0||theid>27){
809 cout<<"in AbsCor EcorrFunctionMC error::thetaId is out of the region [0,27]"<<endl;
810 }
811 double Energy5x5=eg;
812 double x=Energy5x5;
813 //if(Energy5x5>E25max(int(theid))) x=E25max(int(theid));
814
815 int ith=int(theid);
816 double ecor;
817 ecor = m_corFunPar[ith][0]/(m_corFunPar[ith][1]+exp(-x))
818 + m_corFunPar[ith][2]
819 + m_corFunPar[ith][3]*x
820 + m_corFunPar[ith][4]*x*x
821 + m_corFunPar[ith][5]*x*x*x;
822
823 double EnergyCor=Energy5x5/ecor;
824
825 return EnergyCor;
826
827}
828
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: AbsCor.cxx:31
ai
double ai[4]
Definition: AbsCor.cxx:54
cbeam
double cbeam[56]
Definition: AbsCor.cxx:53
tan
double tan(const BesAngle a)
Definition: BesAngle.h:216
runNo
int runNo
Definition: DQA_TO_DB.cxx:12
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
energy
************Class m_ypar INTEGER m_KeyWgt INTEGER m_KeyIHVP INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
Definition: KK2f.h:50
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
exp
EvtComplex exp(const EvtComplex &c)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtComplex.hh:252
AbsCor::execute
StatusCode execute()
Definition: AbsCor.cxx:339
AbsCor::AbsCor
AbsCor(const std::string &name, ISvcLocator *pSvcLocator)
Definition: AbsCor.cxx:55
AbsCor::finalize
StatusCode finalize()
Definition: AbsCor.cxx:751
AbsCor::initialize
StatusCode initialize()
Definition: AbsCor.cxx:84
DstEmcShower::cellId
int cellId() const
Definition: DstEmcShower.h:32
DstEmcShower::position
HepPoint3D position() const
Definition: DstEmcShower.h:34
DstEmcShower::e5x5
double e5x5() const
Definition: DstEmcShower.h:49
DstEmcShower::setStatus
void setStatus(int st)
Definition: DstEmcShower.h:59
DstEmcShower::setEnergy
void setEnergy(double e)
Definition: DstEmcShower.h:63
EmcID::crystal_id
static Identifier crystal_id(const unsigned int barrel_ec, const unsigned int theta_module, const unsigned int phi_module)
For a single crystal.
Definition: EmcID.cxx:71
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: EmcID.cxx:38
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition: EmcID.cxx:43
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition: EmcID.cxx:48
IEmcRecGeoSvc::GetCrystalPoint
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) const =0
EventModel::Dst::DstEmcShowerCol
_EXTERN_ std::string DstEmcShowerCol
Definition: EventModel.h:134
EventModel::EvtRec::EvtRecEvent
_EXTERN_ std::string EvtRecEvent
Definition: EventModel.h:116
EventModel::EvtRec::EvtRecTrackCol
_EXTERN_ std::string EvtRecTrackCol
Definition: EventModel.h:117
EventModel::Recon::RecEmcShowerCol
_EXTERN_ std::string RecEmcShowerCol
Definition: EventModel.h:103