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