CGEM BOSS 6.6.5.f
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("UseTof", usetof=1);
68 declareProperty("DoDataCor", dodatacor = 1);
69}
70
71
72TGraphErrors *dt = new TGraphErrors();
73// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
74StatusCode AbsCor::initialize(){
75 MsgStream log(msgSvc(), name());
76 log << MSG::INFO << "in initialize()" << endmsg;
77 StatusCode status;
78 if(ntOut == true){
79 NTuplePtr nt1(ntupleSvc(), "FILE1/ec");
80 if ( nt1 ) m_tuple1 = nt1;
81 else {
82 m_tuple1 = ntupleSvc()->book ("FILE1/ec", CLID_ColumnWiseTuple, "ks N-Tuple example");
83 if ( m_tuple1 ) {
84 status = m_tuple1->addItem ("ef", m_ef);
85 status = m_tuple1->addItem ("e5", m_e5);
86 status = m_tuple1->addItem ("ec", m_ec);
87 status = m_tuple1->addItem ("ct", m_ct);
88 status = m_tuple1->addItem ("cedge", m_cedge);
89 }
90 else {
91 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
92 return StatusCode::FAILURE;
93 }
94 }
95 }
96
97 m_index = new int*[56];
98 for (int j=0;j<56;j++ ) {
99 m_index[j] = new int[120];
100 for ( int k=0; k<120; k++) {
101 m_index[j][k]=-1;
102 }
103 }
104
105 m_par = new double*[22];
106 for (int j=0;j<22;j++ ) {
107 m_par[j] = new double[11];
108 for ( int k=0; k<11; k++) {
109 m_par[j][k]=-99;
110 }
111 }
112
113 m_parphi = new double*[22];
114 for (int j=0;j<22;j++ ) {
115 m_parphi[j] = new double[5];
116 for ( int k=0; k<5; k++) {
117 m_parphi[j][k]=-99;
118 }
119 }
120
121 ifstream EnParIn;
122 //EnParIn.exceptions( ifstream::failbit | ifstream::badbit );
123 string EnParInFile;
124 EnParInFile = getenv("ABSCORROOT");
125 EnParInFile += "/share/para.dat";
126 EnParIn.open(EnParInFile.c_str());
127 for(int i=0;i<22;i++) {
128
129 EnParIn>>m_par[i][0]>>m_par[i][1]>>m_par[i][2]>>m_par[i][3]>>m_par[i][4]
130 >>m_par[i][5]>>m_par[i][6]>>m_par[i][7]>>m_par[i][8]>>m_par[i][9]
131 >>m_par[i][10];
132 }
133 EnParIn.close();
134
135 ifstream EnParInphi;
136 //EnParInphi.exceptions( ifstream::failbit | ifstream::badbit );
137 string EnParInFilephi = getenv("ABSCORROOT");
138 EnParInFilephi += "/share/paraphi.dat";
139 EnParInphi.open(EnParInFilephi.c_str());
140 for(int i=0;i<22;i++) {
141 EnParInphi>>m_parphi[i][0]>>m_parphi[i][1]>>m_parphi[i][2]>>m_parphi[i][3]>>m_parphi[i][4];
142
143 }
144 EnParInphi.close();
145
146
147
148 string DataPathevse;
149 DataPathevse = getenv("ABSCORROOT");
150 DataPathevse += "/share/barmccorevse10bin.txt";
151 ifstream inpi0;
152 inpi0.exceptions( ifstream::failbit | ifstream::badbit );
153 inpi0.open(DataPathevse.c_str(),ios::in);
154
155 double epoint[11],peak[11],peakerr[11];
156 for(Int_t i=0;i<11;i++){
157 inpi0>>epoint[i];
158 inpi0>>peak[i];
159 inpi0>>peakerr[i];
160 }
161 for(Int_t i=0;i<11;i++){
162 dt->SetPoint(i, epoint[i],peak[i]);
163 dt->SetPointError(i,0,peakerr[i]);
164 }
165
166 string DataPathcbeam;
167 DataPathcbeam = getenv("ABSCORROOT");
168 DataPathcbeam += "/share/cbeam.txt";
169 ifstream incbeam;
170 incbeam.exceptions( ifstream::failbit | ifstream::badbit );
171 incbeam.open(DataPathcbeam.c_str(),ios::in);
172 for(int i=0; i<56; i++){
173 double tid,peak,peake,sig,sige;
174 incbeam>>tid;
175 incbeam>>peak;
176 incbeam>>peake;
177 incbeam>>sig;
178 incbeam>>sige;
179 cbeam[i]=1.0/peak;
180 }
181
182 string DataPathc3p;
183 DataPathc3p = getenv("ABSCORROOT");
184
185 string DataPathc3ptof;
186 DataPathc3ptof = getenv("ABSCORROOT");
187
188 cout<<"mccor = "<<mccor<<endl;
189
190 DataPathc3p += "/share/c3p.txt";
191 DataPathc3ptof += "/share/c3ptof.txt";
192
193 ifstream inc3p;
194 inc3p.exceptions( ifstream::failbit | ifstream::badbit );
195 if(usetof==0)inc3p.open(DataPathc3p.c_str(),ios::in);
196 if(usetof==1)inc3p.open(DataPathc3ptof.c_str(),ios::in);
197 for(int i=0; i<4; i++){
198 double am,ame;
199 inc3p>>am;
200 inc3p>>ame;
201 ai[i]=am;
202 }
203
204
205
206
207 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
208 return StatusCode::SUCCESS;
209
210}
211
212// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
213StatusCode AbsCor::execute() {
214 MsgStream log(msgSvc(), name());
215 log << MSG::INFO << "in execute()" << endreq;
216// SmartDataPtr<Event::EventH> eventHeader(eventSvc(),"/Event");
217 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
218 int runNo=eventHeader->runNumber();
219
220 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
221 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
222 if(evtRecEvent->totalTracks()>evtRecTrkCol->size())return SUCCESS;
223 if(evtRecEvent->totalTracks()>50)return SUCCESS;
224
225 IEmcRecGeoSvc* iGeoSvc;
226 ISvcLocator* svcLocator = Gaudi::svcLocator();
227 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
228 if(sc!=StatusCode::SUCCESS) {
229 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
230 }
231
232 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
233 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
234 if(!(*itTrk)->isEmcShowerValid())continue;
235 RecEmcShower *emcTrk = (*itTrk)->emcShower();
236// if(emcTrk->energy()<0.0005)continue;
237 int st = emcTrk->status();
238 if(st>10000)continue;
239 emcTrk->setStatus(st+20000);
240// cout<<"REC after calib EMC status = "<<emcTrk->status()<<endl;
241 if(emcTrk->e5x5()<0.015)continue;
242// if(emcTrk->getTime()<10||emcTrk->getTime()>30)continue;
243/*
244 if(emcTrk->e5x5()<0.02){
245 emcTrk->setEnergy(emcTrk->e5x5()*1.1);
246 continue;
247 }
248*/
249/*
250 int intid = emcTrk->cellId();
251 Identifier id=Identifier(intid);
252 int getthetaid = EmcID::theta_module(id);
253 int getmodule = EmcID::barrel_ec(id);
254 if(getmodule==1)getthetaid=getthetaid+6;
255 if(getmodule==2)getthetaid=55-getthetaid;
256*/
257 double etof=0;
258
259 if(usetof==1 && (*itTrk)->isTofTrackValid()){
260 SmartRefVector<RecTofTrack> recTofTrackVec=(*itTrk)->tofTrack();
261 if(!recTofTrackVec.empty())etof=recTofTrackVec[0]->energy();
262 if(etof>100.)etof=0;
263 }
264 double energyC=emcTrk->energy()+etof;
265/*
266 if(energyC<=0||energyC>4.99)continue;
267 double cor = 1.0;
268 if(runNo>0)cor = dt->Eval(energyC);
269 if(cor<0.001)continue;
270// cout<<cor<<endl;
271 double energyCC= energyC/cor;
272 emcTrk->setEnergy(energyCC);
273*/
274
275 double lnE = std::log(energyC);
276 double lnEcor=1.0;
277 if(runNo>0 && dodatacor==1) lnEcor = ai[0]+ai[1]*lnE+ai[2]*lnE*lnE+ai[3]*lnE*lnE*lnE;
278 if(lnEcor<0.5)continue;
279// cout<<lnEcor<<" "<<etof<<endl;
280
281 RecEmcID id = Identifier(emcTrk->cellId());
282
283 unsigned int module, ntheta, nphi;
284 module = EmcID::barrel_ec(id);
285 ntheta = EmcID::theta_module(id);
286 nphi = EmcID::phi_module(id);
287
288
289
290 id=EmcID::crystal_id(module,ntheta,nphi);
291
292 unsigned int thetaModule = EmcID::theta_module(id);
293 unsigned int phiModule = EmcID::phi_module(id);
294
295 HepPoint3D pos(emcTrk->position());
296
297 double enecor=1.;
298 if(edgecor==1){
299
300 if(module==1) { //barrel
301 unsigned int theModule;
302 if(thetaModule>21) {
303 theModule = 43 - thetaModule;
304 id = EmcID::crystal_id(module,theModule,phiModule);
305 pos.setZ(-pos.z());
306 } else {
307 theModule = thetaModule;
308 }
309
310 double IRShift;
311 if(theModule==21) {
312 IRShift = 0;
313
314 } else if(theModule==20) {
315 IRShift = 2.5;
316
317 } else {
318 IRShift = 5.;
319
320 }
321
322
323 HepPoint3D IR(0,0,IRShift);
324 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
325 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
326 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
327
328 double theta01,theta23,length,d;
329 theta01 = (center01-IR).theta();
330 theta23 = (center23-IR).theta();
331 length = (center01-IR).mag();
332 d = length*tan(theta01-theta23); //length in x direction
333
334 double xIn;
335 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
336 if (xIn < (-d/2) && theModule!=21){
337
338 theModule=theModule+1 ;
339
340 id = EmcID::crystal_id(module,theModule,phiModule);
341 if(theModule==21) {
342 IRShift = 0;
343
344 } else if(theModule==20) {
345 IRShift = 2.5;
346
347 } else {
348 IRShift = 5.;
349
350 }
351 HepPoint3D IR1(0,0,IRShift);
352 IR=IR1;
353 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
354 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
355 theta01 = (center01-IR).theta();
356 theta23 = (center23-IR).theta();
357 length = (center01-IR).mag();
358 d = length*tan(theta01-theta23); //length in x direction
359
360 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
361
362 } else if (xIn > (d/2) && theModule!=0 ){
363
364 theModule=theModule-1 ;
365 id = EmcID::crystal_id(module,theModule,phiModule);
366 if(theModule==21) {
367 IRShift = 0;
368
369 } else if(theModule==20) {
370 IRShift = 2.5;
371
372 } else {
373 IRShift = 5.;
374
375 }
376
377 HepPoint3D IR1(0,0,IRShift);
378 IR=IR1;
379 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
380 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
381 theta01 = (center01-IR).theta();
382 theta23 = (center23-IR).theta();
383 length = (center01-IR).mag();
384 d = length*tan(theta01-theta23); //length in x direction
385
386 xIn = length*tan(theta01-(pos-IR).theta())-d/2;
387 }
388
389 double yIn;
390 // yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
391 // : pos.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
392
393 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
394 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
395
396 if(nphi==0&&yIn>100) yIn=yIn-360;
397 if(nphi==119&&yIn<-100) yIn=yIn+360;
398
399 if(yIn<-1.5-0.15){
400
401 if(nphi!=0) phiModule=phiModule-1;
402 if(nphi==0) phiModule=119;
403 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
404 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
405 if(phiModule==0&&yIn>100) yIn=yIn-360;
406 if(phiModule==119&&yIn<-100) yIn=yIn+360;
407
408 }
409
410 if(yIn>1.5-0.15){
411
412
413 if(nphi!=119) phiModule=phiModule+1;
414 if(nphi==119) phiModule=0;
415
416 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
417 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
418 if(phiModule==0&&yIn>100) yIn=yIn-360;
419 if(phiModule==119&&yIn<-100) yIn=yIn+360;
420 }
421
422 enecor=m_par[theModule][0]
423 +m_par[theModule][1]*xIn
424 +m_par[theModule][2]*xIn*xIn
425 +m_par[theModule][3]*xIn*xIn*xIn
426 +m_par[theModule][4]*xIn*xIn*xIn*xIn
427 +m_par[theModule][5]*xIn*xIn*xIn*xIn*xIn
428 +m_par[theModule][6]*xIn*xIn*xIn*xIn*xIn*xIn
429 // +m_par[theModule][7]*xIn*xIn*xIn*xIn*xIn*xIn*xIn
430 // +m_par[theModule][8]*xIn*xIn*xIn*xIn*xIn*xIn*xIn*xIn
431 +m_par[theModule][7]*yIn
432 +m_par[theModule][8]*yIn*yIn
433 +m_par[theModule][9]*yIn*yIn*yIn
434 +m_par[theModule][10]*yIn*yIn*yIn*yIn;
435 //////////////////////
436 } else{
437 enecor=1;
438 }
439 }
440 double energyCC= energyC/(lnEcor*enecor);
441// cout<<"Trk Level Corr. and Orig. "<<energyCC<<" "<<emcTrk->energy()<<endl;
442 emcTrk->setEnergy(energyCC);
443 if(ntOut==true){
444 m_ef=energyCC;
445 m_e5=emcTrk->e5x5();
446 m_ec=energyC;
447 m_ct=lnEcor ;
448 m_cedge=enecor;
449 m_tuple1->write();
450 }
451
452
453 }
454//==============Modify Dst===============================================================
455 SmartDataPtr<RecEmcShowerCol> recEmcShowerCol(eventSvc(), EventModel::Recon::RecEmcShowerCol);
456 if(!recEmcShowerCol){
457 return( StatusCode::SUCCESS);
458 }
459 SmartDataPtr<DstEmcShowerCol> dstEmcShowerCol(eventSvc(), EventModel::Dst::DstEmcShowerCol);
460 if(!dstEmcShowerCol){
461 return( StatusCode::SUCCESS);
462 }
463
464// cout<<"Rec and Dst Size "<<recEmcShowerCol->size()<<" "<<dstEmcShowerCol->size()<<endl;
465 if(recEmcShowerCol->size()!=dstEmcShowerCol->size())return SUCCESS;
466 for(int i=0;i<recEmcShowerCol->size();i++){
467 RecEmcShowerCol::iterator iter_rec = recEmcShowerCol->begin()+i;
468 DstEmcShowerCol::iterator iter_dst = dstEmcShowerCol->begin()+i;
469// cout<<"Rec and Dst energy "<<(*iter_rec)->energy()<<" "<<(*iter_dst)->energy()<<endl;
470 (*iter_dst)->setEnergy((*iter_rec)->energy());
471 (*iter_dst)->setStatus((*iter_rec)->status());
472// cout<<"DST after calib EMC status = "<<(*iter_dst)->status()<<endl;
473// cout<<"Rec == Dst? "<<(*iter_rec)->energy()<<" "<<(*iter_dst)->energy()<<endl;
474 }
475 return( StatusCode::SUCCESS);
476}
477
478// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
479StatusCode AbsCor::finalize() {
480
481 MsgStream log(msgSvc(), name());
482 log << MSG::INFO << "in finalize()" << endmsg;
483 return StatusCode::SUCCESS;
484}
485
486
487
488
489
490
491
492
493
494
495
496
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: AbsCor.cxx:31
dt
TGraphErrors * dt
Definition: AbsCor.cxx:72
ai
double ai[4]
Definition: AbsCor.cxx:58
cbeam
double cbeam[56]
Definition: AbsCor.cxx:57
runNo
int runNo
Definition: DedxCalibRunByRun.cxx:14
tan
double tan(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:216
AbsCor::execute
StatusCode execute()
Definition: AbsCor.cxx:213
AbsCor::AbsCor
AbsCor(const std::string &name, ISvcLocator *pSvcLocator)
Definition: AbsCor.cxx:59
AbsCor::finalize
StatusCode finalize()
Definition: AbsCor.cxx:479
AbsCor::initialize
StatusCode initialize()
Definition: AbsCor.cxx:74
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: bak_Identifier-01-02-14/src/EmcID.cxx:71
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: bak_Identifier-01-02-14/src/EmcID.cxx:38
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition: bak_Identifier-01-02-14/src/EmcID.cxx:43
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition: bak_Identifier-01-02-14/src/EmcID.cxx:48
IEmcRecGeoSvc::GetCrystalPoint
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) const =0