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