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