BOSS 7.0.3
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
41
42#include "AbsCor/AbsCor.h"
43#include "Identifier/Identifier.h"
44#include "TGraphErrors.h"
45#include "TGraph2D.h"
46#include "TGraph2DErrors.h"
47#include <vector>
48#include <string>
49#include <fstream>
50#include <strstream>
51
52
53
54
55/////////////////////////////////////////////////////////////////////////////
56//
57double cbeam[56];
58double ai[4];
59AbsCor::AbsCor(const std::string& name, ISvcLocator* pSvcLocator) :
60 Algorithm(name, pSvcLocator) {
61
62 //Declare the properties
63 ntOut = true;
64 declareProperty("NTupleOut",ntOut=false);
65 declareProperty("McCor", mccor=0);
66 declareProperty("EdgeCor", edgecor=0);
67 declareProperty("HotCellMask", hotcellmask=0);
68 declareProperty("UseTof", usetof=1);
69 declareProperty("DoDataCor", dodatacor = 1);
70 declareProperty("DoPi0Cor",dopi0Cor=1);
71 declareProperty("MCUseTof", MCuseTof=1);
72}
73
74
75//TGraphErrors *dt = new TGraphErrors();
76// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
77StatusCode AbsCor::initialize(){
78 MsgStream log(msgSvc(), name());
79 log << MSG::INFO << "in initialize()" << endmsg;
80 StatusCode status;
81 if(ntOut == true){
82 NTuplePtr nt1(ntupleSvc(), "FILE1/ec");
83 if ( nt1 ) m_tuple1 = nt1;
84 else {
85 m_tuple1 = ntupleSvc()->book ("FILE1/ec", CLID_ColumnWiseTuple, "ks N-Tuple example");
86 if ( m_tuple1 ) {
87 status = m_tuple1->addItem ("ef", m_ef);
88 status = m_tuple1->addItem ("e5", m_e5);
89 status = m_tuple1->addItem ("ec", m_ec);
90 status = m_tuple1->addItem ("ct", m_ct);
91 status = m_tuple1->addItem ("cedge", m_cedge);
92 }
93 else {
94 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
95 return StatusCode::FAILURE;
96 }
97 }
98 }
99
100 m_index = new int*[56];
101 for (int j=0;j<56;j++ ) {
102 m_index[j] = new int[120];
103 for ( int k=0; k<120; k++) {
104 m_index[j][k]=-1;
105 }
106 }
107
108 m_par = new double*[22];
109 for (int j=0;j<22;j++ ) {
110 m_par[j] = new double[11];
111 for ( int k=0; k<11; k++) {
112 m_par[j][k]=-99;
113 }
114 }
115
116 m_parphi = new double*[22];
117 for (int j=0;j<22;j++ ) {
118 m_parphi[j] = new double[5];
119 for ( int k=0; k<5; k++) {
120 m_parphi[j][k]=-99;
121 }
122 }
123
124 ifstream EnParIn;
125 //EnParIn.exceptions( ifstream::failbit | ifstream::badbit );
126 string EnParInFile;
127 EnParInFile = getenv("ABSCORROOT");
128 EnParInFile += "/share/para.dat";
129 EnParIn.open(EnParInFile.c_str());
130 for(int i=0;i<22;i++) {
131
132 EnParIn>>m_par[i][0]>>m_par[i][1]>>m_par[i][2]>>m_par[i][3]>>m_par[i][4]
133 >>m_par[i][5]>>m_par[i][6]>>m_par[i][7]>>m_par[i][8]>>m_par[i][9]
134 >>m_par[i][10];
135 }
136 EnParIn.close();
137
138 ifstream EnParInphi;
139 //EnParInphi.exceptions( ifstream::failbit | ifstream::badbit );
140 string EnParInFilephi = getenv("ABSCORROOT");
141 EnParInFilephi += "/share/paraphi.dat";
142 EnParInphi.open(EnParInFilephi.c_str());
143 for(int i=0;i<22;i++) {
144 EnParInphi>>m_parphi[i][0]>>m_parphi[i][1]>>m_parphi[i][2]>>m_parphi[i][3]>>m_parphi[i][4];
145
146 }
147 EnParInphi.close();
148
149
150 /* liucx
151 string DataPathevse;
152 DataPathevse = getenv("ABSCORROOT");
153 DataPathevse += "/share/barmccorevse10bin.txt";
154 ifstream inpi0;
155 inpi0.exceptions( ifstream::failbit | ifstream::badbit );
156 inpi0.open(DataPathevse.c_str(),ios::in);
157
158 double epoint[11],peak[11],peakerr[11];
159 for(Int_t i=0;i<11;i++){
160 inpi0>>epoint[i];
161 inpi0>>peak[i];
162 inpi0>>peakerr[i];
163 }
164 for(Int_t i=0;i<11;i++){
165 dt->SetPoint(i, epoint[i],peak[i]);
166 dt->SetPointError(i,0,peakerr[i]);
167 }
168 */
169
170 string DataPathcbeam;
171 DataPathcbeam = getenv("ABSCORROOT");
172 DataPathcbeam += "/share/cbeam.txt";
173 ifstream incbeam;
174 incbeam.exceptions( ifstream::failbit | ifstream::badbit );
175 incbeam.open(DataPathcbeam.c_str(),ios::in);
176 for(int i=0; i<56; i++){
177 double tid,peak,peake,sig,sige;
178 incbeam>>tid;
179 incbeam>>peak;
180 incbeam>>peake;
181 incbeam>>sig;
182 incbeam>>sige;
183 cbeam[i]=1.0/peak;
184 }
185
186 /////////////////////////////
187 string DataPathc3p;
188 DataPathc3p = getenv("ABSCORROOT");
189
190 string DataPathc3ptof;
191 DataPathc3ptof = getenv("ABSCORROOT");
192
193 cout<<"mccor = "<<mccor<<endl;
194
195 DataPathc3p += "/share/c3p.txt";
196 DataPathc3ptof += "/share/c3ptof.txt";
197
198 ifstream inc3p;
199 inc3p.exceptions( ifstream::failbit | ifstream::badbit );
200 if(usetof==0)inc3p.open(DataPathc3p.c_str(),ios::in);
201 if(usetof==1)inc3p.open(DataPathc3ptof.c_str(),ios::in);
202 for(int i=0; i<4; i++){
203 double am,ame;
204 inc3p>>am;
205 inc3p>>ame;
206 ai[i]=am;
207 }
208
209 string paraPath = getenv("ABSCORROOT");
210 if(paraPath==""){} //cout<<"EmcRec environment not set!";
211
212 // Read energy correction parameters from PhotonCor/McCor
213 if(MCuseTof==1){
214 paraPath += "/share/evsetTof.txt";
215 }
216 if(MCuseTof==0){
217 paraPath += "/share/evset.txt";
218 }
219
220 ifstream in2;
221 in2.open(paraPath.c_str());
222 assert(in2);
223 double energy,thetaid,peak1,peakerr1,res,reserr;
224 dt = new TGraph2DErrors();
225 dtErr = new TGraph2DErrors();
226 //for(int i=0;i<560;i++){
227 for(int i=0;i<1484;i++){ //53*28
228 in2>>energy;
229 in2>>thetaid;
230 in2>>peak1;
231 in2>>peakerr1;
232 in2>>res;
233 in2>>reserr;
234 dt->SetPoint(i,energy,thetaid,peak1);
235 dt->SetPointError(i,0,0,peakerr1);
236 dtErr->SetPoint(i,energy,thetaid,res);
237 dtErr->SetPointError(i,0,0,reserr);
238 if(i<28) e25min[int(thetaid)]=energy;
239 if(i>=1484-28) e25max[int(thetaid)]=energy;
240 // if(i>=560-28) e25max[int(thetaid)]=energy;
241 }
242 in2.close();
243 //-------------------------
244 // Suppression of hot crystals
245 // Reading the map from hotcry.txt (Hajime, Jul 2013)
246 for(int ih=0;ih<10;ih++){
247 hrunstart[ih]=-1;
248 hrunend[ih]=-1;
249 hcell[ih]=-1;
250 }
251 int numhots=4; // numbers of hot crystals
252 int dumring,dumphi,dummod,dumid;
253 string HotList;
254 HotList = getenv("ABSCORROOT");
255 HotList += "/share/hotcry.txt";
256 ifstream hotcrys;
257 hotcrys.exceptions( ifstream::failbit | ifstream::badbit );
258 hotcrys.open(HotList.c_str(),ios::in);
259 for(int il=0; il<numhots; il++){
260 hotcrys>>hrunstart[il];
261 hotcrys>>hrunend[il];
262 hotcrys>>hcell[il];
263 hotcrys>>dumring;
264 hotcrys>>dumphi;
265 hotcrys>>dummod;
266 hotcrys>>dumid;
267 }
268 hotcrys.close();
269 //-------------------------
270
271
272 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
273 return StatusCode::SUCCESS;
274
275
276
277
278}
279
280// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
281StatusCode AbsCor::execute() {
282 MsgStream log(msgSvc(), name());
283 log << MSG::INFO << "in execute()" << endreq;
284// SmartDataPtr<Event::EventH> eventHeader(eventSvc(),"/Event");
285 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
286 int runNo=eventHeader->runNumber();
287
288 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
289 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
290 if(evtRecEvent->totalTracks()>evtRecTrkCol->size())return SUCCESS;
291 if(evtRecEvent->totalTracks()>50)return SUCCESS;
292
293 IEmcRecGeoSvc* iGeoSvc;
294 ISvcLocator* svcLocator = Gaudi::svcLocator();
295 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
296 if(sc!=StatusCode::SUCCESS) {
297 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
298 }
299
300 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
301 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
302 if(!(*itTrk)->isEmcShowerValid())continue;
303 RecEmcShower *emcTrk = (*itTrk)->emcShower();
304// if(emcTrk->energy()<0.0005)continue;
305
306 /* by liucx
307 int st = emcTrk->status();
308 if(st>10000)continue;
309 emcTrk->setStatus(st+20000);
310 */
311
312
313
314// cout<<"REC after calib EMC status = "<<emcTrk->status()<<endl;
315
316 //if(emcTrk->e5x5()<0.015)continue;
317
318// if(emcTrk->getTime()<10||emcTrk->getTime()>30)continue;
319/*
320 if(emcTrk->e5x5()<0.02){
321 emcTrk->setEnergy(emcTrk->e5x5()*1.1);
322 continue;
323 }
324*/
325/*
326 int intid = emcTrk->cellId();
327 Identifier id=Identifier(intid);
328 int getthetaid = EmcID::theta_module(id);
329 int getmodule = EmcID::barrel_ec(id);
330 if(getmodule==1)getthetaid=getthetaid+6;
331 if(getmodule==2)getthetaid=55-getthetaid;
332*/
333
334 RecEmcID id = Identifier(emcTrk->cellId());
335
336 unsigned int module, ntheta, nphi;
337 module = EmcID::barrel_ec(id);
338 ntheta = EmcID::theta_module(id);
339 nphi = EmcID::phi_module(id);
340
341 // id=EmcID::crystal_id(module,ntheta,nphi);
342
343 unsigned int thetaModule = EmcID::theta_module(id);
344 unsigned int phiModule = EmcID::phi_module(id);
345
346
347 double e5x5=emcTrk->e5x5();
348 double etof=0;
349
350 if(usetof==1 && (*itTrk)->isTofTrackValid()){
351 SmartRefVector<RecTofTrack> recTofTrackVec=(*itTrk)->tofTrack();
352 if(!recTofTrackVec.empty())etof=recTofTrackVec[0]->energy();
353 if(etof>100.)etof=0;
354 }
355
356 double energyC;
357
358 int thetaId;
359 if (module==0||module==2) thetaId = thetaModule;
360 if (module==1 && thetaModule<=21) thetaId = thetaModule + 6;
361 if (module==1 && thetaModule>21) thetaId = 43 - thetaModule + 6;
362
363 if(MCuseTof==1){
364 energyC=ECorrMC(e5x5+etof,thetaId);
365
366 }
367 if (MCuseTof==0) {
368 energyC=ECorrMC(e5x5,thetaId);
369
370 }
371
372
373 /*
374 if(usetof==1){
375 energyC=ECorrMC(e5x5+etof,thetaId);
376 }
377 if (usetof==0) {
378 energyC=emcTrk->energy();
379 }
380 */
381 // double energyC=emcTrk->energy()+etof;
382/*
383 if(energyC<=0||energyC>4.99)continue;
384 double cor = 1.0;
385 if(runNo>0)cor = dt->Eval(energyC);
386 if(cor<0.001)continue;
387// cout<<cor<<endl;
388 double energyCC= energyC/cor;
389 emcTrk->setEnergy(energyCC);
390*/
391
392 double lnE = std::log(energyC);
393
394 if (energyC>1.0) lnE=std::log(1.0);
395 if (energyC<0.07) lnE=std::log(0.07);
396
397 double lnEcor=1.0;
398 if(dopi0Cor==1){
399 if(runNo>0 && dodatacor==1) {
400 lnEcor = ai[0]+ai[1]*lnE+ai[2]*lnE*lnE+ai[3]*lnE*lnE*lnE;
401 }
402 }
403 if(lnEcor<0.5)continue;
404// cout<<lnEcor<<" "<<etof<<endl;
405
406
407//////////////////////////////////////now no using in the following
408 HepPoint3D pos(emcTrk->position());
409
410 // If it is "hot", return "9999" (Hajime, Jul 2013)
411 double enecor=1.;
412 if(hotcellmask==1){
413 for(int ih=0;ih<10;ih++){
414 if(hrunstart[ih]==-1||hrunend[ih]==-1||hcell[ih]==-1)continue;
415 if(abs(runNo)<hrunstart[ih]||abs(runNo)>hrunend[ih])continue;
416 if ( emcTrk->cellId()==hcell[ih] ) {emcTrk->setStatus(9999);}
417 }
418 }
419
420 if(edgecor==1){
421
422 if(module==1) { //barrel
423 unsigned int theModule;
424 if(thetaModule>21) {
425 theModule = 43 - thetaModule;
426 id = EmcID::crystal_id(module,theModule,phiModule);
427 pos.setZ(-pos.z());
428 } else {
429 theModule = thetaModule;
430 }
431
432 double IRShift;
433 if(theModule==21) {
434 IRShift = 0;
435
436 } else if(theModule==20) {
437 IRShift = 2.5;
438
439 } else {
440 IRShift = 5.;
441
442 }
443
444
445 HepPoint3D IR(0,0,IRShift);
446 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
447 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
448 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
449
450 double theta01,theta23,length,d;
451 theta01 = (center01-IR).theta();
452 theta23 = (center23-IR).theta();
453 length = (center01-IR).mag();
454 d = length*tan(theta01-theta23); //length in x direction
455
456 double xIn;
457 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
458 if (xIn < (-d/2) && theModule!=21){
459
460 theModule=theModule+1 ;
461
462 id = EmcID::crystal_id(module,theModule,phiModule);
463 if(theModule==21) {
464 IRShift = 0;
465
466 } else if(theModule==20) {
467 IRShift = 2.5;
468
469 } else {
470 IRShift = 5.;
471
472 }
473 HepPoint3D IR1(0,0,IRShift);
474 IR=IR1;
475 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
476 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
477 theta01 = (center01-IR).theta();
478 theta23 = (center23-IR).theta();
479 length = (center01-IR).mag();
480 d = length*tan(theta01-theta23); //length in x direction
481
482 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
483
484 } else if (xIn > (d/2) && theModule!=0 ){
485
486 theModule=theModule-1 ;
487 id = EmcID::crystal_id(module,theModule,phiModule);
488 if(theModule==21) {
489 IRShift = 0;
490
491 } else if(theModule==20) {
492 IRShift = 2.5;
493
494 } else {
495 IRShift = 5.;
496
497 }
498
499 HepPoint3D IR1(0,0,IRShift);
500 IR=IR1;
501 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
502 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
503 theta01 = (center01-IR).theta();
504 theta23 = (center23-IR).theta();
505 length = (center01-IR).mag();
506 d = length*tan(theta01-theta23); //length in x direction
507
508 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
509 }
510
511 double yIn;
512 // yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
513 // : pos.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
514
515 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
516 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
517
518 if(nphi==0&&yIn>100) yIn=yIn-360;
519 if(nphi==119&&yIn<-100) yIn=yIn+360;
520
521 if(yIn<-1.5-0.15){
522
523 if(nphi!=0) phiModule=phiModule-1;
524 if(nphi==0) phiModule=119;
525 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
526 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
527 if(phiModule==0&&yIn>100) yIn=yIn-360;
528 if(phiModule==119&&yIn<-100) yIn=yIn+360;
529
530 }
531
532 if(yIn>1.5-0.15){
533
534
535 if(nphi!=119) phiModule=phiModule+1;
536 if(nphi==119) phiModule=0;
537
538 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
539 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
540 if(phiModule==0&&yIn>100) yIn=yIn-360;
541 if(phiModule==119&&yIn<-100) yIn=yIn+360;
542 }
543
544 enecor=m_par[theModule][0]
545 +m_par[theModule][1]*xIn
546 +m_par[theModule][2]*xIn*xIn
547 +m_par[theModule][3]*xIn*xIn*xIn
548 +m_par[theModule][4]*xIn*xIn*xIn*xIn
549 +m_par[theModule][5]*xIn*xIn*xIn*xIn*xIn
550 +m_par[theModule][6]*xIn*xIn*xIn*xIn*xIn*xIn
551 // +m_par[theModule][7]*xIn*xIn*xIn*xIn*xIn*xIn*xIn
552 // +m_par[theModule][8]*xIn*xIn*xIn*xIn*xIn*xIn*xIn*xIn
553 +m_par[theModule][7]*yIn
554 +m_par[theModule][8]*yIn*yIn
555 +m_par[theModule][9]*yIn*yIn*yIn
556 +m_par[theModule][10]*yIn*yIn*yIn*yIn;
557 //////////////////////
558 } else{
559 enecor=1;
560 }
561 }
562//////////////////////////////////////now no using in the above
563
564 double energyCC= energyC/(lnEcor*enecor);
565 // cout<<"Trk Level Corr. and Orig. "<<energyCC<<" "<<emcTrk->energy()<<endl;
566 emcTrk->setEnergy(energyCC);
567
568 if(ntOut==true){
569 m_ef=energyCC;
570 m_e5=emcTrk->e5x5();
571 m_ec=energyC;
572 m_ct=lnEcor ;
573 m_cedge=enecor;
574 m_tuple1->write();
575 }
576
577
578 }
579//==============Modify Dst===============================================================
580 SmartDataPtr<RecEmcShowerCol> recEmcShowerCol(eventSvc(), EventModel::Recon::RecEmcShowerCol);
581 if(!recEmcShowerCol){
582 return( StatusCode::SUCCESS);
583 }
584 SmartDataPtr<DstEmcShowerCol> dstEmcShowerCol(eventSvc(), EventModel::Dst::DstEmcShowerCol);
585 if(!dstEmcShowerCol){
586 return( StatusCode::SUCCESS);
587 }
588
589// cout<<"Rec and Dst Size "<<recEmcShowerCol->size()<<" "<<dstEmcShowerCol->size()<<endl;
590 if(recEmcShowerCol->size()!=dstEmcShowerCol->size())return SUCCESS;
591 for(int i=0;i<recEmcShowerCol->size();i++){
592 RecEmcShowerCol::iterator iter_rec = recEmcShowerCol->begin()+i;
593 DstEmcShowerCol::iterator iter_dst = dstEmcShowerCol->begin()+i;
594// cout<<"Rec and Dst energy "<<(*iter_rec)->energy()<<" "<<(*iter_dst)->energy()<<endl;
595 (*iter_dst)->setEnergy((*iter_rec)->energy());
596 (*iter_dst)->setStatus((*iter_rec)->status());
597// cout<<"DST after calib EMC status = "<<(*iter_dst)->status()<<endl;
598// cout<<"Rec == Dst? "<<(*iter_rec)->energy()<<" "<<(*iter_dst)->energy()<<endl;
599 }
600 return( StatusCode::SUCCESS);
601}
602
603// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
604StatusCode AbsCor::finalize() {
605
606 MsgStream log(msgSvc(), name());
607 log << MSG::INFO << "in finalize()" << endmsg;
608 return StatusCode::SUCCESS;
609}
610
611
612
613double AbsCor::E25min(int n) const
614{
615 return e25min[n];
616}
617double AbsCor::E25max(int n) const
618{
619 return e25max[n];
620}
621
622
623//The following function is copied from PhotonCor/McCor
624double AbsCor::ECorrMC(double eg, double theid) const
625{
626 double Energy5x5=eg;
627 if(eg<E25min(int(theid))) eg=E25min(int(theid));
628 if(eg>E25max(int(theid))) eg=E25max(int(theid))-0.001;
629
630 if(theid<=0)theid=0.001;
631 if(theid>=27)theid=26.999;
632 Float_t einter = eg + 0.00001;
633 Float_t tinter = theid+0.0001;
634 //cout<<"inter="<< einter<<" "<<tinter<<endl;
635 double ecor=dt->Interpolate(einter,tinter);
636 // cout<<"ecor="<<ecor<<endl;
637 if(!(ecor))return Energy5x5;
638 if(ecor<0.5)return Energy5x5;
639 double EnergyCor=Energy5x5/ecor;
640 return EnergyCor;
641}
642
643// Get energy error
644double AbsCor::ErrMC(double eg, double theid) const
645{
646
647 if(eg<E25min(int(theid))) eg=E25min(int(theid));
648 if(eg>E25max(int(theid))) eg=E25max(int(theid))-0.001;
649 if(theid<=0)theid=0.001;
650 if(theid>=27)theid=26.999;
651 Float_t einter = eg + 0.00001;
652 Float_t tinter = theid+0.0001;
653 double err=dtErr->Interpolate(einter,tinter);
654 return err;
655}
656
657
658
659
660
661
662
663
664
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: AbsCor.cxx:31
ai
double ai[4]
Definition: AbsCor.cxx:58
cbeam
double cbeam[56]
Definition: AbsCor.cxx:57
runNo
int runNo
Definition: DQA_TO_DB.cxx:12
n
const Int_t n
Definition: DataBase/tau_mode.c:65
tan
double tan(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:216
energy
************Class m_ypar INTEGER m_KeyWgt INTEGER m_nphot 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
AbsCor::execute
StatusCode execute()
Definition: AbsCor.cxx:281
AbsCor::AbsCor
AbsCor(const std::string &name, ISvcLocator *pSvcLocator)
Definition: AbsCor.cxx:59
AbsCor::finalize
StatusCode finalize()
Definition: AbsCor.cxx:604
AbsCor::initialize
StatusCode initialize()
Definition: AbsCor.cxx:77
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