BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
BesTrigL1.cxx
Go to the documentation of this file.
1#include "GaudiKernel/AlgFactory.h"
2#include "GaudiKernel/Bootstrap.h"
3#include "GaudiKernel/IDataProviderSvc.h"
4#include "GaudiKernel/IHistogramSvc.h"
5#include "GaudiKernel/ISvcLocator.h"
6#include "GaudiKernel/MsgStream.h"
7#include "GaudiKernel/PropertyMgr.h"
8#include "GaudiKernel/SmartDataPtr.h"
9
10#include "CLHEP/Units/PhysicalConstants.h"
11#include <CLHEP/Geometry/Point3D.h>
12
13#include "EmcRawEvent/EmcDigi.h"
14#include "EmcWaveform/EmcWaveform.h"
15#include "EvTimeEvent/RecEsTime.h"
16#include "EventModel/Event.h"
17#include "EventModel/EventHeader.h"
18#include "EventModel/EventModel.h"
19#include "Identifier/Identifier.h"
20#include "Identifier/MdcID.h"
21#include "Identifier/MucID.h"
22#include "Identifier/TofID.h"
23#include "McTruth/McParticle.h"
24#include "McTruth/MucMcHit.h"
25#include "MdcRawEvent/MdcDigi.h"
26#include "MucRawEvent/MucDigi.h"
27#include "RawDataProviderSvc/TofData.h"
28#include "TofRawEvent/TofDigi.h"
29#include "TrigEvent/TrigData.h"
30#include "TrigEvent/TrigEACC.h"
31#include "TrigEvent/TrigEvent.h"
32#include "TrigEvent/TrigGTD.h"
33#include "TrigEvent/TrigTOFT.h"
34
35#include "EmcCalibConstSvc/EmcCalibConstSvc.h"
36
37#include "CLHEP/Random/RandFlat.h"
38#include "CLHEP/Random/RandGauss.h"
39#include "CLHEP/Random/Random.h"
40#include "Trigger/AsciiData.h"
41#include "Trigger/BesGlobalTrigSvc.h"
42#include "Trigger/BesTrigL1.h"
43#include "Trigger/IBesGlobalTrigSvc.h"
44#include "Trigger/TrigPara.h"
45#include <TCanvas.h>
46#include <TGraph.h>
47#include <map>
48#include <math.h>
49
50using namespace CLHEP;
51using namespace std;
52using namespace Event;
53
54// Declaration of the Algorithm Factory
55// static const AlgFactory<BesTrigL1> s_factory ;
56// const IAlgFactory& BesTrigL1Factory = s_factory ;
57BesTrigL1::BesTrigL1( const std::string& name, ISvcLocator* pSvcLocator )
58 : Algorithm( name, pSvcLocator ), m_pIBGT( 0 ), passNo( 0 ) {
59 declareProperty( "WriteFile", writeFile = 0 );
60 declareProperty( "IfOutEvtId", ifoutEvtId = 0 );
61 declareProperty( "Input", input = "boost.dat" );
62 declareProperty( "Output", output = "boostOut.dat" );
63 declareProperty( "OutEvtIdFile", outEvtId = "evtId.dat" );
64 declareProperty( "TrigRootFlag", mTrigRootFlag = false );
65 declareProperty( "RunMode", m_runMode = 1 );
66 declareProperty( "ClockShift", clock_shift = 0 );
67 nEvent = 0;
68}
69
70// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
71StatusCode BesTrigL1::initialize() {
72 MsgStream log( msgSvc(), name() );
73 log << MSG::INFO << "in initialize()" << endmsg;
74
75 //--------------------------------------
76 // define a pointer of trigger service
77 //--------------------------------------
78 ISvcLocator* svcLocator = Gaudi::svcLocator();
79 StatusCode sc = svcLocator->service( "BesGlobalTrigSvc", m_tmpSvc );
80 m_pIBGT = dynamic_cast<BesGlobalTrigSvc*>( m_tmpSvc );
81 if ( sc != StatusCode::SUCCESS )
82 { log << MSG::DEBUG << "Unable to open trigger service" << endmsg; }
83
84 // set run mode, 0: online, 1: offline
85 m_pIBGT->setRunMode( m_runMode );
86
87 //--------------------------------------------------------------
88 // define a pointer of RawDataProviderSvc, used in tof trigger
89 //--------------------------------------------------------------
90 static const bool CREATEIFNOTTHERE( true );
91 sc = service( "RawDataProviderSvc", m_rawDataProviderSvc, CREATEIFNOTTHERE );
92 if ( !sc.isSuccess() )
93 {
94 log << MSG::ERROR << "Could not load RawDataProviderSvc!" << m_rawDataProviderSvc
95 << endmsg;
96 return sc;
97 }
98
99 //--------------------------------------------------------------
100 // use realization service to get trigger configure parameters
101 //--------------------------------------------------------------
102 IRealizationSvc* tmpReal;
103 sc = svcLocator->service( "RealizationSvc", tmpReal );
104 if ( !sc.isSuccess() )
105 { cout << "FATAL: Could not initialize Realization Service" << endl; }
106 else { m_RealizationSvc = dynamic_cast<RealizationSvc*>( tmpReal ); }
107
108 //-----------------------------------------------------------------------
109 // use EmcCalibConstSvc to convert crystal id(theta, phi) to global id.
110 //-----------------------------------------------------------------------
111 sc = svcLocator->service( "EmcCalibConstSvc", emcCalibConstSvc );
112 if ( sc != StatusCode::SUCCESS )
113 { cout << "EmcRecDigit2Hit Error: Can't get EmcCalibConstSvc." << endl; }
114
115 if ( mTrigRootFlag )
116 {
117 //-----------------------------------------
118 // define ntuples for performance check
119 //-----------------------------------------
120 NTuplePtr nt( ntupleSvc(), "FILE302/trig1" );
121 if ( nt ) m_tuple = nt;
122 else
123 {
124 m_tuple = ntupleSvc()->book( "FILE302/trig1", CLID_ColumnWiseTuple, "TrigL1" );
125 if ( m_tuple )
126 {
127 sc = m_tuple->addItem( "x", m_wire_x );
128 sc = m_tuple->addItem( "y", m_wire_y );
129 sc = m_tuple->addItem( "evtId", m_wire_evtId );
130 sc = m_tuple->addItem( "delta_t", m_delta_tdc );
131 }
132 else
133 {
134 log << MSG::ERROR << "Cannot book N-tuple:" << long( m_tuple ) << endmsg;
135 return StatusCode::FAILURE;
136 }
137 }
138
139 NTuplePtr nt1( ntupleSvc(), "FILE302/trig2" );
140 if ( nt1 ) m_tuple1 = nt1;
141 else
142 {
143 m_tuple1 = ntupleSvc()->book( "FILE302/trig2", CLID_ColumnWiseTuple, "TrigL1" );
144 if ( m_tuple1 )
145 {
146 sc = m_tuple1->addItem( "RunId", m_RunId );
147 sc = m_tuple1->addItem( "EventId", m_EventId );
148 sc = m_tuple1->addItem( "mc_totE_all", m_mc_totE_all );
149 sc = m_tuple1->addItem( "data_totE_all", m_data_totE_all );
150 sc = m_tuple1->addItem( "mc_wetotE", m_wetotE );
151 sc = m_tuple1->addItem( "data_wetotE", m_data_wetotE );
152 sc = m_tuple1->addItem( "mc_eetotE", m_eetotE );
153 sc = m_tuple1->addItem( "data_eetotE", m_data_eetotE );
154 sc = m_tuple1->addItem( "index_btc", m_index_btc, 0, 330 );
155 sc = m_tuple1->addIndexedItem( "btc_e", m_index_btc, m_btc_e );
156 sc = m_tuple1->addIndexedItem( "data_btc", m_index_btc, m_data_btc );
157 sc = m_tuple1->addItem( "cond_id", m_cond_id, 0, 48 );
158 sc = m_tuple1->addIndexedItem( "mc_cond", m_cond_id, m_mc_cond );
159 sc = m_tuple1->addIndexedItem( "data_cond", m_cond_id, m_data_cond );
160 sc = m_tuple1->addItem( "block_id", m_block_id, 0, 12 );
161 sc = m_tuple1->addIndexedItem( "mc_blockE", m_block_id, m_mc_blockE );
162 sc = m_tuple1->addIndexedItem( "data_blockE", m_block_id, m_data_blockE );
163 sc = m_tuple1->addIndexedItem( "R_blockE", m_block_id, m_R_blockE );
164 }
165 else
166 { // did not manage to book the N tuple....
167 log << MSG::ERROR << "Cannot book N-tuple1:" << long( m_tuple1 ) << endmsg;
168 return StatusCode::FAILURE;
169 }
170 }
171
172 NTuplePtr nt2( ntupleSvc(), "FILE302/muc" );
173 if ( nt2 ) m_tuple2 = nt2;
174 else
175 {
176 m_tuple2 = ntupleSvc()->book( "FILE302/muc", CLID_ColumnWiseTuple, "TrigL1" );
177 if ( m_tuple2 )
178 {
179 sc = m_tuple2->addItem( "indexlayerSeg", m_index2, 0, 5000 );
180 sc = m_tuple2->addIndexedItem( "nlayerSeg", m_index2, m_fireLayer, 0, 5000 );
181 sc = m_tuple2->addItem( "indexhitLayer", m_index3, 0, 5000 );
182 sc = m_tuple2->addIndexedItem( "nhitLayer", m_index3, m_hitLayer, 0, 5000 );
183 sc = m_tuple2->addItem( "indexhitSeg", m_index4, 0, 5000 );
184 sc = m_tuple2->addIndexedItem( "nhitSeg", m_index4, m_hitSeg, 0, 5000 );
185 sc = m_tuple2->addItem( "indexpara", m_index5, 0, 5000 );
186 sc = m_tuple2->addIndexedItem( "costheta", m_index5, m_costheta, 0, 5000 );
187 sc = m_tuple2->addIndexedItem( "phi", m_index5, m_phi, 0, 5000 );
188 sc = m_tuple2->addIndexedItem( "p", m_index5, m_p, 0, 5000 );
189 sc = m_tuple2->addIndexedItem( "pdgcode", m_index5, m_pdgcode, 0, 5000 );
190 sc = m_tuple2->addItem( "indexhitphi", m_index6, 0, 5000 );
191 sc = m_tuple2->addIndexedItem( "hitphi", m_index6, m_hitphi, 0, 5000 );
192
193 sc = m_tuple2->addItem( "nlayerEE", m_nlayerEE );
194 sc = m_tuple2->addItem( "nlayerBR", m_nlayerBR );
195 sc = m_tuple2->addItem( "nlayerWE", m_nlayerWE );
196 sc = m_tuple2->addItem( "nlayerTotal", m_nlayerTotal );
197 sc = m_tuple2->addItem( "nhitEE", m_nhitEE );
198 sc = m_tuple2->addItem( "nhitBR", m_nhitBR );
199 sc = m_tuple2->addItem( "nhitWE", m_nhitWE );
200 sc = m_tuple2->addItem( "nhitTotal", m_nhitTotal );
201
202 sc = m_tuple2->addItem( "mumcostheta", m_mumcostheta );
203 sc = m_tuple2->addItem( "mumphi", m_mumphi );
204
205 sc = m_tuple2->addItem( "trackfindall", m_trackfindall );
206 sc = m_tuple2->addItem( "trackfind3l", m_trackfind3l );
207 sc = m_tuple2->addItem( "trackb", m_trackb );
208 sc = m_tuple2->addItem( "tracke", m_tracke );
209 sc = m_tuple2->addItem( "ntrack1", m_ntrack1 );
210 sc = m_tuple2->addItem( "ntrack2", m_ntrack2 );
211 sc = m_tuple2->addItem( "ntrack3", m_ntrack3 );
212
213 sc = m_tuple2->addItem( "ngoodevent", m_ngoodevent );
214 sc = m_tuple2->addItem( "ngoodtrack", m_ngoodtrack );
215 }
216 else
217 { // did not manage to book the N tuple....
218 log << MSG::ERROR << "Cannot book N-tuple2:" << long( m_tuple2 ) << endmsg;
219 return StatusCode::FAILURE;
220 }
221 }
222
223 NTuplePtr nt3( ntupleSvc(), "FILE302/trig3" );
224 if ( nt3 ) m_tuple3 = nt3;
225 else
226 {
227 m_tuple3 = ntupleSvc()->book( "FILE302/trig3", CLID_ColumnWiseTuple, "TrigL1" );
228 if ( m_tuple3 )
229 {
230 sc = m_tuple3->addItem( "evtId", m_evtId );
231 for ( int index = 0; index < 48; index++ )
232 {
233 std::ostringstream osname1;
234 osname1 << "cond" << index << "_1";
235 std::string name1 = osname1.str();
236
237 std::ostringstream osname2;
238 osname2 << "cond" << index << "_0";
239 std::string name2 = osname2.str();
240 m_tuple3->addItem( name1.c_str(), m_condNOne[index] );
241 m_tuple3->addItem( name2.c_str(), m_condNZero[index] );
242 }
243 }
244 else
245 { // did not manage to book the N tuple....
246 log << MSG::ERROR << "Cannot book N-tuple3:" << long( m_tuple3 ) << endmsg;
247 return StatusCode::FAILURE;
248 }
249 }
250 }
251
252 // pointer of mdc trigger
253 m_MdcTSF = MdcTSF::get_Mdc();
254
255 // pointer of tof trigger
256 m_TofHitCount = TofHitCount::get_Tof();
257
258 // pointer of emc trigger
259 m_emcDigi = EmcTCFinder::get_Emc();
260
261 // pointer of muc trigger
262 m_mucDigi = MucTrigHit::get_Muc();
263
264 //-------------------------------------
265 // reset total track and event number
266 //-------------------------------------
267 totalEvent = 0;
268 totalTracks = 0;
269
270 if ( mTrigRootFlag )
271 {
272 sc = service( "THistSvc", m_thistsvc );
273 if ( sc.isFailure() )
274 {
275 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endmsg;
276 return sc;
277 }
278 m_trigCondi_MC = new TH1F( "trgCond_MC", "trgCond_MC", 48, 0, 48 );
279 sc = m_thistsvc->regHist( "/TRG/trgCond_MC", m_trigCondi_MC );
280 m_trigCondi_Data = new TH1F( "trgCond_Data", "trgCond_Data", 48, 0, 48 );
281 sc = m_thistsvc->regHist( "/TRG/trgCond_Data", m_trigCondi_Data );
282 }
283
284 //------------------------------------------------------------------
285 // a pointer used to read emc trigger information from eacc board
286 //------------------------------------------------------------------
287 // eacctrig = new TrigEACC("eacc_trig");
288
289 return StatusCode::SUCCESS;
290}
291StatusCode BesTrigL1::execute() {
292 MsgStream log( msgSvc(), name() );
293
294 log << MSG::DEBUG << "in execute()" << endmsg;
295
296 int event, run;
297 ifpass = false;
298
299 //-------------------
300 // get event header
301 //-------------------
302 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
303 if ( !eventHeader )
304 {
305 log << MSG::FATAL << "Could not find Event Header" << endmsg;
306 return ( StatusCode::FAILURE );
307 }
308 run = eventHeader->runNumber();
309 event = eventHeader->eventNumber();
310
311 if ( mTrigRootFlag )
312 {
313 // fill run id and event id into ntuple
314 m_RunId = run;
315 m_EventId = event;
316 }
317
318 //-------------------------------------------------------------------
319 // using ascii file for output, but an ascii input file is needed.
320 //-------------------------------------------------------------------
321 if ( writeFile == 1 && event == 0 )
322 {
323 readin.open( input.c_str(), ios_base::in );
324 if ( readin ) cout << "Input File is ok " << input << endl;
325 readout.open( output.c_str(), ios_base::out );
326 if ( readout ) cout << "Output File is ok " << output << endl;
327 VERSIONNUM version;
328 readin >> version;
329 readout << version;
330 }
331
332 //---------------------------------
333 // define a map to store mdc hits
334 //---------------------------------
335 multimap<int, uint32_t, less<int>> mdc_hitmap;
336 mdc_hitmap.clear();
337
338 //---------------------------------
339 // define a map to store tof hits
340 //---------------------------------
341 multimap<int, int, less<int>> tof_hitmap;
342 tof_hitmap.clear();
343
344 //----------------------------------------------------
345 // get mdc digits from TDS and save them into a map
346 //----------------------------------------------------
347 SmartDataPtr<MdcDigiCol> mdcDigiCol( eventSvc(), "/Event/Digi/MdcDigiCol" );
348 if ( !mdcDigiCol )
349 {
350 log << MSG::FATAL << "Could not find MDC digi" << endmsg;
351 return ( StatusCode::FAILURE );
352 }
353 for ( MdcDigiCol::iterator iter3 = mdcDigiCol->begin(); iter3 != mdcDigiCol->end(); iter3++ )
354 {
355 Identifier id = ( *iter3 )->identify();
356 int layer = MdcID::layer( id );
357 int wire = MdcID::wire( id );
358 int tdc = ( *iter3 )->getTimeChannel();
359 if ( tdc < 0x7FFFFFFF && tdc > 0 )
360 {
361 if ( layer <= 19 )
362 {
363 typedef pair<int, uint32_t> vpair;
364 uint32_t mdc_Id = ( layer & 0xFFFF ) << 16;
365 mdc_Id = mdc_Id | ( wire & 0xFFFF );
366 mdc_hitmap.insert( vpair( tdc, mdc_Id ) );
367 }
368 if ( layer >= 36 && layer <= 39 )
369 {
370 typedef pair<int, uint32_t> vpair;
371 uint32_t mdc_Id = ( layer & 0xFFFF ) << 16;
372 mdc_Id = mdc_Id | ( wire & 0xFFFF );
373 mdc_hitmap.insert( vpair( tdc, mdc_Id ) );
374 }
375 }
376 }
377
378 //------------------------------------------------------------------
379 // get tof digits from rawDataProviderSvc ant save them into a map
380 //------------------------------------------------------------------
381 TofDataMap tofDigiMap = m_rawDataProviderSvc->tofDataMapEstime();
382 for ( TofDataMap::iterator iter3 = tofDigiMap.begin(); iter3 != tofDigiMap.end(); iter3++ )
383 {
384 Identifier idd = Identifier( iter3->first );
385 TofData* p_tofDigi = iter3->second;
386 int barrel_ec = TofID::barrel_ec( idd );
387 int layer = TofID::layer( idd );
388 int im = TofID::phi_module( idd );
389 if ( barrel_ec == 1 )
390 {
391 if ( ( ( p_tofDigi->quality() ) & 0x5 ) != 0x5 ) continue;
392 double tdc1 = p_tofDigi->tdc1();
393 double tdc2 = p_tofDigi->tdc2();
394 int tdc;
395 if ( tdc1 > tdc2 ) tdc = (int)tdc1;
396 else tdc = (int)tdc2;
397 typedef pair<int, int> vpair;
398 tdc = tdc;
399 tof_hitmap.insert( vpair( tdc, 10000 * barrel_ec + 1000 * layer + im * 10 ) );
400 }
401 else
402 {
403 int tdc1 = (int)p_tofDigi->tdc1();
404 typedef pair<int, int> vpair;
405 tdc1 = tdc1;
406 tof_hitmap.insert( vpair( tdc1, 10000 * barrel_ec + 1000 * layer + im * 10 ) );
407 }
408 }
409
410 //--------------------------
411 // get emc digits from TDS
412 //--------------------------
413 SmartDataPtr<EmcDigiCol> emcDigiCol( eventSvc(), "/Event/Digi/EmcDigiCol" );
414 if ( !emcDigiCol )
415 {
416 log << MSG::FATAL << "Could not find EMC digi" << endmsg;
417 return ( StatusCode::FAILURE );
418 }
419
420 //----------------------------------------------------------
421 // define initialize waveform object for each energy block
422 //----------------------------------------------------------
423 EmcWaveform blockWave[16];
424
425 //------------------------------------------------------------
426 // define a map of trigger cell, contains time and id infor.
427 //------------------------------------------------------------
428 multimap<int, uint32_t, less<int>> emc_TC;
429 emc_TC.clear();
430
431 //---------------------------------------------------------------------------------
432 // get emc analog signal, including trigger cell, energy block and cluster infor.
433 //---------------------------------------------------------------------------------
434 getEmcAnalogSig( emcDigiCol, blockWave, emc_TC );
435
436 //-----------------------------------
437 // find peak time of energy block
438 //-----------------------------------
439 double peak_time = 0;
440 findEmcPeakTime( peak_time, blockWave );
441
442 //--------------------------
443 // get muc digits from TDS
444 //--------------------------
445 SmartDataPtr<MucDigiCol> mucDigiCol( eventSvc(), "/Event/Digi/MucDigiCol" );
446 if ( !mucDigiCol )
447 {
448 log << MSG::FATAL << "Could not find MUC digi" << endmsg;
449 return ( StatusCode::FAILURE );
450 }
451 if ( m_mucDigi ) m_mucDigi->getMucDigi( mucDigiCol );
452
453 //----------------------------------------------
454 // output for debugging and count event number
455 //----------------------------------------------
456 if ( event % 10000 == 0 ) std::cout << "---> EventNo is " << event << std::endl;
457 nEvent++;
458
459 //********************************************************************
460 // start time clock
461 //********************************************************************
462
463 //--------------------------
464 // find start and end time
465 //--------------------------
466 double stime = -1, etime = -1;
467 findSETime( mdc_hitmap, tof_hitmap, emc_TC, stime, etime );
468
469 // calculate total clock number
470 int nclock = 0;
471 if ( stime >= 0 ) { nclock = int( ( etime - stime ) / 24 ) + 1; }
472 else { nclock = 0; }
473
474 //------------------------------------------------------------
475 // define an array to store trigger conditions in each clock
476 //------------------------------------------------------------
477 int** trgcond = new int*[48];
478 for ( int condId = 0; condId < 48; condId++ ) { trgcond[condId] = new int[nclock]; }
479
480 // used for tof status machine
481 int idle_status = -1;
482
483 for ( int iclock = 0; iclock < nclock; iclock++ )
484 {
485 //---------------------------
486 // start mdc trigger logic
487 //---------------------------
488 runAclock_mdc( iclock, stime, mdc_hitmap );
489
490 //---------------------------
491 // start tof trigger logic
492 //---------------------------
493 runAclock_tof( iclock, stime, idle_status, tof_hitmap );
494
495 //--------------------------
496 // start emc trigger logic
497 //--------------------------
498 runAclock_emc( iclock, stime, emc_TC, blockWave );
499
500 //----------------------------------
501 // start track match trigger logic
502 //----------------------------------
503 // m_pIBGT->startTMTrig();
504
505 //--------------------------
506 // set trigger conditions
507 //--------------------------
508 StatusCode status = m_pIBGT->setTrigCondition();
509 if ( status != StatusCode::SUCCESS )
510 {
511 log << MSG::FATAL << "Could not set trigger condition index" << endmsg;
512 return StatusCode::FAILURE;
513 }
514
515 //--------------------------------------------
516 // get each trigger condition in each clock
517 //--------------------------------------------
518 for ( int condId = 0; condId < 48; condId++ )
519 { trgcond[condId][iclock] = m_pIBGT->getTrigCond( condId ); }
520 }
521
522 //------------------------------
523 // stretch trigger conditions
524 //------------------------------
525 stretchTrgCond( nclock, trgcond );
526
527 //-------------------------
528 // SAF delay
529 //-------------------------
530 // trgSAFDelay(nclock, trgcond);
531
532 //-------------------------
533 // GTL delay
534 //-------------------------
535 // trgGTLDelay(nclock, trgcond);
536
537 //********************************************************************
538 // end time clock
539 //********************************************************************
540
541 //-------------------------------------------------------------------------------------------------------------------
542 // deal with emc trigger conditions, in principle, if NClus>=1 is true between peaktime
543 // - 1.6us and peak time, other emc conditions can be true, but not used now.
544 //-------------------------------------------------------------------------------------------------------------------
545 bool ifClus1 = false;
546 for ( int i = 0; i < nclock; i++ )
547 {
548 if ( trgcond[0][i] > 0 ) ifClus1 = true;
549 }
550
551 if ( ifClus1 == false )
552 {
553 for ( int i = 0; i < nclock; i++ )
554 {
555 for ( int j = 0; j < 16; j++ ) { trgcond[j][i] = 0; }
556 }
557 }
558
559 //-----------------------------------------------------------
560 // do logic 'or' for each trigger condition in all clocks.
561 //-----------------------------------------------------------
562 for ( int i = 0; i < nclock; i++ )
563 {
564 for ( int j = 0; j < 48; j++ )
565 {
566 if ( trgcond[j][i] ) m_pIBGT->setTrigCond( j, 1 );
567 }
568 }
569
570 //----------------------------
571 // match with trigger table
572 //----------------------------
573 m_pIBGT->GlobalTrig();
574
575 //--------------------------------------
576 // this event can pass trigger or not
577 //--------------------------------------
578 ifpass = m_pIBGT->getIfpass();
579 if ( ifpass )
580 {
581 passNo++;
582 log << MSG::INFO << "pass event number is " << passNo << endl;
583 }
584
585 //-------------------------------------------
586 // write out events which can pass trigger.
587 //-------------------------------------------
588 if ( writeFile == 2 )
589 {
590 if ( ifpass ) { setFilterPassed( true ); }
591 else { setFilterPassed( false ); }
592 }
593
594 if ( mTrigRootFlag )
595 {
596 //--------------------------------------------
597 // fill histograms, trigger conditions of MC
598 //--------------------------------------------
599 for ( int i = 0; i < 48; i++ )
600 {
601 bool edge = false;
602 int NOne = 0;
603 m_condNOne[i] = -9;
604 m_condNZero[i] = -9;
605 for ( int j = 0; j < nclock; j++ )
606 {
607 m_mc_cond[i] += trgcond[i][j];
608 if ( trgcond[i][j] != 0 )
609 {
610 if ( NOne == 0 )
611 {
612 m_condNZero[i] = j;
613 m_trigCondi_MC->Fill( i );
614 }
615 edge = true;
616 NOne++;
617 }
618 else { edge = false; }
619 if ( edge == false && NOne != 0 ) break;
620 }
621 m_condNOne[i] = NOne;
622 }
623 m_cond_id = 48;
624
625 //-----------------------------------------------
626 // fill histograms, trigger conditions of data
627 //-----------------------------------------------
628 if ( m_runMode == 0 )
629 {
630 SmartDataPtr<TrigData> trigData( eventSvc(), "/Event/Trig/TrigData" );
631 if ( !trigData )
632 {
633 log << MSG::FATAL << "Could not find Trigger Data for physics analysis" << endmsg;
634 return StatusCode::FAILURE;
635 }
636
637 for ( int id = 0; id < 48; id++ )
638 {
639 if ( trigData->getTrigCondition( id ) != 0 ) { m_trigCondi_Data->Fill( id ); }
640 m_data_cond[id] = trigData->getTrigCondition( id );
641 }
642 m_cond_id = 48;
643 }
644 }
645
646 //----------------------------
647 // release memory
648 //----------------------------
649 for ( int condId = 0; condId < 48; condId++ ) { delete trgcond[condId]; }
650 delete trgcond;
651
652 if ( mTrigRootFlag )
653 {
654 m_evtId = event;
655 m_tuple3->write();
656
657 m_mc_totE_all = m_pIBGT->getEmcTotE();
658 m_wetotE = m_pIBGT->getEmcWETotE();
659 m_eetotE = m_pIBGT->getEmcEETotE();
660
661 map<int, vector<complex<int>>, greater<int>> mymap;
662 mymap = m_pIBGT->getEmcClusId();
663 log << MSG::INFO << "EMC test: " << endmsg;
664 int emc_btc_id = 0;
665 for ( map<int, vector<complex<int>>, greater<int>>::iterator iter = mymap.begin();
666 iter != mymap.end(); iter++ )
667 {
668 if ( ( iter->first ) == 1 )
669 {
670 for ( unsigned int i = 0; i < ( iter->second ).size(); i++ )
671 {
672 log << MSG::INFO << "barrel theta is " << ( iter->second[i] ).real() << " phi is "
673 << ( iter->second[i] ).imag() << endmsg;
674 emc_btc_id++;
675 }
676 }
677 if ( ( iter->first ) == 0 )
678 {
679 for ( unsigned int i = 0; i < ( iter->second ).size(); i++ )
680 log << MSG::INFO << "east theta is " << ( iter->second[i] ).real() << " phi is "
681 << ( iter->second[i] ).imag() << endmsg;
682 }
683 if ( ( iter->first ) == 2 )
684 {
685 for ( unsigned int i = 0; i < ( iter->second ).size(); i++ )
686 log << MSG::INFO << "west theta is " << ( iter->second[i] ).real() << " phi is "
687 << ( iter->second[i] ).imag() << endmsg;
688 }
689 }
690
691 // retrieve EMC trigger information from EACC
692 /* SmartDataPtr<TrigGTDCol> trigGTDCol(eventSvc(), "/Event/Trig/TrigGTDCol");
693 if (!trigGTDCol) {
694 log << MSG::FATAL << "Could not find Global Trigger Data" << endmsg;
695 return StatusCode::FAILURE;
696 }
697 eacctrig->initialize();
698 TrigGTDCol::iterator iter5 = trigGTDCol->begin();
699 for (; iter5 != trigGTDCol->end(); iter5++) {
700 uint32_t size = (*iter5)->getDataSize();
701 const uint32_t* ptr = (*iter5)->getDataPtr();
702 //set EACC trigger data
703 if((*iter5)->getId() == 0xD7) {
704 eacctrig->setEACCTrigData((*iter5)->getId(), (*iter5)->getTimeWindow(), size, ptr);
705 }
706 }
707
708 double bmean[12] = {8.02,10.1,12.3,7.43,14.8,13.0,12.5,13.2,10.9,12.3,14.7,15.7};
709 double bsigma[12] = {0.88,0.52,0.9,0.72,0.7,0.82,0.64,0.78,0.72,0.76,0.54,0.64};
710 vector<double> vblockE = m_pIBGT->getEmcBlockE();
711 for(int blockId = 0; blockId < vblockE.size(); blockId++) {
712 //m_mc_blockE[blockId] = vblockE[blockId];
713 int block_time;
714 m_mc_blockE[blockId] = blockWave[blockId+2].max(block_time)*0.333 - 0xa +
715 RandGauss::shoot(bmean[blockId],bsigma[blockId]); m_data_blockE[blockId] =
716 eacctrig->getBBLKCharge(blockId); float r_blockE; if((eacctrig->getBBLKCharge(blockId) -
717 bmean[blockId]) == 0.) r_blockE = 0; else r_blockE =
718 vblockE[blockId]/(eacctrig->getBBLKCharge(blockId) - bmean[blockId]); if(!(r_blockE >=0.
719 || r_blockE <= 0.)) r_blockE = 0; m_R_blockE[blockId] = r_blockE;
720 }
721 m_block_id = vblockE.size();
722
723 m_data_totE_all = eacctrig->getEMCTotalCharge();
724 //endcap energy
725 int ee_endcap = 0, we_endcap = 0;
726 for(int i = 0; i < 2; i++) {
727 ee_endcap += eacctrig->getEBLKCharge(i);
728 we_endcap += eacctrig->getWBLKCharge(i);
729 }
730 m_data_wetotE = we_endcap;
731 m_data_eetotE = ee_endcap;
732
733 m_data_totE_all = eacctrig->getEMCTotalCharge();
734
735 //fill trigger cell energy
736 int window = eacctrig->getTimeWindow();
737 int index_tc = 0;
738 for(int i=0;i<TrigConf::TCTHETANO_B;i++)
739 for(int j=0;j<TrigConf::TCPHINO_B;j++)
740 {
741 m_btc_e[index_tc] = m_pIBGT->getBTCEnergy(i,j);
742 int if_clus = 0;
743 for(int k = 0; k < window; k++) {
744 if(eacctrig->getBTC(i,j,k) == 1) {
745 if_clus = 1;
746 break;
747 }
748 }
749 m_data_btc[index_tc] = if_clus;
750 index_tc++;
751 }
752 m_index_btc = index_tc;
753
754 index_tc = 0;
755 for(int i =0;i<TrigConf::TCTHETANO_E;i++)
756 for(int j =0;j<TrigConf::TCPHINO_E;j++)
757 {
758 //m_wetc_e[index_tc] = m_pIBGT->getWETCEnergy(i,j);
759 //m_eetc_e[index_tc] = m_pIBGT->getEETCEnergy(i,j);
760 index_tc++;
761 }
762 //m_index_etc = index_tc;
763 */
764 m_tuple1->write();
765
766 //----------------------------------------------
767 // check information of MDC, TOF, EMC output
768 //----------------------------------------------
769 vector<int> vstrkId;
770 vector<int> vltrkId;
771 vstrkId = m_pIBGT->getMdcStrkId();
772 vltrkId = m_pIBGT->getMdcLtrkId();
773 log << MSG::INFO << "Mdc test: " << endmsg;
774 for ( unsigned int i = 0; i < vstrkId.size(); i++ )
775 log << MSG::INFO << "short is " << vstrkId[i] << endmsg;
776 for ( unsigned int j = 0; j < vltrkId.size(); j++ )
777 { log << MSG::INFO << "long is " << vltrkId[j] << endmsg; }
778
779 map<int, vector<int>, greater<int>> tofmap;
780 tofmap = m_pIBGT->getTofHitPos();
781 log << MSG::INFO << "TOF test: " << endmsg;
782 for ( map<int, vector<int>, greater<int>>::iterator iter = tofmap.begin();
783 iter != tofmap.end(); iter++ )
784 {
785 if ( iter->first == 0 )
786 {
787 for ( unsigned int i = 0; i < iter->second.size(); i++ )
788 { log << MSG::INFO << "east tof Id is " << iter->second[i] << endmsg; }
789 }
790 if ( iter->first == 1 )
791 {
792 for ( unsigned int i = 0; i < iter->second.size(); i++ )
793 { log << MSG::INFO << " barrel tof Id is " << iter->second[i] << endmsg; }
794 }
795 if ( iter->first == 2 )
796 {
797 for ( unsigned int i = 0; i < iter->second.size(); i++ )
798 { log << MSG::INFO << "west tof Id is " << iter->second[i] << endmsg; }
799 }
800 }
801
802 // Fill ntuple for MUC
803 std::vector<int> vtmp;
804
805 vtmp = m_pIBGT->getMuclayerSeg();
806 m_index2 = 0;
807 for ( std::vector<int>::iterator iter = vtmp.begin(); iter != vtmp.end(); iter++ )
808 {
809 m_fireLayer[m_index2] = (long)*iter;
810 m_index2++;
811 if ( m_index2 > m_index2->range().distance() )
812 {
813 break;
814 cerr << "*********** too many index ************" << endl;
815 }
816 }
817 // find tracks by count the fired layer number
818 long trackb3 = 0, tracke3 = 0, trackb2 = 0, tracke2 = 0, trackb1 = 0, tracke1 = 0;
819 int trackwe = 0, trackee = 0;
820 for ( unsigned int i = 0; i < vtmp.size(); i++ )
821 {
822 if ( 0 <= vtmp[i] && vtmp[i] < 100 )
823 {
824 if ( ( vtmp[i] % 10 ) >= 3 )
825 {
826 tracke3++;
827 trackee++;
828 }
829 }
830 if ( 200 <= vtmp[i] )
831 {
832 if ( ( ( vtmp[i] - 200 ) % 10 ) >= 3 )
833 {
834 tracke3++;
835 trackwe++;
836 }
837 }
838 if ( 100 <= vtmp[i] && vtmp[i] < 200 )
839 {
840 if ( ( ( vtmp[i] - 100 ) % 10 ) >= 3 ) trackb3++;
841 }
842 }
843 m_ntrack3 = trackb3 + tracke3;
844
845 for ( unsigned int i = 0; i < vtmp.size(); i++ )
846 {
847 if ( 0 <= vtmp[i] && vtmp[i] < 100 )
848 {
849 if ( ( vtmp[i] % 10 ) >= 2 ) tracke2++;
850 }
851 if ( 200 <= vtmp[i] )
852 {
853 if ( ( ( vtmp[i] - 200 ) % 10 ) >= 2 ) tracke2++;
854 }
855 if ( 100 <= vtmp[i] && vtmp[i] < 200 )
856 {
857 if ( ( ( vtmp[i] - 100 ) % 10 ) >= 2 ) trackb2++;
858 }
859 }
860 m_ntrack2 = trackb2 + tracke2;
861
862 for ( unsigned int i = 0; i < vtmp.size(); i++ )
863 {
864 if ( 0 <= vtmp[i] && vtmp[i] < 100 )
865 {
866 if ( ( vtmp[i] % 10 ) >= 1 ) tracke1++;
867 }
868 if ( 200 <= vtmp[i] )
869 {
870 if ( ( ( vtmp[i] - 200 ) % 10 ) >= 1 ) tracke1++;
871 }
872 if ( 100 <= vtmp[i] && vtmp[i] < 200 )
873 {
874 if ( ( ( vtmp[i] - 100 ) % 10 ) >= 1 ) trackb1++;
875 }
876 }
877 m_ntrack1 = trackb1 + tracke1;
878 // end of finding tracks by count the fired layer number
879
880 vtmp = m_pIBGT->getMuchitLayer();
881 m_index3 = 0;
882 for ( std::vector<int>::iterator iter = vtmp.begin(); iter != vtmp.end(); iter++ )
883 {
884 m_hitLayer[m_index3] = (long)*iter;
885 m_index3++;
886 if ( m_index3 > m_index3->range().distance() )
887 {
888 break;
889 cerr << "*********** too many index ************" << endl;
890 }
891 }
892
893 vtmp = m_pIBGT->getMuchitSeg();
894 m_index4 = 0;
895 for ( std::vector<int>::iterator iter = vtmp.begin(); iter != vtmp.end(); iter++ )
896 {
897 m_hitSeg[m_index4] = *( iter );
898 m_index4++;
899 if ( m_index4 > m_index4->range().distance() )
900 {
901 break;
902 cerr << "*********** too many index ************" << endl;
903 }
904 }
905 } // end fill ntuple
906
907 //---------------------------------------------------
908 // write out event number which not passed trigger.
909 //---------------------------------------------------
910 if ( ifoutEvtId == 1 )
911 {
912 ofstream eventnum( outEvtId.c_str(), ios_base::app );
913 if ( !ifpass ) eventnum << event << endl;
914 eventnum.close();
915 }
916
917 //-------------------------------------------------
918 // write out event number which passed trigger.
919 //-------------------------------------------------
920 if ( ifoutEvtId == 2 )
921 {
922 ofstream eventnum( outEvtId.c_str(), ios_base::app );
923 if ( ifpass ) eventnum << event << endl;
924 eventnum.close();
925 }
926
927 //--------------------------------------------------------
928 // write out events (passed trigger) into an ascii file
929 //--------------------------------------------------------
930 if ( writeFile == 1 )
931 {
932 EVENT asciiEvt;
933 readin >> asciiEvt;
934 if ( asciiEvt.header.eventNo == event )
935 {
936 if ( ifpass == true ) readout << asciiEvt << endl;
937 }
938 else
939 cout << "********* Event No. from AsciiFile do not equal Event No. from TDS "
940 << asciiEvt.header.eventNo << " " << event << endl;
941 }
942
943 //--------------------------------------------------------------
944 // if it is offline mode, register trigger information into TDS
945 //--------------------------------------------------------------
946 if ( m_runMode == 1 )
947 {
948 const int* trigcond = m_pIBGT->getTrigCond();
949 const int* trigchan = m_pIBGT->getTrigChan();
950 int window = 0;
951 int timing = 0;
952 bool preScale = false;
953
954 StatusCode sc = StatusCode::SUCCESS;
955 TrigEvent* aTrigEvent = new TrigEvent;
956 sc = eventSvc()->registerObject( "/Event/Trig", aTrigEvent );
957 if ( sc != StatusCode::SUCCESS )
958 { log << MSG::DEBUG << "Could not register TrigEvent, you can ignore." << endmsg; }
959
960 TrigData* aTrigData = new TrigData( window, timing, trigcond, trigchan, preScale );
961 sc = eventSvc()->registerObject( "/Event/Trig/TrigData", aTrigData );
962 if ( sc != StatusCode::SUCCESS )
963 { log << MSG::ERROR << "Could not register TrigData!!!!!" << endmsg; }
964 }
965
966 return StatusCode::SUCCESS;
967}
968
969StatusCode BesTrigL1::finalize() {
970
971 MsgStream msg( msgSvc(), name() );
972 msg << MSG::DEBUG << "==> Finalize BesTrigL1" << endmsg;
973
974 if ( writeFile == 1 )
975 {
976 readin.close();
977 readout.close();
978 }
979 cout << "There are total " << passNo << " event pass trigger" << endl;
980 return StatusCode::SUCCESS;
981}
982
983void BesTrigL1::findSETime( multimap<int, uint32_t, less<int>> mdc_hitmap,
984 multimap<int, int, less<int>> tof_hitmap,
985 multimap<int, uint32_t, less<int>> emc_TC, double& stime,
986 double& etime ) {
987 std::multimap<int, uint32_t, less<int>>::iterator mdc_iter = mdc_hitmap.begin();
988 double smdctime = -1, emdctime = -1;
989 if ( mdc_hitmap.size() != 0 )
990 {
991 smdctime = ( mdc_iter->first ) * 0.09375;
992 mdc_iter = mdc_hitmap.end();
993 mdc_iter--;
994 emdctime = ( mdc_iter->first ) * 0.09375;
995 }
996
997 std::multimap<int, int, less<int>>::iterator tof_iter = tof_hitmap.begin();
998 double stoftime = -1, etoftime = -1;
999 if ( tof_hitmap.size() != 0 )
1000 {
1001 stoftime = ( tof_iter->first );
1002 tof_iter = tof_hitmap.end();
1003 tof_iter--;
1004 etoftime = ( tof_iter->first );
1005 }
1006
1007 std::multimap<int, uint32_t, less<int>>::iterator emc_iter = emc_TC.begin();
1008 double semctime = -1, eemctime = -1;
1009 if ( emc_TC.size() != 0 )
1010 {
1011 semctime = ( emc_iter->first );
1012 emc_iter = emc_TC.end();
1013 emc_iter--;
1014 eemctime = ( emc_iter->first );
1015 }
1016
1017 stime = -1, etime = -1;
1018 if ( smdctime >= 0 && stoftime >= 0 )
1019 {
1020 if ( smdctime > stoftime ) stime = stoftime;
1021 else stime = smdctime;
1022
1023 if ( ( emdctime + 800 ) > ( etoftime + 24 ) ) etime = emdctime + 800;
1024 else etime = etoftime + 24;
1025 }
1026 else if ( smdctime < 0 && stoftime >= 0 )
1027 {
1028 stime = stoftime;
1029 etime = etoftime + 24;
1030 }
1031 else if ( smdctime >= 0 && stoftime < 0 )
1032 {
1033 stime = smdctime;
1034 etime = emdctime + 800;
1035 }
1036 else
1037 {
1038 stime = -1;
1039 etime = -1;
1040 }
1041 // compare with emc time
1042 if ( semctime >= 0 && stime >= 0 )
1043 {
1044 if ( semctime > stime ) stime = stime;
1045 else stime = semctime;
1046
1047 if ( ( eemctime + 16 * 24 ) > etime ) etime = eemctime + 16 * 24;
1048 else etime = etime;
1049 }
1050 else if ( semctime < 0 && stime >= 0 )
1051 {
1052 stime = stime;
1053 etime = etime;
1054 }
1055 else if ( semctime >= 0 && stime < 0 )
1056 {
1057 stime = semctime;
1058 etime = eemctime + 16 * 24;
1059 }
1060 else
1061 {
1062 stime = -1;
1063 etime = -1;
1064 }
1065}
1066
1067void BesTrigL1::runAclock_mdc( int iclock, double stime,
1068 multimap<int, uint32_t, less<int>> mdc_hitmap ) {
1069 std::vector<int> vmdcHit;
1070 vmdcHit.clear();
1071
1072 std::multimap<int, uint32_t, less<int>>::iterator mdc_iter = mdc_hitmap.begin();
1073 // int beginclock = int ((mdc_iter->first)*0.09375/24);
1074
1075 //--------------------------
1076 // consider mdc noise
1077 //--------------------------
1078 /*
1079 if((iclock - beginclock) >= 0 && (iclock - beginclock) <= 33) {
1080 for(int i = 0; i < 16; i++) {
1081 for(int hit_id = 0; hit_id < 256; hit_id++) {
1082 int layer, wire;
1083 double ratio = -1;
1084 if(i == 0) layer = 8;
1085 if(i == 1) layer = 9;
1086 if(i == 2) layer = 10;
1087 if(i == 3) layer = 11;
1088 if(i == 4) layer = 12;
1089 if(i == 5) layer = 13;
1090 if(i == 6) layer = 14;
1091 if(i == 7) layer = 15;
1092 if(i == 8) layer = 16;
1093 if(i == 9) layer = 17;
1094 if(i == 10) layer = 18;
1095 if(i == 11) layer = 19;
1096 if(i == 12) layer = 36;
1097 if(i == 13) layer = 37;
1098 if(i == 14) layer = 38;
1099 if(i == 15) layer = 39;
1100
1101 if(hit_id < 76) {
1102 if(i == 0) ratio = hit9[hit_id];
1103 if(i == 1) ratio = hit10[hit_id];
1104 }
1105 if(hit_id < 88) {
1106 if(i == 2) ratio = hit11[hit_id];
1107 if(i == 3) ratio = hit12[hit_id];
1108 }
1109 if(hit_id < 100) {
1110 if(i == 4) ratio = hit13[hit_id];
1111 if(i == 5) ratio = hit14[hit_id];
1112 }
1113 if(hit_id < 112) {
1114 if(i == 6) ratio = hit15[hit_id];
1115 if(i == 7) ratio = hit16[hit_id];
1116 }
1117 if(hit_id < 128) {
1118 if(i == 8) ratio = hit17[hit_id];
1119 if(i == 9) ratio = hit18[hit_id];
1120 }
1121 if(hit_id < 140) {
1122 if(i == 10) ratio = hit19[hit_id];
1123 if(i == 11) ratio = hit20[hit_id];
1124 }
1125 if(i == 12) ratio = hit37[hit_id];
1126 if(i == 13) ratio = hit38[hit_id];
1127 if(i == 14) ratio = hit39[hit_id];
1128 if(i == 15) ratio = hit40[hit_id];
1129
1130 wire = hit_id;
1131
1132 if(RandFlat::shoot() < ratio*(33 - iclock)*24/2000.) {
1133 vmdcHit.push_back(layer);
1134 vmdcHit.push_back(wire);
1135 }
1136 }
1137 }
1138 }
1139 */
1140
1141 for ( std::multimap<int, uint32_t, less<int>>::iterator mdc_iter = mdc_hitmap.begin();
1142 mdc_iter != mdc_hitmap.end(); mdc_iter++ )
1143 {
1144 double time = ( mdc_iter->first ) * 0.09375;
1145 if ( ( time < ( stime + ( iclock + 1 ) * 24. ) ) &&
1146 ( time + 800. ) > ( stime + iclock * 24. ) )
1147 {
1148 uint32_t mdcId = mdc_iter->second;
1149 int layer = ( mdcId & 0xFFFF0000 ) >> 16;
1150 int cell = mdcId & 0xFFFF;
1151 bool firstdc = true;
1152 // for(std::multimap<int,int,less<int> >::iterator tmp_mdc = mdc_hitmap.begin(); tmp_mdc
1153 // != mdc_iter; tmp_mdc++) {
1154 // if(mdcId == (tmp_mdc->second)) firstdc = false;
1155 // }
1156 if ( firstdc == true )
1157 {
1158 vmdcHit.push_back( layer );
1159 vmdcHit.push_back( cell );
1160 }
1161 }
1162 }
1163
1164 // set mdc vector hit
1165 m_MdcTSF->setMdcDigi( vmdcHit );
1166 m_pIBGT->startMdcTrig();
1167}
1168
1169void BesTrigL1::runAclock_tof( int iclock, double stime, int& idle_status,
1170 std::multimap<int, int, less<int>> tof_hitmap ) {
1171 std::vector<int> vtofHit;
1172 vtofHit.clear();
1173
1174 // tof trigger
1175 if ( idle_status != -1 && ( iclock - idle_status ) == 3 ) idle_status = -1;
1176 for ( std::multimap<int, int, less<int>>::iterator tof_iter = tof_hitmap.begin();
1177 tof_iter != tof_hitmap.end(); tof_iter++ )
1178 {
1179 double time = ( tof_iter->first ); // ns
1180 if ( idle_status == -1 )
1181 {
1182 if ( time < ( stime + ( iclock + 1 ) * 24 ) && time >= ( stime + iclock * 24 ) )
1183 {
1184 // if(time < (stime + (iclock + 1)*24) && (time + 24) > (stime + iclock*24)) {
1185 // //stretch signal
1186 vtofHit.push_back( tof_iter->second );
1187 }
1188 }
1189 else
1190 {
1191 if ( ( iclock - idle_status ) == 1 )
1192 {
1193 if ( ( time < ( stime + ( iclock + 1 ) * 24 ) && time >= ( stime + iclock * 24 ) ) ||
1194 ( time < ( stime + iclock * 24 ) && time >= ( stime + ( iclock - 1 ) * 24 ) ) )
1195 { vtofHit.push_back( tof_iter->second ); }
1196 }
1197 if ( ( iclock - idle_status ) == 2 )
1198 {
1199 if ( ( time < ( stime + ( iclock + 1 ) * 24 ) && time >= ( stime + iclock * 24 ) ) ||
1200 ( time < ( stime + iclock * 24 ) && time >= ( stime + ( iclock - 1 ) * 24 ) ) ||
1201 ( time < ( stime + ( iclock - 1 ) * 24 ) &&
1202 time >= ( stime + ( iclock - 2 ) * 24 ) ) )
1203 { vtofHit.push_back( tof_iter->second ); }
1204 }
1205 }
1206 }
1207 if ( idle_status == -1 && vtofHit.size() != 0 ) idle_status = iclock;
1208
1209 // set tof vector hit
1210 m_TofHitCount->setTofDigi( vtofHit );
1211 m_pIBGT->startTofTrig();
1212}
1213
1214void BesTrigL1::runAclock_emc( int iclock, double stime,
1215 std::multimap<int, uint32_t, less<int>> emc_TC,
1216 EmcWaveform* blockWave ) {
1217 std::vector<uint32_t> vemcClus;
1218 std::vector<double> vemcBlkE;
1219
1220 vemcClus.clear();
1221 vemcBlkE.clear();
1222 // std::cout << "iclock, emc_TC size: " << iclock << ", " << emc_TC.size() << std::endl;
1223 // cluster finding in emc trigger
1224 for ( std::multimap<int, uint32_t, less<int>>::iterator emc_iter = emc_TC.begin();
1225 emc_iter != emc_TC.end(); emc_iter++ )
1226 {
1227 double time = ( emc_iter->first );
1228 if ( ( time < ( stime + ( iclock + 1 ) * 24. ) ) &&
1229 ( time + 16 * 24 ) > ( stime + iclock * 24. ) )
1230 { vemcClus.push_back( emc_iter->second ); }
1231 }
1232
1233 // energy adding in emc trigger
1234 for ( int blockId = 0; blockId < 16; blockId++ )
1235 {
1236 double block_ADC = ( blockWave[blockId] ).getADCTrg( (int)stime + iclock * 24 );
1237 vemcBlkE.push_back( block_ADC );
1238 // std::cout << " block_ADC: " << block_ADC << std::endl;
1239 }
1240 // std::cout << "iclock,stime,vemcClus size: " << iclock << "," << stime << ", " <<
1241 // vemcClus.size() << std::endl;
1242 m_emcDigi->setEmcTC( vemcClus );
1243 m_emcDigi->setEmcBE( vemcBlkE ); // set block energy
1244 // start EMC trigger logic
1245 m_pIBGT->startEmcTrig();
1246}
1247
1248void BesTrigL1::getEmcAnalogSig( EmcDigiCol* emcDigiCol, EmcWaveform ( &blockWave )[16],
1249 multimap<int, uint32_t, less<int>>& emc_TC ) {
1250 EmcWaveform eewave[32];
1251 EmcWaveform wewave[32];
1252 EmcWaveform bwave[11][30];
1253
1254 for ( int i = 0; i < 11; i++ )
1255 {
1256 for ( int j = 0; j < 30; j++ ) { bwave[i][j].makeWaveformTrg( 0, 0 ); }
1257 }
1258 for ( int i = 0; i < 32; i++ )
1259 {
1260 if ( i < 16 ) blockWave[i].makeWaveformTrg( 0, 0 );
1261 eewave[i].makeWaveformTrg( 0, 0 );
1262 wewave[i].makeWaveformTrg( 0, 0 );
1263 }
1264
1265 for ( EmcDigiCol::iterator iter3 = emcDigiCol->begin(); iter3 != emcDigiCol->end(); iter3++ )
1266 {
1267 Identifier id = ( *iter3 )->identify();
1268 unsigned int module = EmcID::barrel_ec( id );
1269 unsigned int theta = EmcID::theta_module( id );
1270 unsigned int phi = EmcID::phi_module( id );
1271
1272 int index = emcCalibConstSvc->getIndex( module, theta, phi );
1273 double trgGain = m_RealizationSvc->getTrgGain( index );
1274 double adc = (double)( *iter3 )->getChargeChannel();
1275 double mv = RandGauss::shoot( 978., 14. );
1276
1277 if ( ( *iter3 )->getMeasure() == 0 ) adc = adc * 2 * mv * 2 / 65535. * ( trgGain );
1278 else if ( ( *iter3 )->getMeasure() == 1 ) adc = adc * 16 * mv * 2 / 65535 * ( trgGain );
1279 else adc = adc * 64 * mv * 2 / 65535 * ( trgGain );
1280
1281 unsigned int tdc = ( *iter3 )->getTimeChannel();
1282 int theTC = m_emcDigi->getTCThetaId( module, theta, phi );
1283 int phiTC = m_emcDigi->getTCPhiId( module, theta, phi );
1284 EmcWaveform wave1;
1285 if ( module == 0 )
1286 {
1287 wave1.makeWaveformTrg( adc, tdc + 80 );
1288 eewave[phiTC] += wave1;
1289 }
1290 if ( module == 1 )
1291 {
1292 wave1.makeWaveformTrg( adc, tdc + 80 );
1293 bwave[theTC][phiTC] += wave1;
1294 }
1295 if ( module == 2 )
1296 {
1297 wave1.makeWaveformTrg( adc, tdc + 80 );
1298 wewave[phiTC] += wave1;
1299 }
1300 }
1301
1302 // find barrel cluster
1303 for ( int i = 0; i < 11; i++ )
1304 {
1305 for ( int j = 0; j < 30; j++ )
1306 {
1307 int time_low = bwave[i][j].frontEdgeTrg( m_pIBGT->getL1TC_GATE() );
1308 int time_high = bwave[i][j].frontEdgeTrg( m_pIBGT->getL1TC_THRESH() );
1309 int time = -1;
1310
1311 if ( time_high >= 0 )
1312 {
1313 if ( time_low * 50 + 1500 > time_high * 50 ) time = time_low * 50 + 1500;
1314 else time = time_high * 50;
1315 uint32_t TCID = ( 1 & 0xFF ) << 16;
1316 TCID = TCID | ( ( i & 0xFF ) << 8 );
1317 TCID = TCID | ( j & 0xFF );
1318 typedef pair<int, uint32_t> vpair;
1319 emc_TC.insert( vpair( time, TCID ) );
1320 // std::cout <<"i, j: " << i << ", " << j << " time: " << time << std::endl;
1321 }
1322 if ( time_low >= 0 )
1323 {
1324 int blockId = m_emcDigi->getBLKId( i, j );
1325 blockWave[blockId + 2] += bwave[i][j];
1326 }
1327 }
1328 }
1329 // find end cap cluster
1330 for ( int i = 0; i < 32; i++ )
1331 {
1332 // east end cap
1333 int time_low1 = eewave[i].frontEdgeTrg( m_pIBGT->getL1TC_GATE() );
1334 int time_high1 = eewave[i].frontEdgeTrg( m_pIBGT->getL1TC_THRESH() );
1335 int time1 = -1;
1336 if ( time_high1 >= 0 )
1337 {
1338 if ( time_low1 * 50 + 1500 > time_high1 * 50 ) time1 = time_low1 * 50 + 1500;
1339 else time1 = time_high1 * 50;
1340 uint32_t TCID1 = ( 0 & 0xFF ) << 16;
1341 TCID1 = TCID1 | ( ( 0 & 0xFF ) << 8 );
1342 TCID1 = TCID1 | ( i & 0xFF );
1343 typedef pair<int, uint32_t> vpair;
1344 emc_TC.insert( vpair( time1, TCID1 ) );
1345 }
1346 if ( time_low1 >= 0 )
1347 {
1348 if ( i < 16 ) blockWave[0] += eewave[i];
1349 else blockWave[1] += eewave[i];
1350 }
1351 // west end cap
1352 int time_low2 = wewave[i].frontEdgeTrg( m_pIBGT->getL1TC_GATE() );
1353 int time_high2 = wewave[i].frontEdgeTrg( m_pIBGT->getL1TC_THRESH() );
1354 int time2 = -1;
1355 if ( time_high2 >= 0 )
1356 {
1357 if ( time_low2 * 50 + 1500 > time_high2 * 50 ) time2 = time_low2 * 50 + 1500;
1358 else time2 = time_high2 * 50;
1359 uint32_t TCID2 = ( 2 & 0xFF ) << 16;
1360 TCID2 = TCID2 | ( ( 0 & 0xFF ) << 8 );
1361 TCID2 = TCID2 | ( i & 0xFF );
1362 typedef pair<int, uint32_t> vpair;
1363 emc_TC.insert( vpair( time2, TCID2 ) );
1364 }
1365 if ( time_low2 >= 0 )
1366 {
1367 if ( i < 16 ) blockWave[14] += wewave[i];
1368 else blockWave[15] += wewave[i];
1369 }
1370 }
1371}
1372
1373void BesTrigL1::findEmcPeakTime( double& peak_time, EmcWaveform* blockWave ) {
1374 double tot_block_max = 0;
1375 for ( int i = 0; i < 16; i++ )
1376 {
1377 int block_time;
1378 double block_max = blockWave[i].max( block_time );
1379 tot_block_max += block_max;
1380 }
1381
1382 for ( int i = 0; i < 16; i++ )
1383 {
1384 if ( tot_block_max == 0 ) break;
1385 int block_time;
1386 double block_max = blockWave[i].max( block_time );
1387 block_time = block_time * 50;
1388 peak_time += block_max / tot_block_max * block_time;
1389 }
1390}
1391
1392void BesTrigL1::stretchTrgCond( int nclock, int**& trgcond ) {
1393 int emc_clus = 34;
1394 int emc_ener = 50;
1395 int mdc = 34;
1396 int mdc_n = 68;
1397 int tof = 4;
1398 for ( int icond = 0; icond < 48; icond++ )
1399 {
1400 int sclock = -1;
1401 bool retrig = false;
1402 for ( int iclock = 0; iclock < nclock; iclock++ )
1403 {
1404 if ( icond < 16 )
1405 { // stretch emc trigger conditions
1406 if ( icond < 7 || icond == 12 || icond == 13 )
1407 { // stretch cluster trigger conditions
1408 if ( sclock != -1 && iclock - sclock == emc_clus ) sclock = -1;
1409 if ( sclock == -1 && trgcond[icond][iclock] > 0 )
1410 {
1411 if ( iclock == 0 ) sclock = iclock;
1412 else
1413 {
1414 if ( trgcond[icond][iclock] * trgcond[icond][iclock - 1] == 0 ) sclock = iclock;
1415 }
1416 }
1417 if ( sclock != -1 && iclock - sclock < emc_clus ) trgcond[icond][iclock] = 1;
1418 }
1419 else
1420 { // stretch emc energy trigger conditions, re-triggering is available
1421 if ( sclock != -1 && iclock - sclock == emc_ener ) sclock = -1;
1422 if ( sclock == -1 && trgcond[icond][iclock] > 0 )
1423 {
1424 if ( iclock == 0 ) sclock = iclock;
1425 else
1426 {
1427 if ( trgcond[icond][iclock] * trgcond[icond][iclock - 1] == 0 ) sclock = iclock;
1428 }
1429 }
1430 if ( sclock != -1 && iclock - sclock < emc_ener && trgcond[icond][iclock] == 0 )
1431 retrig = true;
1432 if ( retrig == true )
1433 {
1434 if ( trgcond[icond][iclock] > 0 )
1435 {
1436 sclock = iclock;
1437 retrig = false;
1438 }
1439 }
1440 if ( sclock != -1 && iclock - sclock < emc_ener ) trgcond[icond][iclock] = 1;
1441 }
1442 }
1443 else if ( icond >= 16 && icond < 23 )
1444 { // stretch tof trigger conditions
1445 if ( sclock != -1 && iclock - sclock == tof ) sclock = -1;
1446 if ( sclock == -1 && trgcond[icond][iclock] > 0 )
1447 {
1448 if ( iclock == 0 ) sclock = iclock;
1449 else
1450 {
1451 if ( trgcond[icond][iclock] * trgcond[icond][iclock - 1] == 0 ) sclock = iclock;
1452 }
1453 }
1454 if ( sclock != -1 && iclock - sclock < tof ) trgcond[icond][iclock] = 1;
1455 }
1456 else if ( icond >= 38 )
1457 { // stretch mdc trigger conditions
1458 if ( icond == 39 || icond == 43 )
1459 {
1460 if ( sclock != -1 && iclock - sclock == mdc_n ) sclock = -1;
1461 if ( sclock == -1 && trgcond[icond][iclock] > 0 )
1462 {
1463 if ( iclock == 0 ) sclock = iclock;
1464 else
1465 {
1466 if ( trgcond[icond][iclock] * trgcond[icond][iclock - 1] == 0 ) sclock = iclock;
1467 }
1468 }
1469 if ( sclock != -1 && iclock - sclock < mdc_n ) trgcond[icond][iclock] = 1;
1470 }
1471 else
1472 {
1473 if ( sclock != -1 && iclock - sclock == mdc ) sclock = -1;
1474 if ( sclock == -1 && trgcond[icond][iclock] > 0 )
1475 {
1476 if ( iclock == 0 ) sclock = iclock;
1477 else
1478 {
1479 if ( trgcond[icond][iclock] * trgcond[icond][iclock - 1] == 0 ) sclock = iclock;
1480 }
1481 }
1482 if ( sclock != -1 && iclock - sclock < mdc ) trgcond[icond][iclock] = 1;
1483 }
1484 }
1485 else
1486 { // stretch other trigger conditions, including track match and muc
1487 }
1488 }
1489 }
1490}
1491
1492void BesTrigL1::trgSAFDelay( int nclock, int**& trgcond ) {
1493 // SAF delay time
1494 // int delay[48] = {31,31,31,31,31,31,31,7,7,7,7,7,31,31,7,7,
1495 // 135,135,135,135,135,135,135,83,83,83,6,6,6,83,83,83,
1496 // 97,97,97,97,97,97,86,87,85,87,83,85,83,85,122,122};
1497 int delay[48] = { 24, 24, 24, 24, 24, 24, 24, 7, 7, 7, 7, 7, 24, 24, 7, 7,
1498 0, 0, 0, 0, 0, 0, 0, 83, 83, 83, 6, 6, 6, 83, 83, 83,
1499 97, 97, 97, 97, 97, 97, 0, 0, 0, 0, 0, 0, 0, 0, 122, 122 };
1500
1501 for ( int icond = 0; icond < 48; icond++ )
1502 {
1503 for ( int iclock = nclock - 1; iclock >= 0; iclock-- )
1504 {
1505 if ( iclock < delay[icond] ) trgcond[icond][iclock] = 0;
1506 else { trgcond[icond][iclock] = trgcond[icond][iclock - delay[icond]]; }
1507 }
1508 }
1509}
1510
1511void BesTrigL1::trgGTLDelay( int nclock, int**& trgcond ) {
1512 // GTL delay time
1513 int delay[48] = { 1, 1, 1, 1, 1, 1, 1, 18, 18, 18, 18, 18, 1, 1, 18, 18,
1514 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18,
1515 14, 14, 14, 14, 14, 14, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10 };
1516
1517 for ( int icond = 0; icond < 48; icond++ )
1518 {
1519 for ( int iclock = nclock - 1; iclock >= 0; iclock-- )
1520 {
1521 if ( iclock < delay[icond] ) trgcond[icond][iclock] = 0;
1522 else trgcond[icond][iclock] = trgcond[icond][iclock - delay[icond]];
1523 }
1524 }
1525}
mdc
std::string mdc
Definition CalibModel.cxx:46
time
Double_t time
Definition DataBase/tau_mode.c:7
EmcDigiCol
ObjectVector< EmcDigi > EmcDigiCol
Definition Emc/EmcRawEvent/include/EmcRawEvent/EmcDigi.h:43
TofDataMap
std::multimap< unsigned int, TofData * > TofDataMap
Definition Event/RawDataProviderSvc/include/RawDataProviderSvc/TofData.h:237
imag
double imag(const EvtComplex &c)
Definition EvtComplex.hh:207
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
BesGlobalTrigSvc::setRunMode
void setRunMode(int mode)
Definition BesGlobalTrigSvc.h:109
BesTrigL1::execute
virtual StatusCode execute()
Algorithm execution.
Definition BesTrigL1.cxx:291
BesTrigL1::runAclock_emc
void runAclock_emc(int iclock, double stime, std::multimap< int, uint32_t, less< int > > emc_TC, EmcWaveform *blockWave)
Definition BesTrigL1.cxx:1214
BesTrigL1::findEmcPeakTime
void findEmcPeakTime(double &peak_time, EmcWaveform *blockWave)
Definition BesTrigL1.cxx:1373
BesTrigL1::BesTrigL1
BesTrigL1(const std::string &name, ISvcLocator *pSvcLocator)
Standard constructor.
Definition BesTrigL1.cxx:57
BesTrigL1::runAclock_tof
void runAclock_tof(int iclock, double stime, int &idle_status, std::multimap< int, int, less< int > > tof_hitmap)
Definition BesTrigL1.cxx:1169
BesTrigL1::getEmcAnalogSig
void getEmcAnalogSig(EmcDigiCol *emcDigiCol, EmcWaveform(&blockWave)[16], multimap< int, uint32_t, less< int > > &emc_TC)
Definition BesTrigL1.cxx:1248
BesTrigL1::trgGTLDelay
void trgGTLDelay(int nclock, int **&trgcond)
Definition BesTrigL1.cxx:1511
BesTrigL1::trgSAFDelay
void trgSAFDelay(int nclock, int **&trgcond)
Definition BesTrigL1.cxx:1492
BesTrigL1::finalize
virtual StatusCode finalize()
Algorithm finalization.
Definition BesTrigL1.cxx:969
BesTrigL1::findSETime
void findSETime(multimap< int, uint32_t, less< int > > mdc_hitmap, multimap< int, int, less< int > > tof_hitmap, multimap< int, uint32_t, less< int > > emc_TC, double &stime, double &etime)
Definition BesTrigL1.cxx:983
BesTrigL1::initialize
virtual StatusCode initialize()
Destructor.
Definition BesTrigL1.cxx:71
BesTrigL1::stretchTrgCond
void stretchTrgCond(int nclock, int **&trgcond)
Definition BesTrigL1.cxx:1392
BesTrigL1::runAclock_mdc
void runAclock_mdc(int iclock, double stime, multimap< int, uint32_t, less< int > > mdc_hitmap)
Definition BesTrigL1.cxx:1067
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
Definition EmcID.cxx:36
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition EmcID.cxx:41
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition EmcID.cxx:46
EmcTCFinder::get_Emc
static EmcTCFinder * get_Emc(void)
Definition EmcTCFinder.cxx:32
EmcWaveform::max
double max(int &binOfMax) const
Definition EmcWaveform.cxx:72
EmcWaveform::makeWaveformTrg
void makeWaveformTrg(double energy, double time)
Definition EmcWaveform.cxx:108
EmcWaveform::frontEdgeTrg
int frontEdgeTrg(double thres)
Definition EmcWaveform.cxx:95
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0).
Definition MdcID.cxx:47
MdcID::wire
static int wire(const Identifier &id)
Definition MdcID.cxx:52
MdcTSF::get_Mdc
static MdcTSF * get_Mdc(void)
Definition MdcTSF.cxx:28
MucTrigHit::get_Muc
static MucTrigHit * get_Muc(void)
Definition MucTrigHit.cxx:24
TofData::tdc2
double tdc2()
Definition TofData.cxx:573
TofData::tdc1
double tdc1()
Definition TofData.cxx:561
TofHitCount::get_Tof
static TofHitCount * get_Tof(void)
Definition TofHitCount.cxx:26
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