BOSS 7.0.2
BESIII Offline Software System
Loading...
Searching...
No Matches
TofShower.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2#include "GaudiKernel/Bootstrap.h"
3#include "GaudiKernel/ISvcLocator.h"
4#include "GaudiKernel/PropertyMgr.h"
5#include "GaudiKernel/IJobOptionsSvc.h"
6#include "GaudiKernel/INTupleSvc.h"
7#include "TofEnergyRec/TofShower.h"
8#include "Identifier/Identifier.h"
9#include "Identifier/TofID.h"
10#include "RawDataProviderSvc/TofData.h"
11#include "TofCaliSvc/ITofCaliSvc.h"
12#include "TofQCorrSvc/ITofQCorrSvc.h"
13#include "TofQElecSvc/ITofQElecSvc.h"
14#include "TofEnergyCalibSvc/ITofEnergyCalibSvc.h"
15
16
17#include "DstEvent/TofHitStatus.h"
18#include "globals.hh"
19#include <G4String.hh>
20#include <iostream>
21#include <fstream>
22#include <math.h>
23#include <string>
24#include "TofSim/BesTofDigitizer.hh"
25#include "EvTimeEvent/RecEsTime.h"
26#include "RawEvent/RawDataUtil.h"
27
28
29
30
31TofShower::TofShower():m_output(false),m_isData(true) // m_output wird auf false und m_isData auf true gesetzt
32{
33 IJobOptionsSvc* jobSvc;
34 Gaudi::svcLocator()->service("JobOptionsSvc", jobSvc);
35 jobSvc->setMyProperties("TofShower", &m_propMgr);
36
37}
38
39void TofShower::BookNtuple(NTuple::Tuple*& tuple, NTuple::Tuple*& tuple1, NTuple::Tuple*& tuple2)
40{
41 m_output = true;
42 m_tuple = tuple;
43 if(!m_tuple) {
44 std::cerr << "Invalid ntuple in TofEnergyRec!" << std::endl;
45 } else {
46 m_tuple->addItem ("part", m_part);
47 m_tuple->addItem ("layer", m_layer);
48 m_tuple->addItem ("im", m_im);
49 m_tuple->addItem ("end", m_end);
50 m_tuple->addItem ("zpos", m_zpos);
51 m_tuple->addItem ("adc1", m_adc1);
52 m_tuple->addItem ("adc2", m_adc2);
53 m_tuple->addItem ("tdc1", m_tdc1);
54 m_tuple->addItem ("tdc2", m_tdc2);
55 m_tuple->addItem ("energy", m_energy);
56 }
57
58 m_tuple1 = tuple1;
59 if(!m_tuple1) {
60 std::cerr << "Invalid ntuple1 in TofEnergyRec!" << std::endl;
61 } else {
62 m_tuple1->addItem ("part", m_shower_part);
63 m_tuple1->addItem ("layer", m_shower_layer);
64 m_tuple1->addItem ("im", m_shower_im);
65 m_tuple1->addItem ("zpos", m_shower_zpos);
66 m_tuple1->addItem ("energy", m_shower_energy);
67 }
68
69 m_tuple2 = tuple2;
70 if(!m_tuple2) {
71 std::cerr << "Invalid ntuple2 in TofEnergyRec!" << std::endl;
72 } else {
73 m_tuple2->addItem ("dist", m_seed_dist);
74 }
75}
76
77void TofShower::energyCalib(vector<TofData*>& tofDataVec, RecTofTrackCol* recTofTrackCol)
78{
79 //Get TOF Calibtration Service
80 ISvcLocator* svcLocator = Gaudi::svcLocator();
81 ITofCaliSvc* tofCaliSvc;
82 StatusCode sc = svcLocator->service("TofCaliSvc", tofCaliSvc);
83 if (sc != StatusCode::SUCCESS) {
84 cout << "TofEnergyRec Get Calibration Service Failed !! " << endl;
85 }
86
87 ITofQCorrSvc* tofQCorrSvc;
88 sc = svcLocator->service("TofQCorrSvc", tofQCorrSvc);
89 if (sc != StatusCode::SUCCESS) {
90 cout << "TofEnergyRec Get QCorr Service Failed !! " << endl;
91 }
92
93 ITofQElecSvc* tofQElecSvc;
94 sc = svcLocator->service("TofQElecSvc", tofQElecSvc);
95 if (sc != StatusCode::SUCCESS) {
96 cout << "TofEnergyRec Get QElec Service Failed !! " << endl;
97 }
98
99 ITofEnergyCalibSvc* m_TofEnergyCalibSvc=0;
100 sc = svcLocator->service("TofEnergyCalibSvc", m_TofEnergyCalibSvc);
101 if( sc != StatusCode::SUCCESS){
102 cout << "can not use TofEnergyCalibSvc" << endreq;
103 }
104 else {
105 /*
106 std::cout << "Test TofEnergyCalibSvc getCalibConst= "
107 << m_TofEnergyCalibSvc -> getCalibConst()
108 <<"para1="<<m_TofEnergyCalibSvc -> getPara1()
109 <<"para2="<<m_TofEnergyCalibSvc -> getPara2()
110 <<"para3="<<m_TofEnergyCalibSvc -> getPara3()
111 <<"para4="<<m_TofEnergyCalibSvc -> getPara4()
112 <<"para5="<<m_TofEnergyCalibSvc -> getPara5()
113 << std::endl;
114 */
115 }
116
117
118
119 vector<TofData*>::iterator it;
120 for(it=tofDataVec.begin();
121 it!=tofDataVec.end();
122 it++) {
123
124 Identifier id((*it)->identify());
125 if( TofID::is_scin(id) ) {
126 int barrel_ec = TofID::barrel_ec(id);
127 int layer = TofID::layer(id);
128 int im = TofID::phi_module(id);
129 int end = TofID::end(id);
130
131 if(m_output) {
132 m_part = barrel_ec;
133 m_layer = layer;
134 m_im = im;
135 m_end = end;
136 }
137
138 if((*it)->barrel()) {
139 TofData* bTofData = (*it);
140 bTofData->setZpos(99.);
141 bTofData->setEnergy(0.);
142 if(bTofData->tdc1()<=0||bTofData->tdc1()>8000||bTofData->tdc2()<=0||bTofData->tdc2()>8000) continue;
143
144 double adc1,adc2,tdc1,tdc2;
145 tdc1 = bTofData->tdc1();
146 tdc2 = bTofData->tdc2();
147 adc1 = bTofData->adc1();
148 adc2 = bTofData->adc2();
149
150 //from data CalibSvc
151 double zpos = tofCaliSvc->ZTDC( tdc1, tdc2, bTofData->tofId() );
152 if(fabs(zpos)>115) continue;
153 double tofq = tofCaliSvc->BPh( adc1, adc2, zpos, bTofData->tofId());
154 if(tofq<100||tofq>10000) continue;
155 //double energy = q*0.0036;
156 //double energy = tofq*m_calibConst; //new calibration result in 2009.9.27
157 //cout<< "TofShower Barrel energy " << energy << endl;
158
159 zpos /= 100.; //cm->m
160
161 ///////////////////
162
163
164 double CalibPar[5];
165
166 CalibPar[0]=m_TofEnergyCalibSvc -> getPara1();
167 CalibPar[1]=m_TofEnergyCalibSvc -> getPara2();
168 CalibPar[2]=m_TofEnergyCalibSvc -> getPara3();
169 CalibPar[3]=m_TofEnergyCalibSvc -> getPara4();
170 CalibPar[4]=m_TofEnergyCalibSvc -> getPara5();
171
172 double TofEnCalibConst=CalibPar[0]+CalibPar[1]*zpos+CalibPar[2]*pow(zpos,2)+CalibPar[3]*pow(zpos,3)+CalibPar[4]*pow(zpos,4);
173 double energy = tofq*TofEnCalibConst;
174
175 //cout<<"tofenergy="<<energy<<endl;
176 //cout<<"tofenergy-old="<<tofq*m_calibConst<<endl;
177 ///////////////////
178 bTofData->setZpos(zpos);
179 bTofData->setEnergy(energy);
180
181 if(m_output) {
182 m_part = barrel_ec;
183 m_layer = layer;
184 m_im = im;
185 m_end = end;
186 m_adc1 = bTofData->adc1();
187 m_adc2 = bTofData->adc2();
188 m_tdc1 = bTofData->tdc1();
189 m_tdc2 = bTofData->tdc2();
190 m_zpos = zpos;
191 m_energy = energy;
192 m_tuple->write();
193 }
194
195 }
196 else {
197 //cout<<"endcap"<<endl;
198 //ETofData* eTofData = dynamic_cast<ETofData*>(*it);
199 //TofData* bTofData = (*it);
200 //double energy = 2*eTofData->adcChannel()/140;
201 //eTofData->setEnergy(energy);
202 }
203 }
204 else {
205 int barrel_ec = TofID::barrel_ec(id);
206 int endcap = TofID::endcap(id);
207 int module = TofID::module(id);
208 int strip = TofID::strip(id);
209 int end = TofID::end(id);
210
211 if( barrel_ec!=3 || ( endcap!=0 && endcap!=1 ) ) {
212 cout << "TofEnergyRec Get ETF(MRPC) data ERROR !! barrel_ec:" << barrel_ec << " endcap:" << endcap << endl;
213 }
214
215 if(m_output) {
216 m_part = barrel_ec;
217 m_layer = endcap;
218 m_im = module;
219 m_end = strip;
220 }
221
222 TofData* etfData = (*it);
223 etfData->setZpos(99.);
224 etfData->setEnergy(0.);
225 if(etfData->tdc1()<=0||etfData->tdc1()>8000||etfData->tdc2()<=0||etfData->tdc2()>8000) continue;
226
227 double adc1,adc2,tdc1,tdc2;
228 tdc1 = etfData->tdc1();
229 tdc2 = etfData->tdc2();
230 adc1 = etfData->adc1();
231 adc2 = etfData->adc2();
232
233 //from data CalibSvc
234 double zpos = tofCaliSvc->ZTDC( tdc1, tdc2, etfData->tofId() );
235 if(fabs(zpos)>115) continue;
236 double tofq = tofCaliSvc->BPh( adc1, adc2, zpos, etfData->tofId());
237 if(tofq<100||tofq>10000) continue;
238 //double energy = q*0.0036;
239 double energy = tofq*m_calibConst; //new calibration result in 2009.9.27
240 //cout<< "TofShower Barrel energy " << energy << endl;
241 zpos /= 100.; //cm->m
242
243 etfData->setZpos(zpos);
244 etfData->setEnergy(energy);
245
246 if(m_output) {
247 m_part = barrel_ec;
248 m_layer = endcap;
249 m_im = module;
250 m_end = strip;
251 m_adc1 = etfData->adc1();
252 m_adc2 = etfData->adc2();
253 m_tdc1 = etfData->tdc1();
254 m_tdc2 = etfData->tdc2();
255 m_zpos = zpos;
256 m_energy = energy;
257 m_tuple->write();
258 }
259 }
260 }
261 return;
262}
263
264vector<Identifier> TofShower::getNeighbors(const Identifier& id)
265{
266 vector<int> NeighborVec;
267 vector<Identifier> NeighborIdVec;
268 NeighborVec.clear();
269 NeighborIdVec.clear();
270
271 if( TofID::is_scin(id) ) {
272 int barrel_ec = TofID::barrel_ec(id);
273 int layer = TofID::layer(id);
274 int im = TofID::phi_module(id);
275 int end = TofID::end(id);
276
277 if(barrel_ec==1) { //barrel
278 int num = im+layer*88;
279 if(num<88) { //layer1
280 if(num==0) {
281 NeighborVec.push_back(1);
282 NeighborVec.push_back(87);
283 NeighborVec.push_back(88);
284 NeighborVec.push_back(89);
285 } else if(num==87) {
286 NeighborVec.push_back(0);
287 NeighborVec.push_back(86);
288 NeighborVec.push_back(88);
289 NeighborVec.push_back(175);
290 } else {
291 NeighborVec.push_back(num+1);
292 NeighborVec.push_back(num-1);
293 NeighborVec.push_back(num+88);
294 NeighborVec.push_back(num+88+1);
295 }
296 } else {
297 if(num==88) {
298 NeighborVec.push_back(89);
299 NeighborVec.push_back(175);
300 NeighborVec.push_back(0);
301 NeighborVec.push_back(87);
302 } else if(num==175) {
303 NeighborVec.push_back(88);
304 NeighborVec.push_back(174);
305 NeighborVec.push_back(86);
306 NeighborVec.push_back(87);
307 } else {
308 NeighborVec.push_back(num+1);
309 NeighborVec.push_back(num-1);
310 NeighborVec.push_back(num-88);
311 NeighborVec.push_back(num-88-1);
312 }
313 }
314
315 int size=NeighborVec.size();
316 for(int i=0;i<size;i++) {
317 layer = NeighborVec[i]/88;
318 im = NeighborVec[i]%88;
319 NeighborIdVec.push_back(TofID::cell_id(barrel_ec,layer,im,end));
320 }
321 }
322
323 else { //endcap
324 if(im==0) {
325 NeighborVec.push_back(1);
326 NeighborVec.push_back(47);
327 } else if(im==47) {
328 NeighborVec.push_back(0);
329 NeighborVec.push_back(46);
330 } else {
331 NeighborVec.push_back(im-1);
332 NeighborVec.push_back(im+1);
333 }
334
335 int size=NeighborVec.size();
336 for(int i=0;i<size;i++) {
337 im = NeighborVec[i];
338 NeighborIdVec.push_back(TofID::cell_id(barrel_ec,layer,im,end));
339 }
340 }
341 }
342 else { // mrpc
343 int barrel_ec = TofID::barrel_ec(id);
344 int endcap = TofID::endcap(id);
345 int module = TofID::module(id);
346 int strip = TofID::strip(id);
347 int end = TofID::end(id);
348
349 if( module==0 ) {
350 NeighborVec.push_back(1);
351 NeighborVec.push_back(35);
352 } else if( module==35 ) {
353 NeighborVec.push_back(0);
354 NeighborVec.push_back(34);
355 } {
356 NeighborVec.push_back(module-1);
357 NeighborVec.push_back(module+1);
358 }
359
360 int size=NeighborVec.size();
361 for(int i=0;i<size;i++) {
362 int im = NeighborVec[i];
363 for(int j=-2; j<3; j++ ) {
364 int ist = strip + j;
365 if( ist<0 || ist>11 ) continue;
366 NeighborIdVec.push_back(TofID::cell_id(barrel_ec,endcap,module,ist,end));
367 }
368 }
369
370 }
371
372 return NeighborIdVec;
373}
374
375vector<Identifier> TofShower::getNextNeighbors(const Identifier& id)
376{
377 vector<Identifier> NeighborVec,tmpNeighborVec,tmpNextNeighborVec;
378 vector<Identifier>::iterator ci_NV,ci_tmpNV,ci_tmpNNV;
379 NeighborVec=getNeighbors(id);
380 tmpNeighborVec=getNeighbors(id);
381 bool flag=false; //whether NeighborVec already includes this crystal
382 bool flagNeighbor=false; //whether this crystal belongs to NeighborVec
383
384 //------------------------------------------------------------------
385 for(ci_tmpNV=tmpNeighborVec.begin();
386 ci_tmpNV!=tmpNeighborVec.end();
387 ci_tmpNV++){
388 tmpNextNeighborVec=getNeighbors(*ci_tmpNV);
389 //================================================================
390 for(ci_tmpNNV=tmpNextNeighborVec.begin();
391 ci_tmpNNV!=tmpNextNeighborVec.end();
392 ci_tmpNNV++){
393
394 for(ci_NV=NeighborVec.begin();
395 ci_NV!=NeighborVec.end();
396 ci_NV++){
397 if(*ci_tmpNNV==*ci_NV){ //this crystal is already included
398 flag=true;
399 break;
400 }
401 }
402
403 if(!flag){ //find a new crystal
404 //for(ci_tmpNV1=tmpNeighborVec.begin();
405 // ci_tmpNV1!=tmpNeighborVec.end();
406 // ci_tmpNV1++){
407 // if(*ci_tmpNNV==*ci_tmpNV1){ //this crystal belongs to NeighborVec
408 // flagNeighbor=true;
409 // break;
410 // }
411 //}
412
413 if(*ci_tmpNNV==id)
414 flagNeighbor=true;
415
416 if(!flagNeighbor)
417 NeighborVec.push_back(*ci_tmpNNV);
418 else
419 flagNeighbor=false;
420 }
421 else
422 flag=false;
423 }
424 //================================================================
425 }
426 //------------------------------------------------------------------
427
428 return NeighborVec;
429}
430
431
432
433
434void TofShower::findSeed(vector<TofData*>& tofDataVec)
435{
436 bool max=false;
437 m_seedVec.clear();
438
439 vector<TofData*>::iterator it;
440 for(it=tofDataVec.begin();
441 it!=tofDataVec.end();
442 it++) {
443
444 if( !((*it)->is_mrpc()) ) {
445 if((*it)->barrel()) { //barrel
446 TofData* bTofData = (*it);
447
448 //cout << "Identifier(bTofData->identify()) " << Identifier(bTofData->identify()) << endl;
449 //std::cout << "TofShower bTofData->energy() = " << bTofData->energy() << std::endl;
450 if(bTofData->energy()<5.) continue; //seed energy cut = 6MeV
451
452 max=true;
453 vector<Identifier> NeighborVec=getNextNeighbors(Identifier(bTofData->identify()));
454
455 //const Identifier & help = Identifier(bTofData->identify());
456 //cout << "TofShower::findSeed Barrel Track base "<<TofID::layer(help) << " " << TofID::phi_module(help) << endl;
457
458 vector<Identifier>::iterator iNeigh;
459 for(iNeigh=NeighborVec.begin();
460 iNeigh!=NeighborVec.end();
461 iNeigh++) {
462
463 //const Identifier & help2 = Identifier(*iNeigh);
464 //cout << "TofShower::findSeed Barrel Track neigh "<<TofID::layer(help2) << " " << TofID::phi_module(help2) << endl;
465
466 vector<TofData*>::iterator it2;
467 for(it2=tofDataVec.begin();
468 it2!=tofDataVec.end();
469 it2++) {
470 if((*it2)->identify()==*iNeigh) {
471 TofData* bTofData2 = (*it2);
472 if(bTofData2->energy()>bTofData->energy()) {
473 max=false;
474 }
475 break;
476 }
477 }
478
479 }
480 }
481 else { //endcap
482
483 //cout << "TofShower::findSeed Endcap Track" << endl;
484 TofData* eTofData = (*it);
485 if(eTofData->energy()<5.) continue; //seed energy cut = 5MeV
486 max=true;
487 vector<Identifier> NeighborVec=getNextNeighbors(Identifier(eTofData->identify()));
488 vector<Identifier>::iterator iNeigh;
489 for(iNeigh=NeighborVec.begin(); iNeigh!=NeighborVec.end(); iNeigh++)
490 {
491 vector<TofData*>::iterator it2;
492 for(it2=tofDataVec.begin(); it2!=tofDataVec.end(); it2++)
493 {
494 if((*it2)->identify()==*iNeigh) {
495 TofData* eTofData2 = (*it2);
496 if(eTofData2->energy()>eTofData->energy()) {
497 max=false;
498 }
499 break;
500 }
501 }//close for
502 }//Close for
503 }//close if(!is_mrpc)
504
505 }
506 else
507 {
508
509 TofData* etfData = (*it);
510 if(etfData->energy()<2.) continue; //Seed energy cut = 2 MeV
511 max=true;
512 vector<Identifier> NeighborVec=getNextNeighbors(Identifier(etfData->identify()));
513 vector<Identifier>::iterator iNeigh;
514 for(iNeigh=NeighborVec.begin(); iNeigh!=NeighborVec.end(); iNeigh++)
515 {
516 //cout << "TofShower::findSeed Phimodule MRPC " << TofID::phi_module(*iNeigh) <<endl;
517 vector<TofData*>::iterator it2;
518 for(it2=tofDataVec.begin(); it2!=tofDataVec.end(); it2++)
519 {
520 if((*it2)->identify()==*iNeigh) {
521 TofData* etfData2 = (*it2);
522 if(etfData2->energy()>etfData->energy()) {
523 max=false;
524 }
525 break;
526 }
527 }//close for
528 }//Close for
529 //cout << "TofShower::findSeed Both forloops done" << endl;
530 }//Close else if(is_mrpc)
531
532 if(max) {
533 m_seedVec.push_back(Identifier((*it)->identify()));
534 }
535
536 }//close loop over tofdata
537}//close find seed
538
539
540void TofShower::findShower(vector<TofData*>& tofDataVec, RecTofTrackCol* recTofTrackCol, double t0)
541{
542
543 energyCalib(tofDataVec, recTofTrackCol);
544 //cout << "TofShower::findShower energycalib OK" << endl;
545 findSeed(tofDataVec);
546 //cout << "TofShower::findShower findseed OK" << endl;
547
548 vector<Identifier>::iterator iSeed;
549
550 for(iSeed=m_seedVec.begin(); iSeed!=m_seedVec.end(); iSeed++)
551 {
552
553 bool is_mrpc = TofID::is_mrpc(*iSeed);
554
555 if(!is_mrpc) //Old TofDetector + barrel
556 {
557 int barrel_ec = TofID::barrel_ec(*iSeed);
558 int layer = TofID::layer(*iSeed);
559 int im = TofID::phi_module(*iSeed);
560 im += layer * 88;
561
562 bool neutral=true;
563 //match with Tof charged track
564 int dphi=999;
565 RecTofTrackCol::iterator iTrack, iMatch;
566 for(iTrack=recTofTrackCol->begin(); iTrack!=recTofTrackCol->end(); iTrack++)
567 {
568 TofHitStatus *status = new TofHitStatus;
569 status->setStatus((*iTrack)->status());
570 if( barrel_ec==1 && status->is_barrel() ) {
571 dphi=abs(im-(*iTrack)->tofID());
572 dphi = dphi>=44 ? 88-dphi : dphi;
573 } else if( barrel_ec==2 && !(status->is_barrel()) && ((*iTrack)->tofID()>47) ) {
574 dphi=abs(im-(*iTrack)->tofID()+48);
575 dphi = dphi>=24 ? 48-dphi : dphi;
576 } else if( barrel_ec==0 && !(status->is_barrel()) && ((*iTrack)->tofID()<48) ) {
577 dphi=abs(im-(*iTrack)->tofID());
578 dphi = dphi>=24 ? 48-dphi : dphi;
579 }
580 if(abs(dphi)<=2) {
581 iMatch = iTrack;
582 neutral=false;
583 break;
584 }
585 }
586
587 //energy sum of seed and its neighbors
588 //use avarage mean to calculation position
589 double zpos=0;
590 double energy=0;
591 double seedPos=0;
592
593 //cout << "Energhy =0 " << endl;
594 vector<TofData*>::iterator it;
595 for(it=tofDataVec.begin(); it!=tofDataVec.end(); it++)
596 {
597 if((*it)->identify()==*iSeed) {
598 //cout<<"iSeed="<<*iSeed<<endl;
599 if((*it)->barrel()) {
600 TofData* bTofData = (*it);
601 zpos+=bTofData->zpos()*bTofData->energy();
602 energy+=bTofData->energy();
603 seedPos=bTofData->zpos();
604
605 //cout << "Add energy barrel seed "<< TofID::layer(*iSeed) <<" " << TofID::phi_module(*iSeed) << " " << bTofData->energy() << " ---> "<< energy <<endl;
606 } else {
607 TofData* eTofData = (*it);
608 energy+=eTofData->energy();
609
610 //cout << "Add energy" << endl;
611 }
612 break;
613 }
614 }
615
616 vector<Identifier> NeighborVec=getNextNeighbors(*iSeed);
617 vector<Identifier>::iterator iNeigh;
618 for(iNeigh=NeighborVec.begin(); iNeigh!=NeighborVec.end(); iNeigh++)
619 {
620
621 vector<TofData*>::iterator it2;
622 for(it2=tofDataVec.begin(); it2!=tofDataVec.end(); it2++)
623 {
624 if((*it2)->identify()==*iNeigh) {
625 //cout<<"iNeigh="<<*iNeigh<<endl;
626 if((*it)->barrel()) {
627 TofData* bTofData2 = (*it2);
628
629 if(fabs(bTofData2->zpos())>2) continue;
630 if(m_output) {
631 m_seed_dist = seedPos-bTofData2->zpos();
632 m_tuple2->write();
633 }
634 if(fabs(seedPos-bTofData2->zpos())>0.3) continue;
635 zpos+=bTofData2->zpos()*bTofData2->energy();
636 energy+=bTofData2->energy();
637 //cout << "Add energy barrel neig " << TofID::layer(*iNeigh) <<" " << TofID::phi_module(*iNeigh) << " " << bTofData2->energy() << " ---> "<< energy <<endl;
638 } else {
639 TofData* eTofData2 = (*it2);
640 energy+=eTofData2->energy();
641 }
642 break;
643 }
644 }
645 }
646 if(energy>0) zpos/=energy;
647
648 //for charged track, set energy
649 if(neutral==false) {
650 if(fabs(zpos)<1.15&&energy>5.&&energy<1000) {
651 (*iMatch)->setEnergy(energy/1000);
652 }
653 continue;
654 }
655
656 //for neutral track
657 if(fabs(zpos)<1.15&&energy>5.&&energy<1000) { //shower energy cut = 10MeV
658 RecTofTrack* tof = new RecTofTrack;
659 tof->setTofID(*iSeed);
660 TofHitStatus* hitStatus = new TofHitStatus;
661 hitStatus->init();
662 tof->setStatus( hitStatus->value());
663 tof->setZrHit(zpos);
664 tof->setEnergy(energy/1000); //MeV-->GeV
665 recTofTrackCol->push_back(tof);
666
667 if(m_output) {
668 m_shower_part = barrel_ec;
669 m_shower_layer = layer;
670 m_shower_im = im;
671 m_shower_zpos = zpos;
672 m_shower_energy = energy;
673 m_tuple1->write();
674 }
675 }
676 }//close if(!is_mrpc)
677 else // mrpc_case
678 {
679 //Determine wether the track is a neutral one or not:(Compare with charged tracks)
680 //cout << "TofShower::findShower Seed is mrpc" << endl;
681
682 int barrel_ec = TofID::barrel_ec(*iSeed);
683 int endcap = TofID::endcap(*iSeed);
684 int im = TofID::module(*iSeed);
685 int strip = TofID::strip(*iSeed);
686 int end = TofID::end(*iSeed);
687 im += endcap * 36;
688
689 bool neutral=true;
690 //match with Tof charged track
691 int dphi=999, dtheta=999;
692 RecTofTrackCol::iterator iTrack, iMatch;
693 for(iTrack=recTofTrackCol->begin(); iTrack!=recTofTrackCol->end(); iTrack++)
694 {
695 TofHitStatus *status = new TofHitStatus;
696 status->setStatus((*iTrack)->status());
697 if( barrel_ec==3 && endcap==0 && ((*iTrack)->tofID()<36) ) {
698 dphi = abs(im-(*iTrack)->tofID());
699 dphi = dphi>=18 ? 36-dphi : dphi;
700 dtheta = abs( strip - status->layer() );
701 } else if( barrel_ec==3 && endcap==1 && ((*iTrack)->tofID()>35) ) {
702 dphi=abs(im-(*iTrack)->tofID()+36);
703 dphi = dphi>=18 ? 36-dphi : dphi;
704 dtheta = abs( strip - status->layer() );
705 }
706 if( ( abs(dphi)==0 && abs(dtheta)<=2 ) || ( abs(dphi)==1 && abs(dtheta)<=1 ) ) {
707 iMatch = iTrack;
708 neutral = false;
709 break;
710 }
711 }
712
713 //energy sum of seed and its neighbors
714 //use avarage mean to calculation position
715 double zpos=0;
716 double energy=0;
717 double seedPos=0;
718
719 //cout << "Energhy =0 " << endl;
720 vector<TofData*>::iterator it;
721 for(it=tofDataVec.begin(); it!=tofDataVec.end(); it++)
722 {
723 if((*it)->identify()==*iSeed) {
724 TofData* etfData = (*it);
725 zpos+=etfData->zpos()*etfData->energy();
726 energy+=etfData->energy();
727 seedPos=etfData->zpos();
728 break;
729 }
730 }
731
732 vector<Identifier> NeighborVec=getNextNeighbors(*iSeed);
733 vector<Identifier>::iterator iNeigh;
734 for(iNeigh=NeighborVec.begin(); iNeigh!=NeighborVec.end(); iNeigh++)
735 {
736 vector<TofData*>::iterator it2;
737 for(it2=tofDataVec.begin(); it2!=tofDataVec.end(); it2++)
738 {
739 if((*it)->barrel()) continue;
740 if((*it2)->identify()==*iNeigh) {
741 TofData* etfData2 = (*it2);
742 if(fabs(etfData2->zpos())>2) continue;
743 if(m_output) {
744 m_seed_dist = seedPos-etfData2->zpos();
745 m_tuple2->write();
746 }
747 if(fabs(seedPos-etfData2->zpos())>0.3) continue;
748 zpos+=etfData2->zpos()*etfData2->energy();
749 energy+=etfData2->energy();
750 //cout << "Add energy barrel neig " << TofID::layer(*iNeigh) <<" " << TofID::phi_module(*iNeigh) << " " << bTofData2->energy() << " ---> "<< energy <<endl;
751 break;
752 }
753 }
754 }
755 if(energy>0) zpos/=energy;
756
757 //for charged track, set energy
758 if(neutral==false) {
759 if(fabs(zpos)<1.15&&energy>5.&&energy<1000) {
760 (*iMatch)->setEnergy(energy/1000);
761 }
762 continue;
763 }
764
765 //for neutral track
766 if(fabs(zpos)<1.15&&energy>5.&&energy<1000) { //shower energy cut = 10MeV
767 RecTofTrack* tof = new RecTofTrack;
768 tof->setTofID(*iSeed);
769 TofHitStatus* hitStatus = new TofHitStatus;
770 hitStatus->init();
771 tof->setStatus( hitStatus->value());
772 tof->setZrHit(zpos);
773 tof->setEnergy(energy/1000); //MeV-->GeV
774 recTofTrackCol->push_back(tof);
775
776 if(m_output) {
777 m_shower_part = barrel_ec;
778 m_shower_layer = strip;
779 m_shower_im = im;
780 m_shower_zpos = zpos;
781 m_shower_energy = energy;
782 m_tuple1->write();
783 }
784 }
785 }
786
787 }
788}
789
790void TofShower::readCalibPar()
791{
792 string paraPath = getenv("TOFENERGYRECROOT");
793 paraPath += "/share/peak.dat";
794 ifstream in;
795 in.open(paraPath.c_str());
796 assert(in);
797 for(int i=0;i<176;i++) {
798 in>>m_ecalib[i];
799 }
800 in.close();
801
802 paraPath = getenv("TOFENERGYRECROOT");
803 paraPath += "/share/calib.dat";
804 ifstream in1;
805 in1.open(paraPath.c_str());
806 assert(in1);
807 for(int i=0;i<176;i++) {
808 in1>>m_calib[i][0]>>m_calib[i][1]>>m_calib[i][2]>>m_calib[i][3];
809 }
810 in1.close();
811}
812
813double TofShower::ecalib(const int nsci) const
814{
815 if(nsci<176) {
816 return m_ecalib[nsci];
817 } else {
818 return 0;
819 }
820}
821
822void TofShower::setEcalib(const int nsci, const double ecalib)
823{
824 if(nsci<176) {
825 m_ecalib[nsci]=ecalib;
826 }
827}
828
829double TofShower::calib(const int n, const int m) const
830{
831 if(n<176&&m<4) {
832 return m_calib[n][m];
833 } else {
834 return 0;
835 }
836}
837
838void TofShower::setCalib(const int n, const int m, const double ecalib)
839{
840 if(n<176&&m<4) {
841 m_calib[n][m]=ecalib;
842 }
843}
844
n
const Int_t n
Definition: DataBase/tau_mode.c:65
RecTofTrackCol
ObjectVector< RecTofTrack > RecTofTrackCol
Definition: InstallArea/include/TofRecEvent/TofRecEvent/RecTofTrack.h:33
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
tofQElecSvc
ITofQElecSvc * tofQElecSvc
Definition: TofRawDataProvider.cxx:34
tofQCorrSvc
ITofQCorrSvc * tofQCorrSvc
Definition: TofRawDataProvider.cxx:33
tofCaliSvc
ITofCaliSvc * tofCaliSvc
Definition: TofRawDataProvider.cxx:32
ITofCaliSvc::BPh
virtual const double BPh(double ADC1, double ADC2, double zHit, unsigned int id)=0
ITofCaliSvc::ZTDC
virtual const double ZTDC(double tleft, double tright, unsigned id)=0
TofData::adc2
double adc2()
Definition: TofData.cxx:656
TofData::tdc2
double tdc2()
Definition: TofData.cxx:665
TofData::tdc1
double tdc1()
Definition: TofData.cxx:647
TofData::adc1
double adc1()
Definition: TofData.cxx:638
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:134
TofID::endcap
static int endcap(const Identifier &id)
Definition: TofID.cxx:124
TofID::strip
static int strip(const Identifier &id)
Definition: TofID.cxx:136
TofID::cell_id
static Identifier cell_id(int barrel_ec, int layer, int phi_module, int end)
For a single crystal.
Definition: TofID.cxx:143
TofID::is_scin
static bool is_scin(const Identifier &id)
Definition: TofID.cxx:102
TofID::end
static int end(const Identifier &id)
Definition: TofID.cxx:79
TofID::is_mrpc
static bool is_mrpc(const Identifier &id)
Definition: TofID.cxx:113
TofID::phi_module
static int phi_module(const Identifier &id)
Definition: TofID.cxx:73
TofID::barrel_ec
static int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: TofID.cxx:61
TofID::layer
static int layer(const Identifier &id)
Definition: TofID.cxx:66
TofID::module
static int module(const Identifier &id)
Definition: TofID.cxx:130
TofShower::calib
double calib(const int n, const int m) const
Definition: TofShower.cxx:829
TofShower::ecalib
double ecalib(const int nsci) const
Definition: TofShower.cxx:813
TofShower::readCalibPar
void readCalibPar()
Definition: TofShower.cxx:790
TofShower::findSeed
void findSeed(vector< TofData * > &tofDataVec)
TofShower::BookNtuple
void BookNtuple(NTuple::Tuple *&tuple, NTuple::Tuple *&tuple1, NTuple::Tuple *&tuple2)
Definition: TofShower.cxx:39
TofShower::setEcalib
void setEcalib(const int nsci, const double ecalib)
Definition: TofShower.cxx:822
TofShower::energyCalib
void energyCalib(vector< TofData * > &tofDataVec, RecTofTrackCol *recTofTrackCol)
TofShower::getNextNeighbors
vector< Identifier > getNextNeighbors(const Identifier &id)
Definition: TofShower.cxx:375
TofShower::setCalib
void setCalib(const int n, const int m, const double ecalib)
Definition: TofShower.cxx:838
TofShower::getNeighbors
vector< Identifier > getNeighbors(const Identifier &id)
Definition: TofShower.cxx:264
TofShower::findShower
void findShower(vector< TofData * > &tofDataVec, RecTofTrackCol *recTofTrackCol, double)