BOSS 7.0.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
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 in2>>energy;
228 in2>>thetaid;
229 in2>>peak1;
230 in2>>peakerr1;
231 in2>>res;
232 in2>>reserr;
233 dt->SetPoint(i,energy,thetaid,peak1);
234 dt->SetPointError(i,0,0,peakerr1);
235 dtErr->SetPoint(i,energy,thetaid,res);
236 dtErr->SetPointError(i,0,0,reserr);
237 if(i<28) e25min[int(thetaid)]=energy;
238 if(i>=560-28) e25max[int(thetaid)]=energy;
239 }
240 in2.close();
241 //-------------------------
242 // Suppression of hot crystals
243 // Reading the map from hotcry.txt (Hajime, Jul 2013)
244 for(int ih=0;ih<10;ih++){
245 hrunstart[ih]=-1;
246 hrunend[ih]=-1;
247 hcell[ih]=-1;
248 }
249 int numhots=4; // numbers of hot crystals
250 int dumring,dumphi,dummod,dumid;
251 string HotList;
252 HotList = getenv("ABSCORROOT");
253 HotList += "/share/hotcry.txt";
254 ifstream hotcrys;
255 hotcrys.exceptions( ifstream::failbit | ifstream::badbit );
256 hotcrys.open(HotList.c_str(),ios::in);
257 for(int il=0; il<numhots; il++){
258 hotcrys>>hrunstart[il];
259 hotcrys>>hrunend[il];
260 hotcrys>>hcell[il];
261 hotcrys>>dumring;
262 hotcrys>>dumphi;
263 hotcrys>>dummod;
264 hotcrys>>dumid;
265 }
266 hotcrys.close();
267 //-------------------------
268
269
270 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
271 return StatusCode::SUCCESS;
272
273
274
275
276}
277
278// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
279StatusCode AbsCor::execute() {
280 MsgStream log(msgSvc(), name());
281 log << MSG::INFO << "in execute()" << endreq;
282// SmartDataPtr<Event::EventH> eventHeader(eventSvc(),"/Event");
283 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
284 int runNo=eventHeader->runNumber();
285
286 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
287 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
288 if(evtRecEvent->totalTracks()>evtRecTrkCol->size())return SUCCESS;
289 if(evtRecEvent->totalTracks()>50)return SUCCESS;
290
291 IEmcRecGeoSvc* iGeoSvc;
292 ISvcLocator* svcLocator = Gaudi::svcLocator();
293 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
294 if(sc!=StatusCode::SUCCESS) {
295 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
296 }
297
298 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
299 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
300 if(!(*itTrk)->isEmcShowerValid())continue;
301 RecEmcShower *emcTrk = (*itTrk)->emcShower();
302// if(emcTrk->energy()<0.0005)continue;
303
304 int st = emcTrk->status();
305 if(st>10000)continue;
306 emcTrk->setStatus(st+20000);
307
308// cout<<"REC after calib EMC status = "<<emcTrk->status()<<endl;
309
310 //if(emcTrk->e5x5()<0.015)continue;
311
312// if(emcTrk->getTime()<10||emcTrk->getTime()>30)continue;
313/*
314 if(emcTrk->e5x5()<0.02){
315 emcTrk->setEnergy(emcTrk->e5x5()*1.1);
316 continue;
317 }
318*/
319/*
320 int intid = emcTrk->cellId();
321 Identifier id=Identifier(intid);
322 int getthetaid = EmcID::theta_module(id);
323 int getmodule = EmcID::barrel_ec(id);
324 if(getmodule==1)getthetaid=getthetaid+6;
325 if(getmodule==2)getthetaid=55-getthetaid;
326*/
327
328 RecEmcID id = Identifier(emcTrk->cellId());
329
330 unsigned int module, ntheta, nphi;
331 module = EmcID::barrel_ec(id);
332 ntheta = EmcID::theta_module(id);
333 nphi = EmcID::phi_module(id);
334
335 // id=EmcID::crystal_id(module,ntheta,nphi);
336
337 unsigned int thetaModule = EmcID::theta_module(id);
338 unsigned int phiModule = EmcID::phi_module(id);
339
340
341 double e5x5=emcTrk->e5x5();
342 double etof=0;
343
344 if(usetof==1 && (*itTrk)->isTofTrackValid()){
345 SmartRefVector<RecTofTrack> recTofTrackVec=(*itTrk)->tofTrack();
346 if(!recTofTrackVec.empty())etof=recTofTrackVec[0]->energy();
347 if(etof>100.)etof=0;
348 }
349
350 double energyC;
351
352 int thetaId;
353 if (module==0||module==2) thetaId = thetaModule;
354 if (module==1 && thetaModule<=21) thetaId = thetaModule + 6;
355 if (module==1 && thetaModule>21) thetaId = 43 - thetaModule + 6;
356
357 if(MCuseTof==1){
358 energyC=ECorrMC(e5x5+etof,thetaId);
359
360 }
361 if (MCuseTof==0) {
362 energyC=ECorrMC(e5x5,thetaId);
363
364 }
365
366
367 /*
368 if(usetof==1){
369 energyC=ECorrMC(e5x5+etof,thetaId);
370 }
371 if (usetof==0) {
372 energyC=emcTrk->energy();
373 }
374 */
375 // double energyC=emcTrk->energy()+etof;
376/*
377 if(energyC<=0||energyC>4.99)continue;
378 double cor = 1.0;
379 if(runNo>0)cor = dt->Eval(energyC);
380 if(cor<0.001)continue;
381// cout<<cor<<endl;
382 double energyCC= energyC/cor;
383 emcTrk->setEnergy(energyCC);
384*/
385
386 double lnE = std::log(energyC);
387
388 if (energyC>1.0) lnE=std::log(1.0);
389 if (energyC<0.07) lnE=std::log(0.07);
390
391 double lnEcor=1.0;
392 if(dopi0Cor==1){
393 if(runNo>0 && dodatacor==1) {
394 lnEcor = ai[0]+ai[1]*lnE+ai[2]*lnE*lnE+ai[3]*lnE*lnE*lnE;
395 }
396 }
397 if(lnEcor<0.5)continue;
398// cout<<lnEcor<<" "<<etof<<endl;
399
400
401//////////////////////////////////////now no using in the following
402 HepPoint3D pos(emcTrk->position());
403
404 // If it is "hot", return "9999" (Hajime, Jul 2013)
405 double enecor=1.;
406 if(hotcellmask==1){
407 for(int ih=0;ih<10;ih++){
408 if(hrunstart[ih]==-1||hrunend[ih]==-1||hcell[ih]==-1)continue;
409 if(abs(runNo)<hrunstart[ih]||abs(runNo)>hrunend[ih])continue;
410 if ( emcTrk->cellId()==hcell[ih] ) {emcTrk->setStatus(9999);}
411 }
412 }
413
414 if(edgecor==1){
415
416 if(module==1) { //barrel
417 unsigned int theModule;
418 if(thetaModule>21) {
419 theModule = 43 - thetaModule;
420 id = EmcID::crystal_id(module,theModule,phiModule);
421 pos.setZ(-pos.z());
422 } else {
423 theModule = thetaModule;
424 }
425
426 double IRShift;
427 if(theModule==21) {
428 IRShift = 0;
429
430 } else if(theModule==20) {
431 IRShift = 2.5;
432
433 } else {
434 IRShift = 5.;
435
436 }
437
438
439 HepPoint3D IR(0,0,IRShift);
440 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
441 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
442 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
443
444 double theta01,theta23,length,d;
445 theta01 = (center01-IR).theta();
446 theta23 = (center23-IR).theta();
447 length = (center01-IR).mag();
448 d = length*tan(theta01-theta23); //length in x direction
449
450 double xIn;
451 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
452 if (xIn < (-d/2) && theModule!=21){
453
454 theModule=theModule+1 ;
455
456 id = EmcID::crystal_id(module,theModule,phiModule);
457 if(theModule==21) {
458 IRShift = 0;
459
460 } else if(theModule==20) {
461 IRShift = 2.5;
462
463 } else {
464 IRShift = 5.;
465
466 }
467 HepPoint3D IR1(0,0,IRShift);
468 IR=IR1;
469 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
470 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
471 theta01 = (center01-IR).theta();
472 theta23 = (center23-IR).theta();
473 length = (center01-IR).mag();
474 d = length*tan(theta01-theta23); //length in x direction
475
476 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
477
478 } else if (xIn > (d/2) && theModule!=0 ){
479
480 theModule=theModule-1 ;
481 id = EmcID::crystal_id(module,theModule,phiModule);
482 if(theModule==21) {
483 IRShift = 0;
484
485 } else if(theModule==20) {
486 IRShift = 2.5;
487
488 } else {
489 IRShift = 5.;
490
491 }
492
493 HepPoint3D IR1(0,0,IRShift);
494 IR=IR1;
495 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
496 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
497 theta01 = (center01-IR).theta();
498 theta23 = (center23-IR).theta();
499 length = (center01-IR).mag();
500 d = length*tan(theta01-theta23); //length in x direction
501
502 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
503 }
504
505 double yIn;
506 // yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
507 // : pos.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
508
509 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
510 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
511
512 if(nphi==0&&yIn>100) yIn=yIn-360;
513 if(nphi==119&&yIn<-100) yIn=yIn+360;
514
515 if(yIn<-1.5-0.15){
516
517 if(nphi!=0) phiModule=phiModule-1;
518 if(nphi==0) phiModule=119;
519 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
520 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
521 if(phiModule==0&&yIn>100) yIn=yIn-360;
522 if(phiModule==119&&yIn<-100) yIn=yIn+360;
523
524 }
525
526 if(yIn>1.5-0.15){
527
528
529 if(nphi!=119) phiModule=phiModule+1;
530 if(nphi==119) phiModule=0;
531
532 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
533 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
534 if(phiModule==0&&yIn>100) yIn=yIn-360;
535 if(phiModule==119&&yIn<-100) yIn=yIn+360;
536 }
537
538 enecor=m_par[theModule][0]
539 +m_par[theModule][1]*xIn
540 +m_par[theModule][2]*xIn*xIn
541 +m_par[theModule][3]*xIn*xIn*xIn
542 +m_par[theModule][4]*xIn*xIn*xIn*xIn
543 +m_par[theModule][5]*xIn*xIn*xIn*xIn*xIn
544 +m_par[theModule][6]*xIn*xIn*xIn*xIn*xIn*xIn
545 // +m_par[theModule][7]*xIn*xIn*xIn*xIn*xIn*xIn*xIn
546 // +m_par[theModule][8]*xIn*xIn*xIn*xIn*xIn*xIn*xIn*xIn
547 +m_par[theModule][7]*yIn
548 +m_par[theModule][8]*yIn*yIn
549 +m_par[theModule][9]*yIn*yIn*yIn
550 +m_par[theModule][10]*yIn*yIn*yIn*yIn;
551 //////////////////////
552 } else{
553 enecor=1;
554 }
555 }
556//////////////////////////////////////now no using in the above
557
558 double energyCC= energyC/(lnEcor*enecor);
559 // cout<<"Trk Level Corr. and Orig. "<<energyCC<<" "<<emcTrk->energy()<<endl;
560 emcTrk->setEnergy(energyCC);
561
562 if(ntOut==true){
563 m_ef=energyCC;
564 m_e5=emcTrk->e5x5();
565 m_ec=energyC;
566 m_ct=lnEcor ;
567 m_cedge=enecor;
568 m_tuple1->write();
569 }
570
571
572 }
573//==============Modify Dst===============================================================
574 SmartDataPtr<RecEmcShowerCol> recEmcShowerCol(eventSvc(), EventModel::Recon::RecEmcShowerCol);
575 if(!recEmcShowerCol){
576 return( StatusCode::SUCCESS);
577 }
578 SmartDataPtr<DstEmcShowerCol> dstEmcShowerCol(eventSvc(), EventModel::Dst::DstEmcShowerCol);
579 if(!dstEmcShowerCol){
580 return( StatusCode::SUCCESS);
581 }
582
583// cout<<"Rec and Dst Size "<<recEmcShowerCol->size()<<" "<<dstEmcShowerCol->size()<<endl;
584 if(recEmcShowerCol->size()!=dstEmcShowerCol->size())return SUCCESS;
585 for(int i=0;i<recEmcShowerCol->size();i++){
586 RecEmcShowerCol::iterator iter_rec = recEmcShowerCol->begin()+i;
587 DstEmcShowerCol::iterator iter_dst = dstEmcShowerCol->begin()+i;
588// cout<<"Rec and Dst energy "<<(*iter_rec)->energy()<<" "<<(*iter_dst)->energy()<<endl;
589 (*iter_dst)->setEnergy((*iter_rec)->energy());
590 (*iter_dst)->setStatus((*iter_rec)->status());
591// cout<<"DST after calib EMC status = "<<(*iter_dst)->status()<<endl;
592// cout<<"Rec == Dst? "<<(*iter_rec)->energy()<<" "<<(*iter_dst)->energy()<<endl;
593 }
594 return( StatusCode::SUCCESS);
595}
596
597// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
598StatusCode AbsCor::finalize() {
599
600 MsgStream log(msgSvc(), name());
601 log << MSG::INFO << "in finalize()" << endmsg;
602 return StatusCode::SUCCESS;
603}
604
605
606
607double AbsCor::E25min(int n) const
608{
609 return e25min[n];
610}
611double AbsCor::E25max(int n) const
612{
613 return e25max[n];
614}
615
616
617//The following function is copied from PhotonCor/McCor
618double AbsCor::ECorrMC(double eg, double theid) const
619{
620 double Energy5x5=eg;
621 if(eg<E25min(int(theid))) eg=E25min(int(theid));
622 if(eg>E25max(int(theid))) eg=E25max(int(theid))-0.001;
623
624 if(theid<=0)theid=0.001;
625 if(theid>=27)theid=26.999;
626 Float_t einter = eg + 0.00001;
627 Float_t tinter = theid+0.0001;
628 //cout<<"inter="<< einter<<" "<<tinter<<endl;
629 double ecor=dt->Interpolate(einter,tinter);
630 // cout<<"ecor="<<ecor<<endl;
631 if(!(ecor))return Energy5x5;
632 if(ecor<0.5)return Energy5x5;
633 double EnergyCor=Energy5x5/ecor;
634 return EnergyCor;
635}
636
637// Get energy error
638double AbsCor::ErrMC(double eg, double theid) const
639{
640
641 if(eg<E25min(int(theid))) eg=E25min(int(theid));
642 if(eg>E25max(int(theid))) eg=E25max(int(theid))-0.001;
643 if(theid<=0)theid=0.001;
644 if(theid>=27)theid=26.999;
645 Float_t einter = eg + 0.00001;
646 Float_t tinter = theid+0.0001;
647 double err=dtErr->Interpolate(einter,tinter);
648 return err;
649}
650
651
652
653
654
655
656
657
658
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:279
AbsCor::AbsCor
AbsCor(const std::string &name, ISvcLocator *pSvcLocator)
Definition: AbsCor.cxx:59
AbsCor::finalize
StatusCode finalize()
Definition: AbsCor.cxx:598
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